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