emcsfc_snow2mdl 1.14.0
Loading...
Searching...
No Matches
model_grid.F90
Go to the documentation of this file.
1
4
19 module model_grid
20
21 use program_setup, only : model_lsmask_file, &
22 model_lon_file, &
23 model_lat_file, &
24 gfs_lpl_file
25
26 integer :: grid_id_mdl
27 integer :: imdl
28 integer :: jmdl
29 integer :: ijmdl
30 integer, allocatable :: ipts_mdl(:)
31 integer, allocatable :: jpts_mdl(:)
32
33 integer :: kgds_mdl(200)
34 integer, allocatable :: lonsperlat_mdl (:)
37
38 logical :: thinned
40
41 real, allocatable :: lats_mdl(:)
42 real :: lat11
43 real :: latlast
44 real :: lon11
45 real :: lonlast
46 real, allocatable :: lons_mdl(:)
47 real, allocatable :: lsmask_mdl(:,:)
50 real, allocatable :: lsmask_mdl_sav(:,:)
52 real :: resol_mdl
53
54 contains
84 subroutine read_mdl_grid_info
85 use grib_mod ! grib 2 library
86
87 implicit none
88
89 character*256 :: fngrib
90
91 integer :: i, j, ij, jj
92 integer :: ii, iii, istart, iend, imid
93 integer :: iret
94 integer :: isgrib
95 integer, parameter :: iunit = 14 ! unit of input grib file
96 integer, parameter :: iunit2 = 34 ! unit of input lonsperlat file
97 integer :: jgds(200)
98 integer :: jpds(200)
99 integer :: jdisc, jgdtn, jpdtn, k
100 integer :: jids(200), jgdt(200), jpdt(200)
101 integer :: lskip
102 integer, parameter :: lugi = 0 ! unit of grib index file - not used
103 integer :: kgds(200)
104 integer :: kpds(200)
105 integer :: message_num
106 integer :: numbytes
107 integer :: numpts
108
109 logical*1, allocatable :: lbms(:)
110 logical :: unpack
111
112 real :: gridis, gridie, fraction, x1, r
113 real, allocatable :: lats_mdl_temp (:,:)
114 real, allocatable :: lons_mdl_temp (:,:)
115
116 type(gribfield) :: gfld
117
118 print*,"- READ MODEL GRID INFORMATION"
119
120!-----------------------------------------------------------------------
121! read latitudes on the model grid. first check if file is grib1 or 2.
122!-----------------------------------------------------------------------
123
124 fngrib = model_lat_file
125
126 call grib_check(fngrib, isgrib)
127
128 if (isgrib==0) then
129 print*,'- FATAL ERROR: MODEL LAT FILE MUST BE GRIB1 OR GRIB2 FORMAT'
130 call w3tage('SNOW2MDL')
131 call errexit(90)
132 end if
133
134 print*,"- OPEN MODEL LAT FILE ", trim(fngrib)
135 call baopenr (iunit, fngrib, iret)
136
137 if (iret /= 0) then
138 print*,'- FATAL ERROR: BAD OPEN, IRET IS ', iret
139 call w3tage('SNOW2MDL')
140 call errexit(80)
141 end if
142
143 if (isgrib==1) then ! grib 1 file
144
145!-----------------------------------------------------------------------
146! tell degribber to search for latitudes
147!-----------------------------------------------------------------------
148
149 lskip = -1 ! read beginning of file
150 jgds = -1
151 jpds = -1
152 jpds(5) = 176 ! latitude
153 kgds = -1
154 kpds = -1
155
156 print*,"- GET GRIB HEADER"
157 call getgbh(iunit, lugi, lskip, jpds, jgds, numbytes, &
158 numpts, message_num, kpds, kgds, iret)
159
160 if (iret /= 0) then
161 print*,'- FATAL ERROR: BAD READ OF GRIB HEADER. IRET IS ', iret
162 call w3tage('SNOW2MDL')
163 call errexit(81)
164 end if
165
166!-----------------------------------------------------------------------
167! save gds for gribbing the interpolated data later. also required
168! by ncep ipolates library.
169!-----------------------------------------------------------------------
170
171 kgds_mdl = kgds
172
173!-----------------------------------------------------------------------
174! get model grid specs from header. model resolution (km) is used
175! to determine the interpolation method.
176!-----------------------------------------------------------------------
177
178 grid_id_mdl = kpds(3) ! grib 1 grid id number. sect 1, oct 7
179
180 if (kgds(1) == 4) then ! gaussian grid
181 imdl = kgds(2) ! i-dimension of model grid
182 jmdl = kgds(3) ! j-dimension of model grid
183 resol_mdl = float(kgds(9)) / 1000.0 * 111.0
184 else if (kgds(1) == 203) then ! e-grid
185 imdl = kgds(2) ! i-dimension of model grid
186 jmdl = kgds(3) ! j-dimension of model grid
187 resol_mdl = sqrt( (float(kgds(9)) / 1000.0)**2 + &
188 (float(kgds(10)) / 1000.0)**2 )
189 resol_mdl = resol_mdl * 111.0
190 else if (kgds(1) == 205) then ! b-grid
191 imdl = kgds(2) ! i-dimension of model grid
192 jmdl = kgds(3) ! j-dimension of model grid
193 resol_mdl = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) &
194 * 0.5 * 111.0
195 else
196 print*,'- FATAL ERROR: UNRECOGNIZED MODEL GRID.'
197 call w3tage('SNOW2MDL')
198 call errexit(79)
199 end if
200
201 allocate(lats_mdl_temp(imdl,jmdl))
202 allocate(lbms(imdl*jmdl))
203
204 print*,"- DEGRIB DATA"
205 call getgb(iunit, lugi, (imdl*jmdl), lskip, jpds, jgds, &
206 numpts, message_num, kpds, kgds, lbms, lats_mdl_temp, iret)
207
208 if (iret /= 0) then
209 print*,'- FATAL ERROR: BAD DEGRIB OF FILE. IRET IS ',iret
210 call w3tage('SNOW2MDL')
211 call errexit(82)
212 end if
213
214 deallocate(lbms)
215
216 lat11 = lats_mdl_temp(1,1)
217 latlast = lats_mdl_temp(imdl,jmdl)
218
219 elseif (isgrib==2) then ! grib 2 file
220
221 j = 0 ! search at beginning of file
222 jdisc = 0 ! search for discipline; 0 - meteorological products
223 jpdtn = -1 ! search for any product definition template number
224 jgdtn = -1 ! search for any grid definition template number
225 jids = -9999 ! array of values in identification section, set to wildcard
226 jgdt = -9999 ! array of values in grid definition template 3.m
227 jpdt = -9999 ! array of values in product definition template 4.n
228 jpdt(1) = 191 ! search for parameter category - misc
229 jpdt(2) = 1 ! search for parameter number - latitude
230 unpack = .true. ! unpack data
231
232 call grib2_null(gfld)
233
234 print*,"- DEGRIB DATA"
235 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
236 unpack, k, gfld, iret)
237
238 if (iret /=0) then
239 print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret
240 call w3tage('SNOW2MDL')
241 call errexit(82)
242 endif
243
244!-----------------------------------------------------------------------
245! create the grib 1 gds array from the g2 gdt array. the grib 1
246! gds info is used by ipolates and for gribbing the final snow analysis.
247!-----------------------------------------------------------------------
248
249 call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds_mdl, &
250 imdl, jmdl, resol_mdl)
251
252 grid_id_mdl = 255 ! grib1 grid id number. n/a for grib2.
253 ! set to 'missing'.
254
255 allocate(lats_mdl_temp(imdl,jmdl))
256 lats_mdl_temp = reshape(gfld%fld , (/imdl,jmdl/) )
257
258 lat11 = lats_mdl_temp(1,1)
259 latlast = lats_mdl_temp(imdl,jmdl)
260
261 call grib2_free(gfld)
262
263 endif ! grib1 or grib2?
264
265 call baclose(iunit,iret)
266
267!-----------------------------------------------------------------------
268! read longitudes on the model grid.
269!-----------------------------------------------------------------------
270
271 fngrib = model_lon_file
272
273 call grib_check(fngrib, isgrib)
274
275 if (isgrib==0) then
276 print*,'- FATAL ERROR: MODEL LON FILE MUST BE GRIB1 OR GRIB2 FORMAT'
277 call w3tage('SNOW2MDL')
278 call errexit(91)
279 end if
280
281 print*,"- OPEN MODEL LON FILE ", trim(fngrib)
282 call baopenr (iunit, fngrib, iret)
283
284 if (iret /= 0) then
285 print*,"- FATAL ERROR: BAD OPEN. IRET IS ", iret
286 call w3tage('SNOW2MDL')
287 call errexit(83)
288 end if
289
290 if (isgrib==1) then ! grib 1 file
291
292 lskip = -1
293 kgds = -1
294 kpds = -1
295 jgds = -1
296 jpds = -1
297 jpds(5) = 177 ! longitude
298
299 allocate(lons_mdl_temp(imdl,jmdl))
300 allocate(lbms(imdl*jmdl))
301
302 print*,"- DEGRIB DATA"
303 call getgb(iunit, lugi, (imdl*jmdl), lskip, jpds, jgds, &
304 numpts, message_num, kpds, kgds, lbms, lons_mdl_temp, iret)
305
306 if (iret /= 0) then
307 print*,'- FATAL ERROR: BAD DEGRIB OF DATA. IRET IS ',iret
308 call w3tage('SNOW2MDL')
309 call errexit(84)
310 end if
311
312 deallocate(lbms)
313
314 lon11 = lons_mdl_temp(1,1)
315 lonlast = lons_mdl_temp(imdl,jmdl)
316
317 elseif (isgrib==2) then ! grib2
318
319 j = 0 ! search at beginning of file
320 jdisc = 0 ! search for discipline; 0 - meteorological products
321 jpdtn = -1 ! search for any product definition template number
322 jgdtn = -1 ! search for any grid definition template number
323 jids = -9999 ! array of values in identification section, set to wildcard
324 jgdt = -9999 ! array of values in grid definition template 3.m
325 jpdt = -9999 ! array of values in product definition template 4.n
326 jpdt(1) = 191 ! search for parameter category - misc
327 jpdt(2) = 2 ! search for parameter number - longitude
328 unpack = .true. ! unpack data
329
330 call grib2_null(gfld)
331
332 print*,"- DEGRIB DATA"
333 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
334 unpack, k, gfld, iret)
335
336 if (iret /=0) then
337 print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret
338 call w3tage('SNOW2MDL')
339 call errexit(84)
340 endif
341
342 allocate(lons_mdl_temp(imdl,jmdl))
343 lons_mdl_temp = reshape(gfld%fld , (/imdl,jmdl/) )
344
345 lon11 = lons_mdl_temp(1,1)
346 lonlast = lons_mdl_temp(imdl,jmdl)
347
348 call grib2_free(gfld)
349
350 endif ! grib1 or grib?
351
352 call baclose(iunit, iret)
353
354!-----------------------------------------------------------------------
355! read model land/sea mask.
356!-----------------------------------------------------------------------
357
358 fngrib = model_lsmask_file
359
360 call grib_check(fngrib, isgrib)
361
362 if (isgrib==0) then
363 print*,'- FATAL ERROR: MODEL LANDMASK FILE MUST BE GRIB1 OR GRIB2 FORMAT'
364 call w3tage('SNOW2MDL')
365 call errexit(92)
366 end if
367
368 print*,"- OPEN MODEL LANDMASK FILE ", trim(fngrib)
369 call baopenr (iunit, fngrib, iret)
370
371 if (iret /= 0) then
372 print*,'- FATAL ERROR: BAD OPEN OF FILE. IRET IS ', iret
373 call w3tage('SNOW2MDL')
374 call errexit(85)
375 end if
376
377 if (isgrib==1) then ! grib 1 file
378
379 lskip = -1
380 kgds = -1
381 kpds = -1
382 jpds = -1
383 jgds = -1
384 jpds(5) = 81 ! land-sea mask
385
386 allocate(lsmask_mdl(imdl,jmdl))
387 allocate(lbms(imdl*jmdl))
388
389 print*,"- DEGRIB DATA"
390 call getgb(iunit, lugi, (imdl*jmdl), lskip, jpds, jgds, &
391 numpts, message_num, kpds, kgds, lbms, lsmask_mdl, iret)
392
393 if (iret /= 0) then
394 print*,'- FATAL ERROR: BAD DEGRIB OF DATA. IRET IS ',iret
395 call w3tage('SNOW2MDL')
396 call errexit(86)
397 end if
398
399 deallocate (lbms)
400
401 elseif (isgrib==2) then ! grib2
402
403 j = 0 ! search at beginning of file
404 jdisc = 2 ! search for discipline; 2 - land-sfc products
405 jpdtn = -1 ! search for any product definition template number
406 jgdtn = -1 ! search for any grid definition template number
407 jids = -9999 ! array of values in identification section, set to wildcard
408 jgdt = -9999 ! array of values in grid definition template 3.m
409 jpdt = -9999 ! array of values in product definition template 4.n
410 jpdt(1) = 0 ! search for parameter category - veg_biomass
411 jpdt(2) = 0 ! search for parameter number - landcover
412 unpack = .true. ! unpack data
413
414 call grib2_null(gfld)
415
416 print*,"- DEGRIB DATA"
417 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
418 unpack, k, gfld, iret)
419
420 if (iret /=0) then
421 print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret
422 call w3tage('SNOW2MDL')
423 call errexit(86)
424 endif
425
426 allocate(lsmask_mdl(imdl,jmdl))
427 lsmask_mdl = reshape(gfld%fld , (/imdl,jmdl/) )
428
429 call grib2_free(gfld)
430
431 endif
432
433 call baclose(iunit,iret)
434
435!-----------------------------------------------------------------------
436! global model runs on a thinned grid (# grid points decreases
437! towards the poles). if thinned logical is set, and this is a
438! gaussian grid, modify the land/sea mask to account for the
439! fact that delta x increases toward the poles.
440!-----------------------------------------------------------------------
441
442 thinned=.false.
443 if (kgds(1) == 4 .and. (len_trim(gfs_lpl_file) > 0)) then
444
445 thinned=.true.
446
447 print*,"- RUNNING A THINNED GRID"
448
449 allocate (lonsperlat_mdl(jmdl/2))
450
451 print*,"- OPEN/READ GFS LONSPERLAT FILE: ",trim(gfs_lpl_file)
452 open (iunit2, file=trim(gfs_lpl_file), iostat=iret)
453 if (iret /= 0) then
454 print*,'- FATAL ERROR: BAD OPEN OF LONSPERLAT FILE. ABORT. IRET: ', iret
455 call w3tage('SNOW2MDL')
456 call errexit(76)
457 endif
458
459 read (iunit2,*,iostat=iret) numpts, lonsperlat_mdl
460 close(iunit2)
461 if (iret /= 0) then
462 print*,'- FATAL ERROR: BAD READ OF LONSPERLAT FILE. ABORT. IRET: ', iret
463 call w3tage('SNOW2MDL')
464 call errexit(76)
465 endif
466
467 if (numpts /= (jmdl/2)) then
468 print*,'- FATAL ERROR: WRONG DIMENSIION IN LONSPERLAT FILE. ABORT.'
469 call w3tage('SNOW2MDL')
470 call errexit(76)
471 endif
472
473 allocate (lsmask_mdl_sav(imdl,jmdl))
474 lsmask_mdl_sav = lsmask_mdl
475 lsmask_mdl = 0.0 ! this will identify land points to be processed by
476 ! the ipolates routines.
477
478!-----------------------------------------------------------------------
479! loop over every point on the thinned grid. calculate the start/end
480! bounds with respect to the full grid in the 'i' direction. if
481! the thinned point contains land, set all full grid points within
482! the bounds to be land. this modified mask will identify the
483! points to be processed by ipolates. after the call to ipolates,
484! the thinned points will be set to a linear weighting of the full points
485! located within the thinned point.
486!-----------------------------------------------------------------------
487
488 do j = 1, jmdl
489 jj = j
490 if (j > jmdl/2) jj = jmdl - j + 1
491 r = float(imdl)/ float(lonsperlat_mdl(jj))
492 do i = 1, lonsperlat_mdl(jj)
493 x1=float(i-1)*r
494 imid = nint(x1+1.0) ! for this thinned grid point, this is
495 ! the nearest 'i' index on the full grid.
496 if (lsmask_mdl_sav(imid,j) > 0.0) then
497 gridis = x1+1.0-r/2.
498 istart = nint(gridis)
499 gridie = x1+1.0+r/2.
500 iend = nint(gridie)
501 do ii = istart, iend
502 if (ii == istart) then
503 fraction = 0.5 - (gridis - float(istart))
504 if (fraction < 0.0001) cycle
505 endif
506 if (ii == iend) then
507 fraction = 0.5 + (gridie - float(iend))
508 if (fraction < 0.0001) cycle
509 endif
510 iii = ii
511 if (iii < 1) iii = imdl + iii
512 lsmask_mdl(iii,j) = lsmask_mdl_sav(imid,j)
513 enddo
514 endif
515 enddo
516 enddo
517
518 end if
519
520!-----------------------------------------------------------------------
521! program only worries about land points. save i/j coordinate
522! with respect to 2-d grid.
523!-----------------------------------------------------------------------
524
525 ij = 0
526
527 do j = 1, jmdl
528 do i = 1, imdl
529 if (lsmask_mdl(i,j) > 0.0) then
530 ij = ij+1
531 end if
532 enddo
533 enddo
534
535 ijmdl = ij
536
537 if (ijmdl == 0) then ! grid has only water points, dont run
538 print*,' '
539 print*,'- MODEL GRID ONLY HAS WATER POINTS, DONT CREATE SNOW FILE.'
540 print*,'- NORMAL TERMINATION.'
541 call w3tage('SNOW2MDL')
542 call errexit(0)
543 endif
544
545 allocate (lats_mdl(ijmdl))
546 allocate (lons_mdl(ijmdl))
547 allocate (ipts_mdl(ijmdl))
548 allocate (jpts_mdl(ijmdl))
549
550 ij = 0
551 do j = 1, jmdl
552 do i = 1, imdl
553 if (lsmask_mdl(i,j) > 0.0) then
554 ij = ij+1
555 lats_mdl(ij) = lats_mdl_temp(i,j)
556 lons_mdl(ij) = lons_mdl_temp(i,j)
557 ipts_mdl(ij) = i
558 jpts_mdl(ij) = j
559 end if
560 enddo
561 enddo
562
563 deallocate (lats_mdl_temp, lons_mdl_temp)
564
565 return
566 end subroutine read_mdl_grid_info
567
568
577 subroutine model_grid_cleanup
578
579 implicit none
580
581 if (allocated(lsmask_mdl)) deallocate(lsmask_mdl)
582 if (allocated(lats_mdl)) deallocate(lats_mdl)
583 if (allocated(lons_mdl)) deallocate(lons_mdl)
584 if (allocated(lonsperlat_mdl)) deallocate(lonsperlat_mdl)
585 if (allocated(ipts_mdl)) deallocate(ipts_mdl)
586 if (allocated(jpts_mdl)) deallocate(jpts_mdl)
587
588 return
589
590 end subroutine model_grid_cleanup
591
592 end module model_grid
grib2_free
subroutine grib2_free(gfld)
Deallocate the grib2 gribfield pointers.
Definition grib_utils.F90:636
gdt_to_gds
subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
Convert from the grib2 grid description template array used by the ncep grib2 library,...
Definition grib_utils.F90:140
grib2_null
subroutine grib2_null(gfld)
Nullify the grib2 gribfield pointers.
Definition grib_utils.F90:611
grib_check
subroutine grib_check(file_name, isgrib)
Determine whether file is grib or not.
Definition grib_utils.F90:24
model_grid
Read in data defining the model grid.
Definition model_grid.F90:19
model_grid::jpts_mdl
integer, dimension(:), allocatable jpts_mdl
j index of point on full grid
Definition model_grid.F90:31
model_grid::lonlast
real lonlast
Corner point longitude (imdl,jmdl) of model grid.
Definition model_grid.F90:45
model_grid::jmdl
integer jmdl
j-dimension of model grid
Definition model_grid.F90:28
model_grid::lonsperlat_mdl
integer, dimension(:), allocatable lonsperlat_mdl
Number of longitudes (i-points) for each latitude (row).
Definition model_grid.F90:34
model_grid::imdl
integer imdl
i-dimension of model grid
Definition model_grid.F90:27
model_grid::latlast
real latlast
Corner point latitude (imdl,jmdl) of model grid.
Definition model_grid.F90:43
model_grid::lsmask_mdl
real, dimension(:,:), allocatable lsmask_mdl
land mask of model grid (0 - non land, 1-land) for global grids run thinned, will contain a modified ...
Definition model_grid.F90:47
model_grid::lats_mdl
real, dimension(:), allocatable lats_mdl
Latitudes of model grid points.
Definition model_grid.F90:41
model_grid::thinned
logical thinned
When true, global grids will run thinned (number of i points decrease toward pole)
Definition model_grid.F90:38
model_grid::lsmask_mdl_sav
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.
Definition model_grid.F90:50
model_grid::ijmdl
integer ijmdl
total number of model land points
Definition model_grid.F90:29
model_grid::kgds_mdl
integer, dimension(200) kgds_mdl
holds grib gds info of model grid
Definition model_grid.F90:33
model_grid::resol_mdl
real resol_mdl
approximate model resolution in km.
Definition model_grid.F90:52
model_grid::lat11
real lat11
Corner point latitude (1,1) of model grid.
Definition model_grid.F90:42
model_grid::ipts_mdl
integer, dimension(:), allocatable ipts_mdl
i index of point on full grid
Definition model_grid.F90:30
model_grid::grid_id_mdl
integer grid_id_mdl
grib id of model grid, 4-gaussian, 203-egrid
Definition model_grid.F90:26
model_grid::model_grid_cleanup
subroutine model_grid_cleanup
Clean up allocatable arrays.
Definition model_grid.F90:578
model_grid::lon11
real lon11
Corner point longitude (1,1) of model grid.
Definition model_grid.F90:44
model_grid::lons_mdl
real, dimension(:), allocatable lons_mdl
longitudes of model grid points
Definition model_grid.F90:46
model_grid::read_mdl_grid_info
subroutine read_mdl_grid_info
Read mdl grid.
Definition model_grid.F90:85
program_setup
This module reads in data from the program's configuration namelist.
Definition program_setup.F90:15
program_setup::gfs_lpl_file
character *200, public gfs_lpl_file
GFS gaussian thinned (reduced) grid definition file.
Definition program_setup.F90:29
program_setup::model_lat_file
character *200, public model_lat_file
path/name lats on the model grid
Definition program_setup.F90:32
program_setup::model_lon_file
character *200, public model_lon_file
path/name lons on the model grid
Definition program_setup.F90:33
program_setup::model_lsmask_file
character *200, public model_lsmask_file
path/name nesdis/ims land mask
Definition program_setup.F90:34