emcsfc_snow2mdl  1.3.0
 All Data Structures Files Functions Variables
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_null
subroutine grib2_null(gfld)
Nullify the grib2 gribfield pointers.
Definition: grib_utils.F90:610
model_grid
Read in data defining the model grid.
Definition: model_grid.F90:19
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, to the grib1 grid description section array used by ncep ipolates library.
Definition: grib_utils.F90:139
grib2_free
subroutine grib2_free(gfld)
Deallocate the grib2 gribfield pointers.
Definition: grib_utils.F90:635
model_grid::model_grid_cleanup
subroutine model_grid_cleanup
Clean up allocatable arrays.
Definition: model_grid.F90:577
model_grid::read_mdl_grid_info
subroutine read_mdl_grid_info
Read mdl grid.
Definition: model_grid.F90:84
program_setup
This module reads in data from the program&#39;s configuration namelist.
Definition: program_setup.F90:15
grib_check
subroutine grib_check(file_name, isgrib)
Determine whether file is grib or not.
Definition: grib_utils.F90:23