chgres_cube  1.13.0
 All Data Structures Files Functions Variables
model_grid.F90
Go to the documentation of this file.
1 
4 
9  module model_grid
10 
11  use esmf
12  use esmf_logpublicmod
13 
14  use utilities, only : error_handler, netcdf_err
15  implicit none
16 
17  private
18 
19  character(len=5), allocatable, public :: tiles_target_grid(:)
20 
21  character(len=50), public :: input_grid_type = "latlon"
22 
23 
24  ! Made lsoil_target non-parameter to allow for RAP land surface initiation
25  integer, public :: lsoil_target = 4 ! # soil layers
26 
27  integer, public :: i_input
28 
30  integer, public :: j_input
31 
33  integer, public :: ip1_input
34 
35  integer, public :: jp1_input
36 
37  integer, public :: i_target
38 
40  integer, public :: j_target
41 
43  integer, public :: ip1_target
44 
45  integer, public :: jp1_target
46 
47  integer, public :: num_tiles_input_grid
48 
49  integer, public :: num_tiles_target_grid
50 
51 
52  type(esmf_grid), public :: input_grid
53 
54  type(esmf_grid), public :: target_grid
55 
56 
57  type(esmf_field), public :: latitude_input_grid
58 
59  type(esmf_field), public :: longitude_input_grid
60 
61  type(esmf_field), public :: latitude_s_input_grid
62 
64  type(esmf_field), public :: longitude_s_input_grid
65 
67  type(esmf_field), public :: latitude_w_input_grid
68 
70  type(esmf_field), public :: longitude_w_input_grid
71 
73  type(esmf_field), public :: landmask_target_grid
74 
76  type(esmf_field), public :: land_frac_target_grid
77 
78  type(esmf_field), public :: latitude_target_grid
79 
80  type(esmf_field), public :: latitude_s_target_grid
81 
83  type(esmf_field), public :: latitude_w_target_grid
84 
86  type(esmf_field), public :: longitude_target_grid
87 
88  type(esmf_field), public :: longitude_s_target_grid
89 
91  type(esmf_field), public :: longitude_w_target_grid
92 
94  type(esmf_field), public :: seamask_target_grid
95 
97  type(esmf_field), public :: terrain_target_grid
98 
99 
100  public :: define_target_grid
101  public :: define_input_grid
102  public :: cleanup_input_target_grid_data
103 
104  contains
105 
117  subroutine define_input_grid(localpet, npets)
118 
119  use program_setup, only : input_type
120 
121  implicit none
122 
123  integer, intent(in) :: localpet, npets
124 
125  if (trim(input_type) == "gaussian_nemsio" .or. &
126  trim(input_type) == "gfs_gaussian_nemsio" .or. &
127  trim(input_type) == "gfs_sigio" .or. &
128  trim(input_type) == "gaussian_netcdf") then
129  call define_input_grid_gaussian(npets)
130  elseif (trim(input_type) == "grib2") then
131  call define_input_grid_grib2(npets)
132  else
133  call define_input_grid_mosaic(localpet, npets)
134  endif
135 
136  end subroutine define_input_grid
137 
148  subroutine define_input_grid_gaussian(npets)
149 
150  use nemsio_module
151 
152  use program_setup, only : data_dir_input_grid, &
153  atm_files_input_grid, &
154  sfc_files_input_grid, &
155  input_type, &
156  convert_atm, convert_sfc
157 
158  use sfcio_module
159  use sigio_module
160  use netcdf
161 
162  implicit none
163 
164  integer, intent(in) :: npets
165 
166  character(len=250) :: the_file
167 
168  integer :: i, j, rc, clb(2), cub(2), ncid, id_grid
169  integer(sfcio_intkind) :: rc2
170  integer(sigio_intkind) :: rc3
171 
172  real(esmf_kind_r8), allocatable :: latitude(:,:)
173  real(esmf_kind_r8), allocatable :: longitude(:,:)
174  real(esmf_kind_r8), pointer :: lat_src_ptr(:,:)
175  real(esmf_kind_r8), pointer :: lon_src_ptr(:,:)
176  real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:)
177  real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:)
178  real(esmf_kind_r8) :: deltalon
179  real(esmf_kind_r8), allocatable :: slat(:), wlat(:)
180 
181  type(nemsio_gfile) :: gfile
182  type(esmf_polekind_flag) :: polekindflag(2)
183  type(sfcio_head) :: sfchead
184  type(sigio_head) :: sighead
185 
186  print*,"- DEFINE INPUT GRID OBJECT FOR GAUSSIAN DATA."
187 
188  num_tiles_input_grid = 1
189 
190  if (convert_sfc) then
191  the_file=trim(data_dir_input_grid) // "/" // trim(sfc_files_input_grid(1))
192  elseif (convert_atm) then
193  the_file=trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
194  endif
195 
196  if (trim(input_type) == "gfs_sigio") then ! sigio/sfcio format, used by
197  ! spectral gfs prior to 7/19/2017.
198 
199  if (convert_sfc) then ! sfcio format
200  print*,"- OPEN AND READ ", trim(the_file)
201  call sfcio_sropen(21, trim(the_file), rc2)
202  if (rc2 /= 0) call error_handler("OPENING FILE", rc2)
203  call sfcio_srhead(21, sfchead, rc2)
204  if (rc2 /= 0) call error_handler("READING FILE", rc2)
205  call sfcio_sclose(21, rc2)
206  i_input = sfchead%lonb
207  j_input = sfchead%latb
208  elseif (convert_atm) then ! sigio format
209  print*,"- OPEN AND READ ", trim(the_file)
210  call sigio_sropen(21, trim(the_file), rc3)
211  if (rc3 /= 0) call error_handler("OPENING FILE", rc3)
212  call sigio_srhead(21, sighead, rc3)
213  if (rc3 /= 0) call error_handler("READING FILE", rc3)
214  call sigio_sclose(21, rc3)
215  i_input = sighead%lonb
216  j_input = sighead%latb
217  endif
218 
219  elseif (trim(input_type) == "gaussian_netcdf") then
220 
221  print*,'- OPEN AND READ: ',trim(the_file)
222  rc=nf90_open(trim(the_file),nf90_nowrite,ncid)
223  call netcdf_err(rc, 'opening file')
224 
225  print*,"- READ grid_xt"
226  rc=nf90_inq_dimid(ncid, 'grid_xt', id_grid)
227  call netcdf_err(rc, 'reading grid_xt id')
228  rc=nf90_inquire_dimension(ncid,id_grid,len=i_input)
229  call netcdf_err(rc, 'reading grid_xt')
230 
231  print*,"- READ grid_yt"
232  rc=nf90_inq_dimid(ncid, 'grid_yt', id_grid)
233  call netcdf_err(rc, 'reading grid_yt id')
234  rc=nf90_inquire_dimension(ncid,id_grid,len=j_input)
235  call netcdf_err(rc, 'reading grid_yt')
236 
237  rc = nf90_close(ncid)
238 
239  else ! nemsio format
240 
241  call nemsio_init(iret=rc)
242 
243  print*,"- OPEN AND READ ", trim(the_file)
244  call nemsio_open(gfile, the_file, "read", iret=rc)
245  if (rc /= 0) call error_handler("OPENING FILE", rc)
246 
247  call nemsio_getfilehead(gfile, iret=rc, dimx=i_input, dimy=j_input)
248  if (rc /= 0) call error_handler("READING FILE", rc)
249 
250  call nemsio_close(gfile)
251 
252  endif
253 
254  ip1_input = i_input + 1
255  jp1_input = j_input + 1
256 
257  polekindflag(1:2) = esmf_polekind_monopole
258 
259  print*,"- CALL GridCreate1PeriDim FOR INPUT GRID."
260  input_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
261  maxindex=(/i_input,j_input/), &
262  polekindflag=polekindflag, &
263  periodicdim=1, &
264  poledim=2, &
265  coordsys=esmf_coordsys_sph_deg, &
266  regdecomp=(/1,npets/), &
267  indexflag=esmf_index_global, rc=rc)
268  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
269  call error_handler("IN GridCreate1PeriDim", rc)
270 
271  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
272  latitude_input_grid = esmf_fieldcreate(input_grid, &
273  typekind=esmf_typekind_r8, &
274  staggerloc=esmf_staggerloc_center, &
275  name="input_grid_latitude", rc=rc)
276  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
277  call error_handler("IN FieldCreate", rc)
278 
279  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
280  longitude_input_grid = esmf_fieldcreate(input_grid, &
281  typekind=esmf_typekind_r8, &
282  staggerloc=esmf_staggerloc_center, &
283  name="input_grid_longitude", rc=rc)
284  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
285  call error_handler("IN FieldCreate", rc)
286 
287  allocate(longitude(i_input,j_input))
288  allocate(latitude(i_input,j_input))
289 
290  deltalon = 360.0_esmf_kind_r8 / real(i_input,kind=esmf_kind_r8)
291  do i = 1, i_input
292  longitude(i,:) = real((i-1),kind=esmf_kind_r8) * deltalon
293  enddo
294 
295  allocate(slat(j_input))
296  allocate(wlat(j_input))
297  call splat(4, j_input, slat, wlat)
298 
299  do i = 1, j_input
300  latitude(:,i) = 90.0_esmf_kind_r8 - (acos(slat(i))* 180.0_esmf_kind_r8 / &
301  (4.0_esmf_kind_r8*atan(1.0_esmf_kind_r8)))
302  enddo
303 
304  deallocate(slat, wlat)
305 
306  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
307  call esmf_fieldscatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
308  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
309  call error_handler("IN FieldScatter", rc)
310 
311  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE."
312  call esmf_fieldscatter(latitude_input_grid, latitude, rootpet=0, rc=rc)
313  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
314  call error_handler("IN FieldScatter", rc)
315 
316  print*,"- CALL GridAddCoord FOR INPUT GRID."
317  call esmf_gridaddcoord(input_grid, &
318  staggerloc=esmf_staggerloc_center, rc=rc)
319  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
320  call error_handler("IN GridAddCoord", rc)
321 
322  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
323  nullify(lon_src_ptr)
324  call esmf_gridgetcoord(input_grid, &
325  staggerloc=esmf_staggerloc_center, &
326  coorddim=1, &
327  farrayptr=lon_src_ptr, rc=rc)
328  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
329  call error_handler("IN GridGetCoord", rc)
330 
331  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
332  nullify(lat_src_ptr)
333  call esmf_gridgetcoord(input_grid, &
334  staggerloc=esmf_staggerloc_center, &
335  coorddim=2, &
336  computationallbound=clb, &
337  computationalubound=cub, &
338  farrayptr=lat_src_ptr, rc=rc)
339  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
340  call error_handler("IN GridGetCoord", rc)
341 
342  do j = clb(2), cub(2)
343  do i = clb(1), cub(1)
344  lon_src_ptr(i,j) = longitude(i,j)
345  if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8
346  lat_src_ptr(i,j) = latitude(i,j)
347  enddo
348  enddo
349 
350  print*,"- CALL GridAddCoord FOR INPUT GRID."
351  call esmf_gridaddcoord(input_grid, &
352  staggerloc=esmf_staggerloc_corner, rc=rc)
353  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
354  call error_handler("IN GridAddCoord", rc)
355 
356  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
357  nullify(lon_corner_src_ptr)
358  call esmf_gridgetcoord(input_grid, &
359  staggerloc=esmf_staggerloc_corner, &
360  coorddim=1, &
361  farrayptr=lon_corner_src_ptr, rc=rc)
362  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
363  call error_handler("IN GridGetCoord", rc)
364 
365  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
366  nullify(lat_corner_src_ptr)
367  call esmf_gridgetcoord(input_grid, &
368  staggerloc=esmf_staggerloc_corner, &
369  coorddim=2, &
370  computationallbound=clb, &
371  computationalubound=cub, &
372  farrayptr=lat_corner_src_ptr, rc=rc)
373  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
374  call error_handler("IN GridGetCoord", rc)
375 
376  do j = clb(2), cub(2)
377  do i = clb(1), cub(1)
378  lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon)
379  if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8
380  if (j == 1) then
381  lat_corner_src_ptr(i,j) = 90.0_esmf_kind_r8
382  cycle
383  endif
384  if (j == jp1_input) then
385  lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8
386  cycle
387  endif
388  lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j))
389  enddo
390  enddo
391 
392  deallocate(latitude,longitude)
393 
394  end subroutine define_input_grid_gaussian
395 
402  subroutine define_input_grid_mosaic(localpet, npets)
403 
404  use netcdf
405  use program_setup, only : mosaic_file_input_grid, &
406  orog_dir_input_grid, &
407  orog_files_input_grid
408 
409  implicit none
410 
411  character(len=500) :: the_file
412 
413  integer, intent(in) :: localpet, npets
414 
415  integer :: id_tiles, id_dim, tile
416  integer :: extra, error, ncid
417  integer, allocatable :: decomptile(:,:)
418 
419  real(esmf_kind_r8), allocatable :: latitude_one_tile(:,:)
420  real(esmf_kind_r8), allocatable :: latitude_s_one_tile(:,:)
421  real(esmf_kind_r8), allocatable :: latitude_w_one_tile(:,:)
422  real(esmf_kind_r8), allocatable :: longitude_one_tile(:,:)
423  real(esmf_kind_r8), allocatable :: longitude_s_one_tile(:,:)
424  real(esmf_kind_r8), allocatable :: longitude_w_one_tile(:,:)
425 
426  print*,'- OPEN INPUT GRID MOSAIC FILE: ',trim(mosaic_file_input_grid)
427  error=nf90_open(trim(mosaic_file_input_grid),nf90_nowrite,ncid)
428  call netcdf_err(error, 'opening grid mosaic file')
429 
430  print*,"- READ NUMBER OF TILES"
431  error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
432  call netcdf_err(error, 'reading ntiles id')
433  error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles_input_grid)
434  call netcdf_err(error, 'reading ntiles')
435 
436  error = nf90_close(ncid)
437 
438  print*,'- NUMBER OF TILES, INPUT MODEL GRID IS ', num_tiles_input_grid
439 
440  if (mod(npets,num_tiles_input_grid) /= 0) then
441  call error_handler("MUST RUN WITH A TASK COUNT THAT IS A MULTIPLE OF 6.", 1)
442  endif
443 
444 !-----------------------------------------------------------------------
445 ! Create ESMF grid object for the model grid.
446 !-----------------------------------------------------------------------
447 
448  extra = npets / num_tiles_input_grid
449 
450  allocate(decomptile(2,num_tiles_input_grid))
451 
452  do tile = 1, num_tiles_input_grid
453  decomptile(:,tile)=(/1,extra/)
454  enddo
455 
456  print*,"- CALL GridCreateMosaic FOR INPUT MODEL GRID"
457  input_grid = esmf_gridcreatemosaic(filename=trim(mosaic_file_input_grid), &
458  regdecompptile=decomptile, &
459  staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
460  esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
461  indexflag=esmf_index_global, &
462  tilefilepath=trim(orog_dir_input_grid), &
463  rc=error)
464  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
465  call error_handler("IN GridCreateMosaic", error)
466 
467 !-----------------------------------------------------------------------
468 ! Read the mask and lat/lons.
469 !-----------------------------------------------------------------------
470 
471  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
472  latitude_input_grid = esmf_fieldcreate(input_grid, &
473  typekind=esmf_typekind_r8, &
474  staggerloc=esmf_staggerloc_center, &
475  name="input_grid_latitude", &
476  rc=error)
477  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
478  call error_handler("IN FieldCreate", error)
479 
480  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
481  longitude_input_grid = esmf_fieldcreate(input_grid, &
482  typekind=esmf_typekind_r8, &
483  staggerloc=esmf_staggerloc_center, &
484  name="input_grid_longitude", &
485  rc=error)
486  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
487  call error_handler("IN FieldCreate", error)
488 
489  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE_S."
490  latitude_s_input_grid = esmf_fieldcreate(input_grid, &
491  typekind=esmf_typekind_r8, &
492  staggerloc=esmf_staggerloc_edge2, &
493  name="input_grid_latitude_s", &
494  rc=error)
495  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
496  call error_handler("IN FieldCreate", error)
497 
498  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE_S."
499  longitude_s_input_grid = esmf_fieldcreate(input_grid, &
500  typekind=esmf_typekind_r8, &
501  staggerloc=esmf_staggerloc_edge2, &
502  name="input_grid_longitude_s", &
503  rc=error)
504  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
505  call error_handler("IN FieldCreate", error)
506 
507  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE_W."
508  latitude_w_input_grid = esmf_fieldcreate(input_grid, &
509  typekind=esmf_typekind_r8, &
510  staggerloc=esmf_staggerloc_edge1, &
511  name="input_grid_latitude_w", &
512  rc=error)
513  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
514  call error_handler("IN FieldCreate", error)
515 
516  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE_W."
517  longitude_w_input_grid = esmf_fieldcreate(input_grid, &
518  typekind=esmf_typekind_r8, &
519  staggerloc=esmf_staggerloc_edge1, &
520  name="input_grid_longitude_w", &
521  rc=error)
522  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
523  call error_handler("IN FieldCreate", error)
524 
525  the_file = trim(orog_dir_input_grid) // trim(orog_files_input_grid(1))
526 
527  print*,'- OPEN FIRST INPUT GRID OROGRAPHY FILE: ',trim(the_file)
528  error=nf90_open(trim(the_file),nf90_nowrite,ncid)
529  call netcdf_err(error, 'opening ororgraphy file')
530  print*,"- READ GRID DIMENSIONS"
531  error=nf90_inq_dimid(ncid, 'lon', id_dim)
532  call netcdf_err(error, 'reading lon id')
533  error=nf90_inquire_dimension(ncid,id_dim,len=i_input)
534  call netcdf_err(error, 'reading lon')
535  error=nf90_inq_dimid(ncid, 'lat', id_dim)
536  call netcdf_err(error, 'reading lat id')
537  error=nf90_inquire_dimension(ncid,id_dim,len=j_input)
538  call netcdf_err(error, 'reading lat')
539  error = nf90_close(ncid)
540 
541  print*,"- I/J DIMENSIONS OF THE INPUT GRID TILES ", i_input, j_input
542 
543  ip1_input = i_input + 1
544  jp1_input = j_input + 1
545 
546  if (localpet == 0) then
547  allocate(longitude_one_tile(i_input,j_input))
548  allocate(longitude_s_one_tile(i_input,jp1_input))
549  allocate(longitude_w_one_tile(ip1_input,j_input))
550  allocate(latitude_one_tile(i_input,j_input))
551  allocate(latitude_s_one_tile(i_input,jp1_input))
552  allocate(latitude_w_one_tile(ip1_input,j_input))
553  else
554  allocate(longitude_one_tile(0,0))
555  allocate(longitude_s_one_tile(0,0))
556  allocate(longitude_w_one_tile(0,0))
557  allocate(latitude_one_tile(0,0))
558  allocate(latitude_s_one_tile(0,0))
559  allocate(latitude_w_one_tile(0,0))
560  endif
561 
562  do tile = 1, num_tiles_input_grid
563  if (localpet == 0) then
564  call get_model_latlons(mosaic_file_input_grid, orog_dir_input_grid, num_tiles_input_grid, tile, &
565  i_input, j_input, ip1_input, jp1_input, latitude_one_tile, &
566  latitude_s_one_tile, latitude_w_one_tile, longitude_one_tile, &
567  longitude_s_one_tile, longitude_w_one_tile)
568  endif
569  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE. TILE IS: ", tile
570  call esmf_fieldscatter(latitude_input_grid, latitude_one_tile, rootpet=0, tile=tile, rc=error)
571  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
572  call error_handler("IN FieldScatter", error)
573  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE. TILE IS: ", tile
574  call esmf_fieldscatter(longitude_input_grid, longitude_one_tile, rootpet=0, tile=tile, rc=error)
575  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
576  call error_handler("IN FieldScatter", error)
577  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE_S. TILE IS: ", tile
578  call esmf_fieldscatter(latitude_s_input_grid, latitude_s_one_tile, rootpet=0, tile=tile, rc=error)
579  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
580  call error_handler("IN FieldScatter", error)
581  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE_S. TILE IS: ", tile
582  call esmf_fieldscatter(longitude_s_input_grid, longitude_s_one_tile, rootpet=0, tile=tile, rc=error)
583  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
584  call error_handler("IN FieldScatter", error)
585  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE_W. TILE IS: ", tile
586  call esmf_fieldscatter(latitude_w_input_grid, latitude_w_one_tile, rootpet=0, tile=tile, rc=error)
587  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
588  call error_handler("IN FieldScatter", error)
589  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE_W. TILE IS: ", tile
590  call esmf_fieldscatter(longitude_w_input_grid, longitude_w_one_tile, rootpet=0, tile=tile, rc=error)
591  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
592  call error_handler("IN FieldScatter", error)
593  enddo
594 
595  deallocate(longitude_one_tile)
596  deallocate(longitude_s_one_tile)
597  deallocate(longitude_w_one_tile)
598  deallocate(latitude_one_tile)
599  deallocate(latitude_s_one_tile)
600  deallocate(latitude_w_one_tile)
601 
602  end subroutine define_input_grid_mosaic
603 
610  subroutine define_input_grid_grib2(npets)
611 
612  use grib_mod
613  use gdswzd_mod
614  use program_setup, only : grib2_file_input_grid, data_dir_input_grid
615 
616  implicit none
617 
618  integer, intent(in) :: npets
619 
620  character(len=500) :: the_file
621 
622  integer :: i, j, k, jdisc, jgdtn, jpdtn, lugb, lugi
623  integer :: jids(200), jgdt(200), jpdt(200), rc
624  integer :: kgds(200), nret, clb(2), cub(2)
625 
626  logical :: unpack
627 
628  real :: res
629  real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:)
630  real, allocatable :: rlon_corner(:,:),rlat_corner(:,:)
631  real, allocatable :: rlon_diff(:,:),rlat_diff(:,:)
632  real, allocatable :: xpts_corner(:,:),ypts_corner(:,:)
633  real(esmf_kind_r8), allocatable :: latitude(:,:)
634  real(esmf_kind_r8), allocatable :: longitude(:,:)
635  real(esmf_kind_r8), allocatable :: latitude_corner(:,:)
636  real(esmf_kind_r8), allocatable :: longitude_corner(:,:)
637  real(esmf_kind_r8), pointer :: lat_src_ptr(:,:)
638  real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:)
639  real(esmf_kind_r8), pointer :: lon_src_ptr(:,:)
640  real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:)
641 
642  type(esmf_polekind_flag) :: polekindflag(2)
643 
644  type(gribfield) :: gfld
645 
646  the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid
647 
648  lugb=12
649 
650  print*,"- OPEN AND READ INPUT DATA GRIB2 FILE: ", trim(the_file)
651  call baopenr(lugb,the_file,rc)
652  if (rc /= 0) call error_handler("OPENING FILE", rc)
653 
654 ! Read the first record and get the grid definition template.
655 
656  j = 0 ! Search at beginning of file
657  lugi = 0 ! No grib index file
658  jdisc = -1 ! Search for any discipline
659  jpdtn = -1 ! Search for any product definition template number
660  jgdtn = -1 ! Search for any grid definition template number
661  jids = -9999 ! Array of values in identification section, set to wildcard.
662  jgdt = -9999 ! Array of values in grid definition template, set to wildcard.
663  jpdt = -9999 ! Array of values in product definition template, set to wildcard.
664  unpack = .false. ! unpack data
665 
666  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
667  unpack, k, gfld, rc)
668  if (rc /= 0) call error_handler("DEGRIBBING INPUT FILE.", rc)
669 
670  call baclose(lugb,rc)
671 
672  if (gfld%igdtnum == 0) then
673  print*,"- INPUT DATA ON LAT/LON GRID."
674  input_grid_type = 'latlon'
675  elseif (gfld%igdtnum == 30) then
676  print*,"- INPUT DATA ON LAMBERT CONFORMAL GRID."
677  input_grid_type = 'lambert'
678  elseif (gfld%igdtnum == 32769 .or. gfld%igdtnum == 1) then
679  print*,"- INPUT DATA ON ROTATED LAT/LON GRID."
680  input_grid_type = 'rotated_latlon'
681  else
682  call error_handler("INPUT GRID TEMPLATE NOT SUPPORTED.", 2)
683  endif
684 
685  kgds = 0
686  call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, i_input, j_input, res)
687 
688  ip1_input = i_input + 1
689  jp1_input = j_input + 1
690 
691  allocate(rlat(i_input,j_input))
692  allocate(rlon(i_input,j_input))
693  allocate(rlat_diff(i_input,j_input))
694  allocate(rlon_diff(i_input,j_input))
695  allocate(xpts(i_input,j_input))
696  allocate(ypts(i_input,j_input))
697  allocate(rlat_corner(ip1_input,jp1_input))
698  allocate(rlon_corner(ip1_input,jp1_input))
699  allocate(xpts_corner(ip1_input,jp1_input))
700  allocate(ypts_corner(ip1_input,jp1_input))
701 
702  do j = 1, j_input
703  do i = 1, i_input
704  xpts(i,j) = float(i)
705  ypts(i,j) = float(j)
706  enddo
707  enddo
708 
709  print*,"- COMPUTE GRID CELL CENTER COORDINATES."
710  call gdswzd(kgds,1,(i_input*j_input),-9999.,xpts,ypts,rlon,rlat,nret)
711 
712  if (nret /= (i_input*j_input)) then
713  call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2)
714  endif
715 
716  deallocate(xpts, ypts)
717 
718  do j = 1, jp1_input
719  do i = 1, ip1_input
720  xpts_corner(i,j) = float(i) - 0.5
721  ypts_corner(i,j) = float(j) - 0.5
722  enddo
723  enddo
724 
725  print*,"- COMPUTE GRID CELL CORNER COORDINATES."
726  call gdswzd(kgds,1,(ip1_input*jp1_input),-9999.,xpts_corner,ypts_corner,rlon_corner,rlat_corner,nret)
727 
728  if (nret /= (ip1_input*jp1_input)) then
729  call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2)
730  endif
731 
732  deallocate(xpts_corner, ypts_corner)
733 
734  if (gfld%igdtnum == 0) then ! gfs lat/lon data
735 
736  print*,"- CALL GridCreate1PeriDim FOR INPUT GRID."
737 
738  polekindflag(1:2) = esmf_polekind_monopole
739 
740  input_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
741  maxindex=(/i_input,j_input/), &
742  polekindflag=polekindflag, &
743  periodicdim=1, &
744  poledim=2, &
745  coordsys=esmf_coordsys_sph_deg, &
746  regdecomp=(/1,npets/), &
747  indexflag=esmf_index_global, rc=rc)
748 
749  else
750 
751  print*,"- CALL GridCreateNoPeriDim FOR INPUT GRID."
752 
753  input_grid = esmf_gridcreatenoperidim(maxindex=(/i_input,j_input/), &
754  indexflag=esmf_index_global, &
755  rc=rc)
756 
757  endif
758 
759  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
760  call error_handler("IN GridCreate1PeriDim", rc)
761 
762  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
763  latitude_input_grid = esmf_fieldcreate(input_grid, &
764  typekind=esmf_typekind_r8, &
765  staggerloc=esmf_staggerloc_center, &
766  name="input_grid_latitude", rc=rc)
767  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
768  call error_handler("IN FieldCreate", rc)
769 
770  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
771  longitude_input_grid = esmf_fieldcreate(input_grid, &
772  typekind=esmf_typekind_r8, &
773  staggerloc=esmf_staggerloc_center, &
774  name="input_grid_longitude", rc=rc)
775  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
776  call error_handler("IN FieldCreate", rc)
777 
778  allocate(latitude(i_input,j_input))
779  allocate(longitude(i_input,j_input))
780 
781  latitude = rlat
782  longitude = rlon
783 
784  deallocate (rlat, rlon)
785 
786  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
787  call esmf_fieldscatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
788  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
789  call error_handler("IN FieldScatter", rc)
790 
791  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE."
792  call esmf_fieldscatter(latitude_input_grid, latitude, rootpet=0, rc=rc)
793  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
794  call error_handler("IN FieldScatter", rc)
795 
796  print*,"- CALL GridAddCoord FOR INPUT GRID."
797  call esmf_gridaddcoord(input_grid, &
798  staggerloc=esmf_staggerloc_center, rc=rc)
799  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
800  call error_handler("IN GridAddCoord", rc)
801 
802  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
803  nullify(lon_src_ptr)
804  call esmf_gridgetcoord(input_grid, &
805  staggerloc=esmf_staggerloc_center, &
806  coorddim=1, &
807  farrayptr=lon_src_ptr, rc=rc)
808  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
809  call error_handler("IN GridGetCoord", rc)
810 
811  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
812  nullify(lat_src_ptr)
813  call esmf_gridgetcoord(input_grid, &
814  staggerloc=esmf_staggerloc_center, &
815  coorddim=2, &
816  computationallbound=clb, &
817  computationalubound=cub, &
818  farrayptr=lat_src_ptr, rc=rc)
819  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
820  call error_handler("IN GridGetCoord", rc)
821 
822  do j = clb(2), cub(2)
823  do i = clb(1), cub(1)
824  lon_src_ptr(i,j) = longitude(i,j)
825  if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8
826  lat_src_ptr(i,j) = latitude(i,j)
827  enddo
828  enddo
829 
830  deallocate(latitude, longitude)
831 
832  allocate(latitude_corner(ip1_input,jp1_input))
833  allocate(longitude_corner(ip1_input,jp1_input))
834 
835  latitude_corner = rlat_corner
836  longitude_corner = rlon_corner
837 
838  deallocate (rlat_corner, rlon_corner)
839 
840  print*,"- CALL GridAddCoord FOR INPUT GRID."
841  call esmf_gridaddcoord(input_grid, &
842  staggerloc=esmf_staggerloc_corner, rc=rc)
843  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
844  call error_handler("IN GridAddCoord", rc)
845 
846  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
847  nullify(lon_corner_src_ptr)
848  call esmf_gridgetcoord(input_grid, &
849  staggerloc=esmf_staggerloc_corner, &
850  coorddim=1, &
851  farrayptr=lon_corner_src_ptr, rc=rc)
852  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
853  call error_handler("IN GridGetCoord", rc)
854 
855  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
856  nullify(lat_corner_src_ptr)
857  call esmf_gridgetcoord(input_grid, &
858  staggerloc=esmf_staggerloc_corner, &
859  coorddim=2, &
860  computationallbound=clb, &
861  computationalubound=cub, &
862  farrayptr=lat_corner_src_ptr, rc=rc)
863  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
864  call error_handler("IN GridGetCoord", rc)
865 
866  do j = clb(2), cub(2)
867  do i = clb(1), cub(1)
868  lon_corner_src_ptr(i,j) = longitude_corner(i,j)
869  if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8
870  lat_corner_src_ptr(i,j) = latitude_corner(i,j)
871  enddo
872  enddo
873 
874  deallocate(latitude_corner, longitude_corner)
875 
876  end subroutine define_input_grid_grib2
877 
883  subroutine define_target_grid(localpet, npets)
884 
885  use netcdf
886  use program_setup, only : mosaic_file_target_grid, &
887  orog_dir_target_grid, &
888  orog_files_target_grid, &
889  nsoill_out
890 
891  implicit none
892 
893  integer, intent(in) :: localpet, npets
894 
895  character(len=500) :: the_file
896 
897  integer :: error, ncid, extra
898  integer :: id_tiles
899  integer :: id_dim, id_grid_tiles
900  integer :: tile
901  integer, allocatable :: decomptile(:,:)
902  integer(esmf_kind_i8), allocatable :: landmask_one_tile(:,:)
903  integer(esmf_kind_i8), allocatable :: seamask_one_tile(:,:)
904 
905  real(esmf_kind_r8), allocatable :: land_frac_one_tile(:,:)
906  real(esmf_kind_r8), allocatable :: latitude_one_tile(:,:)
907  real(esmf_kind_r8), allocatable :: latitude_s_one_tile(:,:)
908  real(esmf_kind_r8), allocatable :: latitude_w_one_tile(:,:)
909  real(esmf_kind_r8), allocatable :: longitude_one_tile(:,:)
910  real(esmf_kind_r8), allocatable :: longitude_s_one_tile(:,:)
911  real(esmf_kind_r8), allocatable :: longitude_w_one_tile(:,:)
912  real(esmf_kind_r8), allocatable :: terrain_one_tile(:,:)
913 
914  lsoil_target = nsoill_out
915 
916  print*,'- OPEN TARGET GRID MOSAIC FILE: ',trim(mosaic_file_target_grid)
917  error=nf90_open(trim(mosaic_file_target_grid),nf90_nowrite,ncid)
918  call netcdf_err(error, 'opening grid mosaic file')
919 
920  print*,"- READ NUMBER OF TILES"
921  error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
922  call netcdf_err(error, 'reading ntile id')
923  error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles_target_grid)
924  call netcdf_err(error, 'reading ntiles')
925  error=nf90_inq_varid(ncid, 'gridtiles', id_grid_tiles)
926  call netcdf_err(error, 'reading gridtiles id')
927  allocate(tiles_target_grid(num_tiles_target_grid))
928  tiles_target_grid="NULL"
929  print*,"- READ TILE NAMES"
930  error=nf90_get_var(ncid, id_grid_tiles, tiles_target_grid)
931  call netcdf_err(error, 'reading gridtiles')
932 
933  error = nf90_close(ncid)
934 
935  print*,'- NUMBER OF TILES, TARGET MODEL GRID IS ', num_tiles_target_grid
936 
937  if (mod(npets,num_tiles_target_grid) /= 0) then
938  call error_handler("MUST RUN WITH TASK COUNT THAT IS A MULTIPLE OF # OF TILES.", 1)
939  endif
940 
941 !-----------------------------------------------------------------------
942 ! Get the model grid specs and land mask from the orography files.
943 !-----------------------------------------------------------------------
944 
945  the_file = trim(orog_dir_target_grid) // trim(orog_files_target_grid(1))
946 
947  print*,'- OPEN FIRST TARGET GRID OROGRAPHY FILE: ',trim(the_file)
948  error=nf90_open(trim(the_file),nf90_nowrite,ncid)
949  call netcdf_err(error, 'opening orography file')
950  print*,"- READ GRID DIMENSIONS"
951  error=nf90_inq_dimid(ncid, 'lon', id_dim)
952  call netcdf_err(error, 'reading lon id')
953  error=nf90_inquire_dimension(ncid,id_dim,len=i_target)
954  call netcdf_err(error, 'reading lon')
955  error=nf90_inq_dimid(ncid, 'lat', id_dim)
956  call netcdf_err(error, 'reading lat id')
957  error=nf90_inquire_dimension(ncid,id_dim,len=j_target)
958  call netcdf_err(error, 'reading lat')
959  error = nf90_close(ncid)
960 
961  print*,"- I/J DIMENSIONS OF THE TARGET GRID TILES ", i_target, j_target
962 
963  ip1_target = i_target + 1
964  jp1_target = j_target + 1
965 
966 !-----------------------------------------------------------------------
967 ! Create ESMF grid object for the model grid.
968 !-----------------------------------------------------------------------
969 
970  extra = npets / num_tiles_target_grid
971 
972  allocate(decomptile(2,num_tiles_target_grid))
973 
974  do tile = 1, num_tiles_target_grid
975  decomptile(:,tile)=(/1,extra/)
976  enddo
977 
978  print*,"- CALL GridCreateMosaic FOR TARGET GRID"
979  target_grid = esmf_gridcreatemosaic(filename=trim(mosaic_file_target_grid), &
980  regdecompptile=decomptile, &
981  staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
982  esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
983  indexflag=esmf_index_global, &
984  tilefilepath=trim(orog_dir_target_grid), rc=error)
985  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
986  call error_handler("IN GridCreateMosaic", error)
987 
988 !-----------------------------------------------------------------------
989 ! Set target model landmask (1 - land, 0 - not land) and
990 ! seamask (1 - non-land, 0 -land). Read lat/lon on target grid.
991 !-----------------------------------------------------------------------
992 
993  print*,"- CALL FieldCreate FOR TARGET GRID LAND FRACTION."
994  land_frac_target_grid = esmf_fieldcreate(target_grid, &
995  typekind=esmf_typekind_r8, &
996  staggerloc=esmf_staggerloc_center, &
997  name="target_grid_landmask", rc=error)
998  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
999  call error_handler("IN FieldCreate", error)
1000 
1001  print*,"- CALL FieldCreate FOR TARGET GRID LANDMASK."
1002  landmask_target_grid = esmf_fieldcreate(target_grid, &
1003  typekind=esmf_typekind_i8, &
1004  staggerloc=esmf_staggerloc_center, &
1005  name="target_grid_landmask", rc=error)
1006  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1007  call error_handler("IN FieldCreate", error)
1008 
1009  print*,"- CALL FieldCreate FOR TARGET GRID SEAMASK."
1010  seamask_target_grid = esmf_fieldcreate(target_grid, &
1011  typekind=esmf_typekind_i8, &
1012  staggerloc=esmf_staggerloc_center, &
1013  name="target_grid_seamask", rc=error)
1014  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1015  call error_handler("IN FieldCreate", error)
1016 
1017  print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE."
1018  latitude_target_grid = esmf_fieldcreate(target_grid, &
1019  typekind=esmf_typekind_r8, &
1020  staggerloc=esmf_staggerloc_center, &
1021  name="target_grid_latitude", rc=error)
1022  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1023  call error_handler("IN FieldCreate", error)
1024 
1025  print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE_S."
1026  latitude_s_target_grid = esmf_fieldcreate(target_grid, &
1027  typekind=esmf_typekind_r8, &
1028  staggerloc=esmf_staggerloc_edge2, &
1029  name="target_grid_latitude_s", rc=error)
1030  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1031  call error_handler("IN FieldCreate", error)
1032 
1033  print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE_W."
1034  latitude_w_target_grid = esmf_fieldcreate(target_grid, &
1035  typekind=esmf_typekind_r8, &
1036  staggerloc=esmf_staggerloc_edge1, &
1037  name="target_grid_latitude_w", &
1038  rc=error)
1039  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1040  call error_handler("IN FieldCreate", error)
1041 
1042  print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE."
1043  longitude_target_grid = esmf_fieldcreate(target_grid, &
1044  typekind=esmf_typekind_r8, &
1045  staggerloc=esmf_staggerloc_center, &
1046  name="target_grid_longitude", &
1047  rc=error)
1048  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1049  call error_handler("IN FieldCreate", error)
1050 
1051  print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE_S."
1052  longitude_s_target_grid = esmf_fieldcreate(target_grid, &
1053  typekind=esmf_typekind_r8, &
1054  staggerloc=esmf_staggerloc_edge2, &
1055  name="target_grid_longitude_s", &
1056  rc=error)
1057  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1058  call error_handler("IN FieldCreate", error)
1059 
1060  print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE_W."
1061  longitude_w_target_grid = esmf_fieldcreate(target_grid, &
1062  typekind=esmf_typekind_r8, &
1063  staggerloc=esmf_staggerloc_edge1, &
1064  name="target_grid_longitude_w", &
1065  rc=error)
1066  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1067  call error_handler("IN FieldCreate", error)
1068 
1069  print*,"- CALL FieldCreate FOR TARGET GRID TERRAIN."
1070  terrain_target_grid = esmf_fieldcreate(target_grid, &
1071  typekind=esmf_typekind_r8, &
1072  staggerloc=esmf_staggerloc_center, &
1073  name="target_grid_terrain", &
1074  rc=error)
1075  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1076  call error_handler("IN FieldCreate", error)
1077 
1078  if (localpet == 0) then
1079  allocate(land_frac_one_tile(i_target,j_target))
1080  allocate(landmask_one_tile(i_target,j_target))
1081  allocate(seamask_one_tile(i_target,j_target))
1082  allocate(latitude_one_tile(i_target,j_target))
1083  allocate(latitude_s_one_tile(i_target,jp1_target))
1084  allocate(latitude_w_one_tile(ip1_target,j_target))
1085  allocate(longitude_one_tile(i_target,j_target))
1086  allocate(longitude_s_one_tile(i_target,jp1_target))
1087  allocate(longitude_w_one_tile(ip1_target,j_target))
1088  allocate(terrain_one_tile(i_target,j_target))
1089  else
1090  allocate(land_frac_one_tile(0,0))
1091  allocate(landmask_one_tile(0,0))
1092  allocate(seamask_one_tile(0,0))
1093  allocate(longitude_one_tile(0,0))
1094  allocate(longitude_s_one_tile(0,0))
1095  allocate(longitude_w_one_tile(0,0))
1096  allocate(latitude_one_tile(0,0))
1097  allocate(latitude_s_one_tile(0,0))
1098  allocate(latitude_w_one_tile(0,0))
1099  allocate(terrain_one_tile(0,0))
1100  endif
1101 
1102  do tile = 1, num_tiles_target_grid
1103  if (localpet == 0) then
1104  the_file = trim(orog_dir_target_grid) // trim(orog_files_target_grid(tile))
1105  call get_model_mask_terrain(trim(the_file), i_target, j_target, landmask_one_tile, &
1106  terrain_one_tile, land_frac_one_tile)
1107 
1108  seamask_one_tile = 0 ! all land
1109  where(floor(land_frac_one_tile) == 0) seamask_one_tile = 1 ! at least some non-land.
1110  landmask_one_tile = 0 ! all non-land
1111  where(ceiling(land_frac_one_tile) == 1) landmask_one_tile = 1 ! at least some land
1112 
1113  call get_model_latlons(mosaic_file_target_grid, orog_dir_target_grid, num_tiles_target_grid, tile, &
1114  i_target, j_target, ip1_target, jp1_target, latitude_one_tile, &
1115  latitude_s_one_tile, latitude_w_one_tile, longitude_one_tile, &
1116  longitude_s_one_tile, longitude_w_one_tile)
1117  endif
1118  print*,"- CALL FieldScatter FOR TARGET GRID LAND FRACTION. TILE IS: ", tile
1119  call esmf_fieldscatter(land_frac_target_grid, land_frac_one_tile, rootpet=0, tile=tile, rc=error)
1120  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1121  call error_handler("IN FieldScatter", error)
1122  print*,"- CALL FieldScatter FOR TARGET GRID LANDMASK. TILE IS: ", tile
1123  call esmf_fieldscatter(landmask_target_grid, landmask_one_tile, rootpet=0, tile=tile, rc=error)
1124  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1125  call error_handler("IN FieldScatter", error)
1126  print*,"- CALL FieldScatter FOR TARGET GRID SEAMASK. TILE IS: ", tile
1127  call esmf_fieldscatter(seamask_target_grid, seamask_one_tile, rootpet=0, tile=tile, rc=error)
1128  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1129  call error_handler("IN FieldScatter", error)
1130  print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE. TILE IS: ", tile
1131  call esmf_fieldscatter(longitude_target_grid, longitude_one_tile, rootpet=0, tile=tile, rc=error)
1132  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1133  call error_handler("IN FieldScatter", error)
1134  print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE_S. TILE IS: ", tile
1135  call esmf_fieldscatter(longitude_s_target_grid, longitude_s_one_tile, rootpet=0, tile=tile, rc=error)
1136  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1137  call error_handler("IN FieldScatter", error)
1138  print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE_W. TILE IS: ", tile
1139  call esmf_fieldscatter(longitude_w_target_grid, longitude_w_one_tile, rootpet=0, tile=tile, rc=error)
1140  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1141  call error_handler("IN FieldScatter", error)
1142  print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE. TILE IS: ", tile
1143  call esmf_fieldscatter(latitude_target_grid, latitude_one_tile, rootpet=0, tile=tile, rc=error)
1144  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1145  call error_handler("IN FieldScatter", error)
1146  print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE_S. TILE IS: ", tile
1147  call esmf_fieldscatter(latitude_s_target_grid, latitude_s_one_tile, rootpet=0, tile=tile, rc=error)
1148  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1149  call error_handler("IN FieldScatter", error)
1150  print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE_W. TILE IS: ", tile
1151  call esmf_fieldscatter(latitude_w_target_grid, latitude_w_one_tile, rootpet=0, tile=tile, rc=error)
1152  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1153  call error_handler("IN FieldScatter", error)
1154  print*,"- CALL FieldScatter FOR TARGET GRID TERRAIN. TILE IS: ", tile
1155  call esmf_fieldscatter(terrain_target_grid, terrain_one_tile, rootpet=0, tile=tile, rc=error)
1156  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1157  call error_handler("IN FieldScatter", error)
1158  enddo
1159 
1160  deallocate(land_frac_one_tile)
1161  deallocate(landmask_one_tile)
1162  deallocate(seamask_one_tile)
1163  deallocate(longitude_one_tile)
1164  deallocate(longitude_s_one_tile)
1165  deallocate(longitude_w_one_tile)
1166  deallocate(latitude_one_tile)
1167  deallocate(latitude_s_one_tile)
1168  deallocate(latitude_w_one_tile)
1169  deallocate(terrain_one_tile)
1170 
1171  end subroutine define_target_grid
1172 
1191  subroutine get_model_latlons(mosaic_file, orog_dir, num_tiles, tile, &
1192  i_tile, j_tile, ip1_tile, jp1_tile, &
1193  latitude, latitude_s, latitude_w, &
1194  longitude, longitude_s, longitude_w)
1195 
1196  use netcdf
1197 
1198  implicit none
1199 
1200  character(len=*), intent(in) :: mosaic_file, orog_dir
1201 
1202  integer, intent(in) :: num_tiles, tile
1203  integer, intent(in) :: i_tile, j_tile
1204  integer, intent(in) :: ip1_tile, jp1_tile
1205 
1206  real(esmf_kind_r8), intent(out) :: latitude(i_tile, j_tile)
1207  real(esmf_kind_r8), intent(out) :: latitude_s(i_tile, jp1_tile)
1208  real(esmf_kind_r8), intent(out) :: latitude_w(ip1_tile, j_tile)
1209  real(esmf_kind_r8), intent(out) :: longitude(i_tile, j_tile)
1210  real(esmf_kind_r8), intent(out) :: longitude_s(i_tile, jp1_tile)
1211  real(esmf_kind_r8), intent(out) :: longitude_w(ip1_tile, j_tile)
1212 
1213  character(len=50) :: grid_files(num_tiles)
1214  character(len=255) :: grid_file
1215 
1216  integer :: error, id_var, ncid
1217  integer :: id_dim, nxp, nyp, i, j, ii, jj
1218 
1219  real(esmf_kind_r8), allocatable :: tmpvar(:,:)
1220 
1221  print*,"- READ MODEL GRID FILE"
1222 
1223  print*,'- OPEN MOSAIC FILE: ', trim(mosaic_file)
1224  error=nf90_open(trim(mosaic_file), nf90_nowrite, ncid)
1225  call netcdf_err(error, 'opening mosaic file')
1226 
1227  print*,"- READ GRID FILE NAMES"
1228  error=nf90_inq_varid(ncid, 'gridfiles', id_var)
1229  call netcdf_err(error, 'reading gridfiles id')
1230  error=nf90_get_var(ncid, id_var, grid_files)
1231  call netcdf_err(error, 'reading gridfiles')
1232 
1233  error = nf90_close(ncid)
1234 
1235  grid_file = trim(orog_dir) // trim(grid_files(tile))
1236 
1237  print*,'- OPEN GRID FILE: ', trim(grid_file)
1238  error=nf90_open(trim(grid_file), nf90_nowrite, ncid)
1239  call netcdf_err(error, 'opening grid file')
1240 
1241  print*,'- READ NXP ID'
1242  error=nf90_inq_dimid(ncid, 'nxp', id_dim)
1243  call netcdf_err(error, 'reading nxp id')
1244 
1245  print*,'- READ NXP'
1246  error=nf90_inquire_dimension(ncid,id_dim,len=nxp)
1247  call netcdf_err(error, 'reading nxp')
1248 
1249  print*,'- READ NYP ID'
1250  error=nf90_inq_dimid(ncid, 'nyp', id_dim)
1251  call netcdf_err(error, 'reading nyp id')
1252 
1253  print*,'- READ NYP'
1254  error=nf90_inquire_dimension(ncid,id_dim,len=nyp)
1255  call netcdf_err(error, 'reading nyp')
1256 
1257  if ((nxp/2 /= i_tile) .or. (nyp/2 /= j_tile)) then
1258  call error_handler("DIMENSION MISMATCH IN GRID FILE.", 1)
1259  endif
1260 
1261  allocate(tmpvar(nxp,nyp))
1262 
1263  print*,'- READ LONGITUDE ID'
1264  error=nf90_inq_varid(ncid, 'x', id_var)
1265  call netcdf_err(error, 'reading longitude id')
1266 
1267  print*,'- READ LONGITUDE'
1268  error=nf90_get_var(ncid, id_var, tmpvar)
1269  call netcdf_err(error, 'reading longitude')
1270 
1271  do j = 1, j_tile
1272  do i = 1, i_tile
1273  ii = 2*i
1274  jj = 2*j
1275  longitude(i,j) = tmpvar(ii,jj)
1276  enddo
1277  enddo
1278 
1279  do j = 1, jp1_tile
1280  do i = 1, i_tile
1281  ii = 2*i
1282  jj = (2*j) - 1
1283  longitude_s(i,j) = tmpvar(ii,jj)
1284  enddo
1285  enddo
1286 
1287  do j = 1, j_tile
1288  do i = 1, ip1_tile
1289  ii = (2*i) - 1
1290  jj = 2*j
1291  longitude_w(i,j) = tmpvar(ii,jj)
1292  enddo
1293  enddo
1294 
1295  print*,'- READ LATITUDE ID'
1296  error=nf90_inq_varid(ncid, 'y', id_var)
1297  call netcdf_err(error, 'reading latitude id')
1298 
1299  print*,'- READ LATIITUDE'
1300  error=nf90_get_var(ncid, id_var, tmpvar)
1301  call netcdf_err(error, 'reading latitude')
1302 
1303  do j = 1, j_tile
1304  do i = 1, i_tile
1305  ii = 2*i
1306  jj = 2*j
1307  latitude(i,j) = tmpvar(ii,jj)
1308  enddo
1309  enddo
1310 
1311  do j = 1, jp1_tile
1312  do i = 1, i_tile
1313  ii = 2*i
1314  jj = (2*j) - 1
1315  latitude_s(i,j) = tmpvar(ii,jj)
1316  enddo
1317  enddo
1318 
1319  do j = 1, j_tile
1320  do i = 1, ip1_tile
1321  ii = (2*i) - 1
1322  jj = 2*j
1323  latitude_w(i,j) = tmpvar(ii,jj)
1324  enddo
1325  enddo
1326 
1327  deallocate(tmpvar)
1328 
1329  error = nf90_close(ncid)
1330 
1331  end subroutine get_model_latlons
1332 
1343  subroutine get_model_mask_terrain(orog_file, idim, jdim, mask, terrain, land_frac)
1344 
1345  use netcdf
1346 
1347  implicit none
1348 
1349  character(len=*), intent(in) :: orog_file
1350 
1351  integer, intent(in) :: idim, jdim
1352  integer(esmf_kind_i8), intent(out) :: mask(idim,jdim)
1353 
1354  real(esmf_kind_i8), intent(out) :: terrain(idim,jdim)
1355  real(esmf_kind_i8), intent(out) :: land_frac(idim,jdim)
1356 
1357  integer :: error, lat, lon
1358  integer :: ncid, id_dim, id_var
1359 
1360  real(kind=4), allocatable :: dummy(:,:)
1361 
1362  print*,"- READ MODEL LAND MASK FILE"
1363 
1364  print*,'- OPEN LAND MASK FILE: ', orog_file
1365  error=nf90_open(orog_file,nf90_nowrite,ncid)
1366  call netcdf_err(error, 'opening land mask file')
1367 
1368  print*,"- READ I-DIMENSION"
1369  error=nf90_inq_dimid(ncid, 'lon', id_dim)
1370  call netcdf_err(error, 'reading idim id')
1371  error=nf90_inquire_dimension(ncid,id_dim,len=lon)
1372  call netcdf_err(error, 'reading idim')
1373 
1374  print*,"- READ J-DIMENSION"
1375  error=nf90_inq_dimid(ncid, 'lat', id_dim)
1376  call netcdf_err(error, 'reading jdim id')
1377  error=nf90_inquire_dimension(ncid,id_dim,len=lat)
1378  call netcdf_err(error, 'reading jdim')
1379 
1380  print*,"- I/J DIMENSIONS: ", lon, lat
1381 
1382  if ((lon /= idim) .or. (lat /= jdim)) then
1383  call error_handler("MISMATCH IN DIMENSIONS.", 1)
1384  endif
1385 
1386  allocate(dummy(idim,jdim))
1387 
1388  print*,"- READ LAND MASK"
1389  error=nf90_inq_varid(ncid, 'slmsk', id_var)
1390  call netcdf_err(error, 'reading slmsk id')
1391  error=nf90_get_var(ncid, id_var, dummy)
1392  call netcdf_err(error, 'reading slmsk')
1393  mask = nint(dummy)
1394 
1395  print*,"- READ RAW OROGRAPHY."
1396  error=nf90_inq_varid(ncid, 'orog_raw', id_var)
1397  call netcdf_err(error, 'reading orog_raw id')
1398  error=nf90_get_var(ncid, id_var, dummy)
1399  call netcdf_err(error, 'reading orog_raw')
1400  terrain = dummy
1401 
1402  print*,"- READ LAND FRACTION."
1403  error=nf90_inq_varid(ncid, 'land_frac', id_var)
1404  call netcdf_err(error, 'reading land_frac id')
1405  error=nf90_get_var(ncid, id_var, dummy)
1406  call netcdf_err(error, 'reading orog_raw')
1407  land_frac = dummy
1408 
1409  error = nf90_close(ncid)
1410 
1411  deallocate (dummy)
1412 
1413  end subroutine get_model_mask_terrain
1414 
1418  subroutine cleanup_input_target_grid_data
1419 
1420  implicit none
1421 
1422  integer :: rc
1423 
1424  print*,"- DESTROY MODEL DATA."
1425 
1426  call esmf_fielddestroy(latitude_input_grid,rc=rc)
1427  call esmf_fielddestroy(longitude_input_grid,rc=rc)
1428  if (esmf_fieldiscreated(latitude_s_input_grid)) then
1429  call esmf_fielddestroy(latitude_s_input_grid, rc=rc)
1430  endif
1431  if (esmf_fieldiscreated(latitude_w_input_grid)) then
1432  call esmf_fielddestroy(latitude_w_input_grid, rc=rc)
1433  endif
1434  if (esmf_fieldiscreated(longitude_s_input_grid)) then
1435  call esmf_fielddestroy(longitude_s_input_grid, rc=rc)
1436  endif
1437  if (esmf_fieldiscreated(longitude_w_input_grid)) then
1438  call esmf_fielddestroy(longitude_w_input_grid, rc=rc)
1439  endif
1440  call esmf_fielddestroy(land_frac_target_grid, rc=rc)
1441  call esmf_fielddestroy(landmask_target_grid, rc=rc)
1442  call esmf_fielddestroy(latitude_target_grid, rc=rc)
1443  if (esmf_fieldiscreated(latitude_s_target_grid)) then
1444  call esmf_fielddestroy(latitude_s_target_grid, rc=rc)
1445  endif
1446  if (esmf_fieldiscreated(latitude_w_target_grid)) then
1447  call esmf_fielddestroy(latitude_w_target_grid, rc=rc)
1448  endif
1449  call esmf_fielddestroy(longitude_target_grid, rc=rc)
1450  if (esmf_fieldiscreated(longitude_s_target_grid)) then
1451  call esmf_fielddestroy(longitude_s_target_grid, rc=rc)
1452  endif
1453  if (esmf_fieldiscreated(longitude_w_target_grid)) then
1454  call esmf_fielddestroy(longitude_w_target_grid, rc=rc)
1455  endif
1456  call esmf_fielddestroy(seamask_target_grid, rc=rc)
1457  call esmf_fielddestroy(terrain_target_grid, rc=rc)
1458  call esmf_griddestroy(input_grid, rc=rc)
1459  call esmf_griddestroy(target_grid, rc=rc)
1460 
1461  end subroutine cleanup_input_target_grid_data
1462 
1474  subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
1475 
1476  implicit none
1477 
1478  integer, intent(in ) :: igdtnum, igdtlen, igdstmpl(igdtlen)
1479  integer, intent( out) :: kgds(200), ni, nj
1480  integer :: iscale, i
1481 
1482  real, intent( out) :: res
1483 
1484  real :: clatr, slatr, clonr, dpr, slat
1485  real :: slat0, clat0, clat, clon, rlat
1486  real :: rlon0, rlon, hs
1487 
1488  kgds=0
1489 
1490  if (igdtnum.eq.32769) then ! rot lat/lon b grid
1491 
1492  iscale=igdstmpl(10)*igdstmpl(11)
1493  if (iscale == 0) iscale = 1e6
1494  kgds(1)=205 ! oct 6, rotated lat/lon for Non-E
1495  ! Stagger grid
1496  kgds(2)=igdstmpl(8) ! octs 7-8, Ni
1497  ni = kgds(2)
1498  kgds(3)=igdstmpl(9) ! octs 9-10, Nj
1499  nj = kgds(3)
1500  kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of
1501  ! 1st grid point
1502  kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of
1503  ! 1st grid point
1504 
1505  kgds(6)=0 ! oct 17, resolution and component flags
1506  if (igdstmpl(1)==2 ) kgds(6)=64
1507  if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
1508  if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
1509 
1510  kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20,
1511  ! Lat of cent of rotation
1512  kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23,
1513  ! Lon of cent of rotation
1514  kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25,
1515  ! Di
1516  kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27,
1517  ! Dj
1518 
1519  kgds(11) = 0 ! oct 28, scan mode
1520  if (btest(igdstmpl(19),7)) kgds(11) = 128
1521  if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
1522  if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
1523 
1524  kgds(12)=nint(float(igdstmpl(20))/float(iscale)*1000.) ! octs 29-31, Lat of
1525  ! last grid point
1526  kgds(13)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 32-34, Lon of
1527  ! last grid point
1528 
1529  kgds(19)=0 ! oct 4, # vert coordinate parameters
1530  kgds(20)=255 ! oct 5, used for thinned grids, set to 255
1531 
1532  res = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) &
1533  * 0.5 * 111.0
1534 
1535  elseif (igdtnum.eq.1) then ! rot lat/lon b grid using standard wmo
1536  ! template.
1537 
1538  iscale=igdstmpl(10)*igdstmpl(11)
1539  if (iscale == 0) iscale = 1e6
1540  kgds(1)=205 ! oct 6, rotated lat/lon for Non-E
1541  ! Stagger grid
1542  kgds(2)=igdstmpl(8) ! octs 7-8, Ni
1543  ni = kgds(2)
1544  kgds(3)=igdstmpl(9) ! octs 9-10, Nj
1545  nj = kgds(3)
1546 
1547  kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of
1548  ! 1st grid point
1549  kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of
1550  ! 1st grid point
1551 
1552  kgds(6)=0 ! oct 17, resolution and component flags
1553  if (igdstmpl(1)==2 ) kgds(6)=64
1554  if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
1555  if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
1556 
1557  kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.) ! octs 18-20,
1558  ! Lat of cent of rotation
1559  kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 21-23,
1560  ! Lon of cent of rotation
1561  kgds(7) = kgds(7) + 90000.0
1562  print*, "INPUT LAT, LON CENTER ", kgds(7), kgds(8)
1563 
1564  dpr = 180.0/3.1415926
1565  clatr=cos((float(kgds(4))/1000.0)/dpr)
1566  slatr=sin((float(kgds(4))/1000.0)/dpr)
1567  clonr=cos((float(kgds(5))/1000.0)/dpr)
1568  slat0=sin((float(kgds(7))/1000.0)/dpr)
1569  clat0=cos((float(kgds(7))/1000.0)/dpr)
1570 
1571  slat=clat0*slatr+slat0*clatr*clonr
1572  clat=sqrt(1-slat**2)
1573  clon=(clat0*clatr*clonr-slat0*slatr)/clat
1574  clon=min(max(clon,-1.0),1.0)
1575 
1576  rlat=dpr*asin(slat)
1577 
1578  rlon0=float(kgds(8))/1000.0
1579 
1580  if ((kgds(5)-kgds(8)) > 0) then
1581  hs = -1.0
1582  else
1583  hs = 1.0
1584  endif
1585 
1586  rlon = mod(rlon0+hs*dpr*acos(clon)+3600,360.0)
1587 
1588  kgds(4)=nint(rlat*1000.) ! octs 11-13, Lat of
1589  kgds(5)=nint(rlon*1000.) ! octs 14-16, Lon of
1590 
1591  kgds(12)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 29-31, Lat of
1592  ! last grid point
1593  kgds(13)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 32-34, Lon of
1594  ! last grid point
1595 
1596  clatr=cos((float(kgds(12))/1000.0)/dpr)
1597  slatr=sin((float(kgds(12))/1000.0)/dpr)
1598  clonr=cos((float(kgds(13))/1000.0)/dpr)
1599 
1600  slat=clat0*slatr+slat0*clatr*clonr
1601  rlat=dpr*asin(slat)
1602 
1603  clat=sqrt(1-slat**2)
1604  clon=(clat0*clatr*clonr-slat0*slatr)/clat
1605  clon=min(max(clon,-1.0),1.0)
1606 
1607  if ((kgds(13)-kgds(8)) > 0) then
1608  hs = -1.0
1609  else
1610  hs = 1.0
1611  endif
1612 
1613  rlon = mod(rlon0+hs*dpr*acos(clon)+3600,360.0)
1614 
1615  print*,'got here last point ',kgds(12), kgds(13)
1616  print*,'got here last point rotated ', rlat, rlon
1617 
1618  kgds(12)=nint(rlat*1000.) ! octs 11-13, Lat of
1619  kgds(13)=nint(rlon*1000.) ! octs 14-16, Lon of
1620 
1621  kgds(9)=igdstmpl(17)
1622  kgds(10)=igdstmpl(18)
1623 
1624  kgds(11) = 0 ! oct 28, scan mode
1625  if (btest(igdstmpl(19),7)) kgds(11) = 128
1626  if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
1627  if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
1628 
1629  kgds(19)=0 ! oct 4, # vert coordinate parameters
1630  kgds(20)=255 ! oct 5, used for thinned grids, set to 255
1631 
1632  res = ((float(kgds(9)) / 1.e6) + (float(kgds(10)) / 1.e6)) &
1633  * 0.5 * 111.0
1634 
1635 
1636  do i = 1, 25
1637  print*,'final kgds ',i,kgds(i)
1638  enddo
1639 
1640 
1641  elseif(igdtnum==30) then
1642 
1643  kgds(1)=3 ! oct 6, lambert conformal
1644  kgds(2)=igdstmpl(8) ! octs 7-8, Ni
1645  ni = kgds(2)
1646  kgds(3)=igdstmpl(9) ! octs 9-10, Nj
1647  nj = kgds(3)
1648 
1649  iscale = 1e6
1650  kgds(4) = nint(float(igdstmpl(10))/1000.0)
1651  kgds(5) = nint(float(igdstmpl(11))/1000.0)
1652 
1653  kgds(6)=0 ! oct 17, resolution and component flags
1654  if (igdstmpl(1)==2 ) kgds(6)=64
1655  if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) ) kgds(6)=kgds(6)+128
1656  if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8
1657 
1658  kgds(7) = nint(float(igdstmpl(14))/1000.0)
1659  kgds(8) = nint(float(igdstmpl(15))/1000.0)
1660  kgds(9) = nint(float(igdstmpl(16))/1000.0)
1661  kgds(10) = 0
1662 
1663  kgds(11) = 0 ! oct 28, scan mode
1664  if (btest(igdstmpl(18),7)) kgds(11) = 128
1665  if (btest(igdstmpl(18),6)) kgds(11) = kgds(11) + 64
1666  if (btest(igdstmpl(18),5)) kgds(11) = kgds(11) + 32
1667 
1668  kgds(12) = nint(float(igdstmpl(19))/1000.0)
1669  kgds(13) = nint(float(igdstmpl(20))/1000.0)
1670  kgds(14) = -90
1671  kgds(15) = 0
1672 
1673  elseif(igdtnum==0) then ! lat/lon grid
1674 
1675  iscale=igdstmpl(10)*igdstmpl(11)
1676  if (iscale == 0) iscale = 1e6
1677  kgds(1)=0 ! oct 6, data representation type.
1678  kgds(2)=igdstmpl(8) ! octs 7-8, Ni
1679  ni = kgds(2)
1680  kgds(3)=igdstmpl(9) ! octs 9-10, Nj
1681  nj = kgds(3)
1682  kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point
1683  kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point
1684 
1685  kgds(6)=0 ! oct 17, resolution and component flags
1686  if (igdstmpl(1)==2 ) kgds(6)=64
1687  if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
1688  if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
1689 
1690  kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point
1691  kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point
1692  kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, "i" resolution.
1693  kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, "j" resolution.
1694 
1695  kgds(11) = 0 ! oct 28, scan mode
1696  if (btest(igdstmpl(19),7)) kgds(11) = 128
1697  if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
1698  if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
1699 
1700  kgds(12)=0 ! octs 29-32, reserved
1701  kgds(19)=0 ! oct 4, # vert coordinate parameters
1702  kgds(20)=255 ! oct 5, used for thinned grids, set to 255
1703 
1704  else
1705 
1706  call error_handler("UNRECOGNIZED INPUT GRID TYPE ", 1)
1707 
1708  endif
1709 
1710  end subroutine gdt_to_gds
1711 
1712  end module model_grid
model_grid::get_model_mask_terrain
subroutine get_model_mask_terrain(orog_file, idim, jdim, mask, terrain, land_frac)
Read the model land mask and terrain for a single tile from the orography file.
Definition: model_grid.F90:1343
model_grid::cleanup_input_target_grid_data
subroutine, public cleanup_input_target_grid_data
Deallocate all esmf grid objects.
Definition: model_grid.F90:1418
model_grid::define_input_grid_grib2
subroutine define_input_grid_grib2(npets)
Define input grid object for grib2 input data.
Definition: model_grid.F90:610
model_grid::gdt_to_gds
subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
Convert the GRIB2 grid description template to to the GRIB1 grid description section.
Definition: model_grid.F90:1474
model_grid::define_input_grid_gaussian
subroutine define_input_grid_gaussian(npets)
Define grid object for input data on global gaussian grids.
Definition: model_grid.F90:148
utilities::netcdf_err
subroutine netcdf_err(err, string)
Error handler for netcdf.
Definition: utils.F90:34
model_grid
Sets up the ESMF grid objects for the input data grid and target FV3 grid.
Definition: model_grid.F90:9
model_grid::get_model_latlons
subroutine get_model_latlons(mosaic_file, orog_dir, num_tiles, tile, i_tile, j_tile, ip1_tile, jp1_tile, latitude, latitude_s, latitude_w, longitude, longitude_s, longitude_w)
Read model lat/lons for a single tile from the "grid" specificaton file.
Definition: model_grid.F90:1191
utilities::error_handler
subroutine error_handler(string, rc)
General error handler.
Definition: utils.F90:12
model_grid::define_input_grid_mosaic
subroutine define_input_grid_mosaic(localpet, npets)
Define input grid for tiled data using the 'mosaic', 'grid' and orography files.
Definition: model_grid.F90:402
program_setup
This module contains code to read the setup namelist file, handle the varmap file for GRIB2 data...
Definition: program_setup.F90:9
model_grid::define_input_grid
subroutine, public define_input_grid(localpet, npets)
Driver routine to setup the esmf grid object for the input grid.
Definition: model_grid.F90:117
model_grid::define_target_grid
subroutine, public define_target_grid(localpet, npets)
Setup the esmf grid object for the target grid.
Definition: model_grid.F90:883