chgres_cube  1.9.0
•All Data Structures Namespaces Files Functions Variables Pages
atm_input_data.F90
Go to the documentation of this file.
1 
4 
14 
15 module atm_input_data
16  use esmf
17  use netcdf
18  use nemsio_module
19 
20  use program_setup, only : data_dir_input_grid, &
21  atm_files_input_grid, &
22  grib2_file_input_grid, &
23  atm_core_files_input_grid, &
24  atm_tracer_files_input_grid, &
25  tracers_input, num_tracers_input, &
26  tracers, &
27  get_var_cond, &
28  external_model, &
29  read_from_input, &
30  input_type
31  use model_grid, only : input_grid, &
32  i_input, j_input, &
33  ip1_input, jp1_input, &
34  num_tiles_input_grid, &
35  latitude_input_grid, &
36  longitude_input_grid
37  use utilities, only : error_handler, &
38  netcdf_err, &
39  handle_grib_error, &
40  quicksort, &
41  dint2p
42 implicit none
43 
44  private
45 
46 ! Fields associated with the atmospheric model.
47 
48  type(esmf_field), public :: dzdt_input_grid
49  type(esmf_field) :: dpres_input_grid
50  type(esmf_field), public :: pres_input_grid
51  type(esmf_field), public :: ps_input_grid
52  type(esmf_field), public :: terrain_input_grid
53  type(esmf_field), public :: temp_input_grid
54 
55  type(esmf_field), public :: u_input_grid
56  type(esmf_field), public :: v_input_grid
57  type(esmf_field), public :: wind_input_grid
58  type(esmf_field), allocatable, public :: tracers_input_grid(:)
59 
60  integer, public :: lev_input
61  integer, public :: levp1_input
62 
63  character(len=50), private, allocatable :: slevs(:)
64 
65  public :: read_input_atm_data
66  public :: cleanup_input_atm_data
67  public :: convert_winds
68 
69  contains
70 
75  subroutine read_input_atm_data(localpet)
76 
77  implicit none
78 
79  integer, intent(in) :: localpet
80 
81 !-------------------------------------------------------------------------------
82 ! Read the tiled 'warm' restart files.
83 !-------------------------------------------------------------------------------
84 
85  if (trim(input_type) == "restart") then
86 
87  call read_input_atm_restart_file(localpet)
88 
89 !-------------------------------------------------------------------------------
90 ! Read the gaussian history files in netcdf format.
91 !-------------------------------------------------------------------------------
92 
93  elseif (trim(input_type) == "gaussian_netcdf") then
94 
95  call read_input_atm_gaussian_netcdf_file(localpet)
96 
97 !-------------------------------------------------------------------------------
98 ! Read the tiled history files in netcdf format.
99 !-------------------------------------------------------------------------------
100 
101  elseif (trim(input_type) == "history") then
102 
103  call read_input_atm_tiled_history_file(localpet)
104 
105 !-------------------------------------------------------------------------------
106 ! Read the gaussian history files in nemsio format.
107 !-------------------------------------------------------------------------------
108 
109  elseif (trim(input_type) == "gaussian_nemsio") then ! fv3gfs gaussian nemsio
110 
111  call read_input_atm_gaussian_nemsio_file(localpet)
112 
113 !-------------------------------------------------------------------------------
114 ! Read the spectral gfs gaussian history files in nemsio format.
115 !-------------------------------------------------------------------------------
116 
117  elseif (trim(input_type) == "gfs_gaussian_nemsio") then ! spectral gfs gaussian
118  ! nemsio.
119  call read_input_atm_gfs_gaussian_nemsio_file(localpet)
120 
121 !-------------------------------------------------------------------------------
122 ! Read the spectral gfs gaussian history files in sigio format.
123 !-------------------------------------------------------------------------------
124 
125  elseif (trim(input_type) == "gfs_sigio") then ! spectral gfs sigio format.
126 
127  call read_input_atm_gfs_sigio_file(localpet)
128 
129 !-------------------------------------------------------------------------------
130 ! Read fv3gfs data in grib2 format.
131 !-------------------------------------------------------------------------------
132 
133  elseif (trim(input_type) == "grib2") then
134 
135  call read_input_atm_grib2_file(localpet)
136 
137  endif
138 
139  end subroutine read_input_atm_data
140 
141 
145  subroutine init_atm_esmf_fields
146 
147  implicit none
148 
149  integer :: i, rc
150 
151  print*,"- INITIALIZE ATMOSPHERIC ESMF FIELDS."
152 
153  print*,"- CALL FieldCreate FOR INPUT GRID 3-D WIND."
154  wind_input_grid = esmf_fieldcreate(input_grid, &
155  typekind=esmf_typekind_r8, &
156  staggerloc=esmf_staggerloc_center, &
157  ungriddedlbound=(/1,1/), &
158  ungriddedubound=(/lev_input,3/), rc=rc)
159  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
160  call error_handler("IN FieldCreate", rc)
161 
162  print*,"- CALL FieldCreate FOR INPUT GRID SURFACE PRESSURE."
163  ps_input_grid = esmf_fieldcreate(input_grid, &
164  typekind=esmf_typekind_r8, &
165  staggerloc=esmf_staggerloc_center, rc=rc)
166  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
167  call error_handler("IN FieldCreate", rc)
168 
169  print*,"- CALL FieldCreate FOR INPUT GRID TERRAIN."
170  terrain_input_grid = esmf_fieldcreate(input_grid, &
171  typekind=esmf_typekind_r8, &
172  staggerloc=esmf_staggerloc_center, rc=rc)
173  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
174  call error_handler("IN FieldCreate", rc)
175 
176  print*,"- CALL FieldCreate FOR INPUT GRID TEMPERATURE."
177  temp_input_grid = esmf_fieldcreate(input_grid, &
178  typekind=esmf_typekind_r8, &
179  staggerloc=esmf_staggerloc_center, &
180  ungriddedlbound=(/1/), &
181  ungriddedubound=(/lev_input/), rc=rc)
182  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
183  call error_handler("IN FieldCreate", rc)
184 
185  allocate(tracers_input_grid(num_tracers_input))
186 
187  do i = 1, num_tracers_input
188  print*,"- CALL FieldCreate FOR INPUT GRID TRACER ", trim(tracers_input(i))
189  tracers_input_grid(i) = esmf_fieldcreate(input_grid, &
190  typekind=esmf_typekind_r8, &
191  staggerloc=esmf_staggerloc_center, &
192  ungriddedlbound=(/1/), &
193  ungriddedubound=(/lev_input/), rc=rc)
194  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
195  call error_handler("IN FieldCreate", rc)
196  enddo
197 
198  print*,"- CALL FieldCreate FOR INPUT GRID DZDT."
199  dzdt_input_grid = esmf_fieldcreate(input_grid, &
200  typekind=esmf_typekind_r8, &
201  staggerloc=esmf_staggerloc_center, &
202  ungriddedlbound=(/1/), &
203  ungriddedubound=(/lev_input/), rc=rc)
204  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
205  call error_handler("IN FieldCreate", rc)
206 
207  print*,"- CALL FieldCreate FOR INPUT GRID U."
208  u_input_grid = esmf_fieldcreate(input_grid, &
209  typekind=esmf_typekind_r8, &
210  staggerloc=esmf_staggerloc_center, &
211  ungriddedlbound=(/1/), &
212  ungriddedubound=(/lev_input/), rc=rc)
213  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
214  call error_handler("IN FieldCreate", rc)
215 
216  print*,"- CALL FieldCreate FOR INPUT GRID V."
217  v_input_grid = esmf_fieldcreate(input_grid, &
218  typekind=esmf_typekind_r8, &
219  staggerloc=esmf_staggerloc_center, &
220  ungriddedlbound=(/1/), &
221  ungriddedubound=(/lev_input/), rc=rc)
222  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
223  call error_handler("IN FieldCreate", rc)
224 
225  print*,"- CALL FieldCreate FOR INPUT GRID PRESSURE."
226  pres_input_grid = esmf_fieldcreate(input_grid, &
227  typekind=esmf_typekind_r8, &
228  staggerloc=esmf_staggerloc_center, &
229  ungriddedlbound=(/1/), &
230  ungriddedubound=(/lev_input/), rc=rc)
231  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
232  call error_handler("IN FieldCreate", rc)
233 
234  end subroutine init_atm_esmf_fields
235 
241  subroutine read_input_atm_gfs_sigio_file(localpet)
242 
243  use sigio_module
244 
245  implicit none
246 
247  integer, intent(in) :: localpet
248 
249  character(len=300) :: the_file
250 
251  integer(sigio_intkind) :: iret
252  integer :: rc, i, j, k
253  integer :: clb(3), cub(3)
254 
255  real(esmf_kind_r8) :: ak, bk
256  real(esmf_kind_r8), allocatable :: dummy2d(:,:)
257  real(esmf_kind_r8), allocatable :: dummy3d(:,:,:)
258  real(esmf_kind_r8), allocatable :: dummy3d2(:,:,:)
259  real(esmf_kind_r8), pointer :: pptr(:,:,:), psptr(:,:)
260  real(esmf_kind_r8), allocatable :: pi(:,:,:)
261 
262  type(sigio_head) :: sighead
263  type(sigio_dbta) :: sigdata
264 
265  the_file = trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
266 
267  print*,"- ATMOSPHERIC DATA IN SIGIO FORMAT."
268  print*,"- OPEN AND READ: ", trim(the_file)
269 
270  call sigio_sropen(21, trim(the_file), iret)
271  if (iret /= 0) then
272  rc = iret
273  call error_handler("OPENING SPECTRAL GFS SIGIO FILE.", rc)
274  endif
275  call sigio_srhead(21, sighead, iret)
276  if (iret /= 0) then
277  rc = iret
278  call error_handler("READING SPECTRAL GFS SIGIO FILE.", rc)
279  endif
280 
281  lev_input = sighead%levs
282  levp1_input = lev_input + 1
283 
284  if (num_tracers_input /= sighead%ntrac) then
285  call error_handler("WRONG NUMBER OF TRACERS EXPECTED.", 99)
286  endif
287 
288  if (sighead%idvt == 0 .or. sighead%idvt == 21) then
289  if (trim(tracers_input(1)) /= 'spfh' .or. &
290  trim(tracers_input(2)) /= 'o3mr' .or. &
291  trim(tracers_input(3)) /= 'clwmr') then
292  call error_handler("TRACERS SELECTED DO NOT MATCH FILE CONTENTS.", 99)
293  endif
294  else
295  print*,'- UNRECOGNIZED IDVT: ', sighead%idvt
296  call error_handler("UNRECOGNIZED IDVT", 99)
297  endif
298 
299 !---------------------------------------------------------------------------
300 ! Initialize esmf atmospheric fields.
301 !---------------------------------------------------------------------------
302 
303  call init_atm_esmf_fields
304 
305  if (localpet == 0) then
306  allocate(dummy2d(i_input,j_input))
307  allocate(dummy3d(i_input,j_input,lev_input))
308  allocate(dummy3d2(i_input,j_input,lev_input))
309  else
310  allocate(dummy2d(0,0))
311  allocate(dummy3d(0,0,0))
312  allocate(dummy3d2(0,0,0))
313  endif
314 
315  if (localpet == 0) then
316  call sigio_aldbta(sighead, sigdata, iret)
317  if (iret /= 0) then
318  rc = iret
319  call error_handler("ALLOCATING SIGDATA.", rc)
320  endif
321  call sigio_srdbta(21, sighead, sigdata, iret)
322  if (iret /= 0) then
323  rc = iret
324  call error_handler("READING SIGDATA.", rc)
325  endif
326  call sptez(0,sighead%jcap,4,i_input, j_input, sigdata%ps, dummy2d, 1)
327  dummy2d = exp(dummy2d) * 1000.0
328  print*,'surface pres ',maxval(dummy2d),minval(dummy2d)
329  endif
330 
331  print*,"- CALL FieldScatter FOR SURFACE PRESSURE."
332  call esmf_fieldscatter(ps_input_grid, dummy2d, rootpet=0, rc=rc)
333  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
334  call error_handler("IN FieldScatter", rc)
335 
336  if (localpet == 0) then
337  call sptez(0,sighead%jcap,4,i_input, j_input, sigdata%hs, dummy2d, 1)
338  print*,'terrain ',maxval(dummy2d),minval(dummy2d)
339  endif
340 
341  print*,"- CALL FieldScatter FOR TERRAIN."
342  call esmf_fieldscatter(terrain_input_grid, dummy2d, rootpet=0, rc=rc)
343  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
344  call error_handler("IN FieldScatter", rc)
345 
346  do k = 1, num_tracers_input
347 
348  if (localpet == 0) then
349  call sptezm(0,sighead%jcap,4,i_input, j_input, lev_input, sigdata%q(:,:,k), dummy3d, 1)
350  print*,trim(tracers_input(k)),maxval(dummy3d),minval(dummy3d)
351  endif
352 
353  print*,"- CALL FieldScatter FOR INPUT ", trim(tracers_input(k))
354  call esmf_fieldscatter(tracers_input_grid(k), dummy3d, rootpet=0, rc=rc)
355  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
356  call error_handler("IN FieldScatter", rc)
357 
358  enddo
359 
360  if (localpet == 0) then
361  call sptezm(0,sighead%jcap,4,i_input, j_input, lev_input, sigdata%t, dummy3d, 1)
362  print*,'temp ',maxval(dummy3d),minval(dummy3d)
363  endif
364 
365  print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE."
366  call esmf_fieldscatter(temp_input_grid, dummy3d, rootpet=0, rc=rc)
367  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
368  call error_handler("IN FieldScatter", rc)
369 
370 !---------------------------------------------------------------------------
371 ! The spectral gfs files have omega, not vertical velocity. Set to
372 ! zero for now. Convert from omega to vv in the future?
373 !---------------------------------------------------------------------------
374 
375  if (localpet == 0) then
376  print*,"- NO VERTICAL VELOCITY RECORD. SET TO ZERO."
377  dummy3d = 0.0
378  endif
379 
380  print*,"- CALL FieldScatter FOR INPUT DZDT."
381  call esmf_fieldscatter(dzdt_input_grid, dummy3d, rootpet=0, rc=rc)
382  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
383  call error_handler("IN FieldScatter", rc)
384 
385  if (localpet == 0) then
386  call sptezmv(0, sighead%jcap, 4, i_input, j_input, lev_input, sigdata%d, sigdata%z, dummy3d, dummy3d2, 1)
387  print*,'u ',maxval(dummy3d),minval(dummy3d)
388  print*,'v ',maxval(dummy3d2),minval(dummy3d2)
389  endif
390 
391  print*,"- CALL FieldScatter FOR INPUT U-WIND."
392  call esmf_fieldscatter(u_input_grid, dummy3d, rootpet=0, rc=rc)
393  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
394  call error_handler("IN FieldScatter", rc)
395 
396  print*,"- CALL FieldScatter FOR INPUT V-WIND."
397  call esmf_fieldscatter(v_input_grid, dummy3d2, rootpet=0, rc=rc)
398  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
399  call error_handler("IN FieldScatter", rc)
400 
401  deallocate(dummy2d, dummy3d, dummy3d2)
402 
403  if (localpet == 0) call sigio_axdbta(sigdata, iret)
404 
405  call sigio_sclose(21, iret)
406 
407 !---------------------------------------------------------------------------
408 ! Convert from 2-d to 3-d component winds.
409 !---------------------------------------------------------------------------
410 
411  call convert_winds
412 
413 !---------------------------------------------------------------------------
414 ! Compute 3-d pressure from 'ak' and 'bk'.
415 !---------------------------------------------------------------------------
416 
417  print*,"- COMPUTE 3-D PRESSURE."
418 
419  print*,"- CALL FieldGet FOR 3-D PRES."
420  nullify(pptr)
421  call esmf_fieldget(pres_input_grid, &
422  computationallbound=clb, &
423  computationalubound=cub, &
424  farrayptr=pptr, rc=rc)
425  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
426  call error_handler("IN FieldGet", rc)
427 
428  print*,"- CALL FieldGet FOR SURFACE PRESSURE."
429  nullify(psptr)
430  call esmf_fieldget(ps_input_grid, &
431  farrayptr=psptr, rc=rc)
432  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
433  call error_handler("IN FieldGet", rc)
434 
435 !---------------------------------------------------------------------------
436 ! First, compute interface pressure.
437 !---------------------------------------------------------------------------
438 
439  allocate(pi(clb(1):cub(1),clb(2):cub(2),1:levp1_input),stat=rc)
440 
441  do k=1,levp1_input
442  ak = sighead%vcoord(k,1)
443  bk = sighead%vcoord(k,2)
444  do i= clb(1), cub(1)
445  do j= clb(2), cub(2)
446  pi(i,j,k) = ak + bk*psptr(i,j)
447  enddo
448  enddo
449  enddo
450 
451  if (localpet == 0) then
452  print*,'pres int ',psptr(clb(1),clb(2)),pi(clb(1),clb(2),:)
453  endif
454 
455 !---------------------------------------------------------------------------
456 ! Now comput mid-layer pressure from interface pressure.
457 !---------------------------------------------------------------------------
458 
459  do k=1,lev_input
460  do i= clb(1), cub(1)
461  do j= clb(2), cub(2)
462  pptr(i,j,k) = (pi(i,j,k)+pi(i,j,k+1))/2.0_esmf_kind_r8
463  enddo
464  enddo
465  enddo
466 
467  deallocate(pi)
468 
469  if (localpet == 0) then
470  print*,'pres ',psptr(clb(1),clb(2)),pptr(clb(1),clb(2),:)
471  endif
472 
473  end subroutine read_input_atm_gfs_sigio_file
474 
480  subroutine read_input_atm_gfs_gaussian_nemsio_file(localpet)
481 
482  implicit none
483 
484  integer, intent(in) :: localpet
485 
486  character(len=300) :: the_file
487  character(len=20) :: vlevtyp, vname
488 
489  integer(nemsio_intkind) :: vlev, iret
490  integer :: i, j, k, n, rc
491  integer :: clb(3), cub(3)
492 
493  real(nemsio_realkind), allocatable :: vcoord(:,:,:)
494  real(nemsio_realkind), allocatable :: dummy(:)
495  real(esmf_kind_r8), allocatable :: dummy2d(:,:)
496  real(esmf_kind_r8), allocatable :: dummy3d(:,:,:)
497  real(esmf_kind_r8) :: ak, bk
498  real(esmf_kind_r8), allocatable :: pi(:,:,:)
499  real(esmf_kind_r8), pointer :: pptr(:,:,:), psptr(:,:)
500 
501  type(nemsio_gfile) :: gfile
502 
503  the_file = trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
504 
505  print*,"- READ ATMOS DATA FROM SPECTRAL GFS NEMSIO FILE: ", trim(the_file)
506 
507  print*,"- OPEN FILE."
508  call nemsio_open(gfile, the_file, "read", iret=iret)
509  if (iret /= 0) call error_handler("OPENING SPECTRAL GFS NEMSIO ATM FILE.", iret)
510 
511  print*,"- READ NUMBER OF VERTICAL LEVELS."
512  call nemsio_getfilehead(gfile, iret=iret, dimz=lev_input)
513  if (iret /= 0) call error_handler("READING NUMBER OF VERTICAL LEVLES.", iret)
514 
515  levp1_input = lev_input + 1
516 
517  allocate(vcoord(levp1_input,3,2))
518 
519  print*,"- READ VERTICAL COORDINATE INFO."
520  call nemsio_getfilehead(gfile, iret=iret, vcoord=vcoord)
521  if (iret /= 0) call error_handler("READING VERTICAL COORDINATE INFO.", iret)
522 
523 !---------------------------------------------------------------------------
524 ! Initialize esmf atmospheric fields.
525 !---------------------------------------------------------------------------
526 
527  call init_atm_esmf_fields
528 
529  if (localpet == 0) then
530  allocate(dummy(i_input*j_input))
531  allocate(dummy2d(i_input,j_input))
532  allocate(dummy3d(i_input,j_input,lev_input))
533  else
534  allocate(dummy(0))
535  allocate(dummy2d(0,0))
536  allocate(dummy3d(0,0,0))
537  endif
538 
539 !-----------------------------------------------------------------------
540 ! 3-d fields in gaussian files increment from bottom to model top.
541 ! That is what is expected by this program, so no need to flip indices.
542 !-----------------------------------------------------------------------
543 
544  if (localpet == 0) then
545  print*,"- READ TEMPERATURE."
546  vname = "tmp"
547  vlevtyp = "mid layer"
548  do vlev = 1, lev_input
549  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
550  if (iret /= 0) call error_handler("READING TEMPERATURE RECORD.", iret)
551  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
552 ! print*,'temp check after read ',vlev, dummy3d(1,1,vlev)
553  enddo
554  endif
555 
556  print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE."
557  call esmf_fieldscatter(temp_input_grid, dummy3d, rootpet=0, rc=rc)
558  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
559  call error_handler("IN FieldScatter", rc)
560 
561  do n = 1, num_tracers_input
562 
563  if (localpet == 0) then
564  print*,"- READ ", trim(tracers_input(n))
565  vname = trim(tracers_input(n))
566  vlevtyp = "mid layer"
567  do vlev = 1, lev_input
568  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
569  if (iret /= 0) call error_handler("READING TRACER RECORD.", iret)
570 ! print*,'tracer ',vlev, maxval(dummy),minval(dummy)
571  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
572  enddo
573  endif
574 
575  print*,"- CALL FieldScatter FOR INPUT ", trim(tracers_input(n))
576  call esmf_fieldscatter(tracers_input_grid(n), dummy3d, rootpet=0, rc=rc)
577  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
578  call error_handler("IN FieldScatter", rc)
579 
580  enddo
581 
582  if (localpet == 0) then
583  print*,"- READ U-WINDS."
584  vname = "ugrd"
585  vlevtyp = "mid layer"
586  do vlev = 1, lev_input
587  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
588  if (iret /= 0) call error_handler("READING U-WIND RECORD.", iret)
589 ! print*,'ugrd ',vlev, maxval(dummy),minval(dummy)
590  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
591  enddo
592  endif
593 
594  print*,"- CALL FieldScatter FOR INPUT U-WIND."
595  call esmf_fieldscatter(u_input_grid, dummy3d, rootpet=0, rc=rc)
596  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
597  call error_handler("IN FieldScatter", rc)
598 
599  if (localpet == 0) then
600  print*,"- READ V-WINDS."
601  vname = "vgrd"
602  vlevtyp = "mid layer"
603  do vlev = 1, lev_input
604  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
605  if (iret /= 0) call error_handler("READING V-WIND RECORD.", iret)
606 ! print*,'vgrd ',vlev, maxval(dummy),minval(dummy)
607  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
608  enddo
609  endif
610 
611  print*,"- CALL FieldScatter FOR INPUT V-WIND."
612  call esmf_fieldscatter(v_input_grid, dummy3d, rootpet=0, rc=rc)
613  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
614  call error_handler("IN FieldScatter", rc)
615 
616 !---------------------------------------------------------------------------
617 ! The spectral gfs nemsio files do not have a vertical velocity or
618 ! omega record. So set to zero for now.
619 !---------------------------------------------------------------------------
620 
621  if (localpet == 0) then
622  print*,"- NO VERTICAL VELOCITY RECORD. SET TO ZERO."
623  dummy3d = 0.0
624  endif
625 
626  print*,"- CALL FieldScatter FOR INPUT DZDT."
627  call esmf_fieldscatter(dzdt_input_grid, dummy3d, rootpet=0, rc=rc)
628  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
629  call error_handler("IN FieldScatter", rc)
630 
631  if (localpet == 0) then
632  print*,"- READ HGT."
633  vname = "hgt"
634  vlevtyp = "sfc"
635  vlev = 1
636  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
637  if (iret /= 0) call error_handler("READING HGT RECORD.", iret)
638 ! print*,'hgt ',vlev, maxval(dummy),minval(dummy)
639  dummy2d = reshape(dummy, (/i_input,j_input/))
640  endif
641 
642  print*,"- CALL FieldScatter FOR TERRAIN."
643  call esmf_fieldscatter(terrain_input_grid, dummy2d, rootpet=0, rc=rc)
644  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
645  call error_handler("IN FieldScatter", rc)
646 
647  if (localpet == 0) then
648  print*,"- READ PRES."
649  vname = "pres"
650  vlevtyp = "sfc"
651  vlev = 1
652  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
653  if (iret /= 0) call error_handler("READING PRES RECORD.", iret)
654 ! print*,'pres ',vlev, maxval(dummy),minval(dummy)
655  dummy2d = reshape(dummy, (/i_input,j_input/))
656  endif
657 
658  print*,"- CALL FieldScatter FOR SURFACE PRESSURE."
659  call esmf_fieldscatter(ps_input_grid, dummy2d, rootpet=0, rc=rc)
660  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
661  call error_handler("IN FieldScatter", rc)
662 
663  call nemsio_close(gfile)
664 
665  deallocate(dummy, dummy2d, dummy3d)
666 
667 !---------------------------------------------------------------------------
668 ! Convert from 2-d to 3-d component winds.
669 !---------------------------------------------------------------------------
670 
671  call convert_winds
672 
673 !---------------------------------------------------------------------------
674 ! Compute 3-d pressure from 'ak' and 'bk'.
675 !---------------------------------------------------------------------------
676 
677  print*,"- COMPUTE 3-D PRESSURE."
678 
679  print*,"- CALL FieldGet FOR 3-D PRES."
680  nullify(pptr)
681  call esmf_fieldget(pres_input_grid, &
682  computationallbound=clb, &
683  computationalubound=cub, &
684  farrayptr=pptr, rc=rc)
685  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
686  call error_handler("IN FieldGet", rc)
687 
688  print*,"- CALL FieldGet FOR SURFACE PRESSURE."
689  nullify(psptr)
690  call esmf_fieldget(ps_input_grid, &
691  farrayptr=psptr, rc=rc)
692  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
693  call error_handler("IN FieldGet", rc)
694 
695 !---------------------------------------------------------------------------
696 ! First, compute interface pressure.
697 !---------------------------------------------------------------------------
698 
699  allocate(pi(clb(1):cub(1),clb(2):cub(2),1:levp1_input))
700 
701  do k=1,levp1_input
702  ak = vcoord(k,1,1)
703  bk = vcoord(k,2,1)
704  do i= clb(1), cub(1)
705  do j= clb(2), cub(2)
706  pi(i,j,k) = ak + bk*psptr(i,j)
707  enddo
708  enddo
709  enddo
710 
711  deallocate(vcoord)
712 
713 !---------------------------------------------------------------------------
714 ! Now comput mid-layer pressure from interface pressure.
715 !---------------------------------------------------------------------------
716 
717  do k=1,lev_input
718  do i= clb(1), cub(1)
719  do j= clb(2), cub(2)
720  pptr(i,j,k) = (pi(i,j,k)+pi(i,j,k+1))/2.0
721  enddo
722  enddo
723  enddo
724 
725  deallocate(pi)
726 
727  end subroutine read_input_atm_gfs_gaussian_nemsio_file
728 
733  subroutine read_input_atm_gaussian_nemsio_file(localpet)
734 
735  implicit none
736 
737  integer, intent(in) :: localpet
738 
739  character(len=300) :: the_file
740  character(len=20) :: vlevtyp, vname
741 
742  integer :: i, j, k, n
743  integer :: rc, clb(3), cub(3)
744  integer(nemsio_intkind) :: vlev, iret
745 
746  real(nemsio_realkind), allocatable :: vcoord(:,:,:)
747  real(nemsio_realkind), allocatable :: dummy(:)
748  real(esmf_kind_r8), allocatable :: dummy2d(:,:)
749  real(esmf_kind_r8), allocatable :: dummy3d(:,:,:)
750  real(esmf_kind_r8), pointer :: presptr(:,:,:), psptr(:,:)
751  real(esmf_kind_r8), pointer :: dpresptr(:,:,:)
752  real(esmf_kind_r8), allocatable :: pres_interface(:)
753 
754  type(nemsio_gfile) :: gfile
755 
756  the_file = trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
757 
758  print*,"- READ ATMOS DATA FROM GAUSSIAN NEMSIO FILE: ", trim(the_file)
759 
760  print*,"- OPEN FILE."
761  call nemsio_open(gfile, the_file, "read", iret=iret)
762  if (iret /= 0) call error_handler("OPENING GAUSSIAN NEMSIO ATM FILE.", iret)
763 
764  print*,"- READ NUMBER OF VERTICAL LEVELS."
765  call nemsio_getfilehead(gfile, iret=iret, dimz=lev_input)
766  if (iret /= 0) call error_handler("READING NUMBER OF VERTICAL LEVLES.", iret)
767 
768  levp1_input = lev_input + 1
769 
770  allocate(vcoord(levp1_input,3,2))
771 
772  print*,"- READ VERTICAL COORDINATE INFO."
773  call nemsio_getfilehead(gfile, iret=iret, vcoord=vcoord)
774  if (iret /= 0) call error_handler("READING VERTICAL COORDINATE INFO.", iret)
775 
776 !---------------------------------------------------------------------------
777 ! Initialize esmf atmospheric fields.
778 !---------------------------------------------------------------------------
779 
780  call init_atm_esmf_fields
781 
782  print*,"- CALL FieldCreate FOR INPUT DPRES."
783  dpres_input_grid = esmf_fieldcreate(input_grid, &
784  typekind=esmf_typekind_r8, &
785  staggerloc=esmf_staggerloc_center, &
786  ungriddedlbound=(/1/), &
787  ungriddedubound=(/lev_input/), rc=rc)
788  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
789  call error_handler("IN FieldCreate", rc)
790 
791  if (localpet == 0) then
792  allocate(dummy(i_input*j_input))
793  allocate(dummy2d(i_input,j_input))
794  allocate(dummy3d(i_input,j_input,lev_input))
795  else
796  allocate(dummy(0))
797  allocate(dummy2d(0,0))
798  allocate(dummy3d(0,0,0))
799  endif
800 
801 !-----------------------------------------------------------------------
802 ! 3-d fields in gaussian files increment from bottom to model top.
803 ! That is what is expected by this program, so no need to flip indices.
804 !-----------------------------------------------------------------------
805 
806  if (localpet == 0) then
807  print*,"- READ TEMPERATURE."
808  vname = "tmp"
809  vlevtyp = "mid layer"
810  do vlev = 1, lev_input
811  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
812  if (iret /= 0) call error_handler("READING TEMPERATURE RECORD.", iret)
813  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
814  print*,'temp check after read ',vlev, dummy3d(1,1,vlev)
815  enddo
816  endif
817 
818  print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE."
819  call esmf_fieldscatter(temp_input_grid, dummy3d, rootpet=0, rc=rc)
820  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
821  call error_handler("IN FieldScatter", rc)
822 
823  do n = 1, num_tracers_input
824 
825  if (localpet == 0) then
826  print*,"- READ ", trim(tracers_input(n))
827  vname = trim(tracers_input(n))
828  vlevtyp = "mid layer"
829  do vlev = 1, lev_input
830  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
831  if (iret /= 0) call error_handler("READING TRACER RECORD.", iret)
832  print*,'tracer ',vlev, maxval(dummy),minval(dummy)
833  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
834  enddo
835  endif
836 
837  print*,"- CALL FieldScatter FOR INPUT ", trim(tracers_input(n))
838  call esmf_fieldscatter(tracers_input_grid(n), dummy3d, rootpet=0, rc=rc)
839  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
840  call error_handler("IN FieldScatter", rc)
841 
842  enddo
843 
844  if (localpet == 0) then
845  print*,"- READ U-WINDS."
846  vname = "ugrd"
847  vlevtyp = "mid layer"
848  do vlev = 1, lev_input
849  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
850  if (iret /= 0) call error_handler("READING U-WIND RECORD.", iret)
851  print*,'ugrd ',vlev, maxval(dummy),minval(dummy)
852  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
853  enddo
854  endif
855 
856  print*,"- CALL FieldScatter FOR INPUT U-WIND."
857  call esmf_fieldscatter(u_input_grid, dummy3d, rootpet=0, rc=rc)
858  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
859  call error_handler("IN FieldScatter", rc)
860 
861  if (localpet == 0) then
862  print*,"- READ V-WINDS."
863  vname = "vgrd"
864  vlevtyp = "mid layer"
865  do vlev = 1, lev_input
866  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
867  if (iret /= 0) call error_handler("READING V-WIND RECORD.", iret)
868  print*,'vgrd ',vlev, maxval(dummy),minval(dummy)
869  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
870  enddo
871  endif
872 
873  print*,"- CALL FieldScatter FOR INPUT V-WIND."
874  call esmf_fieldscatter(v_input_grid, dummy3d, rootpet=0, rc=rc)
875  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
876  call error_handler("IN FieldScatter", rc)
877 
878  if (localpet == 0) then
879  print*,"- READ DPRES."
880  vname = "dpres"
881  vlevtyp = "mid layer"
882  do vlev = 1, lev_input
883  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
884  if (iret /= 0) call error_handler("READING DPRES RECORD.", iret)
885  print*,'dpres ',vlev, maxval(dummy),minval(dummy)
886  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
887  enddo
888  endif
889 
890  print*,"- CALL FieldScatter FOR INPUT DPRES."
891  call esmf_fieldscatter(dpres_input_grid, dummy3d, rootpet=0, rc=rc)
892  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
893  call error_handler("IN FieldScatter", rc)
894 
895  if (localpet == 0) then
896  print*,"- READ DZDT."
897  vname = "dzdt"
898  vlevtyp = "mid layer"
899  do vlev = 1, lev_input
900  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
901  if (iret /= 0) call error_handler("READING DZDT RECORD.", iret)
902  print*,'dzdt ',vlev, maxval(dummy),minval(dummy)
903  dummy3d(:,:,vlev) = reshape(dummy, (/i_input,j_input/))
904  enddo
905  endif
906 
907  print*,"- CALL FieldScatter FOR INPUT DZDT."
908  call esmf_fieldscatter(dzdt_input_grid, dummy3d, rootpet=0, rc=rc)
909  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
910  call error_handler("IN FieldScatter", rc)
911 
912  if (localpet == 0) then
913  print*,"- READ HGT."
914  vname = "hgt"
915  vlevtyp = "sfc"
916  vlev = 1
917  call nemsio_readrecv(gfile, vname, vlevtyp, vlev, dummy, 0, iret)
918  if (iret /= 0) call error_handler("READING HGT RECORD.", iret)
919  print*,'hgt ',vlev, maxval(dummy),minval(dummy)
920  dummy2d = reshape(dummy, (/i_input,j_input/))
921  endif
922 
923  print*,"- CALL FieldScatter FOR TERRAIN."
924  call esmf_fieldscatter(terrain_input_grid, dummy2d, rootpet=0, rc=rc)
925  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
926  call error_handler("IN FieldScatter", rc)
927 
928  call nemsio_close(gfile)
929 
930  deallocate(dummy, dummy2d, dummy3d)
931 
932 !---------------------------------------------------------------------------
933 ! Convert from 2-d to 3-d component winds.
934 !---------------------------------------------------------------------------
935 
936  call convert_winds
937 
938 !---------------------------------------------------------------------------
939 ! Compute 3-d pressure. Mid-layer and surface pressure are computed
940 ! from delta p. The surface pressure in the file is not used. After
941 ! the model's write component interpolates from the cubed-sphere grid
942 ! to the gaussian grid, the surface pressure is no longer consistent
943 ! with the delta p (per Jun Wang).
944 !---------------------------------------------------------------------------
945 
946  print*,"- COMPUTE 3-D PRESSURE."
947 
948  print*,"- CALL FieldGet FOR DELTA PRESSURE."
949  nullify(dpresptr)
950  call esmf_fieldget(dpres_input_grid, &
951  computationallbound=clb, &
952  computationalubound=cub, &
953  farrayptr=dpresptr, rc=rc)
954  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
955  call error_handler("IN FieldGet", rc)
956 
957  print*,"- CALL FieldGet FOR 3-D PRESSURE."
958  nullify(presptr)
959  call esmf_fieldget(pres_input_grid, &
960  farrayptr=presptr, rc=rc)
961  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
962  call error_handler("IN FieldGet", rc)
963 
964  print*,"- CALL FieldGet FOR SURFACE PRESSURE."
965  nullify(psptr)
966  call esmf_fieldget(ps_input_grid, &
967  farrayptr=psptr, rc=rc)
968  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
969  call error_handler("IN FieldGet", rc)
970 
971  allocate(pres_interface(levp1_input))
972 
973  if (localpet == 0) then
974  do k = clb(3), cub(3)
975  print*,'dpres is ',cub(1),cub(2),k, dpresptr(cub(1),cub(2),k)
976  enddo
977  endif
978 
979  do i = clb(1), cub(1)
980  do j = clb(2), cub(2)
981  pres_interface(levp1_input) = vcoord(levp1_input,1,1)
982  do k = lev_input, 1, -1
983  pres_interface(k) = pres_interface(k+1) + dpresptr(i,j,k)
984  enddo
985  psptr(i,j) = pres_interface(1)
986  do k = 1, lev_input
987  presptr(i,j,k) = (pres_interface(k) + pres_interface(k+1)) / 2.0_8
988  enddo
989  enddo
990  enddo
991 
992  deallocate(vcoord)
993 
994  if (localpet == 0) then
995  print*,'psfc is ',clb(1),clb(2),psptr(clb(1),clb(2))
996  print*,'pres is ',clb(1),clb(2),presptr(clb(1),clb(2),:)
997  endif
998 
999  print*,'pres check 1',localpet,maxval(presptr(:,:,1)),minval(presptr(:,:,1))
1000  print*,'pres check lev',localpet,maxval(presptr(:,:,lev_input)),minval(presptr(:,:,lev_input))
1001 
1002  deallocate(pres_interface)
1003 
1004  call esmf_fielddestroy(dpres_input_grid, rc=rc)
1005 
1006  end subroutine read_input_atm_gaussian_nemsio_file
1007 
1016  subroutine read_input_atm_restart_file(localpet)
1017 
1018  implicit none
1019 
1020  integer, intent(in) :: localpet
1021 
1022  character(len=500) :: tilefile
1023 
1024  integer :: i, j, k
1025  integer :: clb(3), cub(3)
1026  integer :: rc, tile, ncid, id_var
1027  integer :: error, id_dim
1028 
1029  real(esmf_kind_r8), allocatable :: ak(:)
1030  real(esmf_kind_r8), pointer :: presptr(:,:,:), psptr(:,:)
1031  real(esmf_kind_r8), pointer :: dpresptr(:,:,:)
1032  real(esmf_kind_r8), allocatable :: data_one_tile(:,:)
1033  real(esmf_kind_r8), allocatable :: data_one_tile_3d(:,:,:)
1034  real(esmf_kind_r8), allocatable :: pres_interface(:)
1035 
1036 !---------------------------------------------------------------------------
1037 ! Get number of vertical levels and model top pressure.
1038 !---------------------------------------------------------------------------
1039 
1040  tilefile = trim(data_dir_input_grid) // "/" // trim(atm_core_files_input_grid(7))
1041  print*,"- READ ATM VERTICAL LEVELS FROM: ", trim(tilefile)
1042  error=nf90_open(trim(tilefile),nf90_nowrite,ncid)
1043  call netcdf_err(error, 'opening: '//trim(tilefile) )
1044 
1045  error=nf90_inq_dimid(ncid, 'xaxis_1', id_dim)
1046  call netcdf_err(error, 'reading xaxis_1 id' )
1047  error=nf90_inquire_dimension(ncid,id_dim,len=levp1_input)
1048  call netcdf_err(error, 'reading xaxis_1 value' )
1049 
1050  lev_input = levp1_input - 1
1051 
1052  allocate(ak(levp1_input))
1053 
1054  error=nf90_inq_varid(ncid, 'ak', id_var)
1055  call netcdf_err(error, 'reading field id' )
1056  error=nf90_get_var(ncid, id_var, ak)
1057  call netcdf_err(error, 'reading ak' )
1058 
1059  error = nf90_close(ncid)
1060 
1061 !---------------------------------------------------------------------------
1062 ! Initialize esmf atmospheric fields.
1063 !---------------------------------------------------------------------------
1064 
1065  call init_atm_esmf_fields
1066 
1067  print*,"- CALL FieldCreate FOR INPUT GRID DELTA PRESSURE."
1068  dpres_input_grid = esmf_fieldcreate(input_grid, &
1069  typekind=esmf_typekind_r8, &
1070  staggerloc=esmf_staggerloc_center, &
1071  ungriddedlbound=(/1/), &
1072  ungriddedubound=(/lev_input/), rc=rc)
1073  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1074  call error_handler("IN FieldCreate", rc)
1075 
1076  if (localpet < num_tiles_input_grid) then
1077  allocate(data_one_tile_3d(i_input,j_input,lev_input))
1078  allocate(data_one_tile(i_input,j_input))
1079  else
1080  allocate(data_one_tile_3d(0,0,0))
1081  allocate(data_one_tile(0,0))
1082  endif
1083 
1084  if (localpet < num_tiles_input_grid) then
1085  tile = localpet+1
1086  tilefile= trim(data_dir_input_grid) // "/" // trim(atm_core_files_input_grid(tile))
1087  print*,"- READ ATMOSPHERIC CORE FILE: ", trim(tilefile)
1088  error=nf90_open(trim(tilefile),nf90_nowrite,ncid)
1089  call netcdf_err(error, 'opening: '//trim(tilefile) )
1090  endif
1091 
1092  if (localpet < num_tiles_input_grid) then
1093  error=nf90_inq_varid(ncid, 'phis', id_var)
1094  call netcdf_err(error, 'reading field id' )
1095  error=nf90_get_var(ncid, id_var, data_one_tile)
1096  call netcdf_err(error, 'reading field' )
1097  data_one_tile = data_one_tile / 9.806_8 ! geopotential height
1098  endif
1099 
1100  do tile = 1, num_tiles_input_grid
1101  print*,"- CALL FieldScatter FOR INPUT GRID TERRAIN for tile ",tile
1102  call esmf_fieldscatter(terrain_input_grid, data_one_tile, rootpet=tile-1, tile=tile, rc=rc)
1103  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1104  call error_handler("IN FieldScatter", rc)
1105  enddo
1106 
1107  if (localpet < num_tiles_input_grid) then
1108 ! error=nf90_inq_varid(ncid, 'W', id_var)
1109 ! call netcdf_err(error, 'reading field id' )
1110 ! error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1111 ! call netcdf_err(error, 'reading field' )
1112 ! data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1113 
1114 ! Using 'w' from restart files has caused problems. Set to zero.
1115  data_one_tile_3d = 0.0_8
1116  endif
1117 
1118  do tile = 1, num_tiles_input_grid
1119  print*,"- CALL FieldScatter FOR INPUT GRID VERTICAL VELOCITY for tile ",tile
1120  call esmf_fieldscatter(dzdt_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1121  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1122  call error_handler("IN FieldScatter", rc)
1123  enddo
1124 
1125  if (localpet < num_tiles_input_grid) then
1126  error=nf90_inq_varid(ncid, 'T', id_var)
1127  call netcdf_err(error, 'reading field id' )
1128  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1129  call netcdf_err(error, 'reading field' )
1130  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1131  endif
1132 
1133  do tile = 1, num_tiles_input_grid
1134  print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE."
1135  call esmf_fieldscatter(temp_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1136  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1137  call error_handler("IN FieldScatter", rc)
1138  enddo
1139 
1140  if (localpet < num_tiles_input_grid) then
1141  error=nf90_inq_varid(ncid, 'delp', id_var)
1142  call netcdf_err(error, 'reading field id' )
1143  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1144  call netcdf_err(error, 'reading field' )
1145  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1146  endif
1147 
1148  do tile = 1, num_tiles_input_grid
1149  print*,"- CALL FieldScatter FOR INPUT DELTA PRESSURE."
1150  call esmf_fieldscatter(dpres_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1151  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1152  call error_handler("IN FieldScatter", rc)
1153  enddo
1154 
1155  if (localpet < num_tiles_input_grid) then
1156  error=nf90_inq_varid(ncid, 'ua', id_var)
1157  call netcdf_err(error, 'reading field id' )
1158  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1159  call netcdf_err(error, 'reading field' )
1160  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1161  endif
1162 
1163  do tile = 1, num_tiles_input_grid
1164  print*,"- CALL FieldScatter FOR INPUT GRID U."
1165  call esmf_fieldscatter(u_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1166  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1167  call error_handler("IN FieldScatter", rc)
1168  enddo
1169 
1170  if (localpet < num_tiles_input_grid) then
1171  error=nf90_inq_varid(ncid, 'va', id_var)
1172  call netcdf_err(error, 'reading field id' )
1173  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1174  call netcdf_err(error, 'reading field' )
1175  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1176  endif
1177 
1178  do tile = 1, num_tiles_input_grid
1179  print*,"- CALL FieldScatter FOR INPUT GRID V."
1180  call esmf_fieldscatter(v_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1181  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1182  call error_handler("IN FieldScatter", rc)
1183  enddo
1184 
1185  if (localpet < num_tiles_input_grid) error = nf90_close(ncid)
1186 
1187  if (localpet < num_tiles_input_grid) then
1188  tile = localpet+1
1189  tilefile= trim(data_dir_input_grid) // "/" // trim(atm_tracer_files_input_grid(tile))
1190  print*,"- READ ATMOSPHERIC TRACER FILE: ", trim(tilefile)
1191  error=nf90_open(trim(tilefile),nf90_nowrite,ncid)
1192  call netcdf_err(error, 'opening: '//trim(tilefile) )
1193  endif
1194 
1195  do i = 1, num_tracers_input
1196 
1197  if (localpet < num_tiles_input_grid) then
1198  error=nf90_inq_varid(ncid, tracers_input(i), id_var)
1199  call netcdf_err(error, 'reading field id' )
1200  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1201  call netcdf_err(error, 'reading field' )
1202  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1203  endif
1204 
1205  do tile = 1, num_tiles_input_grid
1206  print*,"- CALL FieldScatter FOR INPUT ", trim(tracers_input(i))
1207  call esmf_fieldscatter(tracers_input_grid(i), data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1208  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1209  call error_handler("IN FieldScatter", rc)
1210  enddo
1211 
1212  enddo
1213 
1214  if (localpet < num_tiles_input_grid) error=nf90_close(ncid)
1215 
1216 !---------------------------------------------------------------------------
1217 ! Convert from 2-d to 3-d cartesian winds.
1218 !---------------------------------------------------------------------------
1219 
1220  call convert_winds
1221 
1222 !---------------------------------------------------------------------------
1223 ! Compute pressures
1224 !---------------------------------------------------------------------------
1225 
1226  print*,"- CALL FieldGet FOR SURFACE PRESSURE."
1227  call esmf_fieldget(ps_input_grid, &
1228  farrayptr=psptr, rc=rc)
1229  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1230  call error_handler("IN FieldGet", rc)
1231 
1232  print*,"- CALL FieldGet FOR PRESSURE."
1233  call esmf_fieldget(pres_input_grid, &
1234  computationallbound=clb, &
1235  computationalubound=cub, &
1236  farrayptr=presptr, rc=rc)
1237  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1238  call error_handler("IN FieldGet", rc)
1239 
1240  print*,"- CALL FieldGet FOR DELTA PRESSURE."
1241  call esmf_fieldget(dpres_input_grid, &
1242  farrayptr=dpresptr, rc=rc)
1243  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1244  call error_handler("IN FieldGet", rc)
1245 
1246  allocate(pres_interface(levp1_input))
1247 
1248  do i = clb(1), cub(1)
1249  do j = clb(2), cub(2)
1250  pres_interface(levp1_input) = ak(1) ! model top in Pa
1251  do k = (levp1_input-1), 1, -1
1252  pres_interface(k) = pres_interface(k+1) + dpresptr(i,j,k)
1253  enddo
1254  do k = 1, lev_input
1255  presptr(i,j,k) = (pres_interface(k) + pres_interface(k+1)) / 2.0_8
1256  enddo
1257  psptr(i,j) = pres_interface(1)
1258  enddo
1259  enddo
1260 
1261  deallocate(ak)
1262  deallocate(pres_interface)
1263 
1264  call esmf_fielddestroy(dpres_input_grid, rc=rc)
1265 
1266  deallocate(data_one_tile_3d, data_one_tile)
1267 
1268  end subroutine read_input_atm_restart_file
1269 
1275  subroutine read_input_atm_gaussian_netcdf_file(localpet)
1276 
1277  use mpi
1278 
1279  implicit none
1280 
1281  integer, intent(in) :: localpet
1282 
1283  character(len=500) :: tilefile
1284 
1285  integer :: start(3), count(3), iscnt
1286  integer :: error, ncid, num_tracers_file
1287  integer :: id_dim, idim_input, jdim_input
1288  integer :: id_var, rc, nprocs, max_procs
1289  integer :: kdim, remainder, myrank, i, j, k, n
1290  integer :: clb(3), cub(3)
1291  integer, allocatable :: kcount(:), startk(:), displ(:)
1292  integer, allocatable :: ircnt(:)
1293 
1294  real(esmf_kind_r8), allocatable :: phalf(:)
1295  real(esmf_kind_r8), allocatable :: pres_interface(:)
1296  real(kind=4), allocatable :: dummy3d(:,:,:)
1297  real(kind=4), allocatable :: dummy3dall(:,:,:)
1298  real(esmf_kind_r8), allocatable :: dummy3dflip(:,:,:)
1299  real(esmf_kind_r8), allocatable :: dummy(:,:)
1300  real(esmf_kind_r8), pointer :: presptr(:,:,:), dpresptr(:,:,:)
1301  real(esmf_kind_r8), pointer :: psptr(:,:)
1302 
1303  print*,"- READ INPUT ATMOS DATA FROM GAUSSIAN NETCDF FILE."
1304 
1305  tilefile = trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
1306  error=nf90_open(trim(tilefile),nf90_nowrite,ncid)
1307  call netcdf_err(error, 'opening: '//trim(tilefile) )
1308 
1309  error=nf90_inq_dimid(ncid, 'grid_xt', id_dim)
1310  call netcdf_err(error, 'reading grid_xt id' )
1311  error=nf90_inquire_dimension(ncid,id_dim,len=idim_input)
1312  call netcdf_err(error, 'reading grid_xt value' )
1313 
1314  error=nf90_inq_dimid(ncid, 'grid_yt', id_dim)
1315  call netcdf_err(error, 'reading grid_yt id' )
1316  error=nf90_inquire_dimension(ncid,id_dim,len=jdim_input)
1317  call netcdf_err(error, 'reading grid_yt value' )
1318 
1319  if (idim_input /= i_input .or. jdim_input /= j_input) then
1320  call error_handler("DIMENSION MISMATCH BETWEEN SFC AND OROG FILES.", 2)
1321  endif
1322 
1323  error=nf90_inq_dimid(ncid, 'pfull', id_dim)
1324  call netcdf_err(error, 'reading pfull id' )
1325  error=nf90_inquire_dimension(ncid,id_dim,len=lev_input)
1326  call netcdf_err(error, 'reading pfull value' )
1327 
1328  error=nf90_inq_dimid(ncid, 'phalf', id_dim)
1329  call netcdf_err(error, 'reading phalf id' )
1330  error=nf90_inquire_dimension(ncid,id_dim,len=levp1_input)
1331  call netcdf_err(error, 'reading phalf value' )
1332  allocate(phalf(levp1_input))
1333  error=nf90_inq_varid(ncid, 'phalf', id_var)
1334  call netcdf_err(error, 'getting phalf varid' )
1335  error=nf90_get_var(ncid, id_var, phalf)
1336  call netcdf_err(error, 'reading phalf varid' )
1337 
1338  error=nf90_get_att(ncid, nf90_global, 'ncnsto', num_tracers_file)
1339  call netcdf_err(error, 'reading ntracer value' )
1340 
1341  call mpi_comm_size(mpi_comm_world, nprocs, error)
1342  print*,'- Running with ', nprocs, ' processors'
1343 
1344  call mpi_comm_rank(mpi_comm_world, myrank, error)
1345  print*,'- myrank/localpet is ',myrank,localpet
1346 
1347  max_procs = nprocs
1348  if (nprocs > lev_input) then
1349  max_procs = lev_input
1350  endif
1351 
1352  kdim = lev_input / max_procs
1353  remainder = lev_input - (max_procs*kdim)
1354 
1355  allocate(kcount(0:nprocs-1))
1356  kcount=0
1357  allocate(startk(0:nprocs-1))
1358  startk=0
1359  allocate(displ(0:nprocs-1))
1360  displ=0
1361  allocate(ircnt(0:nprocs-1))
1362  ircnt=0
1363 
1364  do k = 0, max_procs-2
1365  kcount(k) = kdim
1366  enddo
1367  kcount(max_procs-1) = kdim + remainder
1368 
1369  startk(0) = 1
1370  do k = 1, max_procs-1
1371  startk(k) = startk(k-1) + kcount(k-1)
1372  enddo
1373 
1374  ircnt(:) = idim_input * jdim_input * kcount(:)
1375 
1376  displ(0) = 0
1377  do k = 1, max_procs-1
1378  displ(k) = displ(k-1) + ircnt(k-1)
1379  enddo
1380 
1381  iscnt=idim_input*jdim_input*kcount(myrank)
1382 
1383 ! Account for case if number of tasks exceeds the number of vert levels.
1384 
1385  if (myrank <= max_procs-1) then
1386  allocate(dummy3d(idim_input,jdim_input,kcount(myrank)))
1387  else
1388  allocate(dummy3d(0,0,0))
1389  endif
1390 
1391  if (myrank == 0) then
1392  allocate(dummy3dall(idim_input,jdim_input,lev_input))
1393  dummy3dall = 0.0
1394  allocate(dummy3dflip(idim_input,jdim_input,lev_input))
1395  dummy3dflip = 0.0
1396  allocate(dummy(idim_input,jdim_input))
1397  dummy = 0.0
1398  else
1399  allocate(dummy3dall(0,0,0))
1400  allocate(dummy3dflip(0,0,0))
1401  allocate(dummy(0,0))
1402  endif
1403 
1404 !---------------------------------------------------------------------------
1405 ! Initialize esmf atmospheric fields.
1406 !---------------------------------------------------------------------------
1407 
1408  call init_atm_esmf_fields
1409 
1410  print*,"- CALL FieldCreate FOR INPUT GRID DELTA PRESSURE."
1411  dpres_input_grid = esmf_fieldcreate(input_grid, &
1412  typekind=esmf_typekind_r8, &
1413  staggerloc=esmf_staggerloc_center, &
1414  ungriddedlbound=(/1/), &
1415  ungriddedubound=(/lev_input/), rc=rc)
1416  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1417  call error_handler("IN FieldCreate", rc)
1418 
1419 ! Temperature
1420 
1421  if (myrank <= max_procs-1) then
1422  start = (/1,1,startk(myrank)/)
1423  count = (/idim_input,jdim_input,kcount(myrank)/)
1424  error=nf90_inq_varid(ncid, 'tmp', id_var)
1425  call netcdf_err(error, 'reading tmp field id' )
1426  error=nf90_get_var(ncid, id_var, dummy3d, start=start, count=count)
1427  call netcdf_err(error, 'reading tmp field' )
1428  endif
1429 
1430  call mpi_gatherv(dummy3d, iscnt, mpi_real, &
1431  dummy3dall, ircnt, displ, mpi_real, &
1432  0, mpi_comm_world, error)
1433  if (error /= 0) call error_handler("IN mpi_gatherv of temperature", error)
1434 
1435  if (myrank == 0) then
1436  dummy3dflip(:,:,1:lev_input) = dummy3dall(:,:,lev_input:1:-1)
1437  endif
1438 
1439  print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE "
1440  call esmf_fieldscatter(temp_input_grid, dummy3dflip, rootpet=0, rc=rc)
1441  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1442  call error_handler("IN FieldScatter", rc)
1443 
1444 ! dpres
1445 
1446  if (myrank <= max_procs-1) then
1447  error=nf90_inq_varid(ncid, 'dpres', id_var)
1448  call netcdf_err(error, 'reading dpres field id' )
1449  error=nf90_get_var(ncid, id_var, dummy3d, start=start, count=count)
1450  call netcdf_err(error, 'reading dpres field' )
1451  endif
1452 
1453  call mpi_gatherv(dummy3d, iscnt, mpi_real, &
1454  dummy3dall, ircnt, displ, mpi_real, &
1455  0, mpi_comm_world, error)
1456  if (error /= 0) call error_handler("IN mpi_gatherv of dpres", error)
1457 
1458  if (myrank == 0) then
1459  dummy3dflip(:,:,1:lev_input) = dummy3dall(:,:,lev_input:1:-1)
1460  endif
1461 
1462  print*,"- CALL FieldScatter FOR INPUT GRID DPRES "
1463  call esmf_fieldscatter(dpres_input_grid, dummy3dflip, rootpet=0, rc=rc)
1464  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1465  call error_handler("IN FieldScatter", rc)
1466 
1467 ! ugrd
1468 
1469  if (myrank <= max_procs-1) then
1470  error=nf90_inq_varid(ncid, 'ugrd', id_var)
1471  call netcdf_err(error, 'reading ugrd field id' )
1472  error=nf90_get_var(ncid, id_var, dummy3d, start=start, count=count)
1473  call netcdf_err(error, 'reading ugrd field' )
1474  endif
1475 
1476  call mpi_gatherv(dummy3d, iscnt, mpi_real, &
1477  dummy3dall, ircnt, displ, mpi_real, &
1478  0, mpi_comm_world, error)
1479  if (error /= 0) call error_handler("IN mpi_gatherv of ugrd", error)
1480 
1481  if (myrank == 0) then
1482  dummy3dflip(:,:,1:lev_input) = dummy3dall(:,:,lev_input:1:-1)
1483  endif
1484 
1485  print*,"- CALL FieldScatter FOR INPUT GRID UGRD "
1486  call esmf_fieldscatter(u_input_grid, dummy3dflip, rootpet=0, rc=rc)
1487  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1488  call error_handler("IN FieldScatter", rc)
1489 
1490 ! vgrd
1491 
1492  if (myrank <= max_procs-1) then
1493  error=nf90_inq_varid(ncid, 'vgrd', id_var)
1494  call netcdf_err(error, 'reading vgrd field id' )
1495  error=nf90_get_var(ncid, id_var, dummy3d, start=start, count=count)
1496  call netcdf_err(error, 'reading vgrd field' )
1497  endif
1498 
1499  call mpi_gatherv(dummy3d, iscnt, mpi_real, &
1500  dummy3dall, ircnt, displ, mpi_real, &
1501  0, mpi_comm_world, error)
1502  if (error /= 0) call error_handler("IN mpi_gatherv of vgrd", error)
1503 
1504  if (myrank == 0) then
1505  dummy3dflip(:,:,1:lev_input) = dummy3dall(:,:,lev_input:1:-1)
1506  endif
1507 
1508  print*,"- CALL FieldScatter FOR INPUT GRID VGRD "
1509  call esmf_fieldscatter(v_input_grid, dummy3dflip, rootpet=0, rc=rc)
1510  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1511  call error_handler("IN FieldScatter", rc)
1512 
1513 ! tracers
1514 
1515  do n = 1, num_tracers_input
1516 
1517  if (myrank <= max_procs-1) then
1518  error=nf90_inq_varid(ncid, tracers_input(n), id_var)
1519  call netcdf_err(error, 'reading tracer field id' )
1520  error=nf90_get_var(ncid, id_var, dummy3d, start=start, count=count)
1521  call netcdf_err(error, 'reading tracer field' )
1522  endif
1523 
1524  call mpi_gatherv(dummy3d, iscnt, mpi_real, &
1525  dummy3dall, ircnt, displ, mpi_real, &
1526  0, mpi_comm_world, error)
1527  if (error /= 0) call error_handler("IN mpi_gatherv of tracer", error)
1528 
1529  if (myrank == 0) then
1530  dummy3dflip(:,:,1:lev_input) = dummy3dall(:,:,lev_input:1:-1)
1531  where(dummy3dflip < 0.0) dummy3dflip = 0.0
1532  endif
1533 
1534  print*,"- CALL FieldScatter FOR INPUT GRID ", tracers_input(n)
1535  call esmf_fieldscatter(tracers_input_grid(n), dummy3dflip, rootpet=0, rc=rc)
1536  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1537  call error_handler("IN FieldScatter", rc)
1538 
1539  enddo
1540 
1541 ! dzdt set to zero for now.
1542 
1543  if (myrank == 0) then
1544  dummy3dflip = 0.0
1545  endif
1546 
1547  print*,"- CALL FieldScatter FOR INPUT GRID DZDT"
1548  call esmf_fieldscatter(dzdt_input_grid, dummy3dflip, rootpet=0, rc=rc)
1549  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1550  call error_handler("IN FieldScatter", rc)
1551 
1552  deallocate(dummy3dflip, dummy3dall, dummy3d)
1553 
1554 ! terrain
1555 
1556  if (myrank==0) then
1557  print*,"- READ TERRAIN."
1558  error=nf90_inq_varid(ncid, 'hgtsfc', id_var)
1559  call netcdf_err(error, 'reading hgtsfc field id' )
1560  error=nf90_get_var(ncid, id_var, dummy)
1561  call netcdf_err(error, 'reading hgtsfc field' )
1562  endif
1563 
1564  print*,"- CALL FieldScatter FOR INPUT GRID TERRAIN."
1565  call esmf_fieldscatter(terrain_input_grid, dummy, rootpet=0, rc=rc)
1566  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1567  call error_handler("IN FieldScatter", rc)
1568 
1569 ! surface pressure
1570 
1571  if (myrank==0) then
1572  print*,"- READ SURFACE P."
1573  error=nf90_inq_varid(ncid, 'pressfc', id_var)
1574  call netcdf_err(error, 'reading pressfc field id' )
1575  error=nf90_get_var(ncid, id_var, dummy)
1576  call netcdf_err(error, 'reading pressfc field' )
1577  endif
1578 
1579  print*,"- CALL FieldScatter FOR INPUT GRID SURFACE P."
1580  call esmf_fieldscatter(ps_input_grid, dummy, rootpet=0, rc=rc)
1581  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1582  call error_handler("IN FieldScatter", rc)
1583 
1584  deallocate(kcount, startk, displ, ircnt, dummy)
1585 
1586 !---------------------------------------------------------------------------
1587 ! Convert from 2-d to 3-d cartesian winds.
1588 !---------------------------------------------------------------------------
1589 
1590  call convert_winds
1591 
1592 !---------------------------------------------------------------------------
1593 ! Compute pressure.
1594 !---------------------------------------------------------------------------
1595 
1596  print*,"- CALL FieldGet FOR PRESSURE."
1597  call esmf_fieldget(pres_input_grid, &
1598  computationallbound=clb, &
1599  computationalubound=cub, &
1600  farrayptr=presptr, rc=rc)
1601  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1602  call error_handler("IN FieldGet", rc)
1603 
1604  print*,"- CALL FieldGet FOR DELTA PRESSURE."
1605  call esmf_fieldget(dpres_input_grid, &
1606  farrayptr=dpresptr, rc=rc)
1607  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1608  call error_handler("IN FieldGet", rc)
1609 
1610  print*,"- CALL FieldGet FOR SURFACE PRESSURE."
1611  call esmf_fieldget(ps_input_grid, &
1612  farrayptr=psptr, rc=rc)
1613  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1614  call error_handler("IN FieldGet", rc)
1615 
1616  allocate(pres_interface(levp1_input))
1617 
1618 !---------------------------------------------------------------------------
1619 ! Compute 3-d pressure.
1620 !---------------------------------------------------------------------------
1621 
1622 !---------------------------------------------------------------------------
1623 ! When ingesting gaussian netcdf files, the mid-layer
1624 ! surface pressure are computed top down from delta-p
1625 ! The surface pressure in the file is not used. According
1626 ! to Jun Wang, after the model's write component interpolates from the
1627 ! cubed-sphere grid to the gaussian grid, the surface pressure is
1628 ! no longer consistent with the delta p.
1629 !---------------------------------------------------------------------------
1630 
1631  do i = clb(1), cub(1)
1632  do j = clb(2), cub(2)
1633  pres_interface(levp1_input) = phalf(1) * 100.0_8
1634  do k = lev_input, 1, -1
1635  pres_interface(k) = pres_interface(k+1) + dpresptr(i,j,k)
1636  enddo
1637  psptr(i,j) = pres_interface(1)
1638  do k = 1, lev_input
1639  presptr(i,j,k) = (pres_interface(k) + pres_interface(k+1)) / 2.0_8
1640  enddo
1641  enddo
1642  enddo
1643 
1644  deallocate(pres_interface, phalf)
1645 
1646  call esmf_fielddestroy(dpres_input_grid, rc=rc)
1647 
1648  end subroutine read_input_atm_gaussian_netcdf_file
1649 
1659  subroutine read_input_atm_tiled_history_file(localpet)
1660 
1661  use mpi
1662 
1663  implicit none
1664 
1665  integer, intent(in) :: localpet
1666 
1667  character(len=500) :: tilefile
1668 
1669  integer :: error, ncid, rc, tile
1670  integer :: id_dim, idim_input, jdim_input
1671  integer :: id_var, i, j, k, n
1672  integer :: clb(3), cub(3), num_tracers_file
1673 
1674  real(esmf_kind_r8), allocatable :: data_one_tile(:,:)
1675  real(esmf_kind_r8), allocatable :: data_one_tile_3d(:,:,:)
1676  real(esmf_kind_r8), pointer :: presptr(:,:,:), dpresptr(:,:,:)
1677  real(esmf_kind_r8), pointer :: psptr(:,:)
1678  real(esmf_kind_r8), allocatable :: pres_interface(:), phalf(:)
1679 
1680  print*,"- READ INPUT ATMOS DATA FROM TILED HISTORY FILES."
1681 
1682  tilefile = trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
1683  error=nf90_open(trim(tilefile),nf90_nowrite,ncid)
1684  call netcdf_err(error, 'opening: '//trim(tilefile) )
1685 
1686  error=nf90_inq_dimid(ncid, 'grid_xt', id_dim)
1687  call netcdf_err(error, 'reading grid_xt id' )
1688  error=nf90_inquire_dimension(ncid,id_dim,len=idim_input)
1689  call netcdf_err(error, 'reading grid_xt value' )
1690 
1691  error=nf90_inq_dimid(ncid, 'grid_yt', id_dim)
1692  call netcdf_err(error, 'reading grid_yt id' )
1693  error=nf90_inquire_dimension(ncid,id_dim,len=jdim_input)
1694  call netcdf_err(error, 'reading grid_yt value' )
1695 
1696  if (idim_input /= i_input .or. jdim_input /= j_input) then
1697  call error_handler("DIMENSION MISMATCH BETWEEN SFC AND OROG FILES.", 2)
1698  endif
1699 
1700  error=nf90_inq_dimid(ncid, 'pfull', id_dim)
1701  call netcdf_err(error, 'reading pfull id' )
1702  error=nf90_inquire_dimension(ncid,id_dim,len=lev_input)
1703  call netcdf_err(error, 'reading pfull value' )
1704 
1705  error=nf90_inq_dimid(ncid, 'phalf', id_dim)
1706  call netcdf_err(error, 'reading phalf id' )
1707  error=nf90_inquire_dimension(ncid,id_dim,len=levp1_input)
1708  call netcdf_err(error, 'reading phalf value' )
1709  allocate(phalf(levp1_input))
1710  error=nf90_inq_varid(ncid, 'phalf', id_var)
1711  call netcdf_err(error, 'getting phalf varid' )
1712  error=nf90_get_var(ncid, id_var, phalf)
1713  call netcdf_err(error, 'reading phalf varid' )
1714 
1715  error=nf90_get_att(ncid, nf90_global, 'ncnsto', num_tracers_file)
1716  call netcdf_err(error, 'reading ntracer value' )
1717 
1718  error = nf90_close(ncid)
1719 
1720  print*,'- FILE HAS ', num_tracers_file, ' TRACERS.'
1721  print*,'- WILL PROCESS ', num_tracers_input, ' TRACERS.'
1722 
1723 !---------------------------------------------------------------------------
1724 ! Initialize esmf atmospheric fields.
1725 !---------------------------------------------------------------------------
1726 
1727  call init_atm_esmf_fields
1728 
1729  print*,"- CALL FieldCreate FOR INPUT GRID DELTA PRESSURE."
1730  dpres_input_grid = esmf_fieldcreate(input_grid, &
1731  typekind=esmf_typekind_r8, &
1732  staggerloc=esmf_staggerloc_center, &
1733  ungriddedlbound=(/1/), &
1734  ungriddedubound=(/lev_input/), rc=rc)
1735  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1736  call error_handler("IN FieldCreate", rc)
1737 
1738  if (localpet < num_tiles_input_grid) then
1739  allocate(data_one_tile(i_input,j_input))
1740  allocate(data_one_tile_3d(i_input,j_input,lev_input))
1741  else
1742  allocate(data_one_tile(0,0))
1743  allocate(data_one_tile_3d(0,0,0))
1744  endif
1745 
1746  if (localpet < num_tiles_input_grid) then
1747  tile = localpet+1
1748  tilefile= trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(tile))
1749  print*,"- READ ATMOSPHERIC DATA FROM: ", trim(tilefile)
1750  error=nf90_open(trim(tilefile),nf90_nowrite,ncid)
1751  call netcdf_err(error, 'opening: '//trim(tilefile) )
1752  endif
1753 
1754  if (localpet < num_tiles_input_grid) then
1755 ! print*,"- READ VERTICAL VELOCITY."
1756 ! error=nf90_inq_varid(ncid, 'dzdt', id_var)
1757 ! call netcdf_err(error, 'reading field id' )
1758 ! error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1759 ! call netcdf_err(error, 'reading field' )
1760 ! data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1761 
1762 ! Using w from the tiled history files has caused problems.
1763 ! Set to zero.
1764  data_one_tile_3d = 0.0_8
1765  endif
1766 
1767  do tile = 1, num_tiles_input_grid
1768  print*,"- CALL FieldScatter FOR INPUT GRID VERTICAL VELOCITY."
1769  call esmf_fieldscatter(dzdt_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1770  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1771  call error_handler("IN FieldScatter", rc)
1772  enddo
1773 
1774  do n = 1, num_tracers_input
1775 
1776  if (localpet < num_tiles_input_grid) then
1777  print*,"- READ ", trim(tracers_input(n))
1778  error=nf90_inq_varid(ncid, tracers_input(n), id_var)
1779  call netcdf_err(error, 'reading field id' )
1780  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1781  call netcdf_err(error, 'reading field' )
1782  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1783  endif
1784 
1785  do tile = 1, num_tiles_input_grid
1786  print*,"- CALL FieldScatter FOR INPUT GRID TRACER ", trim(tracers_input(n))
1787  call esmf_fieldscatter(tracers_input_grid(n), data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1788  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1789  call error_handler("IN FieldScatter", rc)
1790  enddo
1791 
1792  enddo
1793 
1794  if (localpet < num_tiles_input_grid) then
1795  print*,"- READ TEMPERATURE."
1796  error=nf90_inq_varid(ncid, 'tmp', id_var)
1797  call netcdf_err(error, 'reading field id' )
1798  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1799  call netcdf_err(error, 'reading field' )
1800  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1801  endif
1802 
1803  do tile = 1, num_tiles_input_grid
1804  print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE."
1805  call esmf_fieldscatter(temp_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1806  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1807  call error_handler("IN FieldScatter", rc)
1808  enddo
1809 
1810  if (localpet < num_tiles_input_grid) then
1811  print*,"- READ U-WIND."
1812  error=nf90_inq_varid(ncid, 'ugrd', id_var)
1813  call netcdf_err(error, 'reading field id' )
1814  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1815  call netcdf_err(error, 'reading field' )
1816  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1817  endif
1818 
1819  do tile = 1, num_tiles_input_grid
1820  print*,"- CALL FieldScatter FOR INPUT GRID U."
1821  call esmf_fieldscatter(u_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1822  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1823  call error_handler("IN FieldScatter", rc)
1824  enddo
1825 
1826  if (localpet < num_tiles_input_grid) then
1827  print*,"- READ V-WIND."
1828  error=nf90_inq_varid(ncid, 'vgrd', id_var)
1829  call netcdf_err(error, 'reading field id' )
1830  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1831  call netcdf_err(error, 'reading field' )
1832  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1833  endif
1834 
1835  do tile = 1, num_tiles_input_grid
1836  print*,"- CALL FieldScatter FOR INPUT GRID V."
1837  call esmf_fieldscatter(v_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1838  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1839  call error_handler("IN FieldScatter", rc)
1840  enddo
1841 
1842  if (localpet < num_tiles_input_grid) then
1843  print*,"- READ SURFACE PRESSURE."
1844  error=nf90_inq_varid(ncid, 'pressfc', id_var)
1845  call netcdf_err(error, 'reading field id' )
1846  error=nf90_get_var(ncid, id_var, data_one_tile)
1847  call netcdf_err(error, 'reading field' )
1848  endif
1849 
1850  do tile = 1, num_tiles_input_grid
1851  print*,"- CALL FieldScatter FOR INPUT GRID SURFACE PRESSURE."
1852  call esmf_fieldscatter(ps_input_grid, data_one_tile, rootpet=tile-1, tile=tile, rc=rc)
1853  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1854  call error_handler("IN FieldScatter", rc)
1855  enddo
1856 
1857  if (localpet < num_tiles_input_grid) then
1858  print*,"- READ TERRAIN."
1859  error=nf90_inq_varid(ncid, 'hgtsfc', id_var)
1860  call netcdf_err(error, 'reading field id' )
1861  error=nf90_get_var(ncid, id_var, data_one_tile)
1862  call netcdf_err(error, 'reading field' )
1863  endif
1864 
1865  do tile = 1, num_tiles_input_grid
1866  print*,"- CALL FieldScatter FOR INPUT GRID TERRAIN."
1867  call esmf_fieldscatter(terrain_input_grid, data_one_tile, rootpet=tile-1, tile=tile, rc=rc)
1868  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1869  call error_handler("IN FieldScatter", rc)
1870  enddo
1871 
1872  if (localpet < num_tiles_input_grid) then
1873  print*,"- READ DELTA PRESSURE."
1874  error=nf90_inq_varid(ncid, 'dpres', id_var)
1875  call netcdf_err(error, 'reading field id' )
1876  error=nf90_get_var(ncid, id_var, data_one_tile_3d)
1877  call netcdf_err(error, 'reading field' )
1878  data_one_tile_3d(:,:,1:lev_input) = data_one_tile_3d(:,:,lev_input:1:-1)
1879  endif
1880 
1881  do tile = 1, num_tiles_input_grid
1882  print*,"- CALL FieldScatter FOR INPUT DELTA PRESSURE."
1883  call esmf_fieldscatter(dpres_input_grid, data_one_tile_3d, rootpet=tile-1, tile=tile, rc=rc)
1884  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1885  call error_handler("IN FieldScatter", rc)
1886  enddo
1887 
1888  if (localpet < num_tiles_input_grid) error = nf90_close(ncid)
1889 
1890  deallocate(data_one_tile_3d, data_one_tile)
1891 
1892 !---------------------------------------------------------------------------
1893 ! Convert from 2-d to 3-d cartesian winds.
1894 !---------------------------------------------------------------------------
1895 
1896  call convert_winds
1897 
1898 !---------------------------------------------------------------------------
1899 ! Compute pressure.
1900 !---------------------------------------------------------------------------
1901 
1902  print*,"- CALL FieldGet FOR PRESSURE."
1903  call esmf_fieldget(pres_input_grid, &
1904  computationallbound=clb, &
1905  computationalubound=cub, &
1906  farrayptr=presptr, rc=rc)
1907  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1908  call error_handler("IN FieldGet", rc)
1909 
1910  print*,"- CALL FieldGet FOR DELTA PRESSURE."
1911  call esmf_fieldget(dpres_input_grid, &
1912  farrayptr=dpresptr, rc=rc)
1913  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1914  call error_handler("IN FieldGet", rc)
1915 
1916  print*,"- CALL FieldGet FOR SURFACE PRESSURE."
1917  call esmf_fieldget(ps_input_grid, &
1918  farrayptr=psptr, rc=rc)
1919  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1920  call error_handler("IN FieldGet", rc)
1921 
1922  allocate(pres_interface(levp1_input))
1923 
1924 !---------------------------------------------------------------------------
1925 ! Compute 3-d pressure.
1926 !---------------------------------------------------------------------------
1927 
1928  do i = clb(1), cub(1)
1929  do j = clb(2), cub(2)
1930  pres_interface(1) = psptr(i,j)
1931  do k = 2, levp1_input
1932  pres_interface(k) = pres_interface(k-1) - dpresptr(i,j,k-1)
1933  enddo
1934  do k = 1, lev_input
1935  presptr(i,j,k) = (pres_interface(k) + pres_interface(k+1)) / 2.0_8
1936  enddo
1937  enddo
1938  enddo
1939 
1940  deallocate(pres_interface, phalf)
1941 
1942  call esmf_fielddestroy(dpres_input_grid, rc=rc)
1943 
1944  end subroutine read_input_atm_tiled_history_file
1945 
1950  subroutine read_input_atm_grib2_file(localpet)
1951 
1952  use mpi
1953  use grib_mod
1954 
1955  use grib2_util, only : rh2spfh, rh2spfh_gfs, convert_omega
1956 
1957  implicit none
1958 
1959  integer, intent(in) :: localpet
1960 
1961  integer, parameter :: ntrac_max=14
1962  integer, parameter :: max_levs=1000
1963 
1964  character(len=300) :: the_file
1965  character(len=20) :: vname, &
1966  trac_names_vmap(ntrac_max), &
1967  tmpstr, &
1968  method, tracers_input_vmap(num_tracers_input), &
1969  tracers_default(ntrac_max)
1970 
1971  integer :: i, j, k, n
1972  integer :: ii,jj
1973  integer :: rc, clb(3), cub(3)
1974  integer :: vlev, iret,varnum, o3n, pdt_num
1975  integer :: intrp_ier, done_print
1976  integer :: trac_names_oct10(ntrac_max)
1977  integer :: tracers_input_oct10(num_tracers_input)
1978  integer :: trac_names_oct11(ntrac_max)
1979  integer :: tracers_input_oct11(num_tracers_input)
1980  integer :: lugb, lugi, jdisc, jpdt(200), jgdt(200), iscale
1981  integer :: jids(200), jpdtn, jgdtn, octet_23, octet_29
1982  integer :: count_spfh, count_rh, count_icmr, count_scliwc
1983  integer :: count_cice, count_rwmr, count_scllwc, count
1984 
1985  logical :: conv_omega=.false., &
1986  hasspfh=.true., &
1987  isnative=.false., &
1988  use_rh=.false. , unpack, &
1989  all_empty, is_missing
1990 
1991  real(esmf_kind_r8), allocatable :: dum2d_1(:,:)
1992 
1993 
1994  real(esmf_kind_r8) :: rlevs_hold(max_levs)
1995  real(esmf_kind_r8), allocatable :: rlevs(:)
1996  real(esmf_kind_r4), allocatable :: dummy2d(:,:)
1997  real(esmf_kind_r8), allocatable :: dummy3d(:,:,:), dummy2d_8(:,:),&
1998  u_tmp_3d(:,:,:), v_tmp_3d(:,:,:)
1999  real(esmf_kind_r8), pointer :: presptr(:,:,:), psptr(:,:),tptr(:,:,:), &
2000  qptr(:,:,:), wptr(:,:,:), &
2001  uptr(:,:,:), vptr(:,:,:)
2002  real(esmf_kind_r4) :: value
2003  real(esmf_kind_r8), parameter :: p0 = 100000.0
2004  real(esmf_kind_r8), allocatable :: dummy3d_col_in(:),dummy3d_col_out(:)
2005  real(esmf_kind_r8), parameter :: intrp_missing = -999.0
2006  real(esmf_kind_r4), parameter :: lev_no_tr_fill = 20000.0
2007  real(esmf_kind_r4), parameter :: lev_no_o3_fill = 40000.0
2008 
2009  type(gribfield) :: gfld
2010 
2011  tracers(:) = "NULL"
2012 
2013  trac_names_oct10 = (/1, 1, 14, 1, 1, 1, 1, 6, 6, 1, 6, 13, 13, 2 /)
2014  trac_names_oct11 = (/0, 22, 192, 23, 24, 25, 32, 1, 29, 100, 28, 193, 192, 2 /)
2015 
2016  trac_names_vmap = (/"sphum ", "liq_wat ", "o3mr ", "ice_wat ", &
2017  "rainwat ", "snowwat ", "graupel ", "cld_amt ", "ice_nc ", &
2018  "rain_nc ", "water_nc", "liq_aero", "ice_aero", &
2019  "sgs_tke "/)
2020 
2021  tracers_default = (/"sphum ", "liq_wat ", "o3mr ", "ice_wat ", &
2022  "rainwat ", "snowwat ", "graupel ", "cld_amt ", "ice_nc ", &
2023  "rain_nc ", "water_nc", "liq_aero", "ice_aero", &
2024  "sgs_tke "/)
2025 
2026  the_file = trim(data_dir_input_grid) // "/" // trim(grib2_file_input_grid)
2027 
2028  print*,"- READ ATMOS DATA FROM GRIB2 FILE: ", trim(the_file)
2029 
2030  if (localpet == 0) then
2031 
2032  lugb=14
2033  lugi=0
2034  call baopenr(lugb,the_file,iret)
2035  if (iret /= 0) call error_handler("ERROR OPENING GRIB2 FILE.", iret)
2036 
2037  jdisc = 0 ! Search for discipline - meteorological products
2038  j = 0 ! Search at beginning of file.
2039  jpdt = -9999 ! Array of values in product definition template, set to wildcard
2040  jids = -9999 ! Array of values in identification section, set to wildcard
2041  jgdt = -9999 ! Array of values in grid definition template, set to wildcard
2042  jgdtn = -1 ! Search for any grid definition number.
2043  jpdtn = -1 ! Search for any product definition template number.
2044  unpack =.false.
2045 
2046  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2047  unpack, k, gfld, iret)
2048 
2049 !----------------------------------------------------------------------
2050 ! Read first record and check if this is NCEP GEFS data.
2051 ! This will determine what product definition template number to
2052 ! search for (Section 4/Octets 8-9).
2053 !
2054 ! Section 1/Octets 6-7 is '7' (NCEP)
2055 ! Section 1/Octets 8-9 is '2' (NCEP Ensemble products).
2056 !----------------------------------------------------------------------
2057 
2058  if (iret == 0) then
2059  if (gfld%idsect(1) == 7 .and. gfld%idsect(2) == 2) then
2060  print*,'- THIS IS NCEP GEFS DATA.'
2061  pdt_num = 1 ! Search for product definition template number 1.
2062  ! Individual ensember forecast.
2063  else
2064  pdt_num = 0 ! Search for product definition template number 0.
2065  ! Analysis or forecast.
2066  endif
2067  else
2068  call error_handler("READING GRIB2 FILE", iret)
2069  endif
2070 
2071 !----------------------------------------------------------------------
2072 ! First, check for the vertical coordinate. If temperture at the 10 hybrid
2073 ! level is found, hybrid coordinates are assumed. Otherwise, data is on
2074 ! isobaric levels.
2075 !----------------------------------------------------------------------
2076 
2077  j = 0
2078  jpdtn = pdt_num ! Search for the specific product definition template number.
2079  jpdt(1) = 0 ! Sect4/oct 10 - Parameter category - temperature field
2080  jpdt(2) = 0 ! Sect4/oct 11 - Parameter number - temperature
2081  jpdt(10) = 105 ! Sect4/oct 23 - Type of level - hybrid
2082  jpdt(12) = 10 ! Sect4/octs 25/28 - Value of hybrid level
2083  unpack=.false.
2084 
2085  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2086  unpack, k, gfld, iret)
2087 
2088  if (iret == 0) then
2089  print*,'- DATA IS ON HYBRID LEVELS.'
2090  octet_23 = 105 ! Section 4/Oct 23 - type of first fixed surface.
2091  octet_29 = 255 ! Section 4/Oct 29 - type of second fixed surface (N/A).
2092  isnative=.true.
2093  else
2094  print*,'- DATA IS ON ISOBARIC LEVELS.'
2095  octet_23 = 100 ! Section 4/Oct 23 - type of first fixed surface.
2096  octet_29 = 255 ! Section 4/Oct 29 - type of second fixed surface (N/A).
2097  isnative=.false.
2098  endif
2099 
2100 ! Now count the number of vertical levels by searching for u-wind.
2101 ! Store the value of each level.
2102 
2103  rlevs_hold = -999.9
2104  lev_input = 0
2105  iret = 0
2106  j = 0
2107  jpdtn = -1
2108  jpdt = -9999
2109 
2110  do
2111  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2112  unpack, k, gfld, iret)
2113 
2114  if (iret /= 0) exit
2115 
2116  if (gfld%discipline == 0) then ! Discipline - meteorological products
2117  if (gfld%ipdtnum == pdt_num) then ! Product definition template number.
2118  if (gfld%ipdtmpl(1) == 2 .and. gfld%ipdtmpl(2) == 2) then ! u-wind
2119  ! Sect4/octs 10 and 11.
2120  if (gfld%ipdtmpl(10) == octet_23 .and. gfld%ipdtmpl(13) == octet_29) then
2121  ! Sect4 octs 23 and 29.
2122  ! Hybrid or isobaric.
2123  lev_input = lev_input + 1
2124  iscale = 10 ** gfld%ipdtmpl(11)
2125  rlevs_hold(lev_input) = float(gfld%ipdtmpl(12))/float(iscale)
2126  endif
2127  endif
2128  endif
2129  endif
2130 
2131  j = k
2132  enddo
2133 
2134  endif ! read file on task 0.
2135 
2136  call mpi_barrier(mpi_comm_world, iret)
2137  call mpi_bcast(isnative,1,mpi_logical,0,mpi_comm_world,iret)
2138  call mpi_bcast(lev_input,1,mpi_integer,0,mpi_comm_world,iret)
2139  call mpi_bcast(pdt_num,1,mpi_integer,0,mpi_comm_world,iret)
2140  call mpi_bcast(rlevs_hold, max_levs, mpi_integer,0,mpi_comm_world,iret)
2141 
2142  allocate(slevs(lev_input))
2143  allocate(rlevs(lev_input))
2144  allocate(dummy3d_col_in(lev_input))
2145  allocate(dummy3d_col_out(lev_input))
2146 
2147  levp1_input = lev_input + 1
2148 
2149 ! Jili Dong add sort to re-order isobaric levels.
2150 
2151  do i = 1, lev_input
2152  rlevs(i) = rlevs_hold(i)
2153  enddo
2154 
2155  call quicksort(rlevs,1,lev_input)
2156 
2157  do i = 1, lev_input
2158  if (isnative) then
2159  write(slevs(i), '(i6)') nint(rlevs(i))
2160  slevs(i) = trim(slevs(i)) // " hybrid"
2161  else
2162  write(slevs(i), '(f11.2)') rlevs(i)
2163  slevs(i) = trim(slevs(i)) // " Pa"
2164  endif
2165  enddo
2166 
2167  if(localpet == 0) then
2168  do i = 1,lev_input
2169  print*, "- LEVEL AFTER SORT = ",trim(slevs(i))
2170  enddo
2171  endif
2172 
2173 ! Check to see if specfic humidity exists at all the same levels as ugrd.
2174 
2175  if (localpet == 0) then
2176 
2177  jpdtn = pdt_num ! Product definition template number.
2178  jpdt = -9999
2179  jpdt(1) = 1 ! Sect4/oct 10 - Parameter category - moisture
2180  jpdt(2) = 0 ! Sect4/oct 11 - Parameter number - specific humidity
2181  jpdt(10) = octet_23 ! Sect4/oct 23 - type of level.
2182  unpack=.false.
2183 
2184  count_spfh=0
2185 
2186  do vlev = 1, lev_input
2187  j = 0
2188  jpdt(12) = nint(rlevs(vlev))
2189 
2190  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2191  unpack, k, gfld, iret)
2192 
2193  if (iret == 0) then
2194  count_spfh = count_spfh + 1
2195  endif
2196  enddo
2197 
2198  jpdt(1) = 1 ! Sec4/oct 10 - Parameter category - moisture
2199  jpdt(2) = 1 ! Sec4/oct 11 - Parameter number - rel humidity
2200  count_rh=0
2201 
2202  do vlev = 1, lev_input
2203  j = 0
2204  jpdt(12) = nint(rlevs(vlev))
2205 
2206  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2207  unpack, k, gfld, iret)
2208 
2209  if (iret == 0) then
2210  count_rh = count_rh + 1
2211  endif
2212  enddo
2213 
2214  if (count_spfh /= lev_input) then
2215  use_rh = .true.
2216  endif
2217 
2218  if (count_spfh == 0 .or. use_rh) then
2219  if (count_rh == 0) then
2220  call error_handler("READING ATMOSPHERIC WATER VAPOR VARIABLE.", 2)
2221  endif
2222  hasspfh = .false. ! Will read rh and convert to specific humidity.
2223  trac_names_oct10(1) = 1
2224  trac_names_oct11(1) = 1
2225  print*,"- FILE CONTAINS RH."
2226  else
2227  print*,"- FILE CONTAINS SPFH."
2228  endif
2229 
2230  endif
2231 
2232  call mpi_barrier(mpi_comm_world, rc)
2233  call mpi_bcast(hasspfh,1,mpi_logical,0,mpi_comm_world,rc)
2234 
2235 ! Search for and count the number of tracers in the file.
2236 
2237  if (localpet == 0) then
2238 
2239  jpdtn = pdt_num ! Product definition template number.
2240  jpdt = -9999
2241  jpdt(10) = octet_23 ! Sect4/oct 23 - type of level.
2242  unpack=.false.
2243 
2244  count_icmr=0
2245  count_scliwc=0
2246  count_cice=0
2247  count_rwmr=0
2248  count_scllwc=0
2249 
2250  do vlev = 1, lev_input
2251 
2252  j = 0
2253  jpdt(1) = 1 ! Sect4/oct 10 - Parameter category - moisture
2254  jpdt(2) = 23 ! Sect4/oct 11 - Parameter number - ice water mixing ratio
2255  jpdt(12) = nint(rlevs(vlev))
2256 
2257  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2258  unpack, k, gfld, iret)
2259 
2260  if (iret == 0) then
2261  count_icmr = count_icmr + 1
2262  endif
2263 
2264  j = 0
2265  jpdt(1) = 1 ! Sect4/oct 10 - Parameter category - moisture
2266  jpdt(2) = 84 ! Sect4/oct 11 - Parameter number - cloud ice water content.
2267  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2268  unpack, k, gfld, iret)
2269 
2270  if (iret == 0) then
2271  count_scliwc = count_scliwc + 1
2272  endif
2273 
2274  j = 0
2275  jpdt(1) = 6 ! Sect4/oct 10 - Parameter category - clouds
2276  jpdt(2) = 0 ! Sect4/oct 11 - Parameter number - cloud ice
2277  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2278  unpack, k, gfld, iret)
2279 
2280  if (iret == 0) then
2281  count_cice = count_cice + 1
2282  endif
2283 
2284  j = 0
2285  jpdt(1) = 1 ! Sect4/oct 10 - Parameter category - moisture
2286  jpdt(2) = 24 ! Sect4/oct 11 - Parameter number - rain mixing ratio
2287  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2288  unpack, k, gfld, iret)
2289 
2290  if (iret == 0) then
2291  count_rwmr = count_rwmr + 1
2292  endif
2293 
2294  j = 0
2295  jpdt(1) = 1 ! Sect4/oct 10 - Parameter category - moisture
2296  jpdt(2) = 83 ! Sect4/oct 11 - Parameter number - specific cloud liquid
2297  ! water content.
2298  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2299  unpack, k, gfld, iret)
2300 
2301  if (iret == 0) then
2302  count_scllwc = count_scllwc + 1
2303  endif
2304 
2305  enddo
2306 
2307  if (count_icmr == 0) then
2308  if (count_scliwc == 0) then
2309  if (count_cice == 0) then
2310  print*,'- FILE DOES NOT CONTAIN CICE.'
2311  else
2312  trac_names_oct10(4) = 6 ! Sect4/oct 10 - Parameter category - clouds
2313  trac_names_oct11(4) = 0 ! Sect4/oct 11 - Parameter number - cloud ice
2314  print*,"- FILE CONTAINS CICE."
2315  endif
2316  else
2317  trac_names_oct10(4) = 1 ! Sect4/oct 10 - Parameter category - moisture
2318  trac_names_oct11(4) = 84 ! Sect4/oct 11 - Parameter number - cloud ice water content.
2319  print*,"- FILE CONTAINS SCLIWC."
2320  endif
2321  else
2322  print*,"- FILE CONTAINS ICMR."
2323  endif ! count of icmr
2324 
2325  if (count_rwmr == 0) then
2326  if (count_scllwc == 0) then
2327  print*,"- FILE DOES NOT CONTAIN SCLLWC."
2328  else
2329  trac_names_oct10(4) = 1 ! Sect4/oct 10 - Parameter category - moisture
2330  trac_names_oct11(4) = 83 ! Sect4/oct 11 - Parameter number - specific cloud liquid
2331  ! water content.
2332  print*,"- FILE CONTAINS SCLLWC."
2333  endif
2334  else
2335  print*,"- FILE CONTAINS CLWMR."
2336  endif
2337 
2338  endif ! count of tracers/localpet = 0
2339 
2340  call mpi_barrier(mpi_comm_world, rc)
2341  call mpi_bcast(trac_names_oct10,ntrac_max,mpi_integer,0,mpi_comm_world,rc)
2342  call mpi_bcast(trac_names_oct11,ntrac_max,mpi_integer,0,mpi_comm_world,rc)
2343 
2344  print*,"- COUNT NUMBER OF TRACERS TO BE READ IN BASED ON PHYSICS SUITE TABLE"
2345  do n = 1, num_tracers_input
2346 
2347  vname = tracers_input(n)
2348 
2349  i = maxloc(merge(1.,0.,trac_names_vmap == vname),dim=1)
2350 
2351  tracers_input_vmap(n)=trac_names_vmap(i)
2352  tracers(n)=tracers_default(i)
2353  if(trim(tracers(n)) .eq. "o3mr") o3n = n
2354 
2355  tracers_input_oct10(n) = trac_names_oct10(i)
2356  tracers_input_oct11(n) = trac_names_oct11(i)
2357 
2358  enddo
2359 
2360 !---------------------------------------------------------------------------
2361 ! Initialize esmf atmospheric fields.
2362 !---------------------------------------------------------------------------
2363 
2364  call init_atm_esmf_fields
2365 
2366  if (localpet == 0) then
2367  allocate(dummy2d(i_input,j_input))
2368  allocate(dummy2d_8(i_input,j_input))
2369  allocate(dummy3d(i_input,j_input,lev_input))
2370  allocate(dum2d_1(i_input,j_input))
2371  else
2372  allocate(dummy2d(0,0))
2373  allocate(dummy2d_8(0,0))
2374  allocate(dummy3d(0,0,0))
2375  allocate(dum2d_1(0,0))
2376  endif
2377 
2378 !----------------------------------------------------------------------------------
2379 ! This program expects field levels from bottom to top. Fields in non-native
2380 ! files read in from top to bottom. We will flip indices later. Fields on
2381 ! native vertical coordinates read from bottom to top so those need no adjustments.
2382 !----------------------------------------------------------------------------------
2383 
2384  if (localpet == 0) then
2385 
2386  print*,"- READ TEMPERATURE."
2387 
2388  jdisc = 0 ! search for discipline - meteorological products
2389  j = 0 ! search at beginning of file.
2390  jpdt = -9999 ! array of values in product definition template, set to wildcard
2391  jids = -9999 ! array of values in identification section, set to wildcard
2392  jgdt = -9999 ! array of values in grid definition template, set to wildcard
2393  jgdtn = -1 ! search for any grid definition number.
2394  jpdtn = pdt_num ! Search for specific product definition template number.
2395  jpdt(1) = 0 ! Sect 4/oct 10 - parameter category - temperature
2396  jpdt(2) = 0 ! Sect 4/oct 11 - parameter number - temperature
2397  jpdt(10) = octet_23 ! Sect4/oct 23 - type of level.
2398 
2399  unpack=.true.
2400 
2401  do vlev = 1, lev_input
2402 
2403  jpdt(12) = nint(rlevs(vlev))
2404 
2405  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2406  unpack, k, gfld, iret)
2407  if (iret /= 0) then
2408  call error_handler("READING IN TEMPERATURE AT LEVEL "//trim(slevs(vlev)),iret)
2409  endif
2410 
2411  dum2d_1 = reshape(gfld%fld, (/i_input,j_input/) )
2412 
2413  dummy3d(:,:,vlev) = dum2d_1
2414 
2415  enddo
2416 
2417  endif ! Read of temperature
2418 
2419  if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE."
2420  call esmf_fieldscatter(temp_input_grid, dummy3d, rootpet=0, rc=rc)
2421  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2422  call error_handler("IN FieldScatter", rc)
2423 
2424 ! Read tracers
2425 
2426  do n = 1, num_tracers_input
2427 
2428  if (localpet == 0) print*,"- READ ", trim(tracers_input_vmap(n))
2429 
2430  vname = tracers_input_vmap(n)
2431  call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, &
2432  this_field_var_name=tmpstr,loc=varnum)
2433 
2434  if (n==1 .and. .not. hasspfh) then
2435  print*,"- CALL FieldGather TEMPERATURE."
2436  call esmf_fieldgather(temp_input_grid,dummy3d,rootpet=0, tile=1, rc=rc)
2437  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2438  call error_handler("IN FieldGet", rc)
2439  endif
2440 
2441  if (localpet == 0) then
2442 
2443  jdisc = 0 ! search for discipline - meteorological products
2444  jpdt = -9999 ! array of values in product definition template, set to wildcard
2445  jids = -9999 ! array of values in identification section, set to wildcard
2446  jgdt = -9999 ! array of values in grid definition template, set to wildcard
2447  jgdtn = -1 ! search for any grid definition number.
2448  jpdtn = pdt_num ! Search for the product definition template number.
2449  jpdt(10) = octet_23 ! Sect4/oct 23 - type of level.
2450  unpack = .false.
2451 
2452  count = 0
2453 
2454  do vlev = 1, lev_input
2455 
2456  j = 0
2457  jpdt(1) = tracers_input_oct10(n)
2458  jpdt(2) = tracers_input_oct11(n)
2459  jpdt(12) = nint(rlevs(vlev))
2460 
2461  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2462  unpack, k, gfld, iret)
2463 
2464  if (iret == 0) then
2465  count = count + 1
2466  endif
2467 
2468  enddo
2469  iret=count
2470 
2471  ! Check to see if file has any data for this tracer
2472  if (iret == 0) then
2473  all_empty = .true.
2474  else
2475  all_empty = .false.
2476  endif
2477 
2478  is_missing = .false.
2479 
2480  do vlev = 1, lev_input
2481 
2482  unpack=.true.
2483  j = 0
2484  jpdt(1) = tracers_input_oct10(n)
2485  jpdt(2) = tracers_input_oct11(n)
2486  jpdt(12) = nint(rlevs(vlev) )
2487 
2488  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2489  unpack, k, gfld, iret)
2490 
2491  if (iret == 0) then ! found data
2492  dummy2d = real((reshape(gfld%fld, (/i_input,j_input/) )), kind=esmf_kind_r4)
2493  else ! did not find data.
2494  if (trim(method) .eq. 'intrp' .and. .not.all_empty) then
2495  dummy2d = intrp_missing
2496  is_missing = .true.
2497  else
2498  ! Abort if input data has some data for current tracer, but has
2499  ! missing data below 200 mb/ above 400mb
2500  if (.not.all_empty .and. n == o3n) then
2501  if (rlevs(vlev) .lt. lev_no_o3_fill) &
2502  call error_handler("TRACER "//trim(tracers(n))//" HAS MISSING DATA AT "//trim(slevs(vlev))//&
2503  ". SET MISSING VARIABLE CONDITION TO 'INTRP' TO AVOID THIS ERROR", 1)
2504  elseif (.not.all_empty .and. n .ne. o3n) then
2505  if (rlevs(vlev) .gt. lev_no_tr_fill) &
2506  call error_handler("TRACER "//trim(tracers(n))//" HAS MISSING DATA AT "//trim(slevs(vlev))//&
2507  ". SET MISSING VARIABLE CONDITION TO 'INTRP' TO AVOID THIS ERROR.", 1)
2508  endif
2509  ! If entire array is empty and method is set to intrp, switch method to fill
2510  if (trim(method) .eq. 'intrp' .and. all_empty) method='set_to_fill'
2511 
2512  call handle_grib_error(vname, slevs(vlev),method,value,varnum,read_from_input,iret,var=dummy2d)
2513  if (iret==1) then ! missing_var_method == skip or no entry
2514  if ( (tracers_input_oct10(n) == 1 .and. tracers_input_oct11(n) == 0) .or. & ! spec humidity
2515  (tracers_input_oct10(n) == 1 .and. tracers_input_oct11(n) == 1) .or. & ! rel humidity
2516  (tracers_input_oct10(n) == 14 .and. tracers_input_oct11(n) == 192) ) then ! ozone
2517  call error_handler("READING IN "//trim(tracers(n))//" AT LEVEL "//trim(slevs(vlev))&
2518  //". SET A FILL VALUE IN THE VARMAP TABLE IF THIS ERROR IS NOT DESIRABLE.",iret)
2519  endif
2520  endif
2521  endif ! method intrp
2522  endif !iret<=0
2523 
2524  if (n==1 .and. .not. hasspfh) then
2525  if (trim(external_model) .eq. 'GFS') then
2526  print *,'- CALL CALRH GFS'
2527  call rh2spfh_gfs(dummy2d,rlevs(vlev),dummy3d(:,:,vlev))
2528  else
2529  print *,'- CALL CALRH non-GFS'
2530  call rh2spfh(dummy2d,rlevs(vlev),dummy3d(:,:,vlev))
2531  end if
2532  endif
2533 
2534  dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8)
2535 
2536  enddo !vlev
2537 
2538 ! Jili Dong interpolation for missing levels
2539  if (is_missing .and. trim(method) .eq. 'intrp') then
2540  print *,'- INTERPOLATE TRACER '//trim(tracers(n))
2541  done_print = 0
2542  do jj = 1, j_input
2543  do ii = 1, i_input
2544  dummy3d_col_in=dummy3d(ii,jj,:)
2545  call dint2p(rlevs,dummy3d_col_in,lev_input,rlevs,dummy3d_col_out, &
2546  lev_input, 2, intrp_missing, intrp_ier)
2547  if (intrp_ier .gt. 0) call error_handler("Interpolation failed.",intrp_ier)
2548  dummy3d(ii,jj,:)=dummy3d_col_out
2549  enddo
2550  enddo
2551  do vlev=1,lev_input
2552  dummy2d = real(dummy3d(:,:,n) , kind=esmf_kind_r4)
2553  if (any(dummy2d .eq. intrp_missing)) then
2554  ! If we're outside the appropriate region, don't fill but error instead
2555  if (n == o3n .and. rlevs(vlev) .lt. lev_no_o3_fill) then
2556  call error_handler("TRACER "//trim(tracers(n))//" HAS MISSING DATA AT "//trim(slevs(vlev)),1)
2557  elseif (n .ne. o3n .and. rlevs(vlev) .gt. lev_no_tr_fill) then
2558  call error_handler("TRACER "//trim(tracers(n))//" HAS MISSING DATA AT "//trim(slevs(vlev)),1)
2559  else ! we're okay to fill missing data with provided fill value
2560  if (done_print .eq. 0) then
2561  print*, "Pressure out of range of existing data. Defaulting to fill value."
2562  done_print = 1
2563  end if !done print
2564  where(dummy2d .eq. intrp_missing) dummy2d = value
2565  dummy3d(:,:,vlev) = dummy2d
2566  end if !n & lev
2567  endif ! intrp_missing
2568  ! zero out negative tracers from interpolation/extrapolation
2569  where(dummy3d(:,:,vlev) .lt. 0.0) dummy3d(:,:,vlev) = 0.0
2570 ! print*,'tracer af intrp',vlev, maxval(dummy3d(:,:,vlev)),minval(dummy3d(:,:,vlev))
2571  end do !nlevs do
2572  end if !if intrp
2573  endif !localpet == 0
2574 
2575  if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT ", trim(tracers_input_vmap(n))
2576  call esmf_fieldscatter(tracers_input_grid(n), dummy3d, rootpet=0, rc=rc)
2577  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2578  call error_handler("IN FieldScatter", rc)
2579 
2580  enddo
2581 
2582  deallocate(dummy3d_col_in, dummy3d_col_out)
2583 
2584  call read_winds(u_tmp_3d,v_tmp_3d,localpet,octet_23,rlevs,lugb,pdt_num)
2585 
2586  if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT U-WIND."
2587  call esmf_fieldscatter(u_input_grid, u_tmp_3d, rootpet=0, rc=rc)
2588  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2589  call error_handler("IN FieldScatter", rc)
2590 
2591  if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT V-WIND."
2592  call esmf_fieldscatter(v_input_grid, v_tmp_3d, rootpet=0, rc=rc)
2593  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2594  call error_handler("IN FieldScatter", rc)
2595 
2596  if (localpet == 0) then
2597 
2598  print*,"- READ SURFACE PRESSURE."
2599  jdisc = 0 ! search for discipline - meteorological products
2600  j = 0 ! search at beginning of file.
2601  jpdt = -9999 ! array of values in product definition template, set to wildcard
2602  jids = -9999 ! array of values in identification section, set to wildcard
2603  jgdt = -9999 ! array of values in grid definition template, set to wildcard
2604  jgdtn = -1 ! search for any grid definition number.
2605  jpdtn = pdt_num ! Search for the product definition template number.
2606  jpdt(1) = 3 ! Sect4/oct 10 - param category - mass
2607  jpdt(2) = 0 ! Sect4/oct 11 - param number - pressure
2608  jpdt(10) = 1 ! Sect4/oct 23 - type of level - ground surface
2609  unpack=.true.
2610 
2611  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2612  unpack, k, gfld, iret)
2613  if (iret /= 0) call error_handler("READING SURFACE PRESSURE RECORD.", iret)
2614 
2615  dummy2d_8 = reshape(gfld%fld, (/i_input,j_input/) )
2616 
2617  endif ! Read surface pressure
2618 
2619  if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID SURFACE PRESSURE."
2620  call esmf_fieldscatter(ps_input_grid, dummy2d_8, rootpet=0, rc=rc)
2621  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2622  call error_handler("IN FieldScatter", rc)
2623 
2624 ! Read dzdt.
2625 
2626  if (localpet == 0) then
2627 
2628  print*,"- READ DZDT."
2629  vname = "dzdt"
2630  call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, &
2631  loc=varnum)
2632 
2633  jdisc = 0 ! search for discipline - meteorological products
2634  j = 0 ! search at beginning of file.
2635  jpdt = -9999 ! array of values in product definition template, set to wildcard
2636  jids = -9999 ! array of values in identification section, set to wildcard
2637  jgdt = -9999 ! array of values in grid definition template, set to wildcard
2638  jgdtn = -1 ! search for any grid definition number.
2639  jpdtn = pdt_num ! Search for the product definition template number.
2640  jpdt(1) = 2 ! Sect4/oct 10 - param category - momentum
2641  jpdt(2) = 9 ! Sect4/oct 11 - param number - dzdt
2642  jpdt(10) = octet_23 ! Sect4/oct 23 - type of level
2643 
2644  unpack=.true.
2645 
2646  do vlev = 1, lev_input
2647 
2648  jpdt(12) = nint(rlevs(vlev))
2649 
2650  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2651  unpack, k, gfld, iret)
2652 
2653  if (iret /= 0) then ! dzdt not found, look for omega.
2654  print*,"DZDT not available at level ", trim(slevs(vlev)), " so checking for VVEL"
2655  jpdt(2) = 8 ! Sect4/oct 11 - parameter number - omega
2656  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2657  unpack, k, gfld, iret)
2658  if (iret /= 0) then
2659  call handle_grib_error(vname, slevs(vlev),method,value,varnum,read_from_input,iret,var8=dum2d_1)
2660  if (iret==1) then ! missing_var_method == skip
2661  cycle
2662  endif
2663  else
2664  conv_omega = .true.
2665  dum2d_1 = reshape(gfld%fld, (/i_input,j_input/) )
2666  endif
2667  else ! found dzdt
2668  dum2d_1 = reshape(gfld%fld, (/i_input,j_input/) )
2669  endif
2670 
2671  dummy3d(:,:,vlev) = dum2d_1
2672 
2673  enddo
2674 
2675  endif ! Read of dzdt
2676 
2677  call mpi_bcast(conv_omega,1,mpi_logical,0,mpi_comm_world,rc)
2678 
2679  if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT DZDT."
2680  call esmf_fieldscatter(dzdt_input_grid, dummy3d, rootpet=0, rc=rc)
2681  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2682  call error_handler("IN FieldScatter", rc)
2683 
2684 ! Read terrain
2685 
2686  if (localpet == 0) then
2687 
2688  print*,"- READ TERRAIN."
2689  jdisc = 0 ! search for discipline - meteorological products
2690  j = 0 ! search at beginning of file.
2691  jpdt = -9999 ! array of values in product definition template, set to wildcard
2692  jids = -9999 ! array of values in identification section, set to wildcard
2693  jgdt = -9999 ! array of values in grid definition template, set to wildcard
2694  jgdtn = -1 ! search for any grid definition number.
2695  jpdtn = pdt_num ! Search for the product definition template number.
2696  jpdt(1) = 3 ! Sect4/oct 10 - param category - mass
2697  jpdt(2) = 5 ! Sect4/oct 11 - param number - geopotential height
2698  jpdt(10) = 1 ! Sect4/oct 23 - type of level - ground surface
2699  unpack=.true.
2700 
2701  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2702  unpack, k, gfld, iret)
2703  if (iret /= 0) call error_handler("READING TERRAIN HEIGHT RECORD.", iret)
2704 
2705  dummy2d_8 = reshape(gfld%fld, (/i_input,j_input/) )
2706 
2707  endif ! read of terrain.
2708 
2709  if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID TERRAIN."
2710  call esmf_fieldscatter(terrain_input_grid, dummy2d_8, rootpet=0, rc=rc)
2711  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2712  call error_handler("IN FieldScatter", rc)
2713 
2714  deallocate(dummy2d, dummy2d_8)
2715 
2716 if (.not. isnative) then
2717 
2718  !---------------------------------------------------------------------------
2719  ! Flip 'z' indices to all 3-d variables. Data is read in from model
2720  ! top to surface. This program expects surface to model top.
2721  !---------------------------------------------------------------------------
2722 
2723  if (localpet == 0) print*,"- CALL FieldGet FOR SURFACE PRESSURE."
2724  nullify(psptr)
2725  call esmf_fieldget(ps_input_grid, &
2726  farrayptr=psptr, rc=rc)
2727  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2728  call error_handler("IN FieldGet", rc)
2729 
2730  nullify(presptr)
2731  if (localpet == 0) print*,"- CALL FieldGet FOR 3-D PRESSURE."
2732  call esmf_fieldget(pres_input_grid, &
2733  computationallbound=clb, &
2734  computationalubound=cub, &
2735  farrayptr=presptr, rc=rc)
2736  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2737  call error_handler("IN FieldGet", rc)
2738 
2739  nullify(tptr)
2740  if (localpet == 0) print*,"- CALL FieldGet TEMPERATURE."
2741  call esmf_fieldget(temp_input_grid, &
2742  farrayptr=tptr, rc=rc)
2743  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2744  call error_handler("IN FieldGet", rc)
2745 
2746  nullify(uptr)
2747  if (localpet == 0) print*,"- CALL FieldGet FOR U"
2748  call esmf_fieldget(u_input_grid, &
2749  farrayptr=uptr, rc=rc)
2750  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2751  call error_handler("IN FieldGet", rc)
2752 
2753  nullify(vptr)
2754  if (localpet == 0) print*,"- CALL FieldGet FOR V"
2755  call esmf_fieldget(v_input_grid, &
2756  farrayptr=vptr, rc=rc)
2757  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2758  call error_handler("IN FieldGet", rc)
2759 
2760  nullify(wptr)
2761  if (localpet == 0) print*,"- CALL FieldGet FOR W"
2762  call esmf_fieldget(dzdt_input_grid, &
2763  farrayptr=wptr, rc=rc)
2764  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2765  call error_handler("IN FieldGet", rc)
2766 
2767  if (localpet == 0) print*,"- CALL FieldGet FOR TRACERS."
2768  do n=1,num_tracers_input
2769  nullify(qptr)
2770  call esmf_fieldget(tracers_input_grid(n), &
2771  farrayptr=qptr, rc=rc)
2772  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2773  call error_handler("IN FieldGet", rc)
2774  do i = clb(1),cub(1)
2775  do j = clb(2),cub(2)
2776  qptr(i,j,:) = qptr(i,j,lev_input:1:-1)
2777  end do
2778  end do
2779  end do
2780 
2781  do i = clb(1),cub(1)
2782  do j = clb(2),cub(2)
2783  presptr(i,j,:) = rlevs(lev_input:1:-1)
2784  tptr(i,j,:) = tptr(i,j,lev_input:1:-1)
2785  uptr(i,j,:) = uptr(i,j,lev_input:1:-1)
2786  vptr(i,j,:) = vptr(i,j,lev_input:1:-1)
2787  wptr(i,j,:) = wptr(i,j,lev_input:1:-1)
2788  end do
2789  end do
2790 
2791  if (localpet == 0) then
2792  print*,'psfc is ',clb(1),clb(2),psptr(clb(1),clb(2))
2793  print*,'pres is ',cub(1),cub(2),presptr(cub(1),cub(2),:)
2794 
2795  print*,'pres check 1',localpet,maxval(presptr(clb(1):cub(1),clb(2):cub(2),1)), &
2796  minval(presptr(clb(1):cub(1),clb(2):cub(2),1))
2797  print*,'pres check lev',localpet,maxval(presptr(clb(1):cub(1),clb(2):cub(2), &
2798  lev_input)),minval(presptr(clb(1):cub(1),clb(2):cub(2),lev_input))
2799  endif
2800 
2801 else ! is native coordinate (hybrid).
2802 
2803 ! For native files, read in pressure field directly from file but don't flip levels
2804 
2805  if (localpet == 0) then
2806 
2807  print*,"- READ PRESSURE."
2808 
2809  jdisc = 0 ! search for discipline - meteorological products
2810  j = 0 ! search at beginning of file.
2811  jpdt = -9999 ! array of values in product definition template, set to wildcard
2812  jids = -9999 ! array of values in identification section, set to wildcard
2813  jgdt = -9999 ! array of values in grid definition template, set to wildcard
2814  jgdtn = -1 ! search for any grid definition number.
2815  jpdtn = pdt_num ! Search for the product definition template number.
2816  jpdt(1) = 3 ! Sect4/oct 10 - parameter category - mass
2817  jpdt(2) = 0 ! Sect4/oct 11 - parameter number - pressure
2818  jpdt(10) = octet_23 ! Sect4/oct 23 - type of level.
2819  unpack=.true.
2820 
2821  do vlev = 1, lev_input
2822 
2823  jpdt(12) = nint(rlevs(vlev))
2824  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2825  unpack, k, gfld, iret)
2826  if (iret /= 0) then
2827  call error_handler("READING IN PRESSURE AT LEVEL "//trim(slevs(vlev)),iret)
2828  endif
2829 
2830  dum2d_1 = reshape(gfld%fld, (/i_input,j_input/) )
2831 
2832  dummy3d(:,:,vlev) = dum2d_1
2833 
2834  enddo
2835 
2836  endif ! localpet == 0
2837 
2838  if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID PRESSURE."
2839  call esmf_fieldscatter(pres_input_grid, dummy3d, rootpet=0, rc=rc)
2840  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2841  call error_handler("IN FieldScatter", rc)
2842 
2843  endif
2844 
2845  deallocate(dummy3d, dum2d_1)
2846 
2847 !---------------------------------------------------------------------------
2848 ! Convert from 2-d to 3-d component winds.
2849 !---------------------------------------------------------------------------
2850 
2851  call convert_winds
2852 
2853 !---------------------------------------------------------------------------
2854 ! Convert dpdt to dzdt if needed
2855 !---------------------------------------------------------------------------
2856 
2857  if (conv_omega) then
2858 
2859  if (localpet == 0) print*,"- CONVERT FROM OMEGA TO DZDT."
2860 
2861  nullify(tptr)
2862  if (localpet == 0) print*,"- CALL FieldGet TEMPERATURE."
2863  call esmf_fieldget(temp_input_grid, &
2864  farrayptr=tptr, rc=rc)
2865  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2866  call error_handler("IN FieldGet", rc)
2867 
2868  nullify(qptr)
2869  if (localpet == 0) print*,"- CALL FieldGet SPECIFIC HUMIDITY."
2870  call esmf_fieldget(tracers_input_grid(1), &
2871  computationallbound=clb, &
2872  computationalubound=cub, &
2873  farrayptr=qptr, rc=rc)
2874  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2875  call error_handler("IN FieldGet", rc)
2876 
2877  nullify(wptr)
2878  if (localpet == 0) print*,"- CALL FieldGet DZDT."
2879  call esmf_fieldget(dzdt_input_grid, &
2880  computationallbound=clb, &
2881  computationalubound=cub, &
2882  farrayptr=wptr, rc=rc)
2883  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2884  call error_handler("IN FieldGet", rc)
2885 
2886  nullify(presptr)
2887  call esmf_fieldget(pres_input_grid, &
2888  farrayptr=presptr, rc=rc)
2889  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2890  call error_handler("IN FieldGet", rc)
2891 
2892  call convert_omega(wptr,presptr,tptr,qptr,clb,cub)
2893 
2894  endif
2895 
2896  if (localpet == 0) call baclose(lugb, rc)
2897 
2898  end subroutine read_input_atm_grib2_file
2899 
2911  subroutine read_winds(u,v,localpet,octet_23,rlevs,lugb,pdt_num)
2912 
2913  use grib_mod
2914  use program_setup, only : get_var_cond
2915 
2916  implicit none
2917 
2918  integer, intent(in) :: localpet, lugb
2919  integer, intent(in) :: pdt_num, octet_23
2920 
2921  real(esmf_kind_r8), intent(inout), allocatable :: u(:,:,:),v(:,:,:)
2922  real(esmf_kind_r8), intent(in), dimension(lev_input) :: rlevs
2923 
2924  real(esmf_kind_r4), dimension(i_input,j_input) :: alpha
2925  real(esmf_kind_r8), dimension(i_input,j_input) :: lon, lat
2926  real(esmf_kind_r4), allocatable :: u_tmp(:,:),v_tmp(:,:)
2927  real(esmf_kind_r8), allocatable :: dum2d(:,:)
2928  real(esmf_kind_r4), dimension(i_input,j_input) :: ws,wd
2929  real(esmf_kind_r4) :: value_u, value_v,lov,latin1,latin2
2930  real(esmf_kind_r8) :: d2r
2931 
2932  integer :: varnum_u, varnum_v, vlev, &
2933  error, iret
2934  integer :: j, k, lugi, jgdtn, jpdtn
2935  integer :: jdisc, jids(200), jgdt(200), jpdt(200)
2936 
2937  character(len=20) :: vname
2938  character(len=50) :: method_u, method_v
2939 
2940  logical :: unpack
2941 
2942  type(gribfield) :: gfld
2943 
2944  d2r=acos(-1.0_esmf_kind_r8) / 180.0_esmf_kind_r8
2945  if (localpet==0) then
2946  allocate(u(i_input,j_input,lev_input))
2947  allocate(v(i_input,j_input,lev_input))
2948  else
2949  allocate(u(0,0,0))
2950  allocate(v(0,0,0))
2951  endif
2952 
2953  vname = "u"
2954  call get_var_cond(vname,this_miss_var_method=method_u, this_miss_var_value=value_u, &
2955  loc=varnum_u)
2956  vname = "v"
2957  call get_var_cond(vname,this_miss_var_method=method_v, this_miss_var_value=value_v, &
2958  loc=varnum_v)
2959 
2960  print*,"- CALL FieldGather FOR INPUT GRID LONGITUDE"
2961  call esmf_fieldgather(longitude_input_grid, lon, rootpet=0, tile=1, rc=error)
2962  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2963  call error_handler("IN FieldGather", error)
2964 
2965  print*,"- CALL FieldGather FOR INPUT GRID LATITUDE"
2966  call esmf_fieldgather(latitude_input_grid, lat, rootpet=0, tile=1, rc=error)
2967  if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2968  call error_handler("IN FieldGather", error)
2969 
2970  if (localpet==0) then
2971 
2972  lugi = 0 ! index file unit number
2973  jdisc = 0 ! search for discipline - meteorological products
2974  j = 0 ! search at beginning of file.
2975  jpdt = -9999 ! array of values in product definition template, set to wildcard
2976  jids = -9999 ! array of values in identification section, set to wildcard
2977  jgdt = -9999 ! array of values in grid definition template, set to wildcard
2978  jgdtn = -1 ! search for any grid definition number.
2979  jpdtn = pdt_num ! Search for the product definition template number.
2980  unpack=.false.
2981 
2982  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
2983  unpack, k, gfld, iret)
2984 
2985  if (iret /= 0) call error_handler("ERROR READING GRIB2 FILE.", iret)
2986 
2987  if (gfld%igdtnum == 32769) then ! grid definition template number - rotated lat/lon grid
2988 
2989  latin1 = real(float(gfld%igdtmpl(15))/1.0E6, kind=esmf_kind_r4)
2990  lov = real(float(gfld%igdtmpl(16))/1.0E6, kind=esmf_kind_r4)
2991 
2992  print*, "- CALL CALCALPHA_ROTLATLON with center lat,lon = ",latin1,lov
2993  call calcalpha_rotlatlon(lat,lon,latin1,lov,alpha)
2994 
2995  elseif (gfld%igdtnum == 30) then ! grid definition template number - lambert conformal grid.
2996 
2997  lov = real(float(gfld%igdtmpl(14))/1.0E6, kind=esmf_kind_r4)
2998  latin1 = real(float(gfld%igdtmpl(19))/1.0E6, kind=esmf_kind_r4)
2999  latin2 = real(float(gfld%igdtmpl(20))/1.0E6, kind=esmf_kind_r4)
3000 
3001  print*, "- CALL GRIDROT for LC grid with lov,latin1/2 = ",lov,latin1,latin2
3002  call gridrot(lov,latin1,latin2,lon,alpha)
3003 
3004  endif
3005 
3006  jpdt(10) = octet_23 ! Sec4/oct 23 - type of level.
3007 
3008  unpack=.true.
3009 
3010  allocate(dum2d(i_input,j_input))
3011  allocate(u_tmp(i_input,j_input))
3012  allocate(v_tmp(i_input,j_input))
3013 
3014  do vlev = 1, lev_input
3015 
3016  vname = ":UGRD:"
3017 
3018  jpdt(1) = 2 ! Sec4/oct 10 - parameter category - momentum
3019  jpdt(2) = 2 ! Sec4/oct 11 - parameter number - u-wind
3020  jpdt(12) = nint(rlevs(vlev)) ! Sect4/octs 25-28 - scaled value of fixed surface.
3021 
3022  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
3023  unpack, k, gfld, iret)
3024 
3025  if (iret /= 0) then
3026  call handle_grib_error(vname, slevs(vlev),method_u,value_u,varnum_u,read_from_input,iret,var=u_tmp)
3027  if (iret==1) then ! missing_var_method == skip
3028  call error_handler("READING IN U AT LEVEL "//trim(slevs(vlev))//". SET A FILL "// &
3029  "VALUE IN THE VARMAP TABLE IF THIS ERROR IS NOT DESIRABLE.",iret)
3030  endif
3031  else
3032  dum2d = reshape(gfld%fld, (/i_input,j_input/) )
3033  u_tmp(:,:) = real(dum2d, kind=esmf_kind_r4)
3034  endif
3035 
3036  vname = ":VGRD:"
3037 
3038  jpdt(2) = 3 ! Sec4/oct 11 - parameter number - v-wind
3039 
3040  call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
3041  unpack, k, gfld, iret)
3042 
3043  if (iret /= 0) then
3044  call handle_grib_error(vname, slevs(vlev),method_v,value_v,varnum_v,read_from_input,iret,var=v_tmp)
3045  if (iret==1) then ! missing_var_method == skip
3046  call error_handler("READING IN V AT LEVEL "//trim(slevs(vlev))//". SET A FILL "// &
3047  "VALUE IN THE VARMAP TABLE IF THIS ERROR IS NOT DESIRABLE.",iret)
3048  endif
3049  else
3050  dum2d = reshape(gfld%fld, (/i_input,j_input/) )
3051  v_tmp(:,:) = real(dum2d, kind=esmf_kind_r4)
3052  endif
3053 
3054  deallocate(dum2d)
3055 
3056  if (gfld%igdtnum == 0) then ! grid definition template number - lat/lon grid
3057  if (external_model == 'UKMET') then
3058  u(:,:,vlev) = u_tmp
3059  v(:,:,vlev) = (v_tmp(:,2:jp1_input) + v_tmp(:,1:j_input))/2
3060  else
3061  u(:,:,vlev) = u_tmp
3062  v(:,:,vlev) = v_tmp
3063  endif
3064  else if (gfld%igdtnum == 32769) then ! grid definition template number - rotated lat/lon grid
3065  ws = sqrt(u_tmp**2 + v_tmp**2)
3066  wd = real((atan2(-u_tmp,-v_tmp) / d2r), kind=esmf_kind_r4) ! calculate grid-relative wind direction
3067  wd = real((wd + alpha + 180.0), kind=esmf_kind_r4) ! Rotate from grid- to earth-relative direction
3068  wd = real((270.0 - wd), kind=esmf_kind_r4) ! Convert from meteorological (true N) to mathematical direction
3069  u(:,:,vlev) = -ws*cos(wd*d2r)
3070  v(:,:,vlev) = -ws*sin(wd*d2r)
3071  else
3072  u(:,:,vlev) = real(u_tmp * cos(alpha) + v_tmp * sin(alpha),esmf_kind_r8)
3073  v(:,:,vlev) = real(v_tmp * cos(alpha) - u_tmp * sin(alpha),esmf_kind_r8)
3074  endif
3075 
3076  print*, 'max, min U ', minval(u(:,:,vlev)), maxval(u(:,:,vlev))
3077  print*, 'max, min V ', minval(v(:,:,vlev)), maxval(v(:,:,vlev))
3078  enddo
3079  endif
3080 
3081 end subroutine read_winds
3082 
3086  subroutine convert_winds
3087 
3088  implicit none
3089 
3090  integer :: clb(4), cub(4)
3091  integer :: i, j, k, rc
3092 
3093  real(esmf_kind_r8) :: latrad, lonrad
3094  real(esmf_kind_r8), pointer :: windptr(:,:,:,:)
3095  real(esmf_kind_r8), pointer :: uptr(:,:,:)
3096  real(esmf_kind_r8), pointer :: vptr(:,:,:)
3097  real(esmf_kind_r8), pointer :: latptr(:,:)
3098  real(esmf_kind_r8), pointer :: lonptr(:,:)
3099 
3100  print*,"- CALL FieldGet FOR 3-D WIND."
3101  call esmf_fieldget(wind_input_grid, &
3102  computationallbound=clb, &
3103  computationalubound=cub, &
3104  farrayptr=windptr, rc=rc)
3105  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
3106  call error_handler("IN FieldGet", rc)
3107 
3108  print*,"- CALL FieldGet FOR U."
3109  call esmf_fieldget(u_input_grid, &
3110  farrayptr=uptr, rc=rc)
3111  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
3112  call error_handler("IN FieldGet", rc)
3113 
3114  print*,"- CALL FieldGet FOR V."
3115  call esmf_fieldget(v_input_grid, &
3116  farrayptr=vptr, rc=rc)
3117  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
3118  call error_handler("IN FieldGet", rc)
3119 
3120  print*,"- CALL FieldGet FOR LATITUDE."
3121  call esmf_fieldget(latitude_input_grid, &
3122  farrayptr=latptr, rc=rc)
3123  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
3124  call error_handler("IN FieldGet", rc)
3125 
3126  print*,"- CALL FieldGet FOR LONGITUDE."
3127  call esmf_fieldget(longitude_input_grid, &
3128  farrayptr=lonptr, rc=rc)
3129  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
3130  call error_handler("IN FieldGet", rc)
3131 
3132  do i = clb(1), cub(1)
3133  do j = clb(2), cub(2)
3134  latrad = latptr(i,j) * acos(-1.) / 180.0
3135  lonrad = lonptr(i,j) * acos(-1.) / 180.0
3136  do k = clb(3), cub(3)
3137  windptr(i,j,k,1) = uptr(i,j,k) * cos(lonrad) - vptr(i,j,k) * sin(latrad) * sin(lonrad)
3138  windptr(i,j,k,2) = uptr(i,j,k) * sin(lonrad) + vptr(i,j,k) * sin(latrad) * cos(lonrad)
3139  windptr(i,j,k,3) = vptr(i,j,k) * cos(latrad)
3140  enddo
3141  enddo
3142  enddo
3143 
3144  call esmf_fielddestroy(u_input_grid, rc=rc)
3145  call esmf_fielddestroy(v_input_grid, rc=rc)
3146 
3147  end subroutine convert_winds
3148 
3162 subroutine gridrot(lov,latin1,latin2,lon,rot)
3163 
3164  use model_grid, only : i_input,j_input
3165  implicit none
3166 
3167 
3168  real(esmf_kind_r4), intent(in) :: lov,latin1,latin2
3169  real(esmf_kind_r4), intent(inout) :: rot(i_input,j_input)
3170  real(esmf_kind_r8), intent(in) :: lon(i_input,j_input)
3171 
3172  real(esmf_kind_r4) :: trot(i_input,j_input), tlon(i_input,j_input)
3173  real(esmf_kind_r4) :: dtor = 3.14159265359_esmf_kind_r4/180.0_esmf_kind_r4
3174  real(esmf_kind_r4) :: an
3175  !trot_tmp = real(lon,esmf_kind_r4)-lov
3176  !trot = trot_tmp
3177  !where(trot_tmp > 180.0) trot = trot-360.0_esmf_kind_r4
3178  !where(trot_tmp < -180.0) trot = trot-360.0_esmf_kind_r4
3179 
3180  if ( (latin1 - latin2) .lt. 0.000001 ) then
3181  an = sin(latin1*dtor)
3182  else
3183  an = real(log( cos(latin1*dtor) / cos(latin2*dtor) ) / & log( tan(dtor*(90.0-latin1)/2.) / tan(dtor*(90.0-latin2)/2.)), kind=esmf_kind_r4)
3184  end if
3185 
3186  tlon = real((mod(lon - lov + 180. + 3600., 360.) - 180.), kind=esmf_kind_r4)
3187  trot = an * tlon
3188 
3189  rot = trot * dtor
3190 
3191 end subroutine gridrot
3192 
3202 subroutine calcalpha_rotlatlon(latgrid,longrid,cenlat,cenlon,alpha)
3203 
3204  use model_grid, only : i_input,j_input
3205  implicit none
3206 
3207  real(esmf_kind_r8), intent(in) :: latgrid(i_input,j_input), &
3208  longrid(i_input,j_input)
3209  real(esmf_kind_r4), intent(in) :: cenlat, cenlon
3210  real(esmf_kind_r4), intent(out) :: alpha(i_input,j_input)
3211 
3212  ! Variables local to subroutine
3213  real(esmf_kind_r8) :: D2R,lon0_r,lat0_r,sphi0,cphi0
3214  real(esmf_kind_r8), DIMENSION(i_input,j_input) :: tlat,tlon,tph,sinalpha
3215 
3216  d2r = acos(-1.0_esmf_kind_r8) / 180.0_esmf_kind_r8
3217  if (cenlon .lt. 0) then
3218  lon0_r = (cenlon + 360.0)*d2r
3219  else
3220  lon0_r = cenlon*d2r
3221  end if
3222  lat0_r=cenlat*d2r
3223  sphi0=sin(lat0_r)
3224  cphi0=cos(lat0_r)
3225 
3226  ! deal with input lat/lon
3227  tlat = latgrid * d2r
3228  tlon = longrid * d2r
3229 
3230  ! Calculate alpha (rotation angle)
3231  tlon = -tlon + lon0_r
3232  tph = asin(cphi0*sin(tlat) - sphi0*cos(tlat)*cos(tlon))
3233  sinalpha = sphi0 * sin(tlon) / cos(tph)
3234  alpha = real((-asin(sinalpha)/D2R), kind=esmf_kind_r4)
3235  ! returns alpha in degrees
3236 end subroutine calcalpha_rotlatlon
3237 
3241 subroutine cleanup_input_atm_data
3242 
3243  implicit none
3244 
3245  integer :: rc, n
3246 
3247  print*,'- DESTROY ATMOSPHERIC INPUT DATA.'
3248 
3249  call esmf_fielddestroy(terrain_input_grid, rc=rc)
3250  call esmf_fielddestroy(pres_input_grid, rc=rc)
3251  call esmf_fielddestroy(dzdt_input_grid, rc=rc)
3252  call esmf_fielddestroy(temp_input_grid, rc=rc)
3253  call esmf_fielddestroy(wind_input_grid, rc=rc)
3254  call esmf_fielddestroy(ps_input_grid, rc=rc)
3255 
3256  do n = 1, num_tracers_input
3257  call esmf_fielddestroy(tracers_input_grid(n), rc=rc)
3258  enddo
3259  deallocate(tracers_input_grid)
3260 
3261  end subroutine cleanup_input_atm_data
3262 
3263 end module atm_input_data
3264 
grib2_util
Utilities for use when reading grib2 data.
Definition: grib2_util.F90:12
program_setup
This module contains code to read the setup namelist file, handle the varmap file for GRIB2 data...
Definition: program_setup.F90:9
program_setup::atm_tracer_files_input_grid
character(len=500), dimension(6), public atm_tracer_files_input_grid
File names of input atmospheric restart tracer files.
Definition: program_setup.F90:24
program_setup::num_tracers_input
integer, public num_tracers_input
Number of atmospheric tracers in input file.
Definition: program_setup.F90:61
atm_input_data::slevs
character(len=50), dimension(:), allocatable, private slevs
The atmospheric levels in the GRIB2 input file.
Definition: atm_input_data.F90:63
atm_input_data::read_input_atm_gfs_sigio_file
subroutine read_input_atm_gfs_sigio_file(localpet)
Read input atmospheric data from spectral gfs (old sigio format).
Definition: atm_input_data.F90:242
model_grid::ip1_input
integer, public ip1_input
i_input plus 1
Definition: model_grid.F90:33
atm_input_data
Read atmospheric data on the input grid.
Definition: atm_input_data.F90:15
model_grid::j_input
integer, public j_input
j-dimension of input grid (or of each global tile)
Definition: model_grid.F90:30
atm_input_data::dpres_input_grid
type(esmf_field) dpres_input_grid
pressure thickness
Definition: atm_input_data.F90:49
atm_input_data::read_input_atm_restart_file
subroutine read_input_atm_restart_file(localpet)
Read input grid fv3 atmospheric data &#39;warm&#39; restart files.
Definition: atm_input_data.F90:1017
atm_input_data::dzdt_input_grid
type(esmf_field), public dzdt_input_grid
vert velocity
Definition: atm_input_data.F90:48
program_setup::tracers
character(len=20), dimension(max_tracers), public tracers
Name of each atmos tracer to be processed.
Definition: program_setup.F90:66
atm_input_data::read_input_atm_data
subroutine, public read_input_atm_data(localpet)
Read input grid atmospheric data driver.
Definition: atm_input_data.F90:76
model_grid::jp1_input
integer, public jp1_input
j_input plus 1
Definition: model_grid.F90:35
atm_input_data::read_input_atm_tiled_history_file
subroutine read_input_atm_tiled_history_file(localpet)
Read input grid fv3 atmospheric tiled history files in netcdf format.
Definition: atm_input_data.F90:1660
atm_input_data::tracers_input_grid
type(esmf_field), dimension(:), allocatable, public tracers_input_grid
tracers
Definition: atm_input_data.F90:58
model_grid
Sets up the ESMF grid objects for the input data grid and target FV3 grid.
Definition: model_grid.F90:9
program_setup::tracers_input
character(len=20), dimension(max_tracers), public tracers_input
Name of each atmos tracer record in the input file.
Definition: program_setup.F90:72
model_grid::latitude_input_grid
type(esmf_field), public latitude_input_grid
latitude of grid center, input grid
Definition: model_grid.F90:57
atm_input_data::pres_input_grid
type(esmf_field), public pres_input_grid
3-d pressure
Definition: atm_input_data.F90:50
program_setup::get_var_cond
subroutine, public get_var_cond(var_name, this_miss_var_method, this_miss_var_value, this_field_var_name, loc)
Search the variable mapping table to find conditions for handling missing variables.
Definition: program_setup.F90:445
atm_input_data::levp1_input
integer, public levp1_input
number of atmos layer interfaces
Definition: atm_input_data.F90:61
model_grid::num_tiles_input_grid
integer, public num_tiles_input_grid
Number of tiles, input grid.
Definition: model_grid.F90:47
grib2_util::rh2spfh_gfs
subroutine, public rh2spfh_gfs(rh_sphum, p, t)
Convert relative humidity to specific humidity (GFS formula) Calculation of saturation water vapor pr...
Definition: grib2_util.F90:81
atm_input_data::v_input_grid
type(esmf_field), public v_input_grid
box center
Definition: atm_input_data.F90:56
program_setup::external_model
character(len=20), public external_model
The model that the input data is derived from.
Definition: program_setup.F90:57
atm_input_data::init_atm_esmf_fields
subroutine init_atm_esmf_fields
Create atmospheric esmf fields.
Definition: atm_input_data.F90:146
atm_input_data::read_winds
subroutine read_winds(u, v, localpet, octet_23, rlevs, lugb, pdt_num)
Read winds from a grib2 file.
Definition: atm_input_data.F90:2912
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
atm_input_data::wind_input_grid
type(esmf_field), public wind_input_grid
3-component wind
Definition: atm_input_data.F90:57
atm_input_data::convert_winds
subroutine, public convert_winds
Convert winds from 2-d to 3-d components.
Definition: atm_input_data.F90:3087
atm_input_data::read_input_atm_gaussian_nemsio_file
subroutine read_input_atm_gaussian_nemsio_file(localpet)
Read input grid atmospheric fv3 gaussian nemsio files.
Definition: atm_input_data.F90:734
atm_input_data::terrain_input_grid
type(esmf_field), public terrain_input_grid
terrain height
Definition: atm_input_data.F90:52
atm_input_data::ps_input_grid
type(esmf_field), public ps_input_grid
surface pressure
Definition: atm_input_data.F90:51
atm_input_data::calcalpha_rotlatlon
subroutine calcalpha_rotlatlon(latgrid, longrid, cenlat, cenlon, alpha)
Calculate rotation angle for rotated latlon grids.
Definition: atm_input_data.F90:3204
atm_input_data::read_input_atm_gaussian_netcdf_file
subroutine read_input_atm_gaussian_netcdf_file(localpet)
Read fv3 netcdf gaussian history file.
Definition: atm_input_data.F90:1276
atm_input_data::gridrot
subroutine gridrot(lov, latin1, latin2, lon, rot)
Compute grid rotation angle for non-latlon grids.
Definition: atm_input_data.F90:3163
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
grib2_util::convert_omega
subroutine, public convert_omega(omega, p, t, q, clb, cub)
Convert omega to vertical velocity.
Definition: grib2_util.F90:218
atm_input_data::u_input_grid
type(esmf_field), public u_input_grid
u/v wind at grid
Definition: atm_input_data.F90:55
atm_input_data::cleanup_input_atm_data
subroutine, public cleanup_input_atm_data
Free up memory associated with atm data.
Definition: atm_input_data.F90:3243
program_setup::input_type
character(len=25), public input_type
Input data type:
Definition: program_setup.F90:43
program_setup::atm_core_files_input_grid
character(len=500), dimension(7), public atm_core_files_input_grid
File names of input atmospheric restart core files.
Definition: program_setup.F90:23
atm_input_data::temp_input_grid
type(esmf_field), public temp_input_grid
temperature
Definition: atm_input_data.F90:53
model_grid::longitude_input_grid
type(esmf_field), public longitude_input_grid
longitude of grid center, input grid
Definition: model_grid.F90:59
atm_input_data::read_input_atm_grib2_file
subroutine read_input_atm_grib2_file(localpet)
Read input grid atmospheric fv3gfs grib2 files.
Definition: atm_input_data.F90:1951
program_setup::grib2_file_input_grid
character(len=500), public grib2_file_input_grid
REQUIRED.
Definition: program_setup.F90:30
atm_input_data::read_input_atm_gfs_gaussian_nemsio_file
subroutine read_input_atm_gfs_gaussian_nemsio_file(localpet)
Read input atmospheric data from spectral gfs (global gaussian in nemsio format.
Definition: atm_input_data.F90:481
program_setup::read_from_input
logical, dimension(:), allocatable, public read_from_input
When false, variable was not read from GRIB2 input file.
Definition: program_setup.F90:63
atm_input_data::lev_input
integer, public lev_input
number of atmospheric layers
Definition: atm_input_data.F90:60
model_grid::i_input
integer, public i_input
i-dimension of input grid (or of each global tile)
Definition: model_grid.F90:27
grib2_util::rh2spfh
subroutine, public rh2spfh(rh_sphum, p, t)
Convert relative humidity to specific humidity.
Definition: grib2_util.F90:38