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