chgres_cube  1.12.0
 All Data Structures Files Functions Variables
atmosphere.F90
Go to the documentation of this file.
1 
4 
19  module atmosphere
20 
21  use esmf
22 
23  use atmosphere_target_data, only : lev_target, levp1_target, nvcoord_target, &
24  vcoord_target, delp_target_grid, &
25  dzdt_target_grid, ps_target_grid, &
26  temp_target_grid, tracers_target_grid, &
27  u_s_target_grid, v_s_target_grid, &
28  u_w_target_grid, v_w_target_grid, &
29  zh_target_grid, qnwfa_climo_target_grid, &
30  qnifa_climo_target_grid
31 
32  use atm_input_data, only : lev_input, &
33  levp1_input, &
34  tracers_input_grid, &
35  dzdt_input_grid, &
36  ps_input_grid, &
37  xwind_input_grid, &
38  ywind_input_grid, &
39  zwind_input_grid, &
40  temp_input_grid, &
41  pres_input_grid, &
42  terrain_input_grid, &
43  read_input_atm_data, &
44  cleanup_input_atm_data
45 
46  use model_grid, only : target_grid, &
47  latitude_s_target_grid, &
48  longitude_s_target_grid, &
49  latitude_w_target_grid, &
50  longitude_w_target_grid, &
51  terrain_target_grid
52 
53  use program_setup, only : vcoord_file_target_grid, &
54  wam_cold_start, wam_parm_file, &
55  cycle_year, cycle_mon, &
56  cycle_day, cycle_hour, &
57  regional, &
58  tracers, num_tracers, &
59  num_tracers_input, &
60  atm_weight_file, &
61  use_thomp_mp_climo
62 
63  use thompson_mp_climo_data, only : read_thomp_mp_climo_data, &
64  cleanup_thomp_mp_climo_input_data, &
65  qnifa_climo_input_grid, &
66  qnwfa_climo_input_grid, &
67  thomp_pres_climo_input_grid, &
68  lev_thomp_mp_climo
69 
70  use write_data, only : write_fv3_atm_header_netcdf, &
71  write_fv3_atm_bndy_data_netcdf, &
72  write_fv3_atm_data_netcdf
73 
74  use utilities, only : error_handler
75 
76  implicit none
77 
78  private
79 
80  type(esmf_field) :: dzdt_b4adj_target_grid
81  type(esmf_field), allocatable :: tracers_b4adj_target_grid(:)
82  type(esmf_field) :: ps_b4adj_target_grid
83  type(esmf_field) :: pres_target_grid
84  type(esmf_field) :: pres_b4adj_target_grid
85  type(esmf_field) :: temp_b4adj_target_grid
86  type(esmf_field) :: terrain_interp_to_target_grid
87  type(esmf_field) :: xwind_target_grid
88  type(esmf_field) :: ywind_target_grid
89  type(esmf_field) :: zwind_target_grid
90  type(esmf_field) :: xwind_b4adj_target_grid
91  type(esmf_field) :: ywind_b4adj_target_grid
92  type(esmf_field) :: zwind_b4adj_target_grid
93  type(esmf_field) :: xwind_s_target_grid
94  type(esmf_field) :: ywind_s_target_grid
95  type(esmf_field) :: zwind_s_target_grid
96  type(esmf_field) :: xwind_w_target_grid
97  type(esmf_field) :: ywind_w_target_grid
98  type(esmf_field) :: zwind_w_target_grid
99 
100 ! Fields associated with thompson microphysics climatological tracers.
101 
102  type(esmf_field) :: qnifa_climo_b4adj_target_grid
104  type(esmf_field) :: qnwfa_climo_b4adj_target_grid
106  type(esmf_field) :: thomp_pres_climo_b4adj_target_grid
108 
109  public :: atmosphere_driver
110  public :: read_vcoord_info
111 
112  contains
113 
118  subroutine atmosphere_driver(localpet)
119 
120  use mpi_f08
121 
122  implicit none
123 
124  integer, intent(in) :: localpet
125 
126  integer :: isrctermprocessing
127  integer :: rc, n
128 
129  type(esmf_regridmethod_flag) :: method
130  type(esmf_routehandle) :: regrid_bl
131 
132  real(esmf_kind_r8), parameter :: p0=101325.0
133  real(esmf_kind_r8), parameter :: rd = 287.058
134  real(esmf_kind_r8), parameter :: grav = 9.81
135  real(esmf_kind_r8), parameter :: lapse = -6.5e-03
136 
137  real(esmf_kind_r8), parameter :: exponent = rd*lapse/grav
138  real(esmf_kind_r8), parameter :: one_over_exponent = 1.0 / exponent
139 
140  real(esmf_kind_r8), pointer :: psptr(:,:), tempptr(:,:,:)
141 
142 !-----------------------------------------------------------------------------------
143 ! Read atmospheric fields on the input grid.
144 !-----------------------------------------------------------------------------------
145 
146  call read_input_atm_data(localpet)
147 
148 !-----------------------------------------------------------------------------------
149 ! Read vertical coordinate info for target grid.
150 !-----------------------------------------------------------------------------------
151 
152  call read_vcoord_info
153 
154 !-----------------------------------------------------------------------------------
155 ! Create target grid field objects to hold data before vertical adjustment.
156 !-----------------------------------------------------------------------------------
157 
158  call create_atm_b4adj_esmf_fields
159 
160 !-----------------------------------------------------------------------------------
161 ! Horizontally interpolate. If specified, use weights from file.
162 !-----------------------------------------------------------------------------------
163 
164  isrctermprocessing = 1
165 
166  if (trim(atm_weight_file) /= "NULL") then
167 
168  print*,"- CALL FieldSMMStore FOR ATMOSPHERIC FIELDS."
169 
170  call esmf_fieldsmmstore(temp_input_grid, &
171  temp_b4adj_target_grid, &
172  atm_weight_file, &
173  routehandle=regrid_bl, &
174  srctermprocessing=isrctermprocessing, rc=rc)
175  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
176  call error_handler("IN FieldSMMStore", rc)
177 
178  else
179 
180  print*,"- CALL FieldRegridStore FOR ATMOSPHERIC FIELDS."
181 
182  method=esmf_regridmethod_bilinear
183 
184  call esmf_fieldregridstore(temp_input_grid, &
185  temp_b4adj_target_grid, &
186  polemethod=esmf_polemethod_allavg, &
187  srctermprocessing=isrctermprocessing, &
188  routehandle=regrid_bl, &
189  regridmethod=method, rc=rc)
190  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
191  call error_handler("IN FieldRegridStore", rc)
192 
193  endif
194 
195  print*,"- CALL Field_Regrid FOR TEMPERATURE."
196  call esmf_fieldregrid(temp_input_grid, &
197  temp_b4adj_target_grid, &
198  routehandle=regrid_bl, &
199  termorderflag=esmf_termorder_srcseq, &
200  rc=rc)
201  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
202  call error_handler("IN FieldRegrid", rc)
203 
204  print*,"- CALL Field_Regrid FOR PRESSURE."
205  call esmf_fieldregrid(pres_input_grid, &
206  pres_b4adj_target_grid, &
207  routehandle=regrid_bl, &
208  termorderflag=esmf_termorder_srcseq, &
209  rc=rc)
210  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
211  call error_handler("IN FieldRegrid", rc)
212 
213  do n = 1, num_tracers_input
214  print*,"- CALL Field_Regrid FOR TRACER ", trim(tracers(n))
215  call esmf_fieldregrid(tracers_input_grid(n), &
216  tracers_b4adj_target_grid(n), &
217  routehandle=regrid_bl, &
218  termorderflag=esmf_termorder_srcseq, rc=rc)
219  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
220  call error_handler("IN FieldRegrid", rc)
221 
222  enddo
223 
224  print*,"- CALL Field_Regrid FOR VERTICAL VELOCITY."
225  call esmf_fieldregrid(dzdt_input_grid, &
226  dzdt_b4adj_target_grid, &
227  routehandle=regrid_bl, &
228  termorderflag=esmf_termorder_srcseq, rc=rc)
229  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
230  call error_handler("IN FieldRegrid", rc)
231 
232  nullify(tempptr)
233  print*,"- CALL FieldGet FOR INPUT GRID VERTICAL VEL."
234  call esmf_fieldget(dzdt_input_grid, &
235  farrayptr=tempptr, rc=rc)
236  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
237  call error_handler("IN FieldGet", rc)
238 
239  print*, "MIN MAX W INPUT = ", minval(tempptr), maxval(tempptr)
240 
241  nullify(tempptr)
242  print*,"- CALL FieldGet FOR VERTICAL VEL B4ADJ."
243  call esmf_fieldget(dzdt_b4adj_target_grid, &
244  farrayptr=tempptr, rc=rc)
245  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
246  call error_handler("IN FieldGet", rc)
247 
248  print*, "MIN MAX W B4ADJ = ", minval(tempptr), maxval(tempptr)
249 
250  nullify(psptr)
251  print*,"- CALL FieldGet FOR INPUT SURFACE PRESSURE."
252  call esmf_fieldget(ps_input_grid, &
253  farrayptr=psptr, rc=rc)
254  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
255  call error_handler("IN FieldGet", rc)
256 
257 !------------------------------------------------------------------------------------
258 ! Assume standard lapse rate when interpolating pressure (per Phil Pegion).
259 !------------------------------------------------------------------------------------
260 
261  psptr = (psptr/p0)**exponent
262 
263  print*,"- CALL Field_Regrid FOR SURFACE PRESSURE."
264  call esmf_fieldregrid(ps_input_grid, &
265  ps_b4adj_target_grid, &
266  routehandle=regrid_bl, &
267  termorderflag=esmf_termorder_srcseq, &
268  rc=rc)
269  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
270  call error_handler("IN FieldRegrid", rc)
271 
272  nullify(psptr)
273  print*,"- CALL FieldGet FOR INPUT SURFACE PRESSURE B4ADJ."
274  call esmf_fieldget(ps_b4adj_target_grid, &
275  farrayptr=psptr, rc=rc)
276  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
277  call error_handler("IN FieldGet", rc)
278 
279  psptr = p0 * psptr**one_over_exponent
280 
281  print*,"- CALL Field_Regrid FOR TERRAIN."
282  call esmf_fieldregrid(terrain_input_grid, &
283  terrain_interp_to_target_grid, &
284  routehandle=regrid_bl, &
285  termorderflag=esmf_termorder_srcseq, &
286  rc=rc)
287  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
288  call error_handler("IN FieldRegrid", rc)
289 
290  print*,"- CALL Field_Regrid FOR x WIND."
291  call esmf_fieldregrid(xwind_input_grid, &
292  xwind_b4adj_target_grid, &
293  routehandle=regrid_bl, &
294  termorderflag=esmf_termorder_srcseq, rc=rc)
295  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
296  call error_handler("IN FieldRegrid", rc)
297 
298  print*,"- CALL Field_Regrid FOR y WIND."
299  call esmf_fieldregrid(ywind_input_grid, &
300  ywind_b4adj_target_grid, &
301  routehandle=regrid_bl, &
302  termorderflag=esmf_termorder_srcseq, rc=rc)
303  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
304  call error_handler("IN FieldRegrid", rc)
305 
306  print*,"- CALL Field_Regrid FOR z WIND."
307  call esmf_fieldregrid(zwind_input_grid, &
308  zwind_b4adj_target_grid, &
309  routehandle=regrid_bl, &
310  termorderflag=esmf_termorder_srcseq, rc=rc)
311  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
312  call error_handler("IN FieldRegrid", rc)
313 
314 
315 
316  print*,"- CALL FieldRegridRelease."
317  call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
318  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
319  call error_handler("IN FieldRegridRelease", rc)
320 
321 !-----------------------------------------------------------------------------------
322 ! Deallocate input fields.
323 !-----------------------------------------------------------------------------------
324 
325  call cleanup_input_atm_data
326 
327 !-----------------------------------------------------------------------------------
328 ! Create target grid field objects to hold data after vertical interpolation.
329 !-----------------------------------------------------------------------------------
330 
331  call create_atm_esmf_fields
332 
333 !-----------------------------------------------------------------------------------
334 ! Adjust surface pressure for terrain differences.
335 !-----------------------------------------------------------------------------------
336 
337  call newps(localpet)
338 
339 !-----------------------------------------------------------------------------------
340 ! Compute 3-d pressure based on adjusted surface pressure.
341 !-----------------------------------------------------------------------------------
342 
343  call newpr1(localpet)
344 
345 !-----------------------------------------------------------------------------------
346 ! Vertically interpolate.
347 !-----------------------------------------------------------------------------------
348 
349  call vintg
350 
351  if( wam_cold_start ) then
352  call vintg_wam(cycle_year,cycle_mon,cycle_day,cycle_hour,wam_parm_file)
353  endif
354 
355 !-----------------------------------------------------------------------------------
356 ! Compute height.
357 !-----------------------------------------------------------------------------------
358 
359  call compute_zh
360 
361 !-----------------------------------------------------------------------------------
362 ! Free up memory.
363 !-----------------------------------------------------------------------------------
364 
365  call cleanup_target_atm_b4adj_data
366 
367 !-----------------------------------------------------------------------------------
368 ! Interpolate winds to 'd' grid.
369 !-----------------------------------------------------------------------------------
370 
371  isrctermprocessing = 1
372  method=esmf_regridmethod_bilinear
373 
374  print*,"- CALL FieldRegridStore FOR X-WIND WEST EDGE."
375  call esmf_fieldregridstore(xwind_target_grid, &
376  xwind_w_target_grid, &
377  polemethod=esmf_polemethod_allavg, &
378  srctermprocessing=isrctermprocessing, &
379  routehandle=regrid_bl, &
380  extrapmethod=esmf_extrapmethod_nearest_stod, &
381  regridmethod=method, rc=rc)
382  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
383  call error_handler("IN FieldRegridStore", rc)
384 
385  print*,"- CALL Field_Regrid FOR X-WIND WEST EDGE."
386  call esmf_fieldregrid(xwind_target_grid, &
387  xwind_w_target_grid, &
388  routehandle=regrid_bl, &
389  termorderflag=esmf_termorder_srcseq, rc=rc)
390  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
391  call error_handler("IN FieldRegrid", rc)
392 
393  print*,"- CALL Field_Regrid FOR Y-WIND WEST EDGE."
394  call esmf_fieldregrid(ywind_target_grid, &
395  ywind_w_target_grid, &
396  routehandle=regrid_bl, &
397  termorderflag=esmf_termorder_srcseq, rc=rc)
398  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
399  call error_handler("IN FieldRegrid", rc)
400 
401  print*,"- CALL Field_Regrid FOR Z-WIND WEST EDGE."
402  call esmf_fieldregrid(zwind_target_grid, &
403  zwind_w_target_grid, &
404  routehandle=regrid_bl, &
405  termorderflag=esmf_termorder_srcseq, rc=rc)
406  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
407  call error_handler("IN FieldRegrid", rc)
408 
409  print*,"- CALL FieldRegridRelease."
410  call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
411  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
412  call error_handler("IN FieldRegridRelease", rc)
413 
414  isrctermprocessing = 1
415  method=esmf_regridmethod_bilinear
416 
417  print*,"- CALL FieldRegridStore FOR X-WIND SOUTH EDGE."
418  call esmf_fieldregridstore(xwind_target_grid, &
419  xwind_s_target_grid, &
420  polemethod=esmf_polemethod_allavg, &
421  srctermprocessing=isrctermprocessing, &
422  routehandle=regrid_bl, &
423  extrapmethod=esmf_extrapmethod_nearest_stod, &
424  regridmethod=method, rc=rc)
425  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
426  call error_handler("IN FieldRegridStore", rc)
427 
428  print*,"- CALL Field_Regrid FOR X-WIND SOUTH EDGE."
429  call esmf_fieldregrid(xwind_target_grid, &
430  xwind_s_target_grid, &
431  routehandle=regrid_bl, &
432  termorderflag=esmf_termorder_srcseq, rc=rc)
433  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
434  call error_handler("IN FieldRegrid", rc)
435 
436  print*,"- CALL Field_Regrid FOR Y-WIND SOUTH EDGE."
437  call esmf_fieldregrid(ywind_target_grid, &
438  ywind_s_target_grid, &
439  routehandle=regrid_bl, &
440  termorderflag=esmf_termorder_srcseq, rc=rc)
441  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
442  call error_handler("IN FieldRegrid", rc)
443 
444  print*,"- CALL Field_Regrid FOR Z-WIND SOUTH EDGE."
445  call esmf_fieldregrid(zwind_target_grid, &
446  zwind_s_target_grid, &
447  routehandle=regrid_bl, &
448  termorderflag=esmf_termorder_srcseq, rc=rc)
449  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
450  call error_handler("IN FieldRegrid", rc)
451 
452  print*,"- CALL FieldRegridRelease."
453  call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
454  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
455  call error_handler("IN FieldRegridRelease", rc)
456 
457 !-----------------------------------------------------------------------------------
458 ! Convert from 3-d to 2-d cartesian winds.
459 !-----------------------------------------------------------------------------------
460 
461  call convert_winds_to_uv
462 
463 !-----------------------------------------------------------------------------------
464 ! If selected, process thompson microphysics climatological fields.
465 !-----------------------------------------------------------------------------------
466 
467  if (use_thomp_mp_climo) then
468  call read_thomp_mp_climo_data
469  call horiz_interp_thomp_mp_climo
470  call vintg_thomp_mp_climo
471  endif
472 
473 !-----------------------------------------------------------------------------------
474 ! Write target data to file.
475 !-----------------------------------------------------------------------------------
476 
477  call write_fv3_atm_header_netcdf(localpet)
478  if (regional <= 1) call write_fv3_atm_data_netcdf(localpet)
479  if (regional >= 1) call write_fv3_atm_bndy_data_netcdf(localpet)
480 
481 !-----------------------------------------------------------------------------------
482 ! Free up memory.
483 !-----------------------------------------------------------------------------------
484 
485  call cleanup_all_target_atm_data
486 
487  end subroutine atmosphere_driver
488 
494  subroutine create_atm_b4adj_esmf_fields
495 
496  implicit none
497 
498  integer :: rc, n
499 
500  allocate(tracers_b4adj_target_grid(num_tracers_input))
501 
502  do n = 1, num_tracers_input
503  print*,"- CALL FieldCreate FOR TARGET GRID TRACER BEFORE ADJUSTMENT ", trim(tracers(n))
504  tracers_b4adj_target_grid(n) = esmf_fieldcreate(target_grid, &
505  typekind=esmf_typekind_r8, &
506  staggerloc=esmf_staggerloc_center, &
507  ungriddedlbound=(/1/), &
508  ungriddedubound=(/lev_input/), rc=rc)
509  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
510  call error_handler("IN FieldCreate", rc)
511  enddo
512 
513  print*,"- CALL FieldCreate FOR TARGET GRID TEMPERATURE BEFORE ADJUSTMENT."
514  temp_b4adj_target_grid = esmf_fieldcreate(target_grid, &
515  typekind=esmf_typekind_r8, &
516  staggerloc=esmf_staggerloc_center, &
517  ungriddedlbound=(/1/), &
518  ungriddedubound=(/lev_input/), rc=rc)
519  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
520  call error_handler("IN FieldCreate", rc)
521 
522  print*,"- CALL FieldCreate FOR TARGET GRID PRESSURE BEFORE ADJUSTMENT."
523  pres_b4adj_target_grid = esmf_fieldcreate(target_grid, &
524  typekind=esmf_typekind_r8, &
525  staggerloc=esmf_staggerloc_center, &
526  ungriddedlbound=(/1/), &
527  ungriddedubound=(/lev_input/), rc=rc)
528  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
529  call error_handler("IN FieldCreate", rc)
530 
531  print*,"- CALL FieldCreate FOR TARGET GRID VERTICAL VELOCITY BEFORE ADJUSTMENT."
532  dzdt_b4adj_target_grid = esmf_fieldcreate(target_grid, &
533  typekind=esmf_typekind_r8, &
534  staggerloc=esmf_staggerloc_center, &
535  ungriddedlbound=(/1/), &
536  ungriddedubound=(/lev_input/), rc=rc)
537  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
538  call error_handler("IN FieldCreate", rc)
539 
540  print*,"- CALL FieldCreate FOR TARGET GRID xwind."
541  xwind_b4adj_target_grid = esmf_fieldcreate(target_grid, &
542  typekind=esmf_typekind_r8, &
543  staggerloc=esmf_staggerloc_center, &
544  ungriddedlbound=(/1/), &
545  ungriddedubound=(/lev_input/), rc=rc)
546  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
547  call error_handler("IN FieldCreate", rc)
548 
549  print*,"- CALL FieldCreate FOR TARGET GRID ywind."
550  ywind_b4adj_target_grid = esmf_fieldcreate(target_grid, &
551  typekind=esmf_typekind_r8, &
552  staggerloc=esmf_staggerloc_center, &
553  ungriddedlbound=(/1/), &
554  ungriddedubound=(/lev_input/), rc=rc)
555  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
556  call error_handler("IN FieldCreate", rc)
557 
558  print*,"- CALL FieldCreate FOR TARGET GRID zwind."
559  zwind_b4adj_target_grid = esmf_fieldcreate(target_grid, &
560  typekind=esmf_typekind_r8, &
561  staggerloc=esmf_staggerloc_center, &
562  ungriddedlbound=(/1/), &
563  ungriddedubound=(/lev_input/), rc=rc)
564  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
565  call error_handler("IN FieldCreate", rc)
566 
567  print*,"- CALL FieldCreate FOR TARGET TERRAIN."
568  terrain_interp_to_target_grid = esmf_fieldcreate(target_grid, &
569  typekind=esmf_typekind_r8, &
570  staggerloc=esmf_staggerloc_center, rc=rc)
571  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
572  call error_handler("IN FieldCreate", rc)
573 
574  print*,"- CALL FieldCreate FOR TARGET SURFACE PRESSURE BEFORE ADJUSTMENT."
575  ps_b4adj_target_grid = esmf_fieldcreate(target_grid, &
576  typekind=esmf_typekind_r8, &
577  staggerloc=esmf_staggerloc_center, rc=rc)
578  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
579  call error_handler("IN FieldCreate", rc)
580 
581  end subroutine create_atm_b4adj_esmf_fields
582 
586  subroutine create_atm_esmf_fields
587 
588  implicit none
589 
590  integer :: rc, n
591 
592  allocate(tracers_target_grid(num_tracers))
593 
594  do n = 1, num_tracers
595  print*,"- CALL FieldCreate FOR TARGET GRID TRACERS ", trim(tracers(n))
596  tracers_target_grid(n) = esmf_fieldcreate(target_grid, &
597  typekind=esmf_typekind_r8, &
598  staggerloc=esmf_staggerloc_center, &
599  ungriddedlbound=(/1/), &
600  ungriddedubound=(/lev_target/), rc=rc)
601  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
602  call error_handler("IN FieldCreate", rc)
603  enddo
604 
605  print*,"- CALL FieldCreate FOR TARGET GRID TEMPERATURE."
606  temp_target_grid = esmf_fieldcreate(target_grid, &
607  typekind=esmf_typekind_r8, &
608  staggerloc=esmf_staggerloc_center, &
609  ungriddedlbound=(/1/), &
610  ungriddedubound=(/lev_target/), rc=rc)
611  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
612  call error_handler("IN FieldCreate", rc)
613 
614  print*,"- CALL FieldCreate FOR TARGET GRID PRESSURE."
615  pres_target_grid = esmf_fieldcreate(target_grid, &
616  typekind=esmf_typekind_r8, &
617  staggerloc=esmf_staggerloc_center, &
618  ungriddedlbound=(/1/), &
619  ungriddedubound=(/lev_target/), rc=rc)
620  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
621  call error_handler("IN FieldCreate", rc)
622 
623  print*,"- CALL FieldCreate FOR TARGET GRID VERTICAL VELOCITY."
624  dzdt_target_grid = esmf_fieldcreate(target_grid, &
625  typekind=esmf_typekind_r8, &
626  staggerloc=esmf_staggerloc_center, &
627  ungriddedlbound=(/1/), &
628  ungriddedubound=(/lev_target/), rc=rc)
629  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
630  call error_handler("IN FieldCreate", rc)
631 
632  print*,"- CALL FieldCreate FOR TARGET GRID DELP."
633  delp_target_grid = esmf_fieldcreate(target_grid, &
634  typekind=esmf_typekind_r8, &
635  staggerloc=esmf_staggerloc_center, &
636  ungriddedlbound=(/1/), &
637  ungriddedubound=(/lev_target/), rc=rc)
638  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
639  call error_handler("IN FieldCreate", rc)
640 
641  print*,"- CALL FieldCreate FOR TARGET GRID xwind."
642  xwind_target_grid = esmf_fieldcreate(target_grid, &
643  typekind=esmf_typekind_r8, &
644  staggerloc=esmf_staggerloc_center, &
645  ungriddedlbound=(/1/), &
646  ungriddedubound=(/lev_target/), rc=rc)
647  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
648  call error_handler("IN FieldCreate", rc)
649 
650  print*,"- CALL FieldCreate FOR TARGET GRID ywind."
651  ywind_target_grid = esmf_fieldcreate(target_grid, &
652  typekind=esmf_typekind_r8, &
653  staggerloc=esmf_staggerloc_center, &
654  ungriddedlbound=(/1/), &
655  ungriddedubound=(/lev_target/), rc=rc)
656  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
657  call error_handler("IN FieldCreate", rc)
658 
659  print*,"- CALL FieldCreate FOR TARGET GRID zwind."
660  zwind_target_grid = esmf_fieldcreate(target_grid, &
661  typekind=esmf_typekind_r8, &
662  staggerloc=esmf_staggerloc_center, &
663  ungriddedlbound=(/1/), &
664  ungriddedubound=(/lev_target/), rc=rc)
665  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
666  call error_handler("IN FieldCreate", rc)
667 
668  print*,"- CALL FieldCreate FOR TARGET HEIGHT."
669  zh_target_grid = esmf_fieldcreate(target_grid, &
670  typekind=esmf_typekind_r8, &
671  staggerloc=esmf_staggerloc_center, &
672  ungriddedlbound=(/1/), &
673  ungriddedubound=(/levp1_target/), rc=rc)
674  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
675  call error_handler("IN FieldCreate", rc)
676 
677  print*,"- CALL FieldCreate FOR TARGET U_S."
678  u_s_target_grid = esmf_fieldcreate(target_grid, &
679  typekind=esmf_typekind_r8, &
680  staggerloc=esmf_staggerloc_edge2, &
681  ungriddedlbound=(/1/), &
682  ungriddedubound=(/lev_target/), rc=rc)
683  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
684  call error_handler("IN FieldCreate", rc)
685 
686  print*,"- CALL FieldCreate FOR TARGET V_S."
687  v_s_target_grid = esmf_fieldcreate(target_grid, &
688  typekind=esmf_typekind_r8, &
689  staggerloc=esmf_staggerloc_edge2, &
690  ungriddedlbound=(/1/), &
691  ungriddedubound=(/lev_target/), rc=rc)
692  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
693  call error_handler("IN FieldCreate", rc)
694 
695  print*,"- CALL FieldCreate FOR TARGET xwind_S."
696  xwind_s_target_grid = esmf_fieldcreate(target_grid, &
697  typekind=esmf_typekind_r8, &
698  staggerloc=esmf_staggerloc_edge2, &
699  ungriddedlbound=(/1/), &
700  ungriddedubound=(/lev_target/), rc=rc)
701  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
702  call error_handler("IN FieldCreate", rc)
703 
704  print*,"- CALL FieldCreate FOR TARGET ywind_S."
705  ywind_s_target_grid = esmf_fieldcreate(target_grid, &
706  typekind=esmf_typekind_r8, &
707  staggerloc=esmf_staggerloc_edge2, &
708  ungriddedlbound=(/1/), &
709  ungriddedubound=(/lev_target/), rc=rc)
710  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
711  call error_handler("IN FieldCreate", rc)
712 
713  print*,"- CALL FieldCreate FOR TARGET zwind_S."
714  zwind_s_target_grid = esmf_fieldcreate(target_grid, &
715  typekind=esmf_typekind_r8, &
716  staggerloc=esmf_staggerloc_edge2, &
717  ungriddedlbound=(/1/), &
718  ungriddedubound=(/lev_target/), rc=rc)
719  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
720  call error_handler("IN FieldCreate", rc)
721 
722  print*,"- CALL FieldCreate FOR TARGET U_W."
723  u_w_target_grid = esmf_fieldcreate(target_grid, &
724  typekind=esmf_typekind_r8, &
725  staggerloc=esmf_staggerloc_edge1, &
726  ungriddedlbound=(/1/), &
727  ungriddedubound=(/lev_target/), rc=rc)
728  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
729  call error_handler("IN FieldCreate", rc)
730 
731  print*,"- CALL FieldCreate FOR TARGET V_W."
732  v_w_target_grid = esmf_fieldcreate(target_grid, &
733  typekind=esmf_typekind_r8, &
734  staggerloc=esmf_staggerloc_edge1, &
735  ungriddedlbound=(/1/), &
736  ungriddedubound=(/lev_target/), rc=rc)
737  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
738  call error_handler("IN FieldCreate", rc)
739 
740  print*,"- CALL FieldCreate FOR TARGET xwind_W."
741  xwind_w_target_grid = esmf_fieldcreate(target_grid, &
742  typekind=esmf_typekind_r8, &
743  staggerloc=esmf_staggerloc_edge1, &
744  ungriddedlbound=(/1/), &
745  ungriddedubound=(/lev_target/), rc=rc)
746  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
747  call error_handler("IN FieldCreate", rc)
748 
749  print*,"- CALL FieldCreate FOR TARGET ywind_W."
750  ywind_w_target_grid = esmf_fieldcreate(target_grid, &
751  typekind=esmf_typekind_r8, &
752  staggerloc=esmf_staggerloc_edge1, &
753  ungriddedlbound=(/1/), &
754  ungriddedubound=(/lev_target/), rc=rc)
755  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
756  call error_handler("IN FieldCreate", rc)
757 
758  print*,"- CALL FieldCreate FOR TARGET zwind_W."
759  zwind_w_target_grid = esmf_fieldcreate(target_grid, &
760  typekind=esmf_typekind_r8, &
761  staggerloc=esmf_staggerloc_edge1, &
762  ungriddedlbound=(/1/), &
763  ungriddedubound=(/lev_target/), rc=rc)
764  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
765  call error_handler("IN FieldCreate", rc)
766 
767  print*,"- CALL FieldCreate FOR TARGET SURFACE PRESSURE."
768  ps_target_grid = esmf_fieldcreate(target_grid, &
769  typekind=esmf_typekind_r8, &
770  staggerloc=esmf_staggerloc_center, rc=rc)
771  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
772  call error_handler("IN FieldCreate", rc)
773 
774  end subroutine create_atm_esmf_fields
775 
779  subroutine convert_winds_to_uv
780 
781  implicit none
782 
783  integer :: clb(3), cub(3)
784  integer :: i, j, k, rc
785 
786  real(esmf_kind_r8), pointer :: latptr(:,:)
787  real(esmf_kind_r8), pointer :: lonptr(:,:)
788  real(esmf_kind_r8), pointer :: uptr(:,:,:)
789  real(esmf_kind_r8), pointer :: vptr(:,:,:)
790  real(esmf_kind_r8), pointer :: xwindptr(:,:,:)
791  real(esmf_kind_r8), pointer :: ywindptr(:,:,:)
792  real(esmf_kind_r8), pointer :: zwindptr(:,:,:)
793  real(esmf_kind_r8) :: latrad, lonrad
794 
795 !-----------------------------------------------------------------------------------
796 ! Convert from 3-d cartesian to 2-cartesian winds
797 !-----------------------------------------------------------------------------------
798 
799  print*,'- CONVERT WINDS.'
800 
801  print*,"- CALL FieldGet FOR xwind_S."
802  call esmf_fieldget(xwind_s_target_grid, &
803  computationallbound=clb, &
804  computationalubound=cub, &
805  farrayptr=xwindptr, rc=rc)
806  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
807  call error_handler("IN FieldGet", rc)
808 
809  print*,"- CALL FieldGet FOR ywind_S."
810  call esmf_fieldget(ywind_s_target_grid, &
811  farrayptr=ywindptr, rc=rc)
812  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
813  call error_handler("IN FieldGet", rc)
814 
815  print*,"- CALL FieldGet FOR zwind_S."
816  call esmf_fieldget(zwind_s_target_grid, &
817  farrayptr=zwindptr, rc=rc)
818  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
819  call error_handler("IN FieldGet", rc)
820 
821  print*,"- CALL FieldGet FOR U_S."
822  call esmf_fieldget(u_s_target_grid, &
823  farrayptr=uptr, rc=rc)
824  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
825  call error_handler("IN FieldGet", rc)
826 
827  print*,"- CALL FieldGet FOR V_S."
828  call esmf_fieldget(v_s_target_grid, &
829  farrayptr=vptr, rc=rc)
830  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
831  call error_handler("IN FieldGet", rc)
832 
833  print*,"- CALL FieldGet FOR LATITUDE_S."
834  call esmf_fieldget(latitude_s_target_grid, &
835  farrayptr=latptr, rc=rc)
836  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
837  call error_handler("IN FieldGet", rc)
838 
839  print*,"- CALL FieldGet FOR LONGITUDE_S."
840  call esmf_fieldget(longitude_s_target_grid, &
841  farrayptr=lonptr, rc=rc)
842  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
843  call error_handler("IN FieldGet", rc)
844 
845  do i = clb(1), cub(1)
846  do j = clb(2), cub(2)
847  latrad = latptr(i,j) * acos(-1.) / 180.0
848  lonrad = lonptr(i,j) * acos(-1.) / 180.0
849  do k = clb(3), cub(3)
850  uptr(i,j,k) = xwindptr(i,j,k) * cos(lonrad) + ywindptr(i,j,k) * sin(lonrad)
851  vptr(i,j,k) = -xwindptr(i,j,k) * sin(latrad) * sin(lonrad) + &
852  ywindptr(i,j,k) * sin(latrad) * cos(lonrad) + &
853  zwindptr(i,j,k) * cos(latrad)
854  enddo
855  enddo
856  enddo
857 
858  print*,"- CALL FieldGet FOR xwind_w."
859  call esmf_fieldget(xwind_w_target_grid, &
860  computationallbound=clb, &
861  computationalubound=cub, &
862  farrayptr=xwindptr, rc=rc)
863  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
864  call error_handler("IN FieldGet", rc)
865 
866  print*,"- CALL FieldGet FOR ywind_w."
867  call esmf_fieldget(ywind_w_target_grid, &
868  farrayptr=ywindptr, rc=rc)
869  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
870  call error_handler("IN FieldGet", rc)
871 
872  print*,"- CALL FieldGet FOR zwind_w."
873  call esmf_fieldget(zwind_w_target_grid, &
874  farrayptr=zwindptr, rc=rc)
875  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
876  call error_handler("IN FieldGet", rc)
877 
878  print*,"- CALL FieldGet FOR U_W."
879  call esmf_fieldget(u_w_target_grid, &
880  farrayptr=uptr, rc=rc)
881  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
882  call error_handler("IN FieldGet", rc)
883 
884  print*,"- CALL FieldGet FOR V_W."
885  call esmf_fieldget(v_w_target_grid, &
886  farrayptr=vptr, rc=rc)
887  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
888  call error_handler("IN FieldGet", rc)
889 
890  print*,"- CALL FieldGet FOR LATITUDE_W."
891  call esmf_fieldget(latitude_w_target_grid, &
892  farrayptr=latptr, rc=rc)
893  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
894  call error_handler("IN FieldGet", rc)
895 
896  print*,"- CALL FieldGet FOR LONGITUDE_W."
897  call esmf_fieldget(longitude_w_target_grid, &
898  farrayptr=lonptr, rc=rc)
899  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
900  call error_handler("IN FieldGet", rc)
901 
902  do i = clb(1), cub(1)
903  do j = clb(2), cub(2)
904  latrad = latptr(i,j) * acos(-1.) / 180.0
905  lonrad = lonptr(i,j) * acos(-1.) / 180.0
906  do k = clb(3), cub(3)
907  uptr(i,j,k) = xwindptr(i,j,k) * cos(lonrad) + ywindptr(i,j,k) * sin(lonrad)
908  vptr(i,j,k) = -xwindptr(i,j,k) * sin(latrad) * sin(lonrad) + &
909  ywindptr(i,j,k) * sin(latrad) * cos(lonrad) + &
910  zwindptr(i,j,k) * cos(latrad)
911  enddo
912  enddo
913  enddo
914 
915  end subroutine convert_winds_to_uv
916 
948  subroutine newpr1(localpet)
949  implicit none
950 
951  integer, intent(in) :: localpet
952 
953  integer :: idsl, idvc, rc
954  integer :: i, j, k, clb(3), cub(3)
955 
956  real(esmf_kind_r8), parameter :: rd=287.05
957  real(esmf_kind_r8), parameter :: cp=1004.6
958  real(esmf_kind_r8), parameter :: rocp=rd/cp
959  real(esmf_kind_r8), parameter :: rocp1=rocp+1
960  real(esmf_kind_r8), parameter :: rocpr=1/rocp
961 
962  real(esmf_kind_r8), pointer :: delp_ptr(:,:,:)
963  real(esmf_kind_r8), pointer :: pptr(:,:,:) ! adjusted 3-d p.
964  real(esmf_kind_r8), pointer :: psptr(:,:) ! adjusted surface p.
965  real(esmf_kind_r8) :: ak, bk
966  real(esmf_kind_r8), allocatable :: pi(:,:,:)
967 
968  print*,"COMPUTE 3-D PRESSURE FROM ADJUSTED SURFACE PRESSURE."
969 
970  idvc = 2 ! hard wire for now.
971  idsl = 2 ! hard wire for now.
972 
973  print*,"- CALL FieldGet FOR 3-D PRES."
974  call esmf_fieldget(pres_target_grid, &
975  computationallbound=clb, &
976  computationalubound=cub, &
977  farrayptr=pptr, rc=rc)
978  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
979  call error_handler("IN FieldGet", rc)
980 
981  print*,"- CALL FieldGet FOR DELP."
982  call esmf_fieldget(delp_target_grid, &
983  computationallbound=clb, &
984  computationalubound=cub, &
985  farrayptr=delp_ptr, rc=rc)
986  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
987  call error_handler("IN FieldGet", rc)
988 
989  print*,"- CALL FieldGet FOR SURFACE PRESSURE AFTER ADJUSTMENT"
990  call esmf_fieldget(ps_target_grid, &
991  farrayptr=psptr, 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  allocate(pi(clb(1):cub(1),clb(2):cub(2),1:levp1_target))
996 
997  if(idvc.eq.2) then
998  do k=1,levp1_target
999  ak = vcoord_target(k,1)
1000  bk = vcoord_target(k,2)
1001  do i= clb(1), cub(1)
1002  do j= clb(2), cub(2)
1003  pi(i,j,k) = ak + bk*psptr(i,j)
1004  enddo
1005  enddo
1006  enddo
1007  do k=1,lev_target
1008  do i= clb(1), cub(1)
1009  do j= clb(2), cub(2)
1010  delp_ptr(i,j,k) = pi(i,j,k) - pi(i,j,k+1)
1011  enddo
1012  enddo
1013  enddo
1014  else
1015  call error_handler("PROGRAM ONLY WORKS WITH IDVC 2", 1)
1016  endif
1017 
1018  if(idsl.eq.2) then
1019  do k=1,lev_target
1020  do i= clb(1), cub(1)
1021  do j= clb(2), cub(2)
1022  pptr(i,j,k) = (pi(i,j,k)+pi(i,j,k+1))/2.0
1023  enddo
1024  enddo
1025  enddo
1026  else
1027  do k=1,lev_target
1028  do i= clb(1), cub(1)
1029  do j= clb(2), cub(2)
1030  pptr(i,j,k) = ((pi(i,j,k)**rocp1-pi(i,j,k+1)**rocp1)/ &
1031  (rocp1*(pi(i,j,k)-pi(i,j,k+1))))**rocpr
1032  enddo
1033  enddo
1034  enddo
1035  endif
1036 
1037  deallocate(pi)
1038 
1039  if (localpet == 0) then
1040  print*,'new pres ',pptr(clb(1),clb(2),:)
1041  print*,'delp ',delp_ptr(clb(1),clb(2),:)
1042  endif
1043 
1044  end subroutine newpr1
1045 
1059  subroutine newps(localpet)
1060 
1061  implicit none
1062 
1063  integer, intent(in) :: localpet
1064  integer :: i, j, k, ii
1065  integer :: clb(3), cub(3), ls, rc
1066 
1067  real(esmf_kind_r8), pointer :: pptr(:,:,:)
1068  real(esmf_kind_r8), pointer :: psptr(:,:)
1069  real(esmf_kind_r8), pointer :: psnewptr(:,:) ! adjusted surface p.
1070  real(esmf_kind_r8), pointer :: tptr(:,:,:)
1071  real(esmf_kind_r8), pointer :: qptr(:,:,:)
1072  real(esmf_kind_r8), pointer :: zsptr(:,:)
1073  real(esmf_kind_r8), pointer :: zsnewptr(:,:)
1074  real(esmf_kind_r8), allocatable :: zu(:,:)
1075  real(esmf_kind_r8), parameter :: beta=-6.5e-3
1076  real(esmf_kind_r8), parameter :: epsilon=1.e-9
1077  real(esmf_kind_r8), parameter :: g=9.80665
1078  real(esmf_kind_r8), parameter :: rd=287.05
1079  real(esmf_kind_r8), parameter :: rv=461.50
1080  real(esmf_kind_r8), parameter :: gor=g/rd
1081  real(esmf_kind_r8), parameter :: fv=rv/rd-1.
1082  real(esmf_kind_r8) :: ftv, fgam, apu, fz0
1083  real(esmf_kind_r8) :: atvu, atv, fz1, fp0
1084  real(esmf_kind_r8) :: apd, azd, agam, azu
1085  real(esmf_kind_r8) :: atvd, fp1, gamma, pu
1086  real(esmf_kind_r8) :: tvu, pd, tvd
1087  real(esmf_kind_r8) :: at, aq, ap, az
1088 
1089  ftv(at,aq)=at*(1+fv*aq)
1090  fgam(apu,atvu,apd,atvd)=-gor*log(atvd/atvu)/log(apd/apu)
1091  fz0(ap,atv,azd,apd)=azd+atv/gor*log(apd/ap)
1092  fz1(ap,atv,azd,apd,agam)=azd-atv/agam*((apd/ap)**(-agam/gor)-1)
1093  fp0(az,azu,apu,atvu)=apu*exp(-gor/atvu*(az-azu))
1094  fp1(az,azu,apu,atvu,agam)=apu*(1+agam/atvu*(az-azu))**(-gor/agam)
1095 
1096  print*,"- ADJUST SURFACE PRESSURE FOR NEW TERRAIN."
1097 
1098  print*,"- CALL FieldGet FOR 3-D PRES."
1099  call esmf_fieldget(pres_b4adj_target_grid, &
1100  computationallbound=clb, &
1101  computationalubound=cub, &
1102  farrayptr=pptr, rc=rc)
1103  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1104  call error_handler("IN FieldGet", rc)
1105 
1106  if(localpet==0) then
1107  print*,'old pres ',pptr(clb(1),clb(2),:)
1108  endif
1109 
1110  print*,"- CALL FieldGet FOR TEMPERATURE"
1111  call esmf_fieldget(temp_b4adj_target_grid, &
1112  farrayptr=tptr, rc=rc)
1113  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1114  call error_handler("IN FieldGet", rc)
1115 
1116 ! Find specific humidity in the array of tracer fields.
1117 
1118  do ii = 1, num_tracers
1119  if (trim(tracers(ii)) == "sphum") exit
1120  enddo
1121 
1122  print*,"- CALL FieldGet FOR SPECIFIC HUMIDITY"
1123  call esmf_fieldget(tracers_b4adj_target_grid(ii), &
1124  farrayptr=qptr, rc=rc)
1125  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1126  call error_handler("IN FieldGet", rc)
1127 
1128  print*,"- CALL FieldGet FOR SURFACE PRESSURE BEFORE ADJUSTMENT"
1129  call esmf_fieldget(ps_b4adj_target_grid, &
1130  farrayptr=psptr, rc=rc)
1131  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1132  call error_handler("IN FieldGet", rc)
1133 
1134  print*,"- CALL FieldGet FOR SURFACE PRESSURE AFTER ADJUSTMENT"
1135  call esmf_fieldget(ps_target_grid, &
1136  farrayptr=psnewptr, rc=rc)
1137  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1138  call error_handler("IN FieldGet", rc)
1139 
1140  print*,"- CALL FieldGet FOR OLD TERRAIN"
1141  call esmf_fieldget(terrain_interp_to_target_grid, &
1142  farrayptr=zsptr, rc=rc)
1143  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1144  call error_handler("IN FieldGet", rc)
1145 
1146  print*,"- CALL FieldGet FOR NEW TERRAIN"
1147  call esmf_fieldget(terrain_target_grid, &
1148  farrayptr=zsnewptr, rc=rc)
1149  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1150  call error_handler("IN FieldGet", rc)
1151 
1152  allocate(zu(clb(1):cub(1),clb(2):cub(2)))
1153 
1154 !-----------------------------------------------------------------------------------
1155 ! Note, this routine was adapted from the spectral GFS which labeled the lowest
1156 ! model layer as '1'.
1157 !-----------------------------------------------------------------------------------
1158 
1159 !-----------------------------------------------------------------------------------
1160 ! Compute surface pressure below the original ground.
1161 !-----------------------------------------------------------------------------------
1162 
1163  ls=0
1164  k=1
1165  gamma=beta
1166  do i=clb(1), cub(1)
1167  do j=clb(2), cub(2)
1168  pu=pptr(i,j,k)
1169  tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1170  zu(i,j)=fz1(pu,tvu,zsptr(i,j),psptr(i,j),gamma)
1171  if(zsnewptr(i,j).le.zu(i,j)) then
1172  pu=pptr(i,j,k)
1173  tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1174  if(abs(gamma).gt.epsilon) then
1175  psnewptr(i,j)=fp1(zsnewptr(i,j),zu(i,j),pu,tvu,gamma)
1176  else
1177  psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1178  endif
1179  else
1180  psnewptr(i,j)=0
1181  ls=ls+1
1182  endif
1183  enddo
1184  enddo
1185 
1186 !-----------------------------------------------------------------------------------
1187 ! Compute surface pressure above the original ground.
1188 !-----------------------------------------------------------------------------------
1189 
1190  do k=2,cub(3)
1191  if(ls.gt.0) then
1192  do i=clb(1),cub(1)
1193  do j=clb(2),cub(2)
1194  if(psnewptr(i,j).eq.0) then
1195  pu=pptr(i,j,k)
1196  tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1197  pd=pptr(i,j,k-1)
1198  tvd=ftv(tptr(i,j,k-1),qptr(i,j,k-1))
1199  gamma=fgam(pu,tvu,pd,tvd)
1200  if(abs(gamma).gt.epsilon) then
1201  zu(i,j)=fz1(pu,tvu,zu(i,j),pd,gamma)
1202  else
1203  zu(i,j)=fz0(pu,tvu,zu(i,j),pd)
1204  endif
1205  if(zsnewptr(i,j).le.zu(i,j)) then
1206  if(abs(gamma).gt.epsilon) then
1207  psnewptr(i,j)=fp1(zsnewptr(i,j),zu(i,j),pu,tvu,gamma)
1208  else
1209  psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1210  endif
1211  ls=ls-1
1212  endif
1213  endif
1214  enddo
1215  enddo
1216  endif
1217  enddo
1218 
1219 !-----------------------------------------------------------------------------------
1220 ! Compute surface pressure over the top.
1221 !-----------------------------------------------------------------------------------
1222 
1223 
1224  if(ls.gt.0) then
1225  k=cub(3)
1226  gamma=0
1227  do i=clb(1),cub(1)
1228  do j=clb(2),cub(2)
1229  if(psnewptr(i,j).eq.0) then
1230  pu=pptr(i,j,k)
1231  tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1232  psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1233  endif
1234  enddo
1235  enddo
1236  endif
1237 
1238  deallocate(zu)
1239 
1240  if (localpet == 0) then
1241 ! do i=clb(1),cub(1)
1242 ! do j=clb(2),cub(2)
1243  do i=clb(1),clb(1)
1244  do j=clb(2),clb(2)
1245  print*,'sfcp adjust ',(zsnewptr(i,j)-zsptr(i,j)), psptr(i,j),psnewptr(i,j)
1246  enddo
1247  enddo
1248  endif
1249 
1250  end subroutine newps
1251 
1256  subroutine read_vcoord_info
1257  implicit none
1258 
1259  integer :: istat, n, k
1260 
1261  print*
1262  print*,"OPEN VERTICAL COORD FILE: ", trim(vcoord_file_target_grid)
1263  open(14, file=trim(vcoord_file_target_grid), form='formatted', iostat=istat, action='read')
1264  if (istat /= 0) then
1265  call error_handler("OPENING VERTICAL COORD FILE", istat)
1266  endif
1267 
1268  read(14, *, iostat=istat) nvcoord_target, lev_target
1269  if (istat /= 0) then
1270  call error_handler("READING VERTICAL COORD FILE", istat)
1271  endif
1272 
1273  levp1_target = lev_target + 1
1274 
1275  allocate(vcoord_target(levp1_target, nvcoord_target))
1276  read(14, *, iostat=istat) ((vcoord_target(n,k), k=1,nvcoord_target), n=1,levp1_target)
1277  if (istat /= 0) then
1278  call error_handler("READING VERTICAL COORD FILE", istat)
1279  endif
1280 
1281  print*
1282 
1283  close(14)
1284 
1285  end subroutine read_vcoord_info
1286 
1291  subroutine horiz_interp_thomp_mp_climo
1292 
1293  implicit none
1294 
1295  integer :: isrctermprocessing, rc
1296 
1297  type(esmf_regridmethod_flag) :: method
1298  type(esmf_routehandle) :: regrid_bl
1299 
1300  isrctermprocessing=1
1301 
1302  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNIFA BEFORE ADJUSTMENT."
1303  qnifa_climo_b4adj_target_grid = esmf_fieldcreate(target_grid, &
1304  typekind=esmf_typekind_r8, &
1305  staggerloc=esmf_staggerloc_center, &
1306  ungriddedlbound=(/1/), &
1307  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1308  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1309  call error_handler("IN FieldCreate", rc)
1310 
1311  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNWFA BEFORE ADJUSTMENT."
1312  qnwfa_climo_b4adj_target_grid = esmf_fieldcreate(target_grid, &
1313  typekind=esmf_typekind_r8, &
1314  staggerloc=esmf_staggerloc_center, &
1315  ungriddedlbound=(/1/), &
1316  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1317  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1318  call error_handler("IN FieldCreate", rc)
1319 
1320  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO PRESSURE BEFORE ADJUSTMENT."
1321  thomp_pres_climo_b4adj_target_grid = esmf_fieldcreate(target_grid, &
1322  typekind=esmf_typekind_r8, &
1323  staggerloc=esmf_staggerloc_center, &
1324  ungriddedlbound=(/1/), &
1325  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1326  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1327  call error_handler("IN FieldCreate", rc)
1328 
1329  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNIFA."
1330  qnifa_climo_target_grid = esmf_fieldcreate(target_grid, &
1331  typekind=esmf_typekind_r8, &
1332  staggerloc=esmf_staggerloc_center, &
1333  ungriddedlbound=(/1/), &
1334  ungriddedubound=(/lev_target/), rc=rc)
1335  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1336  call error_handler("IN FieldCreate", rc)
1337 
1338  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNWFA."
1339  qnwfa_climo_target_grid = esmf_fieldcreate(target_grid, &
1340  typekind=esmf_typekind_r8, &
1341  staggerloc=esmf_staggerloc_center, &
1342  ungriddedlbound=(/1/), &
1343  ungriddedubound=(/lev_target/), rc=rc)
1344  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1345  call error_handler("IN FieldCreate", rc)
1346 
1347  print*,"- CALL FieldRegridStore FOR THOMPSON CLIMO FIELDS."
1348 
1349  method=esmf_regridmethod_bilinear
1350 
1351  call esmf_fieldregridstore(qnifa_climo_input_grid, &
1352  qnifa_climo_b4adj_target_grid, &
1353  polemethod=esmf_polemethod_allavg, &
1354  srctermprocessing=isrctermprocessing, &
1355  routehandle=regrid_bl, &
1356  regridmethod=method, rc=rc)
1357  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1358  call error_handler("IN FieldRegridStore", rc)
1359 
1360  print*,"- CALL Field_Regrid FOR THOMP CLIMO QNIFA."
1361  call esmf_fieldregrid(qnifa_climo_input_grid, &
1362  qnifa_climo_b4adj_target_grid, &
1363  routehandle=regrid_bl, &
1364  termorderflag=esmf_termorder_srcseq, rc=rc)
1365  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1366  call error_handler("IN FieldRegrid", rc)
1367 
1368  print*,"- CALL Field_Regrid FOR THOMP CLIMO QNWFA."
1369  call esmf_fieldregrid(qnwfa_climo_input_grid, &
1370  qnwfa_climo_b4adj_target_grid, &
1371  routehandle=regrid_bl, &
1372  termorderflag=esmf_termorder_srcseq, rc=rc)
1373  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1374  call error_handler("IN FieldRegrid", rc)
1375 
1376  print*,"- CALL Field_Regrid FOR THOMP PRESSURE."
1377  call esmf_fieldregrid(thomp_pres_climo_input_grid, &
1378  thomp_pres_climo_b4adj_target_grid, &
1379  routehandle=regrid_bl, &
1380  termorderflag=esmf_termorder_srcseq, rc=rc)
1381  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1382  call error_handler("IN FieldRegrid", rc)
1383 
1384  print*,"- CALL FieldRegridRelease."
1385  call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
1386  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1387  call error_handler("IN FieldRegridRelease", rc)
1388 
1389 !-----------------------------------------------------------------------------------
1390 ! Free up input data memory.
1391 !-----------------------------------------------------------------------------------
1392 
1393  call cleanup_thomp_mp_climo_input_data
1394 
1395  end subroutine horiz_interp_thomp_mp_climo
1396 
1403  SUBROUTINE vintg_thomp_mp_climo
1404 
1405  implicit none
1406 
1407  INTEGER :: clb(3), cub(3), rc
1408  INTEGER :: im, km1, km2, nt
1409  INTEGER :: i, j, k
1410 
1411  REAL(ESMF_KIND_R8), ALLOCATABLE :: z1(:,:,:), z2(:,:,:)
1412  REAL(ESMF_KIND_R8), ALLOCATABLE :: c1(:,:,:,:),c2(:,:,:,:)
1413 
1414  REAL(ESMF_KIND_R8), POINTER :: qnifa1ptr(:,:,:) ! input
1415  REAL(ESMF_KIND_R8), POINTER :: qnifa2ptr(:,:,:) ! target
1416  REAL(ESMF_KIND_R8), POINTER :: qnwfa1ptr(:,:,:) ! input
1417  REAL(ESMF_KIND_R8), POINTER :: qnwfa2ptr(:,:,:) ! target
1418  REAL(ESMF_KIND_R8), POINTER :: p1ptr(:,:,:) ! input pressure
1419  REAL(ESMF_KIND_R8), POINTER :: p2ptr(:,:,:) ! target pressure
1420 
1421  print*,"- VERTICALY INTERPOLATE THOMP MP CLIMO TRACERS."
1422 
1423  print*,"- CALL FieldGet FOR 3-D THOMP PRES."
1424  call esmf_fieldget(thomp_pres_climo_b4adj_target_grid, &
1425  computationallbound=clb, &
1426  computationalubound=cub, &
1427  farrayptr=p1ptr, rc=rc)
1428  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1429  call error_handler("IN FieldGet", rc)
1430 
1431 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1432 ! The '1'/'2' arrays hold fields before/after interpolation.
1433 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1434 
1435  nt= 2 ! number of thomp tracers
1436 
1437  ALLOCATE(z1(clb(1):cub(1),clb(2):cub(2),lev_thomp_mp_climo))
1438  ALLOCATE(z2(clb(1):cub(1),clb(2):cub(2),lev_target))
1439  ALLOCATE(c1(clb(1):cub(1),clb(2):cub(2),lev_thomp_mp_climo,nt))
1440  ALLOCATE(c2(clb(1):cub(1),clb(2):cub(2),lev_target,nt))
1441 
1442  z1 = -log(p1ptr)
1443 
1444  print*,"- CALL FieldGet FOR 3-D ADJUSTED PRESS"
1445  call esmf_fieldget(pres_target_grid, &
1446  farrayptr=p2ptr, rc=rc)
1447  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1448  call error_handler("IN FieldGet", rc)
1449 
1450  z2 = -log(p2ptr)
1451 
1452 !print*,'pres check 1 ', p1ptr(clb(1),clb(2),:)
1453 !print*,'pres check 2 ', p2ptr(clb(1),clb(2),:)
1454 
1455  print*,"- CALL FieldGet FOR qnifa before vertical adjustment."
1456  call esmf_fieldget(qnifa_climo_b4adj_target_grid, &
1457  farrayptr=qnifa1ptr, rc=rc)
1458  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1459  call error_handler("IN FieldGet", rc)
1460 
1461  c1(:,:,:,1) = qnifa1ptr(:,:,:)
1462 
1463  print*,"- CALL FieldGet FOR qnwfa before vertical adjustment."
1464  call esmf_fieldget(qnwfa_climo_b4adj_target_grid, &
1465  farrayptr=qnwfa1ptr, rc=rc)
1466  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1467  call error_handler("IN FieldGet", rc)
1468 
1469  c1(:,:,:,2) = qnwfa1ptr(:,:,:)
1470 
1471 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1472 ! PERFORM LAGRANGIAN ONE-DIMENSIONAL INTERPOLATION
1473 ! THAT IS 4TH-ORDER IN INTERIOR, 2ND-ORDER IN OUTSIDE INTERVALS
1474 ! AND 1ST-ORDER FOR EXTRAPOLATION.
1475 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1476 
1477  im = (cub(1)-clb(1)+1) * (cub(2)-clb(2)+1)
1478  km1= lev_thomp_mp_climo
1479  km2= lev_target
1480 
1481  CALL terp3(im,1,1,1,1,nt,(im*km1),(im*km2), &
1482  km1,im,im,z1,c1,km2,im,im,z2,c2)
1483 
1484  print*,"- CALL FieldGet FOR ADJUSTED climo qnifa."
1485  call esmf_fieldget(qnifa_climo_target_grid, &
1486  farrayptr=qnifa2ptr, rc=rc)
1487  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1488  call error_handler("IN FieldGet", rc)
1489 
1490  print*,"- CALL FieldGet FOR ADJUSTED climo qnwfa."
1491  call esmf_fieldget(qnwfa_climo_target_grid, &
1492  farrayptr=qnwfa2ptr, rc=rc)
1493  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1494  call error_handler("IN FieldGet", rc)
1495 
1496  DO k=1,lev_target
1497  DO i=clb(1),cub(1)
1498  DO j=clb(2),cub(2)
1499  qnifa2ptr(i,j,k) = c2(i,j,k,1)
1500  qnwfa2ptr(i,j,k) = c2(i,j,k,2)
1501  ENDDO
1502  ENDDO
1503  ENDDO
1504 
1505  DEALLOCATE (z1, z2, c1, c2)
1506 
1507  call esmf_fielddestroy(qnifa_climo_b4adj_target_grid, rc=rc)
1508  call esmf_fielddestroy(qnwfa_climo_b4adj_target_grid, rc=rc)
1509  call esmf_fielddestroy(thomp_pres_climo_b4adj_target_grid, rc=rc)
1510 
1511  END SUBROUTINE vintg_thomp_mp_climo
1512 
1513 
1527  SUBROUTINE vintg_wam (YEAR,MONTH,DAY,HOUR,PF)
1528 
1529  IMPLICIT NONE
1530 
1531  include 'mpif.h'
1532 
1533  INTEGER, INTENT(IN) :: year,month,day,hour
1534  CHARACTER(*), INTENT(IN) :: pf
1535 
1536  REAL(ESMF_KIND_R8), PARAMETER :: amo = 15.9994 ! molecular weight of o
1537  REAL(ESMF_KIND_R8), PARAMETER :: amo2 = 31.999 !molecular weight of o2
1538  REAL(ESMF_KIND_R8), PARAMETER :: amn2 = 28.013 !molecular weight of n2
1539 
1540  REAL(ESMF_KIND_R8) :: coe,wfun(10),deglat,hold
1541  REAL(ESMF_KIND_R8) :: summass,qvmass,o3mass
1542  INTEGER :: i, j, k, ii, clb(3), cub(3), rc, kref
1543  INTEGER :: idat(8),jdow,jday,icday
1544 
1545  REAL(ESMF_KIND_R8), ALLOCATABLE :: temp(:),on(:),o2n(:),n2n(:),prmb(:)
1546 
1547  REAL(ESMF_KIND_R8), POINTER :: latptr(:,:) ! output latitude
1548  REAL(ESMF_KIND_R8), POINTER :: p1ptr(:,:,:) ! input pressure
1549  REAL(ESMF_KIND_R8), POINTER :: p2ptr(:,:,:) ! output pressure
1550  REAL(ESMF_KIND_R8), POINTER :: dzdt2ptr(:,:,:) ! output vvel
1551  REAL(ESMF_KIND_R8), POINTER :: t2ptr(:,:,:) ! output temperature
1552  REAL(ESMF_KIND_R8), POINTER :: q2ptr(:,:,:) ! output tracer
1553  REAL(ESMF_KIND_R8), POINTER :: qvptr(:,:,:) ! output tracer
1554  REAL(ESMF_KIND_R8), POINTER :: qoptr(:,:,:) ! output tracer
1555  REAL(ESMF_KIND_R8), POINTER :: o2ptr(:,:,:) ! output tracer
1556  REAL(ESMF_KIND_R8), POINTER :: o3ptr(:,:,:) ! output tracer
1557  REAL(ESMF_KIND_R8), POINTER :: xwind2ptr(:,:,:) ! output wind (x component)
1558  REAL(ESMF_KIND_R8), POINTER :: ywind2ptr(:,:,:) ! output wind (y component)
1559  REAL(ESMF_KIND_R8), POINTER :: zwind2ptr(:,:,:) ! output wind (z component)
1560 
1561 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1562 
1563  print*,"VINTG_WAM:- VERTICALY EXTEND FIELDS FOR WAM COLD START."
1564 
1565 ! prepare date
1566  idat = 0
1567  jdow = 0
1568  jday = 0
1569  icday = 0
1570  idat(1)=year
1571  idat(2)=month
1572  idat(3)=day
1573  idat(5)=hour
1574  CALL w3doxdat(idat,jdow,icday,jday)
1575  print *,"VINTG_WAM: WAM START DATE FOR ICDAY=",icday
1576 
1577 ! prepare weighting function
1578  DO k=1,10
1579  wfun(k) = (k-1.0) / 9.0
1580  ENDDO
1581 
1582  ALLOCATE(temp(lev_target))
1583  ALLOCATE(prmb(lev_target))
1584  ALLOCATE( on(lev_target))
1585  ALLOCATE( o2n(lev_target))
1586  ALLOCATE( n2n(lev_target))
1587 
1588 ! p1 (pascal)
1589  print*,"VINTG_WAM:- CALL FieldGet FOR 3-D PRES."
1590  call esmf_fieldget(pres_b4adj_target_grid, &
1591  computationallbound=clb, &
1592  computationalubound=cub, &
1593  farrayptr=p1ptr, rc=rc)
1594  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1595  call error_handler("IN FieldGet", rc)
1596 !print*,"VINTG_WAM: p1ptr ",(p1ptr(1,1,k),k=1,LEV_INPUT)
1597 
1598 ! p2 (pascal)
1599  print*,"VINTG_WAM:- CALL FieldGet FOR 3-D ADJUSTED PRESS"
1600  call esmf_fieldget(pres_target_grid, &
1601  farrayptr=p2ptr, rc=rc)
1602  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1603  call error_handler("IN FieldGet", rc)
1604 !print*,"VINTG_WAM: p2ptr ",(p2ptr(1,1,k),k=1,LEV_TARGET)
1605 
1606 ! latitude in degree
1607  print*,"VINTG_WAM - CALL FieldGet FOR LATITUDE_S."
1608  call esmf_fieldget(latitude_s_target_grid, &
1609  farrayptr=latptr, rc=rc)
1610  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1611  call error_handler("IN FieldGet", rc)
1612 !print*,"VINTG_WAM: latptr ",(latptr(1,j),j=clb(2),cub(2))
1613 
1614 ! temp
1615  print*,"VINTG_WAM:- CALL FieldGet FOR 3-D ADJUSTED TEMP."
1616  call esmf_fieldget(temp_target_grid, &
1617  farrayptr=t2ptr, rc=rc)
1618  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1619  call error_handler("IN FieldGet", rc)
1620 
1621 ! dzdt
1622  print*,"VINTG_WAM:- CALL FieldGet FOR ADJUSTED VERTICAL VELOCITY."
1623  call esmf_fieldget(dzdt_target_grid, &
1624  farrayptr=dzdt2ptr, rc=rc)
1625  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1626  call error_handler("IN FieldGet", rc)
1627 
1628 ! wind
1629  print*,"VINTG_WAM:- CALL FieldGet FOR ADJUSTED WIND COMPONENTS."
1630 
1631  call esmf_fieldget(xwind_target_grid, &
1632  farrayptr=xwind2ptr, rc=rc)
1633  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1634  call error_handler("IN FieldGet", rc)
1635 
1636  call esmf_fieldget(ywind_target_grid, &
1637  farrayptr=ywind2ptr, rc=rc)
1638  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1639  call error_handler("IN FieldGet", rc)
1640 
1641  call esmf_fieldget(zwind_target_grid, &
1642  farrayptr=zwind2ptr, rc=rc)
1643  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1644  call error_handler("IN FieldGet", rc)
1645 
1646 !
1647 ! determine vertical blending point and modified extrapolation values
1648 !
1649  DO i=clb(1),cub(1)
1650  DO j=clb(2),cub(2)
1651 
1652  DO k=1,lev_target
1653  IF(p2ptr(i,j,k).le.p1ptr(i,j,lev_input)) THEN
1654  kref =k-1
1655  EXIT
1656  ENDIF
1657  ENDDO
1658 !
1659  DO k=kref,lev_target
1660  coe = p2ptr(i,j,k) / p2ptr(i,j,kref)
1661  xwind2ptr(i,j,k) = coe*xwind2ptr(i,j,k)
1662  ywind2ptr(i,j,k) = coe*ywind2ptr(i,j,k)
1663  zwind2ptr(i,j,k) = coe*zwind2ptr(i,j,k)
1664  dzdt2ptr(i,j,k) = coe*dzdt2ptr(i,j,k)
1665  ENDDO
1666 
1667  ENDDO
1668  ENDDO
1669 
1670 !
1671 ! point necessary tracers
1672 !
1673  DO ii = 1, num_tracers
1674 
1675  print*,"VINTG_WAM:- CALL FieldGet FOR 3-D TRACER ", trim(tracers(ii))
1676  call esmf_fieldget(tracers_target_grid(ii), &
1677  farrayptr=q2ptr, rc=rc)
1678  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1679  call error_handler("IN FieldGet", rc)
1680 
1681  DO j=clb(2),cub(2)
1682  DO i=clb(1),cub(1)
1683  DO k=1,lev_target
1684  IF(p2ptr(i,j,k).le.p1ptr(i,j,lev_input)) THEN
1685  kref =k-1
1686  EXIT
1687  ENDIF
1688  ENDDO
1689 !
1690  DO k=kref,lev_target
1691  coe = min(1.0, p2ptr(i,j,k) / p2ptr(i,j,kref) )
1692  q2ptr(i,j,k) = coe * q2ptr(i,j,k)
1693  ENDDO
1694  ENDDO
1695  ENDDO
1696 
1697  IF (trim(tracers(ii)) == "sphum") qvptr => q2ptr
1698  IF (trim(tracers(ii)) == "spo" ) qoptr => q2ptr
1699  IF (trim(tracers(ii)) == "spo2" ) o2ptr => q2ptr
1700  IF (trim(tracers(ii)) == "spo3" ) o3ptr => q2ptr
1701 
1702  ENDDO
1703 
1704 !
1705 ! obtained wam gases distribution and temperature profile
1706 !
1707  DO i=clb(1),cub(1)
1708  DO j=clb(2),cub(2)
1709 !
1710  deglat = latptr(i,j)
1711  DO k=1,lev_target
1712  prmb(k) = p2ptr(i,j,k) * 0.01
1713  ENDDO
1714  CALL gettemp(icday,deglat,prmb,lev_target,pf,temp,on,o2n,n2n)
1715 !
1716  DO k=1,lev_target
1717  summass = on(k)*amo+o2n(k)*amo2+n2n(k)*amn2
1718  qvmass = summass*qvptr(i,j,k)/(1.-qvptr(i,j,k))
1719  summass = summass+qvmass
1720  o3mass = summass*o3ptr(i,j,k)
1721  summass = summass+o3mass
1722  hold = 1.0 / summass
1723  qoptr(i,j,k) = on(k)*amo *hold
1724  o2ptr(i,j,k) = o2n(k)*amo2*hold
1725  o3ptr(i,j,k) = o3mass * hold
1726  qvptr(i,j,k) = qvmass * hold
1727  ENDDO
1728 !
1729  DO k=1,lev_target
1730  IF(p2ptr(i,j,k).le.p1ptr(i,j,lev_input)) THEN
1731  kref =k-1
1732  EXIT
1733  ENDIF
1734  ENDDO
1735 !
1736  DO k=kref,lev_target
1737  t2ptr(i,j,k) = temp(k)
1738  ENDDO
1739  DO k=kref-10,kref-1
1740  t2ptr(i,j,k) = wfun(k-kref+11) * temp(k) + &
1741  (1.- wfun(k-kref+11)) * t2ptr(i,j,k)
1742  ENDDO
1743  ENDDO
1744  ENDDO
1745 
1746  DEALLOCATE (temp, prmb, on, o2n, n2n)
1747 
1748  END SUBROUTINE vintg_wam
1749 
1763  SUBROUTINE vintg
1764  use mpi_f08
1765 
1766  IMPLICIT NONE
1767 
1768  REAL(ESMF_KIND_R8), PARAMETER :: dltdz=-6.5e-3*287.05/9.80665
1769  REAL(ESMF_KIND_R8), PARAMETER :: dlpvdrt=-2.5e6/461.50
1770  REAL(ESMF_KIND_R8), PARAMETER :: one = 1.0_esmf_kind_r8
1771 
1772  INTEGER :: i, j, k, clb(3), cub(3), rc
1773  INTEGER :: im, km1, km2, nt, ii
1774 
1775  REAL(ESMF_KIND_R8) :: dz
1776  REAL(ESMF_KIND_R8), ALLOCATABLE :: z1(:,:,:), z2(:,:,:)
1777  REAL(ESMF_KIND_R8), ALLOCATABLE :: c1(:,:,:,:),c2(:,:,:,:)
1778 
1779  REAL(ESMF_KIND_R8), POINTER :: p1ptr(:,:,:) ! input pressure
1780  REAL(ESMF_KIND_R8), POINTER :: p2ptr(:,:,:) ! output pressure
1781  REAL(ESMF_KIND_R8), POINTER :: dzdt1ptr(:,:,:) ! input vvel
1782  REAL(ESMF_KIND_R8), POINTER :: dzdt2ptr(:,:,:) ! output vvel
1783  REAL(ESMF_KIND_R8), POINTER :: t1ptr(:,:,:) ! input temperature
1784  REAL(ESMF_KIND_R8), POINTER :: t2ptr(:,:,:) ! output temperature
1785  REAL(ESMF_KIND_R8), POINTER :: q1ptr(:,:,:) ! input tracer
1786  REAL(ESMF_KIND_R8), POINTER :: q2ptr(:,:,:) ! output tracer
1787  REAL(ESMF_KIND_R8), POINTER :: xwind1ptr(:,:,:) ! input wind (x component)
1788  REAL(ESMF_KIND_R8), POINTER :: ywind1ptr(:,:,:) ! input wind (y component)
1789  REAL(ESMF_KIND_R8), POINTER :: zwind1ptr(:,:,:) ! input wind (z component)
1790  REAL(ESMF_KIND_R8), POINTER :: xwind2ptr(:,:,:) ! output wind (x component)
1791  REAL(ESMF_KIND_R8), POINTER :: ywind2ptr(:,:,:) ! output wind (y component)
1792  REAL(ESMF_KIND_R8), POINTER :: zwind2ptr(:,:,:) ! output wind (z component)
1793 
1794 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1795 ! COMPUTE LOG PRESSURE INTERPOLATING COORDINATE
1796 ! AND COPY INPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS
1797 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1798 
1799  print*,"- VERTICALY INTERPOLATE FIELDS."
1800 
1801  print*,"- CALL FieldGet FOR 3-D PRES."
1802  call esmf_fieldget(pres_b4adj_target_grid, &
1803  computationallbound=clb, &
1804  computationalubound=cub, &
1805  farrayptr=p1ptr, rc=rc)
1806  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1807  call error_handler("IN FieldGet", rc)
1808 
1809 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1810 ! The '1'/'2' arrays hold fields before/after interpolation.
1811 ! Note the 'z' component of the horizontal wind will be treated as a
1812 ! tracer. So add one extra third dimension to these 3-d arrays.
1813 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1814 
1815  ALLOCATE(z1(clb(1):cub(1),clb(2):cub(2),lev_input))
1816  ALLOCATE(z2(clb(1):cub(1),clb(2):cub(2),lev_target))
1817  ALLOCATE(c1(clb(1):cub(1),clb(2):cub(2),lev_input,num_tracers_input+5))
1818  ALLOCATE(c2(clb(1):cub(1),clb(2):cub(2),lev_target,num_tracers_input+5))
1819 
1820  z1 = -log(p1ptr)
1821 
1822  print*,"- CALL FieldGet FOR 3-D ADJUSTED PRESS"
1823  call esmf_fieldget(pres_target_grid, &
1824  farrayptr=p2ptr, rc=rc)
1825  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1826  call error_handler("IN FieldGet", rc)
1827 
1828  z2 = -log(p2ptr)
1829 
1830  print*,"- CALL FieldGet FOR x WIND."
1831  call esmf_fieldget(xwind_b4adj_target_grid, &
1832  farrayptr=xwind1ptr, rc=rc)
1833  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1834  call error_handler("IN FieldGet", rc)
1835 
1836  c1(:,:,:,1) = xwind1ptr(:,:,:)
1837 
1838  print*,"- CALL FieldGet FOR y WIND."
1839  call esmf_fieldget(ywind_b4adj_target_grid, &
1840  farrayptr=ywind1ptr, rc=rc)
1841  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1842  call error_handler("IN FieldGet", rc)
1843 
1844  c1(:,:,:,2) = ywind1ptr(:,:,:)
1845 
1846  print*,"- CALL FieldGet FOR z WIND."
1847  call esmf_fieldget(zwind_b4adj_target_grid, &
1848  farrayptr=zwind1ptr, rc=rc)
1849  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1850  call error_handler("IN FieldGet", rc)
1851 
1852  c1(:,:,:,3) = zwind1ptr(:,:,:)
1853 
1854  print*,"- CALL FieldGet FOR VERTICAL VELOCITY."
1855  call esmf_fieldget(dzdt_b4adj_target_grid, &
1856  farrayptr=dzdt1ptr, rc=rc)
1857  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1858  call error_handler("IN FieldGet", rc)
1859 
1860  c1(:,:,:,4) = dzdt1ptr(:,:,:)
1861  print*,"MIN MAX W TARGETB4 IN VINTG = ", minval(dzdt1ptr(:,:,:)), maxval(dzdt1ptr(:,:,:))
1862 
1863  print*,"- CALL FieldGet FOR 3-D TEMP."
1864  call esmf_fieldget(temp_b4adj_target_grid, &
1865  farrayptr=t1ptr, rc=rc)
1866  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1867  call error_handler("IN FieldGet", rc)
1868 
1869  c1(:,:,:,5) = t1ptr(:,:,:)
1870 
1871  DO i = 1, num_tracers_input
1872 
1873  print*,"- CALL FieldGet FOR 3-D TRACERS ", trim(tracers(i))
1874  call esmf_fieldget(tracers_b4adj_target_grid(i), &
1875  farrayptr=q1ptr, rc=rc)
1876  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1877  call error_handler("IN FieldGet", rc)
1878 
1879  c1(:,:,:,5+i) = q1ptr(:,:,:)
1880 
1881  ENDDO
1882 
1883 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1884 ! PERFORM LAGRANGIAN ONE-DIMENSIONAL INTERPOLATION
1885 ! THAT IS 4TH-ORDER IN INTERIOR, 2ND-ORDER IN OUTSIDE INTERVALS
1886 ! AND 1ST-ORDER FOR EXTRAPOLATION.
1887 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1888 
1889  im = (cub(1)-clb(1)+1) * (cub(2)-clb(2)+1)
1890  km1= lev_input
1891  km2= lev_target
1892  nt= num_tracers_input + 1 ! treat 'z' wind as tracer.
1893 
1894  CALL terp3(im,1,1,1,1,4+nt,(im*km1),(im*km2), &
1895  km1,im,im,z1,c1,km2,im,im,z2,c2)
1896 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1897 ! COPY OUTPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS
1898 ! EXCEPT BELOW THE INPUT DOMAIN, LET TEMPERATURE INCREASE WITH A FIXED
1899 ! LAPSE RATE AND LET THE RELATIVE HUMIDITY REMAIN CONSTANT.
1900 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1901 
1902  print*,"- CALL FieldGet FOR 3-D ADJUSTED TEMP."
1903  call esmf_fieldget(temp_target_grid, &
1904  farrayptr=t2ptr, rc=rc)
1905  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1906  call error_handler("IN FieldGet", rc)
1907 
1908  print*,"- CALL FieldGet FOR ADJUSTED VERTICAL VELOCITY."
1909  call esmf_fieldget(dzdt_target_grid, &
1910  farrayptr=dzdt2ptr, rc=rc)
1911  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1912  call error_handler("IN FieldGet", rc)
1913 
1914  print*,"- CALL FieldGet FOR ADJUSTED xwind."
1915  call esmf_fieldget(xwind_target_grid, &
1916  farrayptr=xwind2ptr, rc=rc)
1917  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1918  call error_handler("IN FieldGet", rc)
1919 
1920  print*,"- CALL FieldGet FOR ADJUSTED ywind."
1921  call esmf_fieldget(ywind_target_grid, &
1922  farrayptr=ywind2ptr, rc=rc)
1923  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1924  call error_handler("IN FieldGet", rc)
1925 
1926  print*,"- CALL FieldGet FOR ADJUSTED zwind."
1927  call esmf_fieldget(zwind_target_grid, &
1928  farrayptr=zwind2ptr, rc=rc)
1929  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1930  call error_handler("IN FieldGet", rc)
1931 
1932  DO k=1,lev_target
1933  DO i=clb(1),cub(1)
1934  DO j=clb(2),cub(2)
1935  xwind2ptr(i,j,k)=c2(i,j,k,1)
1936  ywind2ptr(i,j,k)=c2(i,j,k,2)
1937  zwind2ptr(i,j,k)=c2(i,j,k,3)
1938  dzdt2ptr(i,j,k)=c2(i,j,k,4)
1939  dz=z2(i,j,k)-z1(i,j,1)
1940  IF(dz.GE.0) THEN
1941  t2ptr(i,j,k)=c2(i,j,k,5)
1942  ELSE
1943  t2ptr(i,j,k)=c1(i,j,1,5)*exp(dltdz*dz)
1944  ENDIF
1945  ENDDO
1946  ENDDO
1947  ENDDO
1948 
1949  DO ii = 1, num_tracers_input
1950 
1951  print*,"- CALL FieldGet FOR 3-D TRACER ", trim(tracers(ii))
1952  call esmf_fieldget(tracers_target_grid(ii), &
1953  farrayptr=q2ptr, rc=rc)
1954  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1955  call error_handler("IN FieldGet", rc)
1956 
1957  IF (trim(tracers(ii)) == "sphum") THEN ! specific humidity
1958 
1959  DO k=1,lev_target
1960  DO i=clb(1),cub(1)
1961  DO j=clb(2),cub(2)
1962  dz=z2(i,j,k)-z1(i,j,1)
1963  IF(dz.GE.0) THEN
1964  q2ptr(i,j,k) = c2(i,j,k,5+ii)
1965  ELSE
1966  q2ptr(i,j,k) = c1(i,j,1,5+ii)*exp(dlpvdrt*(one/t2ptr(i,j,k)-one/t1ptr(i,j,1))-dz)
1967  ENDIF
1968  ENDDO
1969  ENDDO
1970  ENDDO
1971 
1972  ELSE ! all other tracers
1973 
1974  DO k=1,lev_target
1975  DO i=clb(1),cub(1)
1976  DO j=clb(2),cub(2)
1977  q2ptr(i,j,k) = c2(i,j,k,5+ii)
1978  ENDDO
1979  ENDDO
1980  ENDDO
1981 
1982  ENDIF
1983 
1984  ENDDO
1985 
1986  DEALLOCATE (z1, z2, c1, c2)
1987 
1988  END SUBROUTINE vintg
1989 
2026  SUBROUTINE terp3(IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2, &
2027  km1,kxz1,kxq1,z1,q1,km2,kxz2,kxq2,z2,q2)
2028  IMPLICIT NONE
2029  INTEGER im,ixz1,ixq1,ixz2,ixq2,nm,nxq1,nxq2
2030  INTEGER km1,kxz1,kxq1,km2,kxz2,kxq2
2031  INTEGER i,k1,k2,n
2032  INTEGER k1s(im,km2)
2033  REAL(ESMF_KIND_R8), PARAMETER :: one = 1.0_esmf_kind_r8
2034  REAL(ESMF_KIND_R8) :: z1(1+(im-1)*ixz1+(km1-1)*kxz1)
2035  REAL(ESMF_KIND_R8) :: q1(1+(im-1)*ixq1+(km1-1)*kxq1+(nm-1)*nxq1)
2036  REAL(ESMF_KIND_R8) :: z2(1+(im-1)*ixz2+(km2-1)*kxz2)
2037  REAL(ESMF_KIND_R8) :: q2(1+(im-1)*ixq2+(km2-1)*kxq2+(nm-1)*nxq2)
2038 ! REAL(ESMF_KIND_R8) :: J2(1+(IM-1)*IXQ2+(KM2-1)*KXQ2+(NM-1)*NXQ2)
2039  REAL(ESMF_KIND_R8) :: ffa(im),ffb(im),ffc(im),ffd(im)
2040  REAL(ESMF_KIND_R8) :: gga(im),ggb(im),ggc(im),ggd(im)
2041  REAL(ESMF_KIND_R8) :: z1a,z1b,z1c,z1d,q1a,q1b,q1c,q1d,z2s,q2s
2042 ! REAL(ESMF_KIND_R8) :: J2S
2043 
2044 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2045 ! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT.
2046  CALL rsearch(im,km1,ixz1,kxz1,z1,km2,ixz2,kxz2,z2,1,im,k1s)
2047 
2048 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2049 ! GENERALLY INTERPOLATE CUBICALLY WITH MONOTONIC CONSTRAINT
2050 ! FROM TWO NEAREST INPUT POINTS ON EITHER SIDE OF THE OUTPUT POINT,
2051 ! BUT WITHIN THE TWO EDGE INTERVALS INTERPOLATE LINEARLY.
2052 ! KEEP THE OUTPUT FIELDS CONSTANT OUTSIDE THE INPUT DOMAIN.
2053 
2054 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(IM,IXZ1,IXQ1,IXZ2), &
2055 !$OMP& SHARED(IXQ2,NM,NXQ1,NXQ2,KM1,KXZ1,KXQ1,Z1,Q1,KM2,KXZ2), &
2056 !$OMP& SHARED(KXQ2,Z2,Q2,K1S)
2057  DO k2=1,km2
2058  DO i=1,im
2059  k1=k1s(i,k2)
2060  IF(k1.EQ.1.OR.k1.EQ.km1-1) THEN
2061  z2s=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
2062  z1a=z1(1+(i-1)*ixz1+(k1-1)*kxz1)
2063  z1b=z1(1+(i-1)*ixz1+(k1+0)*kxz1)
2064  ffa(i)=(z2s-z1b)/(z1a-z1b)
2065  ffb(i)=(z2s-z1a)/(z1b-z1a)
2066  gga(i)=one/(z1a-z1b)
2067  ggb(i)=one/(z1b-z1a)
2068  ELSEIF(k1.GT.1.AND.k1.LT.km1-1) THEN
2069  z2s=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
2070  z1a=z1(1+(i-1)*ixz1+(k1-2)*kxz1)
2071  z1b=z1(1+(i-1)*ixz1+(k1-1)*kxz1)
2072  z1c=z1(1+(i-1)*ixz1+(k1+0)*kxz1)
2073  z1d=z1(1+(i-1)*ixz1+(k1+1)*kxz1)
2074  ffa(i)=(z2s-z1b)/(z1a-z1b)* &
2075  (z2s-z1c)/(z1a-z1c)* &
2076  (z2s-z1d)/(z1a-z1d)
2077  ffb(i)=(z2s-z1a)/(z1b-z1a)* &
2078  (z2s-z1c)/(z1b-z1c)* &
2079  (z2s-z1d)/(z1b-z1d)
2080  ffc(i)=(z2s-z1a)/(z1c-z1a)* &
2081  (z2s-z1b)/(z1c-z1b)* &
2082  (z2s-z1d)/(z1c-z1d)
2083  ffd(i)=(z2s-z1a)/(z1d-z1a)* &
2084  (z2s-z1b)/(z1d-z1b)* &
2085  (z2s-z1c)/(z1d-z1c)
2086  gga(i)= one/(z1a-z1b)* &
2087  (z2s-z1c)/(z1a-z1c)* &
2088  (z2s-z1d)/(z1a-z1d)+ &
2089  (z2s-z1b)/(z1a-z1b)* &
2090  one/(z1a-z1c)* &
2091  (z2s-z1d)/(z1a-z1d)+ &
2092  (z2s-z1b)/(z1a-z1b)* &
2093  (z2s-z1c)/(z1a-z1c)* &
2094  one/(z1a-z1d)
2095  ggb(i)= one/(z1b-z1a)* &
2096  (z2s-z1c)/(z1b-z1c)* &
2097  (z2s-z1d)/(z1b-z1d)+ &
2098  (z2s-z1a)/(z1b-z1a)* &
2099  one/(z1b-z1c)* &
2100  (z2s-z1d)/(z1b-z1d)+ &
2101  (z2s-z1a)/(z1b-z1a)* &
2102  (z2s-z1c)/(z1b-z1c)* &
2103  one/(z1b-z1d)
2104  ggc(i)= one/(z1c-z1a)* &
2105  (z2s-z1b)/(z1c-z1b)* &
2106  (z2s-z1d)/(z1c-z1d)+ &
2107  (z2s-z1a)/(z1c-z1a)* &
2108  one/(z1c-z1b)* &
2109  (z2s-z1d)/(z1c-z1d)+ &
2110  (z2s-z1a)/(z1c-z1a)* &
2111  (z2s-z1b)/(z1c-z1b)* &
2112  one/(z1c-z1d)
2113  ggd(i)= one/(z1d-z1a)* &
2114  (z2s-z1b)/(z1d-z1b)* &
2115  (z2s-z1c)/(z1d-z1c)+ &
2116  (z2s-z1a)/(z1d-z1a)* &
2117  one/(z1d-z1b)* &
2118  (z2s-z1c)/(z1d-z1c)+ &
2119  (z2s-z1a)/(z1d-z1a)* &
2120  (z2s-z1b)/(z1d-z1b)* &
2121  one/(z1d-z1c)
2122  ENDIF
2123  ENDDO
2124 
2125 ! INTERPOLATE.
2126  DO n=1,nm
2127  DO i=1,im
2128  k1=k1s(i,k2)
2129  IF(k1.EQ.0) THEN
2130  q2s=q1(1+(i-1)*ixq1+(n-1)*nxq1)
2131 ! J2S=0
2132  ELSEIF(k1.EQ.km1) THEN
2133  q2s=q1(1+(i-1)*ixq1+(km1-1)*kxq1+(n-1)*nxq1)
2134 ! J2S=0
2135  ELSEIF(k1.EQ.1.OR.k1.EQ.km1-1) THEN
2136  q1a=q1(1+(i-1)*ixq1+(k1-1)*kxq1+(n-1)*nxq1)
2137  q1b=q1(1+(i-1)*ixq1+(k1+0)*kxq1+(n-1)*nxq1)
2138  q2s=ffa(i)*q1a+ffb(i)*q1b
2139 ! J2S=GGA(I)*Q1A+GGB(I)*Q1B
2140  ELSE
2141  q1a=q1(1+(i-1)*ixq1+(k1-2)*kxq1+(n-1)*nxq1)
2142  q1b=q1(1+(i-1)*ixq1+(k1-1)*kxq1+(n-1)*nxq1)
2143  q1c=q1(1+(i-1)*ixq1+(k1+0)*kxq1+(n-1)*nxq1)
2144  q1d=q1(1+(i-1)*ixq1+(k1+1)*kxq1+(n-1)*nxq1)
2145  q2s=ffa(i)*q1a+ffb(i)*q1b+ffc(i)*q1c+ffd(i)*q1d
2146 ! J2S=GGA(I)*Q1A+GGB(I)*Q1B+GGC(I)*Q1C+GGD(I)*Q1D
2147  IF(q2s.LT.min(q1b,q1c)) THEN
2148  q2s=min(q1b,q1c)
2149 ! J2S=0
2150  ELSEIF(q2s.GT.max(q1b,q1c)) THEN
2151  q2s=max(q1b,q1c)
2152 ! J2S=0
2153  ENDIF
2154  ENDIF
2155  q2(1+(i-1)*ixq2+(k2-1)*kxq2+(n-1)*nxq2)=q2s
2156 ! J2(1+(I-1)*IXQ2+(K2-1)*KXQ2+(N-1)*NXQ2)=J2S
2157  ENDDO
2158  ENDDO
2159  ENDDO
2160 !$OMP END PARALLEL DO
2161 
2162  END SUBROUTINE terp3
2163 
2220  SUBROUTINE rsearch(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,IXL2,KXL2,L2)
2221  IMPLICIT NONE
2222 
2223  INTEGER,INTENT(IN) :: im,km1,ixz1,kxz1,km2,ixz2,kxz2,ixl2,kxl2
2224  INTEGER,INTENT(OUT) :: l2(1+(im-1)*ixl2+(km2-1)*kxl2)
2225 
2226  REAL(ESMF_KIND_R8),INTENT(IN) :: z1(1+(im-1)*ixz1+(km1-1)*kxz1)
2227  REAL(ESMF_KIND_R8),INTENT(IN) :: z2(1+(im-1)*ixz2+(km2-1)*kxz2)
2228 
2229  INTEGER :: i,k2,l
2230 
2231  REAL(ESMF_KIND_R8) :: z
2232 
2233 
2234 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2235 ! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT.
2236  DO i=1,im
2237  IF (z1(1+(i-1)*ixz1).LE.z1(1+(i-1)*ixz1+(km1-1)*kxz1)) THEN
2238 ! INPUT COORDINATE IS MONOTONICALLY ASCENDING.
2239  DO k2=1,km2
2240  z=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
2241  l=0
2242  DO
2243  IF(z.LT.z1(1+(i-1)*ixz1+l*kxz1)) EXIT
2244  l=l+1
2245  IF(l.EQ.km1) EXIT
2246  ENDDO
2247  l2(1+(i-1)*ixl2+(k2-1)*kxl2)=l
2248  ENDDO
2249  ELSE
2250 ! INPUT COORDINATE IS MONOTONICALLY DESCENDING.
2251  DO k2=1,km2
2252  z=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
2253  l=0
2254  DO
2255  IF(z.GT.z1(1+(i-1)*ixz1+l*kxz1)) EXIT
2256  l=l+1
2257  IF(l.EQ.km1) EXIT
2258  ENDDO
2259  l2(1+(i-1)*ixl2+(k2-1)*kxl2)=l
2260  ENDDO
2261  ENDIF
2262  ENDDO
2263 
2264  END SUBROUTINE rsearch
2265 
2268  subroutine compute_zh
2269 
2270  implicit none
2271 
2272  integer :: i,ii, j,k, rc, clb(2), cub(2)
2273 
2274  real(esmf_kind_r8), allocatable :: pe0(:), pn0(:)
2275  real(esmf_kind_r8), pointer :: psptr(:,:)
2276  real(esmf_kind_r8), pointer :: zhsfcptr(:,:)
2277  real(esmf_kind_r8), pointer :: zhptr(:,:,:)
2278  real(esmf_kind_r8), pointer :: tptr(:,:,:)
2279  real(esmf_kind_r8), pointer :: qptr(:,:,:)
2280  real(esmf_kind_r8) :: ak, bk, zvir, grd
2281  real(esmf_kind_r8), parameter :: grav = 9.80665
2282  real(esmf_kind_r8), parameter :: rdgas = 287.05
2283  real(esmf_kind_r8), parameter :: rvgas = 461.50
2284 
2285  print*,"- COMPUTE HEIGHT"
2286 
2287  print*,"- CALL FieldGet FOR SURFACE PRESSURE"
2288  call esmf_fieldget(ps_target_grid, &
2289  computationallbound=clb, &
2290  computationalubound=cub, &
2291  farrayptr=psptr, rc=rc)
2292  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2293  call error_handler("IN FieldGet", rc)
2294 
2295  print*,"- CALL FieldGet FOR TERRAIN HEIGHT"
2296  call esmf_fieldget(terrain_target_grid, &
2297  farrayptr=zhsfcptr, rc=rc)
2298  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2299  call error_handler("IN FieldGet", rc)
2300 
2301  print*,"- CALL FieldGet FOR HEIGHT"
2302  call esmf_fieldget(zh_target_grid, &
2303  farrayptr=zhptr, rc=rc)
2304  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2305  call error_handler("IN FieldGet", rc)
2306 
2307  print*,"- CALL FieldGet FOR TEMPERATURE"
2308  call esmf_fieldget(temp_target_grid, &
2309  farrayptr=tptr, rc=rc)
2310  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2311  call error_handler("IN FieldGet", rc)
2312 
2313  do ii = 1, num_tracers
2314  if (trim(tracers(ii)) == "sphum") exit
2315  enddo
2316 
2317  print*,"- CALL FieldGet FOR SPECIFIC HUMIDITY"
2318  call esmf_fieldget(tracers_target_grid(ii), &
2319  farrayptr=qptr, rc=rc)
2320  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2321  call error_handler("IN FieldGet", rc)
2322 
2323  grd = grav/rdgas
2324  zvir = rvgas/rdgas - 1.0_esmf_kind_r8
2325 
2326  allocate(pe0(levp1_target))
2327  allocate(pn0(levp1_target))
2328 
2329  do j = clb(2), cub(2)
2330  do i = clb(1), cub(1)
2331 
2332  do k = 1, levp1_target
2333  ak = vcoord_target(k,1)
2334  ak = max(ak, 1.e-9)
2335  bk = vcoord_target(k,2)
2336 
2337  pe0(k) = ak + bk*psptr(i,j)
2338  pn0(k) = log(pe0(k))
2339  enddo
2340 
2341  zhptr(i,j,1) = zhsfcptr(i,j)
2342 
2343  do k = 2, levp1_target
2344  zhptr(i,j,k) = zhptr(i,j,k-1)+tptr(i,j,k-1)*(1.+zvir*qptr(i,j,k-1))* &
2345  (pn0(k-1)-pn0(k))/grd
2346  enddo
2347 
2348  enddo
2349  enddo
2350 
2351  deallocate(pe0, pn0)
2352 
2353  end subroutine compute_zh
2354 
2358  subroutine cleanup_target_atm_b4adj_data
2359 
2360  implicit none
2361 
2362  integer :: i, rc
2363 
2364  print*,"- DESTROY TARGET GRID ATMOSPHERIC BEFORE ADJUSTMENT FIELDS."
2365 
2366  call esmf_fielddestroy(xwind_b4adj_target_grid, rc=rc)
2367  call esmf_fielddestroy(ywind_b4adj_target_grid, rc=rc)
2368  call esmf_fielddestroy(zwind_b4adj_target_grid, rc=rc)
2369  call esmf_fielddestroy(dzdt_b4adj_target_grid, rc=rc)
2370  call esmf_fielddestroy(ps_b4adj_target_grid, rc=rc)
2371  call esmf_fielddestroy(pres_b4adj_target_grid, rc=rc)
2372  call esmf_fielddestroy(temp_b4adj_target_grid, rc=rc)
2373  call esmf_fielddestroy(terrain_interp_to_target_grid, rc=rc)
2374 
2375  do i = 1, num_tracers_input
2376  call esmf_fielddestroy(tracers_b4adj_target_grid(i), rc=rc)
2377  enddo
2378 
2379  deallocate(tracers_b4adj_target_grid)
2380 
2381  end subroutine cleanup_target_atm_b4adj_data
2382 
2385  subroutine cleanup_all_target_atm_data
2386 
2387  use atmosphere_target_data, only : cleanup_atmosphere_target_data
2388 
2389  implicit none
2390 
2391  integer :: rc
2392 
2393  print*,"- DESTROY LOCAL TARGET GRID ATMOSPHERIC FIELDS."
2394 
2395  call esmf_fielddestroy(pres_target_grid, rc=rc)
2396  call esmf_fielddestroy(xwind_target_grid, rc=rc)
2397  call esmf_fielddestroy(ywind_target_grid, rc=rc)
2398  call esmf_fielddestroy(zwind_target_grid, rc=rc)
2399  call esmf_fielddestroy(xwind_s_target_grid, rc=rc)
2400  call esmf_fielddestroy(ywind_s_target_grid, rc=rc)
2401  call esmf_fielddestroy(zwind_s_target_grid, rc=rc)
2402  call esmf_fielddestroy(xwind_w_target_grid, rc=rc)
2403  call esmf_fielddestroy(ywind_w_target_grid, rc=rc)
2404  call esmf_fielddestroy(zwind_w_target_grid, rc=rc)
2405 
2406  call cleanup_atmosphere_target_data
2407 
2408  end subroutine cleanup_all_target_atm_data
2409 
2410  end module atmosphere
atmosphere::cleanup_target_atm_b4adj_data
subroutine cleanup_target_atm_b4adj_data
Cleanup atmospheric field (before adjustment) objects.
Definition: atmosphere.F90:2358
write_data
Definition: write_data.F90:1
atmosphere::vintg
subroutine vintg
Vertically interpolate upper-air fields.
Definition: atmosphere.F90:1763
write_data::write_fv3_atm_bndy_data_netcdf
subroutine, public write_fv3_atm_bndy_data_netcdf(localpet)
Writes atmospheric fields along the lateral boundary.
Definition: write_data.F90:99
atmosphere::convert_winds_to_uv
subroutine convert_winds_to_uv
Convert 3-d component winds to u and v.
Definition: atmosphere.F90:779
atmosphere_target_data::cleanup_atmosphere_target_data
subroutine, public cleanup_atmosphere_target_data
Free up memory for fields and variables in this module.
Definition: atmosphere_target_data.F90:45
write_data::write_fv3_atm_data_netcdf
subroutine, public write_fv3_atm_data_netcdf(localpet)
Write atmospheric coldstart files (netcdf format).
Definition: write_data.F90:1197
atmosphere::read_vcoord_info
subroutine, public read_vcoord_info
Reads model vertical coordinate definition file (as specified by namelist variable vcoord_file_target...
Definition: atmosphere.F90:1256
atmosphere::vintg_thomp_mp_climo
subroutine vintg_thomp_mp_climo
Vertically interpolate atmospheric fields to target FV3 grid.
Definition: atmosphere.F90:1403
atmosphere::newps
subroutine newps(localpet)
Computes adjusted surface pressure given a new terrain height.
Definition: atmosphere.F90:1059
atmosphere::create_atm_b4adj_esmf_fields
subroutine create_atm_b4adj_esmf_fields
Create target grid field objects to hold data before vertical interpolation.
Definition: atmosphere.F90:494
model_grid
Sets up the ESMF grid objects for the input data grid and target FV3 grid.
Definition: model_grid.F90:9
atmosphere
Process atmospheric fields.
Definition: atmosphere.F90:19
thompson_mp_climo_data::cleanup_thomp_mp_climo_input_data
subroutine, public cleanup_thomp_mp_climo_input_data
Free up memory associated with this module.
Definition: thompson_mp_climo_data.F90:328
atmosphere::terp3
subroutine terp3(IM, IXZ1, IXQ1, IXZ2, IXQ2, NM, NXQ1, NXQ2, KM1, KXZ1, KXQ1, Z1, Q1, KM2, KXZ2, KXQ2, Z2, Q2)
Cubically interpolate in one dimension.
Definition: atmosphere.F90:2026
write_data::write_fv3_atm_header_netcdf
subroutine, public write_fv3_atm_header_netcdf(localpet)
Writes atmospheric header file in netcdf format.
Definition: write_data.F90:28
atmosphere::atmosphere_driver
subroutine, public atmosphere_driver(localpet)
Driver routine to process for atmospheric fields.
Definition: atmosphere.F90:118
thompson_mp_climo_data
Module to read the Thompson climatological MP data file and set up the associated esmf field and grid...
Definition: thompson_mp_climo_data.F90:9
utilities::error_handler
subroutine error_handler(string, rc)
General error handler.
Definition: utils.F90:12
atm_input_data
Read atmospheric data on the input grid.
Definition: atm_input_data.F90:15
thompson_mp_climo_data::read_thomp_mp_climo_data
subroutine, public read_thomp_mp_climo_data
Read Thompson climatological MP data file and time interpolate data to current cycle time...
Definition: thompson_mp_climo_data.F90:50
atmosphere::vintg_wam
subroutine vintg_wam(YEAR, MONTH, DAY, HOUR, PF)
Vertically extend model top into thermosphere for whole atmosphere model.
Definition: atmosphere.F90:1527
atmosphere::create_atm_esmf_fields
subroutine create_atm_esmf_fields
Create target grid field objects.
Definition: atmosphere.F90:586
atmosphere_target_data
Module to hold variables and ESMF fields associated with the target grid atmospheric data...
Definition: atmosphere_target_data.F90:9
gettemp
subroutine gettemp(iday, xlat, pr, np, pf, temp, n_o, n_o2, n_n2)
Routine that computes temperature and neutral density values utilizing MSIS 2.1.
Definition: wam_climo_data.f90:23
atmosphere::rsearch
subroutine rsearch(IM, KM1, IXZ1, KXZ1, Z1, KM2, IXZ2, KXZ2, Z2, IXL2, KXL2, L2)
Search for a surrounding real interval.
Definition: atmosphere.F90:2220
atmosphere::cleanup_all_target_atm_data
subroutine cleanup_all_target_atm_data
Cleanup target grid atmospheric field objects.
Definition: atmosphere.F90:2385
atmosphere::horiz_interp_thomp_mp_climo
subroutine horiz_interp_thomp_mp_climo
Horizontally interpolate thompson microphysics data to the target model grid.
Definition: atmosphere.F90:1291
program_setup
This module contains code to read the setup namelist file, handle the varmap file for GRIB2 data...
Definition: program_setup.F90:9
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:77
atmosphere::compute_zh
subroutine compute_zh
Compute vertical level height.
Definition: atmosphere.F90:2268
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:3282
atmosphere::newpr1
subroutine newpr1(localpet)
Computes 3-D pressure given an adjusted surface pressure.
Definition: atmosphere.F90:948