emcsfc_snow2mdl  1.4.0
 All Data Structures Files Functions Variables
snow2mdl.F90
Go to the documentation of this file.
1 
4 
20  module snow2mdl
21 
22  use program_setup, only : lat_threshold, &
23  model_snow_file, &
24  min_snow_depth, &
25  snow_cvr_threshold, &
26  grib_day, &
27  grib_century, &
28  grib_hour, &
29  grib_month, &
30  grib_year, &
31  output_grib2
32 
33  use model_grid, only : resol_mdl, &
34  imdl, &
35  jmdl, &
36  ijmdl, &
37  ipts_mdl, &
38  jpts_mdl, &
39  lsmask_mdl, &
40  lsmask_mdl_sav, &
41  lats_mdl, &
42  lat11, latlast, &
43  lon11, lonlast, &
44  lons_mdl, &
45  kgds_mdl, &
46  grid_id_mdl, &
47  thinned, &
48  lonsperlat_mdl
49 
50  use snowdat, only : nesdis_res, &
51  afwa_res, &
52  autosnow_res, &
53  inesdis, &
54  jnesdis, &
55  mesh_nesdis, &
56  snow_cvr_nesdis, &
57  sea_ice_nesdis, &
58  bitmap_nesdis, &
59  iafwa, &
60  jafwa, &
61  snow_dep_afwa_global, &
62  snow_dep_afwa_nh, &
63  snow_dep_afwa_sh, &
64  bitmap_afwa_global, &
65  bitmap_afwa_nh, &
66  bitmap_afwa_sh, &
67  iautosnow, &
68  jautosnow, &
69  snow_cvr_autosnow, &
70  bitmap_autosnow, &
71  use_global_afwa, &
72  use_nh_afwa, &
73  use_sh_afwa, &
74  use_nesdis, &
75  use_autosnow, &
76  kgds_nesdis, &
77  kgds_afwa_global, &
78  kgds_afwa_nh, &
79  kgds_afwa_sh, &
80  kgds_autosnow, &
81  bad_afwa_nh, bad_afwa_sh
82 
83  private
84 
85  real, allocatable :: snow_cvr_mdl(:,:)
86  real, allocatable :: snow_dep_mdl(:,:)
87 
88  public :: interp
89 
90  contains
159  subroutine interp
160  use gdswzd_mod
161 
162  implicit none
163 
164  integer :: i, j, ii, jj, ij
165  integer :: ijmdl2, istart, iend, imid, iii
166  integer, allocatable :: idum(:,:)
167  integer :: int_opt, ipopt(20)
168  integer :: kgds_mdl_tmp(200)
169  integer :: no, ibo, iret, nret
170 
171  logical*1, allocatable :: bitmap_mdl(:)
172 
173  real :: gridi(1)
174  real :: gridj(1)
175  real, allocatable :: lsmask_1d(:)
176  real, allocatable :: snow_cvr_mdl_1d(:)
177  real, allocatable :: snow_dep_mdl_tmp(:)
178  real :: sumc, sumd, x1, r, fraction, gridis, gridie
179  real, parameter :: undefined_value = -999.
180 
181 !----------------------------------------------------------------------
182 ! for model grids fully or partially located in the southern
183 ! hemisphere, the user must select sh data (either autosnow or afwa).
184 ! this restriction is relaxed for the ndas domain, which is
185 ! has a small part below the equator.
186 !----------------------------------------------------------------------
187 
188  if (minval(lats_mdl) < -10.0 .and. .not. use_global_afwa .and. .not. use_sh_afwa .and. .not. use_autosnow) then
189  print*,"- FATAL ERROR: MUST SELECT EITHER AFWA OR AUTOSNOW DATA FOR MODEL GRID WITH SH POINTS."
190  call w3tage('SNOW2MDL')
191  call errexit(54)
192  endif
193 
194 !----------------------------------------------------------------------
195 ! if model grid is totally within the southern hemisphere, set flags
196 ! so that nh afwa and nh nesdis/ims is not processed. these flags
197 ! are set to false if user inadvertantly selects these data sources.
198 !----------------------------------------------------------------------
199 
200  if (maxval(lats_mdl) < 0.0) then
201  use_nh_afwa=.false.
202  use_nesdis=.false.
203  endif
204 
205 !----------------------------------------------------------------------
206 ! if model grid is partially or fully located in the northern
207 ! hemisphere, the user must select either nesdis/ims or nh afwa data.
208 !----------------------------------------------------------------------
209 
210  if (maxval(lats_mdl) > 0.0 .and. .not. use_global_afwa .and. .not. use_nh_afwa .and. .not. use_nesdis) then
211  print*,"- FATAL ERROR: MUST SELECT EITHER NESDIS/IMS OR AFWA DATA FOR MODEL GRID WITH NH POINTS."
212  call w3tage('SNOW2MDL')
213  call errexit(54)
214  endif
215 
216 !----------------------------------------------------------------------
217 ! if model grid is totally within the northern hemisphere, set flag
218 ! so that sh data is not processed. these flags are set to false
219 ! if user inadvertantly selects this data.
220 !----------------------------------------------------------------------
221 
222  if (minval(lats_mdl) > 0.0) then ! is model grid totally within northern hemisphere?
223  use_sh_afwa=.false.
224  use_autosnow=.false.
225  endif
226 
227 !----------------------------------------------------------------------
228 ! if selected, interpolate nesdis/ims data to model grid.
229 !----------------------------------------------------------------------
230 
231  nesdis_ims : if (use_nesdis) then
232 
233  ipopt = 0
234  if (nesdis_res < (0.5*resol_mdl)) then
235  print*,"- INTERPOLATE NH NESDIS/IMS DATA TO MODEL GRID USING BUDGET METHOD."
236  ipopt(1)=2 ! break model grid cell into 25 points.
237  ipopt(2:4)=1 ! 25 points are weighted equally.
238  ipopt(5)=10 ! 10% coverage of valid data in box
239  ipopt(20) = nint(100.0 / nesdis_res) + 1 ! search box width of 100 km.
240  kgds_mdl_tmp = kgds_mdl
241  kgds_mdl_tmp(1) = kgds_mdl_tmp(1) - 255 ! subset of grid
242  int_opt = 3
243  no = ijmdl
244  else
245  print*,"- INTERPOLATE NH NESDIS/IMS DATA TO MODEL GRID USING NEIGHBOR METHOD."
246  ipopt(1) = nint(100.0 / nesdis_res) + 1 ! search box width of 100 km.
247  kgds_mdl_tmp = kgds_mdl
248  kgds_mdl_tmp(1) = -1 ! for subsection of model grid.
249  int_opt = 2
250  no = ijmdl ! an input when kgds(1) < 0 (subset of grid)
251  end if
252 
253  allocate (snow_cvr_mdl_1d(ijmdl))
254  snow_cvr_mdl_1d = 0.0
255 
256  allocate (bitmap_mdl(ijmdl))
257  bitmap_mdl=.false. ! if interpolation routine can't find data
258  ! at a point, this flag is false.
259 
260  call ipolates(int_opt, ipopt, kgds_nesdis, kgds_mdl_tmp, &
261  (inesdis*jnesdis), ijmdl, &
262  1, 1, bitmap_nesdis, snow_cvr_nesdis, &
263  no, lats_mdl, lons_mdl, ibo, bitmap_mdl, &
264  snow_cvr_mdl_1d, iret)
265 
266  deallocate (bitmap_nesdis, snow_cvr_nesdis)
267 
268  if (iret /= 0) then
269  print*,"- FATAL ERROR: IN INTERPOLATION ROUTINE. IRET IS: ", iret
270  call w3tage('SNOW2MDL')
271  call errexit(55)
272  endif
273 
274 !-----------------------------------------------------------------------
275 ! if the interpolation routines did not find valid nesdis/ims data
276 ! in the vicinity of the model point, need to set a default value
277 ! of snow cover. south of user-defined latitude threshold, set
278 ! to zero. otherwise, see if the nearest neighbor nesdis/ims point is
279 ! sea ice. if so, assume model point is snow covered.
280 !-----------------------------------------------------------------------
281 
282  do ij = 1, ijmdl
283  if (lats_mdl(ij) < 0.0) cycle ! only consider nh model points
284  if (.not. bitmap_mdl(ij)) then
285  if (lats_mdl(ij) <= lat_threshold) then
286  snow_cvr_mdl_1d(ij) = 0.0
287  else
288  call gdswzd(kgds_nesdis,-1,1,undefined_value,gridi,gridj, &
289  lons_mdl(ij),lats_mdl(ij),nret)
290  if (nret /= 1) then
291  print*,"- WARNING: MODEL POINT OUTSIDE NESDIS/IMS GRID: ", ipts_mdl(ij), jpts_mdl(ij)
292  snow_cvr_mdl_1d(ij) = 0.0
293  else
294  ii = nint(gridi(1))
295  jj = nint(gridj(1))
296  if (sea_ice_nesdis(ii,jj) == 1) then
297  snow_cvr_mdl_1d(ij) = 100.0
298  else
299  snow_cvr_mdl_1d(ij) = 0.0
300  end if
301  end if
302  end if
303  end if
304  enddo
305 
306  deallocate (sea_ice_nesdis)
307  deallocate (bitmap_mdl)
308 
309  endif nesdis_ims
310 
311 !----------------------------------------------------------------------
312 ! now interpolate nh afwa snow depth data.
313 !----------------------------------------------------------------------
314 
315  global_afwa : if (use_global_afwa) then
316 
317 !----------------------------------------------------------------------
318 ! determine interpolation method based on the resolution of
319 ! afwa data and the model grid.
320 !----------------------------------------------------------------------
321 
322  ipopt = 0
323  if (afwa_res < (0.5*resol_mdl)) then
324  print*,"- INTERPOLATE GLOBAL AFWA DATA TO MODEL GRID USING BUDGET METHOD."
325  ipopt(1)=-1 ! break model grid cell into 25 points.
326  ipopt(2)=-1 ! 25 points are weighted equally.
327  ipopt(20) = nint(100.0 / afwa_res) + 1 ! search box width of 100 km.
328  kgds_mdl_tmp = kgds_mdl
329  kgds_mdl_tmp(1) = kgds_mdl_tmp(1) - 255 ! subset of grid
330  no = ijmdl
331  int_opt = 3
332  else
333  print*,"- INTERPOLATE GLOBAL AFWA DATA TO MODEL GRID USING NEIGHBOR METHOD."
334  ipopt(1) = nint(100.0 / afwa_res) + 1 ! search box width of 100 km.
335  kgds_mdl_tmp = kgds_mdl
336  kgds_mdl_tmp(1) = -1 ! for subsection of model grid.
337  int_opt = 2
338  no = ijmdl ! an input when kgds(1) < 0 (subset of grid)
339  end if
340 
341  allocate (snow_dep_mdl_tmp(ijmdl))
342  snow_dep_mdl_tmp = 0.0
343 
344  allocate (bitmap_mdl(ijmdl))
345  bitmap_mdl = .false.
346 
347  call ipolates(int_opt, ipopt, kgds_afwa_global, kgds_mdl_tmp, &
348  (iafwa*jafwa), ijmdl, &
349  1, 1, bitmap_afwa_global, snow_dep_afwa_global, &
350  no, lats_mdl, lons_mdl, ibo, bitmap_mdl, &
351  snow_dep_mdl_tmp, iret)
352 
353  deallocate(bitmap_afwa_global, snow_dep_afwa_global)
354 
355  if (iret /= 0) then
356  print*,"- FATAL ERROR IN INTERPOLATION ROUTINE. IRET IS: ", iret
357  call w3tage('SNOW2MDL')
358  call errexit(55)
359  endif
360 
361 !----------------------------------------------------------------------
362 ! if interpolation did not find afwa data near the model point, then
363 ! use a nominal value based on latitude threshold. this value
364 ! may be overwritten below depending on nesdis/ims cover data (if
365 ! user selects nesdis/ims data).
366 !----------------------------------------------------------------------
367 
368  do ij = 1, ijmdl
369  if (.not. bitmap_mdl(ij)) then
370  if (abs(lats_mdl(ij)) >= lat_threshold) then
371  snow_dep_mdl_tmp(ij) = min_snow_depth
372  else
373  snow_dep_mdl_tmp(ij) = 0.0
374  endif
375  endif
376  enddo
377 
378  deallocate(bitmap_mdl)
379 
380  endif global_afwa
381 
382 !----------------------------------------------------------------------
383 ! now interpolate nh afwa snow depth data.
384 !----------------------------------------------------------------------
385 
386  nh_afwa : if (use_nh_afwa) then
387 
388 !----------------------------------------------------------------------
389 ! determine interpolation method based on the resolution of
390 ! afwa data and the model grid.
391 !----------------------------------------------------------------------
392 
393  ipopt = 0
394  if (afwa_res < (0.5*resol_mdl)) then
395  print*,"- INTERPOLATE NH AFWA DATA TO MODEL GRID USING BUDGET METHOD."
396  ipopt(1)=-1 ! break model grid cell into 25 points.
397  ipopt(2)=-1 ! 25 points are weighted equally.
398  ipopt(20) = nint(100.0 / afwa_res) + 1 ! search box width of 100 km.
399  kgds_mdl_tmp = kgds_mdl
400  kgds_mdl_tmp(1) = kgds_mdl_tmp(1) - 255 ! subset of grid
401  no = ijmdl
402  int_opt = 3
403  else
404  print*,"- INTERPOLATE NH AFWA DATA TO MODEL GRID USING NEIGHBOR METHOD."
405  ipopt(1) = nint(100.0 / afwa_res) + 1 ! search box width of 100 km.
406  kgds_mdl_tmp = kgds_mdl
407  kgds_mdl_tmp(1) = -1 ! for subsection of model grid.
408  int_opt = 2
409  no = ijmdl ! an input when kgds(1) < 0 (subset of grid)
410  end if
411 
412  allocate (snow_dep_mdl_tmp(ijmdl))
413  snow_dep_mdl_tmp = 0.0
414 
415  allocate (bitmap_mdl(ijmdl))
416  bitmap_mdl = .false.
417 
418  call ipolates(int_opt, ipopt, kgds_afwa_nh, kgds_mdl_tmp, &
419  (iafwa*jafwa), ijmdl, &
420  1, 1, bitmap_afwa_nh, snow_dep_afwa_nh, &
421  no, lats_mdl, lons_mdl, ibo, bitmap_mdl, &
422  snow_dep_mdl_tmp, iret)
423 
424  deallocate(bitmap_afwa_nh, snow_dep_afwa_nh)
425 
426  if (iret /= 0) then
427  print*,"- FATAL ERROR IN INTERPOLATION ROUTINE. IRET IS: ", iret
428  call w3tage('SNOW2MDL')
429  call errexit(55)
430  endif
431 
432 !----------------------------------------------------------------------
433 ! if interpolation did not find afwa data near the model point, then
434 ! use a nominal value based on latitude threshold. this value
435 ! may be overwritten below depending on nesdis/ims cover data (if
436 ! user selects nesdis/ims data).
437 !----------------------------------------------------------------------
438 
439  do ij = 1, ijmdl
440  if (lats_mdl(ij) >= 0.) then ! only consider model pts in n hemi.
441  if (.not. bitmap_mdl(ij)) then
442  if (abs(lats_mdl(ij)) >= lat_threshold) then
443  snow_dep_mdl_tmp(ij) = min_snow_depth
444  else
445  snow_dep_mdl_tmp(ij) = 0.0
446  endif
447  endif
448  endif
449  enddo
450 
451  deallocate(bitmap_mdl)
452 
453  endif nh_afwa
454 
455 !----------------------------------------------------------------------
456 ! if nh data selected, use it to determine the cover and depth
457 ! on the model grid.
458 !----------------------------------------------------------------------
459 
460  allocate (snow_dep_mdl(imdl,jmdl))
461  allocate (snow_cvr_mdl(imdl,jmdl))
462  snow_cvr_mdl = 0.0
463  snow_dep_mdl = 0.0
464 
465  if (use_global_afwa .and. use_nesdis) then ! set depth/cover on nesdis ims/afwa blend
466  print*,"- BLEND NESDIS/IMS AND AFWA DATA IN NH."
467  do ij = 1, ijmdl
468  if (lats_mdl(ij) >= 0.0) then
469  if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
470  snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = &
471  max(snow_dep_mdl_tmp(ij), min_snow_depth)
472  endif
473  snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
474  endif
475  enddo
476  deallocate (snow_cvr_mdl_1d)
477  elseif (use_nh_afwa .and. use_nesdis) then ! set depth/cover on nesdis ims/afwa blend
478  print*,"- BLEND NESDIS/IMS AND AFWA DATA IN NH."
479  do ij = 1, ijmdl
480  if (lats_mdl(ij) >= 0.0) then
481  if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
482  snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = &
483  max(snow_dep_mdl_tmp(ij), min_snow_depth)
484  endif
485  snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
486  endif
487  enddo
488  deallocate (snow_cvr_mdl_1d)
489  deallocate (snow_dep_mdl_tmp)
490  elseif (use_global_afwa) then ! set depth/cover on afwa only
491  print*,"- SET DEPTH/COVER FROM AFWA DATA IN NH."
492  do ij = 1, ijmdl
493  if (lats_mdl(ij) >= 0.0) then
494  if (snow_dep_mdl_tmp(ij) > 0.0) then
495  snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = 100.0
496  snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij)
497  endif
498  endif
499  enddo
500  elseif (use_nh_afwa) then ! set depth/cover on afwa only
501  print*,"- SET DEPTH/COVER FROM AFWA DATA IN NH."
502  do ij = 1, ijmdl
503  if (lats_mdl(ij) >= 0.0) then
504  if (snow_dep_mdl_tmp(ij) > 0.0) then
505  snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = 100.0
506  snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij)
507  endif
508  endif
509  enddo
510  deallocate (snow_dep_mdl_tmp)
511  elseif (use_nesdis) then ! set depth/cover on nesdis/ims only
512  print*,"- SET DEPTH/COVER FROM NESDIS/IMS DATA IN NH."
513  do ij = 1, ijmdl
514  if (lats_mdl(ij) >= 0.0) then
515  if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
516  snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = min_snow_depth
517  endif
518  snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
519  endif
520  enddo
521  deallocate (snow_cvr_mdl_1d)
522  end if
523 
524 !----------------------------------------------------------------------
525 ! if selected, interpolate autosnow data to model grid.
526 !----------------------------------------------------------------------
527 
528  autosnow : if (use_autosnow) then
529 
530  ipopt = 0
531  if (autosnow_res < (0.5*resol_mdl)) then
532  print*,"- INTERPOLATE AUTOSNOW DATA TO MODEL GRID USING BUDGET METHOD."
533  ipopt(1)=2 ! break model grid cell into 25 points.
534  ipopt(2:4)=1 ! 25 points are weighted equally.
535  ipopt(5)=10 ! 10% coverage of valid data in box
536  ipopt(20) = nint(100.0 / autosnow_res) + 1 ! search box width of 100 km.
537  kgds_mdl_tmp = kgds_mdl
538  kgds_mdl_tmp(1) = kgds_mdl_tmp(1) - 255 ! subset of grid
539  int_opt = 3
540  no = ijmdl
541  else
542  print*,"- INTERPOLATE AUTOSNOW DATA TO MODEL GRID USING NEIGHBOR METHOD."
543  ipopt(1) = nint(100.0 / autosnow_res) + 1 ! search box width of 100 km.
544  kgds_mdl_tmp = kgds_mdl
545  kgds_mdl_tmp(1) = -1 ! for subsection of model grid.
546  int_opt = 2
547  no = ijmdl ! an input when kgds(1) < 0 (subset of grid)
548  end if
549 
550  allocate (snow_cvr_mdl_1d(ijmdl))
551  snow_cvr_mdl_1d = 0.0
552 
553  allocate (bitmap_mdl(ijmdl))
554  bitmap_mdl=.false. ! if interpolation routine can't find data
555  ! at a point, this flag is false.
556 
557  call ipolates(int_opt, ipopt, kgds_autosnow, kgds_mdl_tmp, &
558  (iautosnow*jautosnow), ijmdl, &
559  1, 1, bitmap_autosnow, snow_cvr_autosnow, &
560  no, lats_mdl, lons_mdl, ibo, bitmap_mdl, &
561  snow_cvr_mdl_1d, iret)
562 
563  deallocate (snow_cvr_autosnow, bitmap_autosnow)
564 
565  if (iret /= 0) then
566  print*,"- FATAL ERROR IN INTERPOLATION ROUTINE. IRET IS: ", iret
567  call w3tage('SNOW2MDL')
568  call errexit(55)
569  endif
570 
571 !----------------------------------------------------------------------
572 ! if interpolation fails to find autosnow data, set the cover
573 ! at the model point to a nomimal value.
574 !----------------------------------------------------------------------
575 
576  do ij = 1, ijmdl
577  if (lats_mdl(ij) < 0.0) then
578  if (.not. bitmap_mdl(ij)) then
579  if (abs(lats_mdl(ij)) <= lat_threshold) then
580  snow_cvr_mdl_1d(ij) = 0.0
581  else
582  snow_cvr_mdl_1d(ij) = 100.0
583  end if
584  end if
585  end if
586  enddo
587 
588  deallocate (bitmap_mdl)
589 
590  endif autosnow
591 
592 !----------------------------------------------------------------------
593 ! now interpolate sh afwa snow depth data.
594 !----------------------------------------------------------------------
595 
596  sh_afwa : if (use_sh_afwa) then
597 
598 !----------------------------------------------------------------------
599 ! determine interpolation method based on the resolution of
600 ! afwa data and the model grid.
601 !----------------------------------------------------------------------
602 
603  ipopt = 0
604  if (afwa_res < (0.5*resol_mdl)) then
605  print*,"- INTERPOLATE SH AFWA DATA TO MODEL GRID USING BUDGET METHOD."
606  ipopt(1)=-1 ! break model grid cell into 25 points.
607  ipopt(2)=-1 ! 25 points are weighted equally.
608  ipopt(20) = nint(100.0 / afwa_res) + 1 ! search box width of 100 km.
609  kgds_mdl_tmp = kgds_mdl
610  kgds_mdl_tmp(1) = kgds_mdl_tmp(1) - 255 ! subset of grid
611  no = ijmdl
612  int_opt = 3
613  else
614  print*,"- INTERPOLATE SH AFWA DATA TO MODEL GRID USING NEIGHBOR METHOD."
615  ipopt(1) = nint(100.0 / afwa_res) + 1 ! search box width of 100 km.
616  kgds_mdl_tmp = kgds_mdl
617  kgds_mdl_tmp(1) = -1 ! for subsection of model grid.
618  int_opt = 2
619  no = ijmdl ! an input when kgds(1) < 0 (subset of grid)
620  end if
621 
622  allocate (snow_dep_mdl_tmp(ijmdl))
623  snow_dep_mdl_tmp = 0.0
624 
625  allocate (bitmap_mdl(ijmdl))
626  bitmap_mdl = .false.
627 
628  call ipolates(int_opt, ipopt, kgds_afwa_sh, kgds_mdl_tmp, &
629  (iafwa*jafwa), ijmdl, &
630  1, 1, bitmap_afwa_sh, snow_dep_afwa_sh, &
631  no, lats_mdl, lons_mdl, ibo, bitmap_mdl, &
632  snow_dep_mdl_tmp, iret)
633 
634  if (iret /= 0) then
635  print*,"- FATAL ERROR IN INTERPOLATION ROUTINE. IRET IS: ", iret
636  call w3tage('SNOW2MDL')
637  call errexit(55)
638  endif
639 
640  deallocate (bitmap_afwa_sh, snow_dep_afwa_sh)
641 
642 !----------------------------------------------------------------------
643 ! if interpolation does not find afwa data, set model point to
644 ! a nominal value.
645 !----------------------------------------------------------------------
646 
647  do ij = 1, ijmdl
648  if (lats_mdl(ij) < 0.) then
649  if (.not. bitmap_mdl(ij)) then
650  if (abs(lats_mdl(ij)) >= lat_threshold) then
651  snow_dep_mdl_tmp(ij) = min_snow_depth
652  else
653  snow_dep_mdl_tmp(ij) = 0.0
654  endif
655  endif
656  endif
657  enddo
658 
659  deallocate(bitmap_mdl)
660 
661  endif sh_afwa
662 
663 !----------------------------------------------------------------------
664 ! if sh data selected, use it to determine the cover and depth
665 ! on the model grid.
666 !----------------------------------------------------------------------
667 
668  if ((use_sh_afwa .or. use_global_afwa) .and. use_autosnow) then ! set depth/cover on autosnow/afwa blend
669  print*,"- BLEND AUTOSNOW AND AFWA DATA IN SH."
670  do ij = 1, ijmdl
671  if (lats_mdl(ij) < 0.0) then
672  if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
673  snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = &
674  max(snow_dep_mdl_tmp(ij), min_snow_depth)
675  endif
676  snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
677  endif
678  enddo
679  deallocate (snow_cvr_mdl_1d)
680  deallocate (snow_dep_mdl_tmp)
681  elseif (use_sh_afwa .or. use_global_afwa) then ! set depth/cover on afwa only
682  print*,"- SET DEPTH/COVER FROM AFWA DATA IN SH."
683  do ij = 1, ijmdl
684  if (lats_mdl(ij) < 0.0) then
685  if (snow_dep_mdl_tmp(ij) > 0.0) then
686  snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = 100.0
687  snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij)
688  endif
689  endif
690  enddo
691  deallocate (snow_dep_mdl_tmp)
692  elseif (use_autosnow) then ! set depth/cover on autosnow only
693  print*,"- SET DEPTH/COVER FROM AUTOSNOW IN SH."
694  do ij = 1, ijmdl
695  if (lats_mdl(ij) < 0.0) then
696  if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
697  snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = min_snow_depth
698  endif
699  snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
700  endif
701  enddo
702  deallocate (snow_cvr_mdl_1d)
703  end if
704 
705 !----------------------------------------------------------------------
706 ! if a global model grid, and if running on thinned grid, then
707 ! take a linear weighting of full points located within the thin points.
708 ! "4" is grid indicator for a gaussian grid.
709 !----------------------------------------------------------------------
710 
711  if (kgds_mdl(1) == 4 .and. thinned) then
712 
713  ijmdl2 = sum(lonsperlat_mdl) * 2
714  allocate (snow_cvr_mdl_1d(ijmdl2))
715  allocate (lsmask_1d(ijmdl2))
716  allocate (snow_dep_mdl_tmp(ijmdl2))
717 
718  lsmask_1d = 0.0
719  snow_cvr_mdl_1d = 0.0
720  snow_dep_mdl_tmp = 0.0
721 
722  ij = 0
723  do j = 1, jmdl
724  jj = j
725  if (jj > jmdl/2) jj = jmdl - j + 1
726  r = float(imdl) / float(lonsperlat_mdl(jj))
727  do i = 1, lonsperlat_mdl(jj)
728  ij = ij + 1
729  x1 = (i-1)*r
730  imid = nint(x1+1.0)
731  lsmask_1d(ij) = lsmask_mdl_sav(imid,j)
732  if (lsmask_mdl_sav(imid,j) == 0.0) cycle
733  gridis=x1+1.0-r/2.
734  istart=nint(gridis)
735  gridie=x1+1.0+r/2.
736  iend=nint(gridie)
737  sumc = 0.0 ! cover
738  sumd = 0.0 ! depth
739  do ii = istart, iend
740  if (ii == istart) then
741  fraction = 0.5 - (gridis - float(istart))
742  elseif (ii == iend) then
743  fraction = 0.5 + (gridie - float(iend))
744  else
745  fraction = 1.0
746  endif
747  if (fraction < 0.0001) cycle
748  iii = ii
749  if (iii < 1) iii = imdl + iii
750  sumc = sumc + fraction * snow_cvr_mdl(iii,j)
751  sumd = sumd + fraction * snow_dep_mdl(iii,j)
752  enddo
753  snow_cvr_mdl_1d(ij) = sumc / r
754  snow_dep_mdl_tmp(ij) = 0.0
755  if (snow_cvr_mdl_1d(ij) > snow_cvr_threshold) then
756  snow_dep_mdl_tmp(ij) = max(sumd / r,min_snow_depth)
757  end if
758  enddo
759  enddo
760 
761  deallocate (lsmask_mdl_sav)
762 
763 !----------------------------------------------------------------------
764 ! now place thinned points into 2-d array for output.
765 !----------------------------------------------------------------------
766 
767  allocate (idum(imdl,jmdl))
768  idum = 0
769  call uninterpred(1, idum, lsmask_1d, lsmask_mdl, imdl, jmdl, ijmdl2, lonsperlat_mdl)
770  call uninterpred(1, idum, snow_cvr_mdl_1d, snow_cvr_mdl, imdl, jmdl, ijmdl2, lonsperlat_mdl)
771  deallocate(snow_cvr_mdl_1d)
772  call uninterpred(1, idum, snow_dep_mdl_tmp, snow_dep_mdl, imdl, jmdl, ijmdl2, lonsperlat_mdl)
773  deallocate(snow_dep_mdl_tmp)
774  deallocate(idum)
775 
776  end if
777 
778 !----------------------------------------------------------------------
779 ! grib the interpolated data.
780 !----------------------------------------------------------------------
781 
782  if (output_grib2) then
783  print*,"- OUTPUT SNOW ANALYSIS DATA IN GRIB2 FORMAT"
784  call write_grib2
785  else
786  print*,"- OUTPUT SNOW ANALYSIS DATA IN GRIB1 FORMAT"
787  call write_grib1
788  endif
789 
790  deallocate (snow_cvr_mdl)
791  deallocate (snow_dep_mdl)
792 
793  return
794 
795  end subroutine interp
796 
810  subroutine write_grib2
811 
812  use grib_mod
813 
814  implicit none
815 
816  character(len=1), allocatable :: cgrib(:)
817 
818  integer, parameter :: numcoord = 0
819 
820  integer :: coordlist(numcoord)
821  integer :: lugb, lcgrib, iret
822  integer :: igds(5)
823  integer :: listsec0(2)
824  integer :: listsec1(13)
825  integer :: ideflist, idefnum, ipdsnum, idrsnum
826  integer :: igdstmplen, ipdstmplen, idrstmplen
827  integer :: ipdstmpl(15)
828  integer, allocatable :: igdstmpl(:), idrstmpl(:)
829  integer :: ngrdpts, ibmap, lengrib
830 
831  logical*1, allocatable :: bmap(:), bmap2d(:,:)
832 
833  real, allocatable :: fld(:)
834 
835 !----------------------------------------------------------------------
836 ! Setup variables and arrays required by grib2 library.
837 !----------------------------------------------------------------------
838 
839  call grib2_check(kgds_mdl, igdstmplen)
840 
841  allocate(igdstmpl(igdstmplen))
842 
843  call init_grib2(grib_century,grib_year, grib_month, grib_day, grib_hour, &
844  kgds_mdl, lat11, latlast, lon11, lonlast, &
845  listsec0, listsec1, igds, ipdstmpl, ipdsnum, igdstmpl, &
846  igdstmplen, idefnum, ideflist, ngrdpts)
847 
848  lcgrib = imdl*jmdl*4
849  allocate(cgrib(lcgrib)) ! this variable holds the grib2 message
850 
851  iret=0
852 
853 !----------------------------------------------------------------------
854 ! Create sections 0 and 1. There is no section 2, local use section.
855 !----------------------------------------------------------------------
856 
857  print*,"- CREATE SECTIONS 0 AND 1"
858  call gribcreate(cgrib,lcgrib,listsec0,listsec1,iret)
859  if (iret /= 0) goto 900
860 
861 !----------------------------------------------------------------------
862 ! Create section 3, the grid description section.
863 !----------------------------------------------------------------------
864 
865  print*,"- CREATE SECTION 3"
866  call addgrid(cgrib,lcgrib,igds,igdstmpl,igdstmplen, &
867  ideflist,idefnum,iret)
868  if (iret /= 0) goto 900
869 
870 !----------------------------------------------------------------------
871 ! Setup arrays for section 5, the data representation section.
872 !----------------------------------------------------------------------
873 
874  idrsnum = 0 ! section 5, use simple packing
875  idrstmplen = 5
876  allocate (idrstmpl(idrstmplen))
877  idrstmpl = 0
878  idrstmpl(3) = 0 ! decimal scaling factor
879 
880  allocate(fld(ngrdpts))
881  fld = reshape(snow_cvr_mdl, (/imdl*jmdl/) )
882 
883  ibmap = 0 ! bitmap applies
884  allocate(bmap2d(imdl,jmdl))
885  bmap2d=.true.
886  where (lsmask_mdl < 0.5) bmap2d=.false.
887  allocate(bmap(ngrdpts))
888  bmap = reshape(bmap2d, (/imdl*jmdl/) )
889  deallocate (bmap2d)
890 
891  coordlist=0 ! not used
892 
893 !----------------------------------------------------------------------
894 ! Create section 4 (product definition section) and 5 (data
895 ! representation section) for cover.
896 !----------------------------------------------------------------------
897 
898  ipdstmpl(1) = 1 ! section 4, oct 10; parameter category, moisture
899  ipdstmpl(2) = 42 ! section 4, oct 11; parameter, snow cover
900 
901  print*,"- CREATE SECTIONS 4 AND 5 FOR SNOW COVER"
902  call addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,ipdstmplen, &
903  coordlist,numcoord,idrsnum,idrstmpl, &
904  idrstmplen,fld,ngrdpts,ibmap,bmap,iret)
905  if (iret /= 0) goto 900
906 
907 !----------------------------------------------------------------------
908 ! for regional model, if afwa data not used, don't output a depth
909 ! record. this tells the sfcupdate code to update the first
910 ! guess snow with the new snow cover analysis.
911 !----------------------------------------------------------------------
912 
913  if (kgds_mdl(1) /= 4) then
914  if (.not. use_global_afwa .and. .not. use_nh_afwa .and. .not. use_sh_afwa) goto 88
915  endif
916 
917 !----------------------------------------------------------------------
918 ! Create section 4 (product definition section) and 5 (data
919 ! representation section) for depth.
920 !----------------------------------------------------------------------
921 
922  fld= reshape(snow_dep_mdl, (/imdl*jmdl/) )
923 
924  ipdstmpl(1) = 1 ! section 4, oct 10; parameter category, moisture
925  ipdstmpl(2) = 11 ! section 4, oct 11; parameter, snow depth
926 
927  idrstmpl = 0 ! section 5
928  idrstmpl(3) = 4 ! section 5, decimal scaling factor
929 
930  print*,"- CREATE SECTIONS 4 AND 5 FOR SNOW DEPTH"
931  call addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,ipdstmplen, &
932  coordlist,numcoord,idrsnum,idrstmpl, &
933  idrstmplen,fld,ngrdpts,ibmap,bmap,iret)
934  if (iret /= 0) goto 900
935 
936  88 continue
937 
938 !----------------------------------------------------------------------
939 ! Create section 8 - end section.
940 !----------------------------------------------------------------------
941 
942  call gribend(cgrib,lcgrib,lengrib,iret)
943  if (iret /= 0) goto 900
944 
945 !----------------------------------------------------------------------
946 ! Now output grib message to file.
947 !----------------------------------------------------------------------
948 
949  lugb=53
950  print*,"- OPEN OUTPUT GRIB FILE ", trim(model_snow_file)
951  call baopenw(lugb, model_snow_file, iret)
952 
953  if (iret /= 0) then
954  print*,'- FATAL ERROR: BAD OPEN OF OUTPUT GRIB FILE. IRET IS ', iret
955  call w3tage('SNOW2MDL')
956  call errexit(49)
957  end if
958 
959  print*,'- WRITE OUTPUT GRIB FILE.'
960  call wryte(lugb, lengrib, cgrib)
961 
962  call baclose(lugb, iret)
963 
964  deallocate(fld, bmap, idrstmpl, igdstmpl, cgrib)
965 
966  return
967 
968  900 continue
969  print*,'- FATAL ERROR CREATING GRIB2 MESSAGE. IRET IS ', iret
970  call w3tage('SNOW2MDL')
971  call errexit(48)
972 
973  end subroutine write_grib2
974 
990  subroutine write_grib1
991 
992  implicit none
993 
994  integer :: iret
995  integer, parameter :: lugb = 64 ! unit number of output grib file
996  integer :: kpds(200)
997 
998  logical*1 :: lbms(imdl,jmdl)
999 
1000 !----------------------------------------------------------------------
1001 ! set up pds section. don't need to set the gds section.
1002 ! since the model grid is not changing, use the kgds array
1003 ! already determined in module model_grid.
1004 !----------------------------------------------------------------------
1005 
1006  kpds = 0
1007 
1008  kpds(1) = 7 ! center id
1009  kpds(2) = 25 ! process id number. this determined from the
1010  ! input data as we are simply interpolating
1011  ! that data to a different grid. should
1012  ! i request a process id for my codes?
1013  kpds(3) = grid_id_mdl ! grid specified in gds
1014  kpds(4) = 192 ! include gds and a bit map section
1015  kpds(5) = 238 ! parameter number for snow cover
1016  kpds(6) = 1 ! level - ground or water surface
1017  kpds(7) = 0 ! height pressure of level
1018  kpds(8) = grib_year ! year of century the time info is determined
1019  kpds(9) = grib_month ! month by operational requirements
1020  kpds(10) = grib_day ! day
1021  kpds(11) = grib_hour ! hour
1022  kpds(12) = 0 ! minute
1023  kpds(13) = 1 ! fcst time unit - hour
1024  kpds(14) = 0 ! period of time, p1. set to '0' for analysis
1025  kpds(15) = 0 ! number of time units, p2.
1026  kpds(16) = 1 ! initialized analysis product
1027  kpds(17) = 0 ! number in average
1028  kpds(18) = 1 ! grib edition 1
1029  kpds(19) = 3 ! parameter table version number
1030  kpds(20) = 0 ! number missing from avg/accum
1031  kpds(21) = grib_century ! century - set as in the input file
1032  kpds(22) = 0 ! decimal scale factor
1033  kpds(23) = 4 ! subcenter - emc
1034  kpds(24) = 0 ! reserved
1035  kpds(25) = 0 ! reserved
1036 
1037  lbms = .false. ! set bitmap section
1038 
1039  where(lsmask_mdl > 0.0) lbms = .true.
1040 
1041  print*,"- OPEN OUTPUT GRIB FILE ", trim(model_snow_file)
1042  call baopenw(lugb, model_snow_file, iret)
1043 
1044  if (iret /= 0) then
1045  print*,'- FATAL ERROR OPENING OUTPUT GRIB FILE. IRET IS ', iret
1046  call w3tage('SNOW2MDL')
1047  call errexit(59)
1048  end if
1049 
1050  print*,"- WRITE OUTPUT GRIB FILE ", trim(model_snow_file)
1051  call putgb(lugb, (imdl*jmdl), kpds, kgds_mdl, lbms, &
1052  snow_cvr_mdl, iret)
1053 
1054  if (iret /= 0) then
1055  print*,'- FATAL ERROR WRITING OUTPUT GRIB FILE. IRET IS ', iret
1056  call w3tage('SNOW2MDL')
1057  call errexit(58)
1058  end if
1059 
1060 ! for regional model, if afwa data not used, don't output a depth
1061 ! record. this tells the sfcupdate code to update the first
1062 ! guess snow with the new snow cover analysis.
1063 
1064  if (kgds_mdl(1) /= 4) then
1065  if (.not. use_global_afwa .and. .not. use_nh_afwa .and. .not. use_sh_afwa) goto 88
1066  endif
1067 
1068  kpds(5) = 66 ! parameter number for snow depth
1069  kpds(22) = 4 ! scaling factor. to nearest mm
1070 
1071  call putgb(lugb, (imdl*jmdl), kpds, kgds_mdl, lbms, &
1072  snow_dep_mdl, iret)
1073 
1074  if (iret /= 0) then
1075  print*,'- FATAL ERROR WRITING OUTPUT GRIB FILE. IRET IS ', iret
1076  call w3tage('SNOW2MDL')
1077  call errexit(57)
1078  end if
1079 
1080  88 call baclose(lugb, iret)
1081 
1082  return
1083 
1084  end subroutine write_grib1
1085 
1102  subroutine uninterpred(iord,kmsk,fi,f,lonl,latd,len,lonsperlat)
1103 
1104  implicit none
1105 
1106  integer, intent(in) :: len
1107  integer, intent(in) :: iord
1108  integer, intent(in) :: lonl
1109  integer, intent(in) :: latd
1110  integer, intent(in) :: lonsperlat(latd/2)
1111  integer, intent(in) :: kmsk(lonl*latd)
1112  integer :: j,lons,jj,latd2,ii,i
1113 
1114  real, intent(in) :: fi(len)
1115  real, intent(out) :: f(lonl,latd)
1116 
1117  latd2 = latd / 2
1118  ii = 1
1119 
1120  do j=1,latd
1121 
1122  jj = j
1123  if (j .gt. latd2) jj = latd - j + 1
1124  lons=lonsperlat(jj)
1125 
1126  if(lons.ne.lonl) then
1127  call intlon(iord,1,1,lons,lonl,kmsk(ii),fi(ii),f(1,j))
1128  else
1129  do i=1,lonl
1130  f(i,j) = fi(ii+i-1)
1131  enddo
1132  endif
1133 
1134  ii = ii + lons
1135 
1136  enddo
1137 
1138  end subroutine uninterpred
1139 
1153  subroutine intlon(iord,imon,imsk,m1,m2,k1,f1,f2)
1154 
1155  implicit none
1156 
1157  integer,intent(in) :: iord,imon,imsk,m1,m2
1158  integer,intent(in) :: k1(m1)
1159  integer :: i2,in,il,ir
1160 
1161  real,intent(in) :: f1(m1)
1162  real,intent(out) :: f2(m2)
1163  real :: r,x1
1164 
1165  r=real(m1)/real(m2)
1166  do i2=1,m2
1167  x1=(i2-1)*r
1168  il=int(x1)+1
1169  ir=mod(il,m1)+1
1170  if(iord.eq.2.and.(imsk.eq.0.or.k1(il).eq.k1(ir))) then
1171  f2(i2)=f1(il)*(il-x1)+f1(ir)*(x1-il+1)
1172  else
1173  in=mod(nint(x1),m1)+1
1174  f2(i2)=f1(in)
1175  endif
1176  enddo
1177 
1178  end subroutine intlon
1179 
1180  end module snow2mdl
init_grib2
subroutine init_grib2(century, year, month, day, hour, kgds, lat11, latlast, lon11, lonlast, listsec0, listsec1, igds, ipdstmpl, ipdsnum, igdstmpl, igdstmplen, idefnum, ideflist, ngrdpts)
Initialize grib2 arrays required by the ncep g2 library according to grib1 gds information.
Definition: grib_utils.F90:397
grib2_check
subroutine grib2_check(kgds, igdstmplen)
Determine length of grib2 gds template array, which is a function of the map projection.
Definition: grib_utils.F90:350
model_grid
Read in data defining the model grid.
Definition: model_grid.F90:19
snow2mdl::interp
subroutine, public interp
Interpolate snow data to model grid.
Definition: snow2mdl.F90:159
snow2mdl::uninterpred
subroutine uninterpred(iord, kmsk, fi, f, lonl, latd, len, lonsperlat)
Fills out full grid using thinned grid data.
Definition: snow2mdl.F90:1102
snow2mdl::write_grib1
subroutine write_grib1
Write grib1 snow cover and depth on the model grid.
Definition: snow2mdl.F90:990
program_setup
This module reads in data from the program&#39;s configuration namelist.
Definition: program_setup.F90:15
snowdat
Read and qc afwa, nesdis/ims and autosnow snow data.
Definition: snowdat.F90:26
snow2mdl::intlon
subroutine intlon(iord, imon, imsk, m1, m2, k1, f1, f2)
Convert data from the thinned (or reduced) to the full grid along a single row.
Definition: snow2mdl.F90:1153
snow2mdl
Interpolate snow data to model grid and grib the result.
Definition: snow2mdl.F90:20
snow2mdl::write_grib2
subroutine write_grib2
Write grib2 snow cover and depth on the model grid.
Definition: snow2mdl.F90:810