emcsfc_snow2mdl 1.14.0
Loading...
Searching...
No Matches
snow2mdl.F90
Go to the documentation of this file.
1
4
20 module snow2mdl
21
22 use program_setup, only : lat_threshold, &
26 grib_day, &
28 grib_hour, &
29 grib_month, &
30 grib_year, &
32
33 use model_grid, only : resol_mdl, &
34 imdl, &
35 jmdl, &
36 ijmdl, &
37 ipts_mdl, &
38 jpts_mdl, &
39 lsmask_mdl, &
41 lats_mdl, &
42 lat11, latlast, &
43 lon11, lonlast, &
44 lons_mdl, &
45 kgds_mdl, &
47 thinned, &
49
50 use snowdat, only : nesdis_res, &
51 afwa_res, &
53 inesdis, &
54 jnesdis, &
59 iafwa, &
60 jafwa, &
67 iautosnow, &
68 jautosnow, &
74 use_nesdis, &
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, &
360 no, lats_mdl, lons_mdl, ibo, bitmap_mdl, &
361 snow_dep_mdl_tmp, iret)
362
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, &
432 no, lats_mdl, lons_mdl, ibo, bitmap_mdl, &
433 snow_dep_mdl_tmp, iret)
434
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
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, &
572 no, lats_mdl, lons_mdl, ibo, bitmap_mdl, &
573 snow_cvr_mdl_1d, iret)
574
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, &
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
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
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
subroutine grib2_check(kgds, igdstmplen)
Determine length of grib2 gds template array, which is a function of the map projection.
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.
Read in data defining the model grid.
integer, dimension(:), allocatable jpts_mdl
j index of point on full grid
real lonlast
Corner point longitude (imdl,jmdl) of model grid.
integer jmdl
j-dimension of model grid
integer, dimension(:), allocatable lonsperlat_mdl
Number of longitudes (i-points) for each latitude (row).
integer imdl
i-dimension of model grid
real latlast
Corner point latitude (imdl,jmdl) of model grid.
real, dimension(:,:), allocatable lsmask_mdl
land mask of model grid (0 - non land, 1-land) for global grids run thinned, will contain a modified ...
real, dimension(:), allocatable lats_mdl
Latitudes of model grid points.
logical thinned
When true, global grids will run thinned (number of i points decrease toward pole)
real, dimension(:,:), allocatable lsmask_mdl_sav
saved copy of land mask of model grid (0 - non land, 1-land) only used for global thinned grids.
integer ijmdl
total number of model land points
integer, dimension(200) kgds_mdl
holds grib gds info of model grid
real resol_mdl
approximate model resolution in km.
real lat11
Corner point latitude (1,1) of model grid.
integer, dimension(:), allocatable ipts_mdl
i index of point on full grid
integer grid_id_mdl
grib id of model grid, 4-gaussian, 203-egrid
real lon11
Corner point longitude (1,1) of model grid.
real, dimension(:), allocatable lons_mdl
longitudes of model grid points
This module reads in data from the program's configuration namelist.
integer, public grib_day
date of the final merged snow product that will be placed in grib header.
real, public lat_threshold
equatorward of this latitude, model points with undefined cover or depth (because the interpolation r...
character *200, public model_snow_file
path/name nesdis/ims snow cover
integer, public grib_year
date of the final merged snow product that will be placed in grib header.
real, public min_snow_depth
minimum snow depth in meters at model points with coverage exceeding threshold.
real, public snow_cvr_threshold
if percent coverage according to nesdis/ims or autosnow exceeds this value, then non-zero snow depth ...
logical, public output_grib2
when true, output model snow analysis is grib 2.
integer, public grib_month
date of the final merged snow product that will be placed in grib header.
integer, public grib_century
date of the final merged snow product that will be placed in grib header.
integer, public grib_hour
date of the final merged snow product that will be placed in grib header.
Interpolate snow data to model grid and grib the result.
Definition snow2mdl.F90:20
subroutine, public interp
Interpolate snow data to model grid.
Definition snow2mdl.F90:158
real, dimension(:,:), allocatable snow_cvr_mdl
snow cover on model grid in percent
Definition snow2mdl.F90:85
real, dimension(:,:), allocatable snow_dep_mdl
snow depth on model grid in meters
Definition snow2mdl.F90:86
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.
subroutine write_grib1
Write grib1 snow cover and depth on the model grid.
subroutine write_grib2
Write grib2 snow cover and depth on the model grid.
Definition snow2mdl.F90:824
subroutine uninterpred(iord, kmsk, fi, f, lonl, latd, len, lonsperlat)
Fills out full grid using thinned grid data.
Read and qc afwa, nesdis/ims and autosnow snow data.
Definition snowdat.F90:26
integer iafwa
i-dimension of afwa grid
Definition snowdat.F90:40
real autosnow_res
Resolution of autosnow in km.
Definition snowdat.F90:83
integer jnesdis
j-dimension of nesdis grid
Definition snowdat.F90:45
logical *1, dimension(:,:), allocatable bitmap_afwa_sh
The southern hemisphere afwa data grib bitmap.
Definition snowdat.F90:72
logical *1, dimension(:,:), allocatable bitmap_afwa_nh
The northern hemisphere afwa data grib bitmap.
Definition snowdat.F90:70
integer mesh_nesdis
nesdis/ims data is 96th mesh (or bediant)
Definition snowdat.F90:58
logical bad_afwa_sh
When true, the southern hemisphere afwa data failed its quality control check.
Definition snowdat.F90:62
logical use_nesdis
True if nesdis/ims data to be used.
Definition snowdat.F90:81
logical *1, dimension(:,:), allocatable bitmap_afwa_global
The global afwa data grib bitmap.
Definition snowdat.F90:68
real, dimension(:,:), allocatable snow_dep_afwa_nh
Northern hemisphere afwa snow depth.
Definition snowdat.F90:89
logical *1, dimension(:,:), allocatable bitmap_nesdis
nesdis data grib bitmap (false-non land, true-land).
Definition snowdat.F90:74
logical use_autosnow
True if autosnow data to be used.
Definition snowdat.F90:80
integer jafwa
j-dimension of afwa grid
Definition snowdat.F90:43
integer iautosnow
i-dimension of autosnow grid
Definition snowdat.F90:41
logical use_sh_afwa
True if southern hemisphere afwa data to be used.
Definition snowdat.F90:78
logical use_global_afwa
True if global hemisphere afwa data to be used.
Definition snowdat.F90:79
real afwa_res
Resolution of afwa data in km.
Definition snowdat.F90:84
real nesdis_res
Resolution of the nesdis data in km.
Definition snowdat.F90:85
integer, dimension(200) kgds_afwa_sh
grib1 grid description section for southern hemisphere 16th mesh afwa data.
Definition snowdat.F90:52
real, dimension(:,:), allocatable snow_dep_afwa_sh
Southern hemisphere afwa snow depth.
Definition snowdat.F90:90
integer *1, dimension(:,:), allocatable sea_ice_nesdis
nesdis/ims sea ice flag (0-open water, 1-ice)
Definition snowdat.F90:59
logical *1, dimension(:,:), allocatable bitmap_autosnow
autosnow data grib bitmap (false-non land, true-land).
Definition snowdat.F90:75
integer jautosnow
j-dimension of autosnow grid
Definition snowdat.F90:44
integer inesdis
i-dimension of nesdis grid
Definition snowdat.F90:42
real, dimension(:,:), allocatable snow_cvr_nesdis
nesdis/ims snow cover flag (0-no, 100-yes)
Definition snowdat.F90:86
integer, dimension(200) kgds_autosnow
autosnow grid description section (grib section 2)
Definition snowdat.F90:56
integer, dimension(200) kgds_afwa_nh
grib1 grid description section for northern hemisphere 16th mesh afwa data.
Definition snowdat.F90:48
integer, dimension(200) kgds_nesdis
nesdis/ims grid description section (grib section 2)
Definition snowdat.F90:57
integer, dimension(200) kgds_afwa_global
grib1 grid description section for global afwa data.
Definition snowdat.F90:46
real, dimension(:,:), allocatable snow_cvr_autosnow
autosnow snow cover flag (0-no, 100-yes)
Definition snowdat.F90:87
real, dimension(:,:), allocatable snow_dep_afwa_global
The global afwa snow depth.
Definition snowdat.F90:88
logical use_nh_afwa
True if northern hemisphere afwa data to be used.
Definition snowdat.F90:77
logical bad_afwa_nh
When true, the northern hemisphere afwa data failed its quality control check.
Definition snowdat.F90:60