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