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