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