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