chgres_cube  1.4.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  implicit none
15 
16  private
17 
18  character(len=5), allocatable, public :: tiles_target_grid(:)
19 
20  character(len=10), public :: inv_file = "chgres.inv"
21 
22  character(len=50), public :: input_grid_type = "latlon"
23 
24 
25  ! Made lsoil_target non-parameter to allow for RAP land surface initiation
26  integer, public :: lsoil_target = 4 ! # soil layers
27 
28  integer, public :: i_input
29 
31  integer, public :: j_input
32 
34  integer, public :: ip1_input
35 
36  integer, public :: jp1_input
37 
38  integer, public :: i_target
39 
41  integer, public :: j_target
42 
44  integer, public :: ip1_target
45 
46  integer, public :: jp1_target
47 
48  integer, public :: num_tiles_input_grid
49 
50  integer, public :: num_tiles_target_grid
51 
52 
53  type(esmf_grid), public :: input_grid
54 
55  type(esmf_grid), public :: target_grid
56 
57 
58  type(esmf_field), public :: latitude_input_grid
59 
60  type(esmf_field), public :: longitude_input_grid
61 
62  type(esmf_field), public :: latitude_s_input_grid
63 
65  type(esmf_field), public :: longitude_s_input_grid
66 
68  type(esmf_field), public :: latitude_w_input_grid
69 
71  type(esmf_field), public :: longitude_w_input_grid
72 
74 
75  type(esmf_field), public :: landmask_target_grid
76 
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, external_model
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(localpet, npets)
130  elseif (trim(external_model) == "GFS" .and. trim(input_type) == "grib2") then
131  call define_input_grid_gfs_grib2(localpet,npets)
132  elseif (trim(input_type) == "grib2") then
133  call define_input_grid_grib2(localpet,npets)
134  else
135  call define_input_grid_mosaic(localpet, npets)
136  endif
137 
138  end subroutine define_input_grid
139 
151  subroutine define_input_grid_gaussian(localpet, npets)
152 
153  use nemsio_module
154 
155  use program_setup, only : data_dir_input_grid, &
156  atm_files_input_grid, &
157  sfc_files_input_grid, &
158  input_type, &
159  convert_atm, convert_sfc
160 
161  use sfcio_module
162  use sigio_module
163  use netcdf
164 
165  implicit none
166 
167  integer, intent(in) :: localpet, npets
168 
169  character(len=250) :: the_file
170 
171  integer :: i, j, rc, clb(2), cub(2), ncid, id_grid
172  integer(sfcio_intkind) :: rc2
173  integer(sigio_intkind) :: rc3
174 
175  real(esmf_kind_r8), allocatable :: latitude(:,:)
176  real(esmf_kind_r8), allocatable :: longitude(:,:)
177  real(esmf_kind_r8), pointer :: lat_src_ptr(:,:)
178  real(esmf_kind_r8), pointer :: lon_src_ptr(:,:)
179  real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:)
180  real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:)
181  real(esmf_kind_r8) :: deltalon
182  real(esmf_kind_r8), allocatable :: slat(:), wlat(:)
183 
184  type(nemsio_gfile) :: gfile
185  type(esmf_polekind_flag) :: polekindflag(2)
186  type(sfcio_head) :: sfchead
187  type(sigio_head) :: sighead
188 
189  print*,"- DEFINE INPUT GRID OBJECT FOR GAUSSIAN DATA."
190 
191  num_tiles_input_grid = 1
192 
193  if (convert_sfc) then
194  the_file=trim(data_dir_input_grid) // "/" // trim(sfc_files_input_grid(1))
195  elseif (convert_atm) then
196  the_file=trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
197  endif
198 
199  if (trim(input_type) == "gfs_sigio") then ! sigio/sfcio format, used by
200  ! spectral gfs prior to 7/19/2017.
201 
202  if (convert_sfc) then ! sfcio format
203  print*,"- OPEN AND READ ", trim(the_file)
204  call sfcio_sropen(21, trim(the_file), rc2)
205  if (rc2 /= 0) call error_handler("OPENING FILE", rc2)
206  call sfcio_srhead(21, sfchead, rc2)
207  if (rc2 /= 0) call error_handler("READING FILE", rc2)
208  call sfcio_sclose(21, rc2)
209  i_input = sfchead%lonb
210  j_input = sfchead%latb
211  elseif (convert_atm) then ! sigio format
212  print*,"- OPEN AND READ ", trim(the_file)
213  call sigio_sropen(21, trim(the_file), rc3)
214  if (rc3 /= 0) call error_handler("OPENING FILE", rc3)
215  call sigio_srhead(21, sighead, rc3)
216  if (rc3 /= 0) call error_handler("READING FILE", rc3)
217  call sigio_sclose(21, rc3)
218  i_input = sighead%lonb
219  j_input = sighead%latb
220  endif
221 
222  elseif (trim(input_type) == "gaussian_netcdf") then
223 
224  print*,'- OPEN AND READ: ',trim(the_file)
225  rc=nf90_open(trim(the_file),nf90_nowrite,ncid)
226  call netcdf_err(rc, 'opening file')
227 
228  print*,"- READ grid_xt"
229  rc=nf90_inq_dimid(ncid, 'grid_xt', id_grid)
230  call netcdf_err(rc, 'reading grid_xt id')
231  rc=nf90_inquire_dimension(ncid,id_grid,len=i_input)
232  call netcdf_err(rc, 'reading grid_xt')
233 
234  print*,"- READ grid_yt"
235  rc=nf90_inq_dimid(ncid, 'grid_yt', id_grid)
236  call netcdf_err(rc, 'reading grid_yt id')
237  rc=nf90_inquire_dimension(ncid,id_grid,len=j_input)
238  call netcdf_err(rc, 'reading grid_yt')
239 
240  rc = nf90_close(ncid)
241 
242  else ! nemsio format
243 
244  call nemsio_init(iret=rc)
245 
246  print*,"- OPEN AND READ ", trim(the_file)
247  call nemsio_open(gfile, the_file, "read", iret=rc)
248  if (rc /= 0) call error_handler("OPENING FILE", rc)
249 
250  call nemsio_getfilehead(gfile, iret=rc, dimx=i_input, dimy=j_input)
251  if (rc /= 0) call error_handler("READING FILE", rc)
252 
253  call nemsio_close(gfile)
254 
255  endif
256 
257  ip1_input = i_input + 1
258  jp1_input = j_input + 1
259 
260  polekindflag(1:2) = esmf_polekind_monopole
261 
262  print*,"- CALL GridCreate1PeriDim FOR INPUT GRID."
263  input_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
264  maxindex=(/i_input,j_input/), &
265  polekindflag=polekindflag, &
266  periodicdim=1, &
267  poledim=2, &
268  coordsys=esmf_coordsys_sph_deg, &
269  regdecomp=(/1,npets/), &
270  indexflag=esmf_index_global, rc=rc)
271  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
272  call error_handler("IN GridCreate1PeriDim", rc)
273 
274  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
275  latitude_input_grid = esmf_fieldcreate(input_grid, &
276  typekind=esmf_typekind_r8, &
277  staggerloc=esmf_staggerloc_center, &
278  name="input_grid_latitude", rc=rc)
279  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
280  call error_handler("IN FieldCreate", rc)
281 
282  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
283  longitude_input_grid = esmf_fieldcreate(input_grid, &
284  typekind=esmf_typekind_r8, &
285  staggerloc=esmf_staggerloc_center, &
286  name="input_grid_longitude", rc=rc)
287  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
288  call error_handler("IN FieldCreate", rc)
289 
290  allocate(longitude(i_input,j_input))
291  allocate(latitude(i_input,j_input))
292 
293  deltalon = 360.0_esmf_kind_r8 / real(i_input,kind=esmf_kind_r8)
294  do i = 1, i_input
295  longitude(i,:) = real((i-1),kind=esmf_kind_r8) * deltalon
296  enddo
297 
298  allocate(slat(j_input))
299  allocate(wlat(j_input))
300  call splat(4, j_input, slat, wlat)
301 
302  do i = 1, j_input
303  latitude(:,i) = 90.0_esmf_kind_r8 - (acos(slat(i))* 180.0_esmf_kind_r8 / &
304  (4.0_esmf_kind_r8*atan(1.0_esmf_kind_r8)))
305  enddo
306 
307  deallocate(slat, wlat)
308 
309  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
310  call esmf_fieldscatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
311  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
312  call error_handler("IN FieldScatter", rc)
313 
314  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE."
315  call esmf_fieldscatter(latitude_input_grid, latitude, rootpet=0, rc=rc)
316  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
317  call error_handler("IN FieldScatter", rc)
318 
319  print*,"- CALL GridAddCoord FOR INPUT GRID."
320  call esmf_gridaddcoord(input_grid, &
321  staggerloc=esmf_staggerloc_center, rc=rc)
322  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
323  call error_handler("IN GridAddCoord", rc)
324 
325  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
326  nullify(lon_src_ptr)
327  call esmf_gridgetcoord(input_grid, &
328  staggerloc=esmf_staggerloc_center, &
329  coorddim=1, &
330  farrayptr=lon_src_ptr, rc=rc)
331  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
332  call error_handler("IN GridGetCoord", rc)
333 
334  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
335  nullify(lat_src_ptr)
336  call esmf_gridgetcoord(input_grid, &
337  staggerloc=esmf_staggerloc_center, &
338  coorddim=2, &
339  computationallbound=clb, &
340  computationalubound=cub, &
341  farrayptr=lat_src_ptr, rc=rc)
342  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
343  call error_handler("IN GridGetCoord", rc)
344 
345  do j = clb(2), cub(2)
346  do i = clb(1), cub(1)
347  lon_src_ptr(i,j) = longitude(i,j)
348  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
349  lat_src_ptr(i,j) = latitude(i,j)
350  enddo
351  enddo
352 
353  print*,"- CALL GridAddCoord FOR INPUT GRID."
354  call esmf_gridaddcoord(input_grid, &
355  staggerloc=esmf_staggerloc_corner, rc=rc)
356  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
357  call error_handler("IN GridAddCoord", rc)
358 
359  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
360  nullify(lon_corner_src_ptr)
361  call esmf_gridgetcoord(input_grid, &
362  staggerloc=esmf_staggerloc_corner, &
363  coorddim=1, &
364  farrayptr=lon_corner_src_ptr, rc=rc)
365  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
366  call error_handler("IN GridGetCoord", rc)
367 
368  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
369  nullify(lat_corner_src_ptr)
370  call esmf_gridgetcoord(input_grid, &
371  staggerloc=esmf_staggerloc_corner, &
372  coorddim=2, &
373  computationallbound=clb, &
374  computationalubound=cub, &
375  farrayptr=lat_corner_src_ptr, rc=rc)
376  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
377  call error_handler("IN GridGetCoord", rc)
378 
379  do j = clb(2), cub(2)
380  do i = clb(1), cub(1)
381  lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon)
382  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
383  if (j == 1) then
384  lat_corner_src_ptr(i,j) = 90.0_esmf_kind_r8
385  cycle
386  endif
387  if (j == jp1_input) then
388  lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8
389  cycle
390  endif
391  lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j))
392  enddo
393  enddo
394 
395  deallocate(latitude,longitude)
396 
397  end subroutine define_input_grid_gaussian
398 
405  subroutine define_input_grid_mosaic(localpet, npets)
406 
407  use netcdf
408  use program_setup, only : mosaic_file_input_grid, &
409  orog_dir_input_grid, &
410  orog_files_input_grid
411 
412  implicit none
413 
414  character(len=500) :: the_file
415 
416  integer, intent(in) :: localpet, npets
417 
418  integer :: id_tiles, id_dim, tile
419  integer :: extra, error, ncid
420  integer, allocatable :: decomptile(:,:)
421 
422  integer(esmf_kind_i8), allocatable :: landmask_one_tile(:,:)
423 
424  real(esmf_kind_r8), allocatable :: latitude_one_tile(:,:)
425  real(esmf_kind_r8), allocatable :: latitude_s_one_tile(:,:)
426  real(esmf_kind_r8), allocatable :: latitude_w_one_tile(:,:)
427  real(esmf_kind_r8), allocatable :: longitude_one_tile(:,:)
428  real(esmf_kind_r8), allocatable :: longitude_s_one_tile(:,:)
429  real(esmf_kind_r8), allocatable :: longitude_w_one_tile(:,:)
430 
431  print*,'- OPEN INPUT GRID MOSAIC FILE: ',trim(mosaic_file_input_grid)
432  error=nf90_open(trim(mosaic_file_input_grid),nf90_nowrite,ncid)
433  call netcdf_err(error, 'opening grid mosaic file')
434 
435  print*,"- READ NUMBER OF TILES"
436  error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
437  call netcdf_err(error, 'reading ntiles id')
438  error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles_input_grid)
439  call netcdf_err(error, 'reading ntiles')
440 
441  error = nf90_close(ncid)
442 
443  print*,'- NUMBER OF TILES, INPUT MODEL GRID IS ', num_tiles_input_grid
444 
445  if (mod(npets,num_tiles_input_grid) /= 0) then
446  call error_handler("MUST RUN WITH A TASK COUNT THAT IS A MULTIPLE OF 6.", 1)
447  endif
448 
449 !-----------------------------------------------------------------------
450 ! Create ESMF grid object for the model grid.
451 !-----------------------------------------------------------------------
452 
453  extra = npets / num_tiles_input_grid
454 
455  allocate(decomptile(2,num_tiles_input_grid))
456 
457  do tile = 1, num_tiles_input_grid
458  decomptile(:,tile)=(/1,extra/)
459  enddo
460 
461  print*,"- CALL GridCreateMosaic FOR INPUT MODEL GRID"
462  input_grid = esmf_gridcreatemosaic(filename=trim(mosaic_file_input_grid), &
463  regdecompptile=decomptile, &
464  staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
465  esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
466  indexflag=esmf_index_global, &
467  tilefilepath=trim(orog_dir_input_grid), &
468  rc=error)
469  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
470  call error_handler("IN GridCreateMosaic", error)
471 
472 !-----------------------------------------------------------------------
473 ! Read the mask and lat/lons.
474 !-----------------------------------------------------------------------
475 
476  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
477  latitude_input_grid = esmf_fieldcreate(input_grid, &
478  typekind=esmf_typekind_r8, &
479  staggerloc=esmf_staggerloc_center, &
480  name="input_grid_latitude", &
481  rc=error)
482  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
483  call error_handler("IN FieldCreate", error)
484 
485  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
486  longitude_input_grid = esmf_fieldcreate(input_grid, &
487  typekind=esmf_typekind_r8, &
488  staggerloc=esmf_staggerloc_center, &
489  name="input_grid_longitude", &
490  rc=error)
491  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
492  call error_handler("IN FieldCreate", error)
493 
494  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE_S."
495  latitude_s_input_grid = esmf_fieldcreate(input_grid, &
496  typekind=esmf_typekind_r8, &
497  staggerloc=esmf_staggerloc_edge2, &
498  name="input_grid_latitude_s", &
499  rc=error)
500  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
501  call error_handler("IN FieldCreate", error)
502 
503  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE_S."
504  longitude_s_input_grid = esmf_fieldcreate(input_grid, &
505  typekind=esmf_typekind_r8, &
506  staggerloc=esmf_staggerloc_edge2, &
507  name="input_grid_longitude_s", &
508  rc=error)
509  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
510  call error_handler("IN FieldCreate", error)
511 
512  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE_W."
513  latitude_w_input_grid = esmf_fieldcreate(input_grid, &
514  typekind=esmf_typekind_r8, &
515  staggerloc=esmf_staggerloc_edge1, &
516  name="input_grid_latitude_w", &
517  rc=error)
518  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
519  call error_handler("IN FieldCreate", error)
520 
521  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE_W."
522  longitude_w_input_grid = esmf_fieldcreate(input_grid, &
523  typekind=esmf_typekind_r8, &
524  staggerloc=esmf_staggerloc_edge1, &
525  name="input_grid_longitude_w", &
526  rc=error)
527  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
528  call error_handler("IN FieldCreate", error)
529 
530  the_file = trim(orog_dir_input_grid) // trim(orog_files_input_grid(1))
531 
532  print*,'- OPEN FIRST INPUT GRID OROGRAPHY FILE: ',trim(the_file)
533  error=nf90_open(trim(the_file),nf90_nowrite,ncid)
534  call netcdf_err(error, 'opening ororgraphy file')
535  print*,"- READ GRID DIMENSIONS"
536  error=nf90_inq_dimid(ncid, 'lon', id_dim)
537  call netcdf_err(error, 'reading lon id')
538  error=nf90_inquire_dimension(ncid,id_dim,len=i_input)
539  call netcdf_err(error, 'reading lon')
540  error=nf90_inq_dimid(ncid, 'lat', id_dim)
541  call netcdf_err(error, 'reading lat id')
542  error=nf90_inquire_dimension(ncid,id_dim,len=j_input)
543  call netcdf_err(error, 'reading lat')
544  error = nf90_close(ncid)
545 
546  print*,"- I/J DIMENSIONS OF THE INPUT GRID TILES ", i_input, j_input
547 
548  ip1_input = i_input + 1
549  jp1_input = j_input + 1
550 
551  if (localpet == 0) then
552  allocate(longitude_one_tile(i_input,j_input))
553  allocate(longitude_s_one_tile(i_input,jp1_input))
554  allocate(longitude_w_one_tile(ip1_input,j_input))
555  allocate(latitude_one_tile(i_input,j_input))
556  allocate(latitude_s_one_tile(i_input,jp1_input))
557  allocate(latitude_w_one_tile(ip1_input,j_input))
558  allocate(landmask_one_tile(i_input,j_input))
559  else
560  allocate(longitude_one_tile(0,0))
561  allocate(longitude_s_one_tile(0,0))
562  allocate(longitude_w_one_tile(0,0))
563  allocate(latitude_one_tile(0,0))
564  allocate(latitude_s_one_tile(0,0))
565  allocate(latitude_w_one_tile(0,0))
566  allocate(landmask_one_tile(0,0))
567  endif
568 
569  do tile = 1, num_tiles_input_grid
570  if (localpet == 0) then
571  call get_model_latlons(mosaic_file_input_grid, orog_dir_input_grid, num_tiles_input_grid, tile, &
572  i_input, j_input, ip1_input, jp1_input, latitude_one_tile, &
573  latitude_s_one_tile, latitude_w_one_tile, longitude_one_tile, &
574  longitude_s_one_tile, longitude_w_one_tile)
575  endif
576  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE. TILE IS: ", tile
577  call esmf_fieldscatter(latitude_input_grid, latitude_one_tile, rootpet=0, tile=tile, rc=error)
578  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
579  call error_handler("IN FieldScatter", error)
580  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE. TILE IS: ", tile
581  call esmf_fieldscatter(longitude_input_grid, longitude_one_tile, rootpet=0, tile=tile, rc=error)
582  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
583  call error_handler("IN FieldScatter", error)
584  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE_S. TILE IS: ", tile
585  call esmf_fieldscatter(latitude_s_input_grid, latitude_s_one_tile, rootpet=0, tile=tile, rc=error)
586  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
587  call error_handler("IN FieldScatter", error)
588  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE_S. TILE IS: ", tile
589  call esmf_fieldscatter(longitude_s_input_grid, longitude_s_one_tile, rootpet=0, tile=tile, rc=error)
590  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
591  call error_handler("IN FieldScatter", error)
592  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE_W. TILE IS: ", tile
593  call esmf_fieldscatter(latitude_w_input_grid, latitude_w_one_tile, rootpet=0, tile=tile, rc=error)
594  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
595  call error_handler("IN FieldScatter", error)
596  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE_W. TILE IS: ", tile
597  call esmf_fieldscatter(longitude_w_input_grid, longitude_w_one_tile, rootpet=0, tile=tile, rc=error)
598  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
599  call error_handler("IN FieldScatter", error)
600  enddo
601 
602  deallocate(longitude_one_tile)
603  deallocate(longitude_s_one_tile)
604  deallocate(longitude_w_one_tile)
605  deallocate(latitude_one_tile)
606  deallocate(latitude_s_one_tile)
607  deallocate(latitude_w_one_tile)
608  deallocate(landmask_one_tile)
609 
610  end subroutine define_input_grid_mosaic
611 
618  subroutine define_input_grid_gfs_grib2(localpet, npets)
619 
620  use wgrib2api
621  use mpi
622  use program_setup, only : data_dir_input_grid, &
623  grib2_file_input_grid
624 
625  implicit none
626 
627  integer, intent(in) :: localpet, npets
628 
629  character(len=250) :: the_file
630 
631  integer :: i, j, rc, clb(2), cub(2)
632  integer :: ierr
633 
634  real(esmf_kind_r8), allocatable :: latitude(:,:)
635  real(esmf_kind_r8), allocatable :: longitude(:,:)
636  real(esmf_kind_r4), allocatable :: lat4(:,:), lon4(:,:)
637  real(esmf_kind_r8), pointer :: lat_src_ptr(:,:)
638  real(esmf_kind_r8), pointer :: lon_src_ptr(:,:)
639  real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:)
640  real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:)
641  real(esmf_kind_r8) :: deltalon
642 
643  type(esmf_polekind_flag) :: polekindflag(2)
644 
645  print*,"- DEFINE INPUT GRID OBJECT FOR INPUT GRIB2 DATA."
646 
647  num_tiles_input_grid = 1
648 
649  the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid
650  if (localpet == 0) then
651  print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file)
652  rc=grb2_mk_inv(the_file,inv_file)
653  if (rc /=0) call error_handler("OPENING GRIB2 FILE",rc)
654  endif
655 
656 ! Wait for localpet 0 to create inventory
657  call mpi_barrier(mpi_comm_world, ierr)
658 
659  rc = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, &
660  lat=lat4, lon=lon4)
661  if (rc /= 1) call error_handler("READING GRIB2 FILE", rc)
662 
663  ip1_input = i_input + 1
664  jp1_input = j_input + 1
665 
666  polekindflag(1:2) = esmf_polekind_monopole
667 
668  print*,"- CALL GridCreate1PeriDim FOR INPUT GRID."
669  input_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
670  maxindex=(/i_input,j_input/), &
671  polekindflag=polekindflag, &
672  periodicdim=1, &
673  poledim=2, &
674  coordsys=esmf_coordsys_sph_deg, &
675  regdecomp=(/1,npets/), &
676  indexflag=esmf_index_global, rc=rc)
677  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
678  call error_handler("IN GridCreate1PeriDim", rc)
679 
680  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
681  latitude_input_grid = esmf_fieldcreate(input_grid, &
682  typekind=esmf_typekind_r8, &
683  staggerloc=esmf_staggerloc_center, &
684  name="input_grid_latitude", rc=rc)
685  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
686  call error_handler("IN FieldCreate", rc)
687 
688  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
689  longitude_input_grid = esmf_fieldcreate(input_grid, &
690  typekind=esmf_typekind_r8, &
691  staggerloc=esmf_staggerloc_center, &
692  name="input_grid_longitude", rc=rc)
693  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
694  call error_handler("IN FieldCreate", rc)
695 
696  allocate(longitude(i_input,j_input))
697  allocate(latitude(i_input,j_input))
698 
699  do i = 1, i_input
700  longitude(i,:) = real(lon4(i,:),kind=esmf_kind_r8)
701  enddo
702 
703  do i = 1, j_input
704  latitude(:,i) = real(lat4(:,i),kind=esmf_kind_r8)
705  enddo
706 
707  deallocate(lat4, lon4)
708 
709  deltalon = abs(longitude(2,1)-longitude(1,1))
710 
711  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
712  call esmf_fieldscatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
713  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
714  call error_handler("IN FieldScatter", rc)
715 
716  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE."
717  call esmf_fieldscatter(latitude_input_grid, latitude, rootpet=0, rc=rc)
718  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
719  call error_handler("IN FieldScatter", rc)
720 
721  print*,"- CALL GridAddCoord FOR INPUT GRID."
722  call esmf_gridaddcoord(input_grid, &
723  staggerloc=esmf_staggerloc_center, rc=rc)
724  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
725  call error_handler("IN GridAddCoord", rc)
726 
727  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
728  nullify(lon_src_ptr)
729  call esmf_gridgetcoord(input_grid, &
730  staggerloc=esmf_staggerloc_center, &
731  coorddim=1, &
732  farrayptr=lon_src_ptr, rc=rc)
733  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
734  call error_handler("IN GridGetCoord", rc)
735 
736  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
737  nullify(lat_src_ptr)
738  call esmf_gridgetcoord(input_grid, &
739  staggerloc=esmf_staggerloc_center, &
740  coorddim=2, &
741  computationallbound=clb, &
742  computationalubound=cub, &
743  farrayptr=lat_src_ptr, rc=rc)
744  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
745  call error_handler("IN GridGetCoord", rc)
746 
747  do j = clb(2), cub(2)
748  do i = clb(1), cub(1)
749  lon_src_ptr(i,j) = longitude(i,j)
750  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
751  lat_src_ptr(i,j) = latitude(i,j)
752  enddo
753  enddo
754 
755  print*,"- CALL GridAddCoord FOR INPUT GRID."
756  call esmf_gridaddcoord(input_grid, &
757  staggerloc=esmf_staggerloc_corner, rc=rc)
758  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
759  call error_handler("IN GridAddCoord", rc)
760 
761  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
762  nullify(lon_corner_src_ptr)
763  call esmf_gridgetcoord(input_grid, &
764  staggerloc=esmf_staggerloc_corner, &
765  coorddim=1, &
766  farrayptr=lon_corner_src_ptr, rc=rc)
767  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
768  call error_handler("IN GridGetCoord", rc)
769 
770  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
771  nullify(lat_corner_src_ptr)
772  call esmf_gridgetcoord(input_grid, &
773  staggerloc=esmf_staggerloc_corner, &
774  coorddim=2, &
775  computationallbound=clb, &
776  computationalubound=cub, &
777  farrayptr=lat_corner_src_ptr, rc=rc)
778  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
779  call error_handler("IN GridGetCoord", rc)
780 
781  do j = clb(2), cub(2)
782  do i = clb(1), cub(1)
783  lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon)
784  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
785  if (j == 1) then
786  lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8
787  cycle
788  endif
789  if (j == jp1_input) then
790  lat_corner_src_ptr(i,j) = +90.0_esmf_kind_r8
791  cycle
792  endif
793  lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j))
794  enddo
795  enddo
796 
797  deallocate(latitude,longitude)
798 
799  end subroutine define_input_grid_gfs_grib2
800 
807  subroutine define_input_grid_grib2(localpet, npets)
808 
809  use mpi
810  use netcdf
811  use wgrib2api
812  use program_setup, only : grib2_file_input_grid, data_dir_input_grid, &
813  fix_dir_input_grid
814  implicit none
815 
816  character(len=500) :: the_file, temp_file
817 
818  integer, intent(in) :: localpet, npets
819 
820  integer :: error, extra, i, j, clb(2), cub(2)
821 
822  real(esmf_kind_r4), allocatable :: latitude_one_tile(:,:), lat_corners(:,:)
823  real(esmf_kind_r4), allocatable :: longitude_one_tile(:,:), lon_corners(:,:)
824  real(esmf_kind_r8) :: lat_target(i_target,j_target), &
825  lon_target(i_target,j_target)
826  real(esmf_kind_r8) :: deltalon, dx
827  integer :: ncid,id_var, id_dim
828  real(esmf_kind_r8), pointer :: lat_src_ptr(:,:), lon_src_ptr(:,:)
829  character(len=10000) :: temp_msg
830  character(len=10) :: temp_num = 'NA'
831 
832  num_tiles_input_grid = 1
833 
834  !inv_file = "chgres.inv"
835  the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid
836  temp_file = trim(fix_dir_input_grid)//"/latlon_grid3.32769.nc"
837 
838  call esmf_fieldgather(latitude_target_grid, lat_target, rootpet=0, tile=1, rc=error)
839  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
840  call error_handler("IN FieldGather", error)
841  call esmf_fieldgather(longitude_target_grid, lon_target, rootpet=0, tile=1, rc=error)
842  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
843  call error_handler("IN FieldGather", error)
844 
845  if (localpet==0) then
846  print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file)
847  error=grb2_mk_inv(the_file,inv_file)
848  if (error /=0) call error_handler("OPENING GRIB2 FILE",error)
849  error = grb2_inq(the_file, inv_file,grid_desc=temp_msg)
850  i = index(temp_msg, "grid_template=") + len("grid_template=")
851  j = index(temp_msg,":winds(")
852  temp_num=temp_msg(i:j-1)
853  endif
854  call mpi_barrier(mpi_comm_world, error)
855  call mpi_bcast(temp_num,10,mpi_char,0,mpi_comm_world,error)
856 
857  ! Wgrib2 can't properly read the lat/lon arrays of data on NCEP rotated lat/lon
858  ! grids, so read in lat/lon from fixed coordinate file
859  if (trim(temp_num)=="3.32769" .or. trim(temp_num)=="32769") then
860 
861  input_grid_type = "rotated_latlon"
862 
863  error=nf90_open(trim(temp_file),nf90_nowrite,ncid)
864  call netcdf_err(error, 'opening: '//trim(temp_file))
865 
866  error=nf90_inq_dimid(ncid, 'nx', id_dim)
867  call netcdf_err(error, 'reading nx id' )
868  error=nf90_inquire_dimension(ncid,id_dim,len=i_input)
869  call netcdf_err(error, 'reading nx value' )
870 
871  error=nf90_inq_dimid(ncid, 'ny', id_dim)
872  call netcdf_err(error, 'reading ny id' )
873  error=nf90_inquire_dimension(ncid,id_dim,len=j_input)
874  call netcdf_err(error, 'reading ny value' )
875 
876  allocate(longitude_one_tile(i_input,j_input))
877  allocate(latitude_one_tile(i_input,j_input))
878 
879  error=nf90_inq_varid(ncid, 'gridlat', id_var)
880  call netcdf_err(error, 'reading field id' )
881  error=nf90_get_var(ncid, id_var, latitude_one_tile)
882  call netcdf_err(error, 'reading field' )
883 
884  error=nf90_inq_varid(ncid, 'gridlon', id_var)
885  call netcdf_err(error, 'reading field id' )
886  error=nf90_get_var(ncid, id_var, longitude_one_tile)
887  call netcdf_err(error, 'reading field' )
888 
889  elseif (temp_num == "3.0" .or. temp_num == "3.30" .or. temp_num=="30" .or. temp_num == "0") then
890 
891  if (temp_num =="3.0" .or. temp_num == "0") input_grid_type = "latlon"
892  if (temp_num =="3.30" .or. temp_num=='30') input_grid_type = "lambert"
893 
894  error = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, &
895  lat=latitude_one_tile, lon=longitude_one_tile)
896  if (error /= 1) call error_handler("READING FILE", error)
897 
898 
899  if (localpet==0) print*, "from file lon(1:10,1) = ", longitude_one_tile(1:10,1)
900  if (localpet==0) print*, "from file lat(1,1:10) = ", latitude_one_tile(1,1:10)
901  elseif (temp_num=="NA") then
902  error = 0
903  call error_handler("Grid template number cannot be read from the input file. Please " //&
904  "check that the wgrib2 executable is in your path.", error)
905  else
906  error = 0
907  call error_handler("Unknown input file grid template number. Must be one of: " //&
908  "3, 3.30, 3.32769", error)
909  endif
910 
911  print*,"- I/J DIMENSIONS OF THE INPUT GRID TILES ", i_input, j_input
912 
913  ip1_input = i_input + 1
914  jp1_input = j_input + 1
915 
916 !-----------------------------------------------------------------------
917 ! Create ESMF grid object for the model grid.
918 !-----------------------------------------------------------------------
919 
920  extra = npets / num_tiles_input_grid
921 
922  print*,"- CALL GridCreateNoPeriDim FOR INPUT MODEL GRID"
923  input_grid = esmf_gridcreatenoperidim(maxindex=(/i_input,j_input/), &
924  indexflag=esmf_index_global, &
925  rc=error)
926  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
927  call error_handler("IN GridCreateNoPeriDim", error)
928 
929 
930 !-----------------------------------------------------------------------
931 ! Read the mask and lat/lons.
932 !-----------------------------------------------------------------------
933 
934  print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
935  latitude_input_grid = esmf_fieldcreate(input_grid, &
936  typekind=esmf_typekind_r8, &
937  staggerloc=esmf_staggerloc_center, &
938  name="input_grid_latitude", &
939  rc=error)
940  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
941  call error_handler("IN FieldCreate", error)
942 
943  print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
944  longitude_input_grid = esmf_fieldcreate(input_grid, &
945  typekind=esmf_typekind_r8, &
946  staggerloc=esmf_staggerloc_center, &
947  name="input_grid_longitude", &
948  rc=error)
949  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
950  call error_handler("IN FieldCreate", error)
951 
952  print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE. "
953  call esmf_fieldscatter(latitude_input_grid, real(latitude_one_tile,esmf_kind_r8), rootpet=0, rc=error)
954  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
955  call error_handler("IN FieldScatter", error)
956 
957 print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
958  call esmf_fieldscatter(longitude_input_grid, real(longitude_one_tile,esmf_kind_r8), rootpet=0, rc=error)
959  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
960  call error_handler("IN FieldScatter", error)
961 
962 
963  print*,"- CALL GridAddCoord FOR INPUT GRID."
964  call esmf_gridaddcoord(input_grid, &
965  staggerloc=esmf_staggerloc_center, rc=error)
966  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
967  call error_handler("IN GridAddCoord", error)
968 
969 
970  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
971  nullify(lon_src_ptr)
972  call esmf_gridgetcoord(input_grid, &
973  staggerloc=esmf_staggerloc_center, &
974  coorddim=1, &
975  farrayptr=lon_src_ptr, rc=error)
976  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
977  call error_handler("IN GridGetCoord", error)
978 
979  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
980  nullify(lat_src_ptr)
981  call esmf_gridgetcoord(input_grid, &
982  staggerloc=esmf_staggerloc_center, &
983  coorddim=2, &
984  computationallbound=clb, &
985  computationalubound=cub, &
986  farrayptr=lat_src_ptr, rc=error)
987  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
988  call error_handler("IN GridGetCoord", error)
989 
990  do j = clb(2),cub(2)
991  do i = clb(1), cub(1)
992  lon_src_ptr(i,j)=real(longitude_one_tile(i,j),esmf_kind_r8)
993  lat_src_ptr(i,j)=real(latitude_one_tile(i,j),esmf_kind_r8)
994  enddo
995  enddo
996 
997  print*,"- CALL GridAddCoord FOR INPUT GRID."
998  call esmf_gridaddcoord(input_grid, &
999  staggerloc=esmf_staggerloc_corner, rc=error)
1000  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1001  call error_handler("IN GridAddCoord", error)
1002 
1003  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
1004  nullify(lon_src_ptr)
1005  call esmf_gridgetcoord(input_grid, &
1006  staggerloc=esmf_staggerloc_corner, &
1007  coorddim=1, &
1008  farrayptr=lon_src_ptr, rc=error)
1009  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1010  call error_handler("IN GridGetCoord", error)
1011 
1012  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
1013  nullify(lat_src_ptr)
1014  call esmf_gridgetcoord(input_grid, &
1015  staggerloc=esmf_staggerloc_corner, &
1016  coorddim=2, &
1017  computationallbound=clb, &
1018  computationalubound=cub, &
1019  farrayptr=lat_src_ptr, rc=error)
1020  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1021  call error_handler("IN GridGetCoord", error)
1022 
1023 
1024  ! If we have data on a lat/lon or lambert grid, create staggered coordinates
1025  if(trim(input_grid_type)=="latlon" .or. trim(input_grid_type) == "lambert") then
1026  if (trim(input_grid_type) == "latlon") then
1027 
1028  deltalon = abs(longitude_one_tile(2,1)-longitude_one_tile(1,1))
1029  do j = clb(2), cub(2)
1030  do i = clb(1), cub(1)
1031 
1032  if (i == ip1_input) then
1033  lon_src_ptr(i,j) = longitude_one_tile(i_input,1) + (0.5_esmf_kind_r8*deltalon)
1034  else
1035  lon_src_ptr(i,j) = longitude_one_tile(i,1) - (0.5_esmf_kind_r8*deltalon)
1036  endif
1037 
1038  if (j == jp1_input) then
1039  lat_src_ptr(i,j) = latitude_one_tile(1,j_input) + (0.5_esmf_kind_r8*deltalon)
1040  else
1041  lat_src_ptr(i,j) = latitude_one_tile(1,j) - (0.5_esmf_kind_r8*deltalon)
1042  endif
1043 
1044  enddo
1045  enddo
1046  else
1047  if (localpet==0) then
1048  !cmdline_msg = "wgrib2 "//trim(the_file)//" -d 1 -grid &> temp2.out"
1049  !call system(cmdline_msg)
1050  !open(4,file="temp2.out")
1051  !do i = 1,6
1052  ! read(4,"(A)") temp_msg2
1053  !enddo
1054  !close(4)
1055  print*, trim(temp_msg)
1056  i = index(temp_msg, "Dx ") + len("Dx ")
1057  j = index(temp_msg," m Dy")
1058  read(temp_msg(i:j-1),"(F9.6)") dx
1059  endif
1060  call mpi_barrier(mpi_comm_world,error)
1061  call mpi_bcast(dx,1,mpi_real8,0,mpi_comm_world,error)
1062 
1063  call get_cell_corners(real(latitude_one_tile,esmf_kind_r8), &
1064  real(longitude_one_tile, esmf_kind_r8), &
1065  lat_src_ptr, lon_src_ptr, dx, clb, cub)
1066  endif
1067  elseif (trim(input_grid_type) == "rotated_latlon") then !Read the corner coords from file
1068 
1069  allocate(lon_corners(ip1_input,jp1_input))
1070  allocate(lat_corners(ip1_input,jp1_input))
1071 
1072  error=nf90_inq_varid(ncid, 'gridlon_corners', id_var)
1073  call netcdf_err(error, 'reading field id' )
1074  error=nf90_get_var(ncid, id_var, lon_corners)
1075  call netcdf_err(error, 'reading field' )
1076 
1077  error=nf90_inq_varid(ncid, 'gridlat_corners', id_var)
1078  call netcdf_err(error, 'reading field id' )
1079  error=nf90_get_var(ncid, id_var, lat_corners)
1080  call netcdf_err(error, 'reading field' )
1081 
1082  do j = clb(2),cub(2)
1083  do i = clb(1), cub(1)
1084  lon_src_ptr(i,j)=real(lon_corners(i,j),esmf_kind_r8)
1085  lat_src_ptr(i,j)=real(lat_corners(i,j),esmf_kind_r8)
1086  enddo
1087  enddo
1088 
1089  error= nf90_close(ncid)
1090  endif
1091 
1092  nullify(lon_src_ptr)
1093  nullify(lat_src_ptr)
1094 
1095  deallocate(longitude_one_tile)
1096  deallocate(latitude_one_tile)
1097 
1098  end subroutine define_input_grid_grib2
1099 
1105  subroutine define_target_grid(localpet, npets)
1106 
1107  use netcdf
1108  use program_setup, only : mosaic_file_target_grid, &
1109  orog_dir_target_grid, &
1110  orog_files_target_grid, &
1111  nsoill_out
1112 
1113  implicit none
1114 
1115  integer, intent(in) :: localpet, npets
1116 
1117  character(len=500) :: the_file
1118 
1119  integer :: error, ncid, extra
1120  integer :: id_tiles
1121  integer :: id_dim, id_grid_tiles
1122  integer :: tile
1123  integer, allocatable :: decomptile(:,:)
1124  integer(esmf_kind_i8), allocatable :: landmask_one_tile(:,:)
1125  integer(esmf_kind_i8), allocatable :: seamask_one_tile(:,:)
1126 
1127  real(esmf_kind_r8), allocatable :: latitude_one_tile(:,:)
1128  real(esmf_kind_r8), allocatable :: latitude_s_one_tile(:,:)
1129  real(esmf_kind_r8), allocatable :: latitude_w_one_tile(:,:)
1130  real(esmf_kind_r8), allocatable :: longitude_one_tile(:,:)
1131  real(esmf_kind_r8), allocatable :: longitude_s_one_tile(:,:)
1132  real(esmf_kind_r8), allocatable :: longitude_w_one_tile(:,:)
1133  real(esmf_kind_r8), allocatable :: terrain_one_tile(:,:)
1134 
1135  lsoil_target = nsoill_out
1136 
1137  print*,'- OPEN TARGET GRID MOSAIC FILE: ',trim(mosaic_file_target_grid)
1138  error=nf90_open(trim(mosaic_file_target_grid),nf90_nowrite,ncid)
1139  call netcdf_err(error, 'opening grid mosaic file')
1140 
1141  print*,"- READ NUMBER OF TILES"
1142  error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
1143  call netcdf_err(error, 'reading ntile id')
1144  error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles_target_grid)
1145  call netcdf_err(error, 'reading ntiles')
1146  error=nf90_inq_varid(ncid, 'gridtiles', id_grid_tiles)
1147  call netcdf_err(error, 'reading gridtiles id')
1148  allocate(tiles_target_grid(num_tiles_target_grid))
1149  tiles_target_grid="NULL"
1150  print*,"- READ TILE NAMES"
1151  error=nf90_get_var(ncid, id_grid_tiles, tiles_target_grid)
1152  call netcdf_err(error, 'reading gridtiles')
1153 
1154  error = nf90_close(ncid)
1155 
1156  print*,'- NUMBER OF TILES, TARGET MODEL GRID IS ', num_tiles_target_grid
1157 
1158  if (mod(npets,num_tiles_target_grid) /= 0) then
1159  call error_handler("MUST RUN WITH TASK COUNT THAT IS A MULTIPLE OF # OF TILES.", 1)
1160  endif
1161 
1162 !-----------------------------------------------------------------------
1163 ! Get the model grid specs and land mask from the orography files.
1164 !-----------------------------------------------------------------------
1165 
1166  the_file = trim(orog_dir_target_grid) // trim(orog_files_target_grid(1))
1167 
1168  print*,'- OPEN FIRST TARGET GRID OROGRAPHY FILE: ',trim(the_file)
1169  error=nf90_open(trim(the_file),nf90_nowrite,ncid)
1170  call netcdf_err(error, 'opening orography file')
1171  print*,"- READ GRID DIMENSIONS"
1172  error=nf90_inq_dimid(ncid, 'lon', id_dim)
1173  call netcdf_err(error, 'reading lon id')
1174  error=nf90_inquire_dimension(ncid,id_dim,len=i_target)
1175  call netcdf_err(error, 'reading lon')
1176  error=nf90_inq_dimid(ncid, 'lat', id_dim)
1177  call netcdf_err(error, 'reading lat id')
1178  error=nf90_inquire_dimension(ncid,id_dim,len=j_target)
1179  call netcdf_err(error, 'reading lat')
1180  error = nf90_close(ncid)
1181 
1182  print*,"- I/J DIMENSIONS OF THE TARGET GRID TILES ", i_target, j_target
1183 
1184  ip1_target = i_target + 1
1185  jp1_target = j_target + 1
1186 
1187 !-----------------------------------------------------------------------
1188 ! Create ESMF grid object for the model grid.
1189 !-----------------------------------------------------------------------
1190 
1191  extra = npets / num_tiles_target_grid
1192 
1193  allocate(decomptile(2,num_tiles_target_grid))
1194 
1195  do tile = 1, num_tiles_target_grid
1196  decomptile(:,tile)=(/1,extra/)
1197  enddo
1198 
1199  print*,"- CALL GridCreateMosaic FOR TARGET GRID"
1200  target_grid = esmf_gridcreatemosaic(filename=trim(mosaic_file_target_grid), &
1201  regdecompptile=decomptile, &
1202  staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
1203  esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
1204  indexflag=esmf_index_global, &
1205  tilefilepath=trim(orog_dir_target_grid), rc=error)
1206  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1207  call error_handler("IN GridCreateMosaic", error)
1208 
1209 !-----------------------------------------------------------------------
1210 ! Set target model landmask (1 - land, 0 - not land) and
1211 ! seamask (1 - non-land, 0 -land). Read lat/lon on target grid.
1212 !-----------------------------------------------------------------------
1213 
1214  print*,"- CALL FieldCreate FOR TARGET GRID LANDMASK."
1215  landmask_target_grid = esmf_fieldcreate(target_grid, &
1216  typekind=esmf_typekind_i8, &
1217  staggerloc=esmf_staggerloc_center, &
1218  name="target_grid_landmask", rc=error)
1219  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1220  call error_handler("IN FieldCreate", error)
1221 
1222  print*,"- CALL FieldCreate FOR TARGET GRID SEAMASK."
1223  seamask_target_grid = esmf_fieldcreate(target_grid, &
1224  typekind=esmf_typekind_i8, &
1225  staggerloc=esmf_staggerloc_center, &
1226  name="target_grid_seamask", rc=error)
1227  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1228  call error_handler("IN FieldCreate", error)
1229 
1230  print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE."
1231  latitude_target_grid = esmf_fieldcreate(target_grid, &
1232  typekind=esmf_typekind_r8, &
1233  staggerloc=esmf_staggerloc_center, &
1234  name="target_grid_latitude", rc=error)
1235  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1236  call error_handler("IN FieldCreate", error)
1237 
1238  print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE_S."
1239  latitude_s_target_grid = esmf_fieldcreate(target_grid, &
1240  typekind=esmf_typekind_r8, &
1241  staggerloc=esmf_staggerloc_edge2, &
1242  name="target_grid_latitude_s", rc=error)
1243  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1244  call error_handler("IN FieldCreate", error)
1245 
1246  print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE_W."
1247  latitude_w_target_grid = esmf_fieldcreate(target_grid, &
1248  typekind=esmf_typekind_r8, &
1249  staggerloc=esmf_staggerloc_edge1, &
1250  name="target_grid_latitude_w", &
1251  rc=error)
1252  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1253  call error_handler("IN FieldCreate", error)
1254 
1255  print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE."
1256  longitude_target_grid = esmf_fieldcreate(target_grid, &
1257  typekind=esmf_typekind_r8, &
1258  staggerloc=esmf_staggerloc_center, &
1259  name="target_grid_longitude", &
1260  rc=error)
1261  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1262  call error_handler("IN FieldCreate", error)
1263 
1264  print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE_S."
1265  longitude_s_target_grid = esmf_fieldcreate(target_grid, &
1266  typekind=esmf_typekind_r8, &
1267  staggerloc=esmf_staggerloc_edge2, &
1268  name="target_grid_longitude_s", &
1269  rc=error)
1270  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1271  call error_handler("IN FieldCreate", error)
1272 
1273  print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE_W."
1274  longitude_w_target_grid = esmf_fieldcreate(target_grid, &
1275  typekind=esmf_typekind_r8, &
1276  staggerloc=esmf_staggerloc_edge1, &
1277  name="target_grid_longitude_w", &
1278  rc=error)
1279  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1280  call error_handler("IN FieldCreate", error)
1281 
1282  print*,"- CALL FieldCreate FOR TARGET GRID TERRAIN."
1283  terrain_target_grid = esmf_fieldcreate(target_grid, &
1284  typekind=esmf_typekind_r8, &
1285  staggerloc=esmf_staggerloc_center, &
1286  name="target_grid_terrain", &
1287  rc=error)
1288  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1289  call error_handler("IN FieldCreate", error)
1290 
1291  if (localpet == 0) then
1292  allocate(landmask_one_tile(i_target,j_target))
1293  allocate(seamask_one_tile(i_target,j_target))
1294  allocate(latitude_one_tile(i_target,j_target))
1295  allocate(latitude_s_one_tile(i_target,jp1_target))
1296  allocate(latitude_w_one_tile(ip1_target,j_target))
1297  allocate(longitude_one_tile(i_target,j_target))
1298  allocate(longitude_s_one_tile(i_target,jp1_target))
1299  allocate(longitude_w_one_tile(ip1_target,j_target))
1300  allocate(terrain_one_tile(i_target,j_target))
1301  else
1302  allocate(landmask_one_tile(0,0))
1303  allocate(seamask_one_tile(0,0))
1304  allocate(longitude_one_tile(0,0))
1305  allocate(longitude_s_one_tile(0,0))
1306  allocate(longitude_w_one_tile(0,0))
1307  allocate(latitude_one_tile(0,0))
1308  allocate(latitude_s_one_tile(0,0))
1309  allocate(latitude_w_one_tile(0,0))
1310  allocate(terrain_one_tile(0,0))
1311  endif
1312 
1313  do tile = 1, num_tiles_target_grid
1314  if (localpet == 0) then
1315  the_file = trim(orog_dir_target_grid) // trim(orog_files_target_grid(tile))
1316  call get_model_mask_terrain(trim(the_file), i_target, j_target, landmask_one_tile, &
1317  terrain_one_tile)
1318  seamask_one_tile = 0
1319  where(landmask_one_tile == 0) seamask_one_tile = 1
1320  call get_model_latlons(mosaic_file_target_grid, orog_dir_target_grid, num_tiles_target_grid, tile, &
1321  i_target, j_target, ip1_target, jp1_target, latitude_one_tile, &
1322  latitude_s_one_tile, latitude_w_one_tile, longitude_one_tile, &
1323  longitude_s_one_tile, longitude_w_one_tile)
1324  endif
1325  print*,"- CALL FieldScatter FOR TARGET GRID LANDMASK. TILE IS: ", tile
1326  call esmf_fieldscatter(landmask_target_grid, landmask_one_tile, rootpet=0, tile=tile, rc=error)
1327  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1328  call error_handler("IN FieldScatter", error)
1329  print*,"- CALL FieldScatter FOR TARGET GRID SEAMASK. TILE IS: ", tile
1330  call esmf_fieldscatter(seamask_target_grid, seamask_one_tile, rootpet=0, tile=tile, rc=error)
1331  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1332  call error_handler("IN FieldScatter", error)
1333  print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE. TILE IS: ", tile
1334  call esmf_fieldscatter(longitude_target_grid, longitude_one_tile, rootpet=0, tile=tile, rc=error)
1335  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1336  call error_handler("IN FieldScatter", error)
1337  print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE_S. TILE IS: ", tile
1338  call esmf_fieldscatter(longitude_s_target_grid, longitude_s_one_tile, rootpet=0, tile=tile, rc=error)
1339  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1340  call error_handler("IN FieldScatter", error)
1341  print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE_W. TILE IS: ", tile
1342  call esmf_fieldscatter(longitude_w_target_grid, longitude_w_one_tile, rootpet=0, tile=tile, rc=error)
1343  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1344  call error_handler("IN FieldScatter", error)
1345  print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE. TILE IS: ", tile
1346  call esmf_fieldscatter(latitude_target_grid, latitude_one_tile, rootpet=0, tile=tile, rc=error)
1347  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1348  call error_handler("IN FieldScatter", error)
1349  print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE_S. TILE IS: ", tile
1350  call esmf_fieldscatter(latitude_s_target_grid, latitude_s_one_tile, rootpet=0, tile=tile, rc=error)
1351  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1352  call error_handler("IN FieldScatter", error)
1353  print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE_W. TILE IS: ", tile
1354  call esmf_fieldscatter(latitude_w_target_grid, latitude_w_one_tile, rootpet=0, tile=tile, rc=error)
1355  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1356  call error_handler("IN FieldScatter", error)
1357  print*,"- CALL FieldScatter FOR TARGET GRID TERRAIN. TILE IS: ", tile
1358  call esmf_fieldscatter(terrain_target_grid, terrain_one_tile, rootpet=0, tile=tile, rc=error)
1359  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1360  call error_handler("IN FieldScatter", error)
1361  enddo
1362 
1363  deallocate(landmask_one_tile)
1364  deallocate(seamask_one_tile)
1365  deallocate(longitude_one_tile)
1366  deallocate(longitude_s_one_tile)
1367  deallocate(longitude_w_one_tile)
1368  deallocate(latitude_one_tile)
1369  deallocate(latitude_s_one_tile)
1370  deallocate(latitude_w_one_tile)
1371  deallocate(terrain_one_tile)
1372 
1373  end subroutine define_target_grid
1374 
1393  subroutine get_model_latlons(mosaic_file, orog_dir, num_tiles, tile, &
1394  i_tile, j_tile, ip1_tile, jp1_tile, &
1395  latitude, latitude_s, latitude_w, &
1396  longitude, longitude_s, longitude_w)
1397 
1398  use netcdf
1399 
1400  implicit none
1401 
1402  character(len=*), intent(in) :: mosaic_file, orog_dir
1403 
1404  integer, intent(in) :: num_tiles, tile
1405  integer, intent(in) :: i_tile, j_tile
1406  integer, intent(in) :: ip1_tile, jp1_tile
1407 
1408  real(esmf_kind_r8), intent(out) :: latitude(i_tile, j_tile)
1409  real(esmf_kind_r8), intent(out) :: latitude_s(i_tile, jp1_tile)
1410  real(esmf_kind_r8), intent(out) :: latitude_w(ip1_tile, j_tile)
1411  real(esmf_kind_r8), intent(out) :: longitude(i_tile, j_tile)
1412  real(esmf_kind_r8), intent(out) :: longitude_s(i_tile, jp1_tile)
1413  real(esmf_kind_r8), intent(out) :: longitude_w(ip1_tile, j_tile)
1414 
1415  character(len=25) :: grid_files(num_tiles)
1416  character(len=255) :: grid_file
1417 
1418  integer :: error, id_var, ncid
1419  integer :: id_dim, nxp, nyp, i, j, ii, jj
1420 
1421  real(esmf_kind_r8), allocatable :: tmpvar(:,:)
1422 
1423  print*,"- READ MODEL GRID FILE"
1424 
1425  print*,'- OPEN MOSAIC FILE: ', trim(mosaic_file)
1426  error=nf90_open(trim(mosaic_file), nf90_nowrite, ncid)
1427  call netcdf_err(error, 'opening mosaic file')
1428 
1429  print*,"- READ GRID FILE NAMES"
1430  error=nf90_inq_varid(ncid, 'gridfiles', id_var)
1431  call netcdf_err(error, 'reading gridfiles id')
1432  error=nf90_get_var(ncid, id_var, grid_files)
1433  call netcdf_err(error, 'reading gridfiles')
1434 
1435  error = nf90_close(ncid)
1436 
1437  grid_file = trim(orog_dir) // trim(grid_files(tile))
1438 
1439  print*,'- OPEN GRID FILE: ', trim(grid_file)
1440  error=nf90_open(trim(grid_file), nf90_nowrite, ncid)
1441  call netcdf_err(error, 'opening grid file')
1442 
1443  print*,'- READ NXP ID'
1444  error=nf90_inq_dimid(ncid, 'nxp', id_dim)
1445  call netcdf_err(error, 'reading nxp id')
1446 
1447  print*,'- READ NXP'
1448  error=nf90_inquire_dimension(ncid,id_dim,len=nxp)
1449  call netcdf_err(error, 'reading nxp')
1450 
1451  print*,'- READ NYP ID'
1452  error=nf90_inq_dimid(ncid, 'nyp', id_dim)
1453  call netcdf_err(error, 'reading nyp id')
1454 
1455  print*,'- READ NYP'
1456  error=nf90_inquire_dimension(ncid,id_dim,len=nyp)
1457  call netcdf_err(error, 'reading nyp')
1458 
1459  if ((nxp/2 /= i_tile) .or. (nyp/2 /= j_tile)) then
1460  call error_handler("DIMENSION MISMATCH IN GRID FILE.", 1)
1461  endif
1462 
1463  allocate(tmpvar(nxp,nyp))
1464 
1465  print*,'- READ LONGITUDE ID'
1466  error=nf90_inq_varid(ncid, 'x', id_var)
1467  call netcdf_err(error, 'reading longitude id')
1468 
1469  print*,'- READ LONGITUDE'
1470  error=nf90_get_var(ncid, id_var, tmpvar)
1471  call netcdf_err(error, 'reading longitude')
1472 
1473  do j = 1, j_tile
1474  do i = 1, i_tile
1475  ii = 2*i
1476  jj = 2*j
1477  longitude(i,j) = tmpvar(ii,jj)
1478  enddo
1479  enddo
1480 
1481  do j = 1, jp1_tile
1482  do i = 1, i_tile
1483  ii = 2*i
1484  jj = (2*j) - 1
1485  longitude_s(i,j) = tmpvar(ii,jj)
1486  enddo
1487  enddo
1488 
1489  do j = 1, j_tile
1490  do i = 1, ip1_tile
1491  ii = (2*i) - 1
1492  jj = 2*j
1493  longitude_w(i,j) = tmpvar(ii,jj)
1494  enddo
1495  enddo
1496 
1497  print*,'- READ LATITUDE ID'
1498  error=nf90_inq_varid(ncid, 'y', id_var)
1499  call netcdf_err(error, 'reading latitude id')
1500 
1501  print*,'- READ LATIITUDE'
1502  error=nf90_get_var(ncid, id_var, tmpvar)
1503  call netcdf_err(error, 'reading latitude')
1504 
1505  do j = 1, j_tile
1506  do i = 1, i_tile
1507  ii = 2*i
1508  jj = 2*j
1509  latitude(i,j) = tmpvar(ii,jj)
1510  enddo
1511  enddo
1512 
1513  do j = 1, jp1_tile
1514  do i = 1, i_tile
1515  ii = 2*i
1516  jj = (2*j) - 1
1517  latitude_s(i,j) = tmpvar(ii,jj)
1518  enddo
1519  enddo
1520 
1521  do j = 1, j_tile
1522  do i = 1, ip1_tile
1523  ii = (2*i) - 1
1524  jj = 2*j
1525  latitude_w(i,j) = tmpvar(ii,jj)
1526  enddo
1527  enddo
1528 
1529  deallocate(tmpvar)
1530 
1531  error = nf90_close(ncid)
1532 
1533  end subroutine get_model_latlons
1534 
1547  subroutine get_cell_corners( latitude, longitude, latitude_sw, longitude_sw, dx,clb,cub)
1548  implicit none
1549 
1550  real(esmf_kind_r8), intent(in) :: latitude(i_input,j_input)
1551  real(esmf_kind_r8), intent(inout), pointer :: latitude_sw(:,:)
1552  real(esmf_kind_r8), intent(in) :: longitude(i_input, j_input)
1553  real(esmf_kind_r8), intent(inout), pointer :: longitude_sw(:,:)
1554  real(esmf_kind_r8), intent(in) :: dx
1555 
1556  integer, intent(in) :: clb(2), cub(2)
1557 
1558  real(esmf_kind_r8) :: lat1, lon1, lat2, lon2, d, brng
1559 
1560 
1561  real(esmf_kind_r8), parameter :: pi = 3.14159265359
1562  real(esmf_kind_r8), parameter :: r = 6370000.0
1563  real(esmf_kind_r8), parameter :: bearingindegrees = 135.0
1564 
1565  integer :: i, j
1566 
1567  d = sqrt((dx**2.0_esmf_kind_r8)/2.0_esmf_kind_r8)
1568 
1569  do j = clb(2),cub(2)
1570  do i = clb(1), cub(1)
1571 
1572  if (j == jp1_input .and. i == ip1_input) then
1573  lat1 = latitude(i_input,j_input) * ( pi / 180.0_esmf_kind_r8 )
1574  lon1 = longitude(i_input,j_input) * ( pi / 180.0_esmf_kind_r8 )
1575  brng = 315.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8
1576  lat2 = asin( sin( lat1 ) * cos( d / r ) + cos( lat1 ) * sin( d / r ) * cos( brng ) );
1577  lon2= lon1 + atan2( sin( brng ) * sin( d / r ) * cos( lat1 ), cos( d / r ) - sin( lat1 ) * sin( lat2 ) );
1578  latitude_sw(ip1_input,jp1_input) = lat2 * 180.0_esmf_kind_r8 / pi
1579  longitude_sw(ip1_input,jp1_input) = lon2 * 180.0_esmf_kind_r8 / pi
1580  cycle
1581  endif
1582 
1583  if (i == ip1_input) then
1584  brng = 225.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8
1585  lat1 = latitude(i_input,j) * ( pi / 180.0_esmf_kind_r8 )
1586  lon1 = longitude(i_input,j) * ( pi / 180.0_esmf_kind_r8 )
1587  lat2 = asin( sin( lat1 ) * cos( d / r ) + cos( lat1 ) * sin( d / r ) * cos( brng ) );
1588  lon2= lon1 + atan2( sin( brng ) * sin( d / r ) * cos( lat1 ), cos( d / r ) - sin( lat1 ) * sin( lat2 ) );
1589  latitude_sw(ip1_input,j) = lat2 * 180.0_esmf_kind_r8 / pi
1590  longitude_sw(ip1_input,j) = lon2 * 180.0_esmf_kind_r8 / pi
1591  cycle
1592  endif
1593 
1594  if (j == jp1_input) then
1595  brng = 45.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8
1596  lat1 = latitude(i,j_input) * ( pi / 180.0_esmf_kind_r8 )
1597  lon1 = longitude(i,j_input) * ( pi / 180.0_esmf_kind_r8 )
1598  lat2 = asin( sin( lat1 ) * cos( d / r ) + cos( lat1 ) * sin( d / r ) * cos( brng ) );
1599  lon2= lon1 + atan2( sin( brng ) * sin( d / r ) * cos( lat1 ), cos( d / r ) - sin( lat1 ) * sin( lat2 ) );
1600  latitude_sw(i,jp1_input) = lat2 * 180.0_esmf_kind_r8 / pi
1601  longitude_sw(i,jp1_input) = lon2 * 180.0_esmf_kind_r8 / pi
1602  cycle
1603  endif
1604 
1605  lat1 = latitude(i,j) * ( pi / 180.0_esmf_kind_r8 )
1606  lon1 = longitude(i,j) * ( pi / 180.0_esmf_kind_r8 )
1607 
1608  brng = bearingindegrees * ( pi / 180.0_esmf_kind_r8 );
1609  lat2 = asin( sin( lat1 ) * cos( d / r ) + cos( lat1 ) * sin( d / r ) * cos( brng ) );
1610  lon2= lon1 + atan2( sin( brng ) * sin( d / r ) * cos( lat1 ), cos( d / r ) - sin( lat1 ) * sin( lat2 ) );
1611 
1612  latitude_sw(i,j) = lat2 * 180.0_esmf_kind_r8 / pi
1613  longitude_sw(i,j) = lon2 * 180.0_esmf_kind_r8 / pi
1614 
1615  enddo
1616  enddo
1617 
1618  end subroutine get_cell_corners
1619 
1629  subroutine get_model_mask_terrain(orog_file, idim, jdim, mask, terrain)
1630 
1631  use netcdf
1632 
1633  implicit none
1634 
1635  character(len=*), intent(in) :: orog_file
1636 
1637  integer, intent(in) :: idim, jdim
1638  integer(esmf_kind_i8), intent(out) :: mask(idim,jdim)
1639 
1640  real(esmf_kind_i8), intent(out) :: terrain(idim,jdim)
1641 
1642  integer :: error, lat, lon
1643  integer :: ncid, id_dim, id_var
1644 
1645  real(kind=4), allocatable :: dummy(:,:)
1646 
1647  print*,"- READ MODEL LAND MASK FILE"
1648 
1649  print*,'- OPEN LAND MASK FILE: ', orog_file
1650  error=nf90_open(orog_file,nf90_nowrite,ncid)
1651  call netcdf_err(error, 'opening land mask file')
1652 
1653  print*,"- READ I-DIMENSION"
1654  error=nf90_inq_dimid(ncid, 'lon', id_dim)
1655  call netcdf_err(error, 'reading idim id')
1656  error=nf90_inquire_dimension(ncid,id_dim,len=lon)
1657  call netcdf_err(error, 'reading idim')
1658 
1659  print*,"- READ J-DIMENSION"
1660  error=nf90_inq_dimid(ncid, 'lat', id_dim)
1661  call netcdf_err(error, 'reading jdim id')
1662  error=nf90_inquire_dimension(ncid,id_dim,len=lat)
1663  call netcdf_err(error, 'reading jdim')
1664 
1665  print*,"- I/J DIMENSIONS: ", lon, lat
1666 
1667  if ((lon /= idim) .or. (lat /= jdim)) then
1668  call error_handler("MISMATCH IN DIMENSIONS.", 1)
1669  endif
1670 
1671  allocate(dummy(idim,jdim))
1672 
1673  print*,"- READ LAND MASK"
1674  error=nf90_inq_varid(ncid, 'slmsk', id_var)
1675  call netcdf_err(error, 'reading slmsk id')
1676  error=nf90_get_var(ncid, id_var, dummy)
1677  call netcdf_err(error, 'reading slmsk')
1678  mask = nint(dummy)
1679 
1680  print*,"- READ RAW OROGRAPHY."
1681  error=nf90_inq_varid(ncid, 'orog_raw', id_var)
1682  call netcdf_err(error, 'reading orog_raw id')
1683  error=nf90_get_var(ncid, id_var, dummy)
1684  call netcdf_err(error, 'reading orog_raw')
1685  terrain = dummy
1686 
1687  error = nf90_close(ncid)
1688 
1689  deallocate (dummy)
1690 
1691  end subroutine get_model_mask_terrain
1692 
1696  subroutine cleanup_input_target_grid_data
1697 
1698  implicit none
1699 
1700  integer :: rc
1701 
1702  print*,"- DESTROY MODEL DATA."
1703 
1704  if (esmf_fieldiscreated(latitude_s_input_grid)) then
1705  call esmf_fielddestroy(latitude_s_input_grid, rc=rc)
1706  endif
1707  if (esmf_fieldiscreated(latitude_w_input_grid)) then
1708  call esmf_fielddestroy(latitude_w_input_grid, rc=rc)
1709  endif
1710  if (esmf_fieldiscreated(longitude_s_input_grid)) then
1711  call esmf_fielddestroy(longitude_s_input_grid, rc=rc)
1712  endif
1713  if (esmf_fieldiscreated(longitude_w_input_grid)) then
1714  call esmf_fielddestroy(longitude_w_input_grid, rc=rc)
1715  endif
1716  call esmf_fielddestroy(landmask_target_grid, rc=rc)
1717  call esmf_fielddestroy(latitude_target_grid, rc=rc)
1718  call esmf_fielddestroy(latitude_s_target_grid, rc=rc)
1719  call esmf_fielddestroy(latitude_w_target_grid, rc=rc)
1720  call esmf_fielddestroy(longitude_target_grid, rc=rc)
1721  call esmf_fielddestroy(longitude_s_target_grid, rc=rc)
1722  call esmf_fielddestroy(longitude_w_target_grid, rc=rc)
1723  call esmf_fielddestroy(seamask_target_grid, rc=rc)
1724  call esmf_fielddestroy(terrain_target_grid, rc=rc)
1725  call esmf_griddestroy(input_grid, rc=rc)
1726  call esmf_griddestroy(target_grid, rc=rc)
1727 
1728  end subroutine cleanup_input_target_grid_data
1729 
1730  end module model_grid
model_grid::get_model_mask_terrain
subroutine get_model_mask_terrain(orog_file, idim, jdim, mask, terrain)
Read the model land mask and terrain for a single tile from the orography file.
Definition: model_grid.F90:1629
model_grid::cleanup_input_target_grid_data
subroutine, public cleanup_input_target_grid_data
Deallocate all esmf grid objects.
Definition: model_grid.F90:1696
model_grid::define_input_grid_gfs_grib2
subroutine define_input_grid_gfs_grib2(localpet, npets)
Define input grid object for GFS grib2 data.
Definition: model_grid.F90:618
model_grid::get_cell_corners
subroutine get_cell_corners(latitude, longitude, latitude_sw, longitude_sw, dx, clb, cub)
For grids with equal cell sizes (e.g., lambert conformal), get latitude and longitude of the grid cel...
Definition: model_grid.F90:1547
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:1393
model_grid::define_input_grid_grib2
subroutine define_input_grid_grib2(localpet, npets)
Define input grid object for non-GFS grib2 data.
Definition: model_grid.F90:807
netcdf_err
subroutine netcdf_err(err, string)
Error handler for netcdf.
Definition: utils.F90:31
model_grid::define_input_grid_gaussian
subroutine define_input_grid_gaussian(localpet, npets)
Define grid object for input data on global gaussian grids.
Definition: model_grid.F90:151
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:405
error_handler
subroutine error_handler(string, rc)
General error handler.
Definition: utils.F90:9
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:1105