chgres_cube  1.3.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  regional, &
44  tracers, num_tracers, &
45  atm_weight_file, &
46  use_thomp_mp_climo
47 
48  use thompson_mp_climo_data, only : read_thomp_mp_climo_data, &
49  cleanup_thomp_mp_climo_input_data, &
50  qnifa_climo_input_grid, &
51  qnwfa_climo_input_grid, &
52  thomp_pres_climo_input_grid, &
53  lev_thomp_mp_climo
54 
55  implicit none
56 
57  private
58 
59  integer, public :: lev_target
60  integer, public :: levp1_target
61  integer, public :: nvcoord_target
62 
63  real(esmf_kind_r8), allocatable, public :: vcoord_target(:,:)
64 
65  type(esmf_field), public :: delp_target_grid
66  type(esmf_field), public :: dzdt_target_grid
67  type(esmf_field) :: dzdt_b4adj_target_grid
68  type(esmf_field), allocatable, public :: tracers_target_grid(:)
69  type(esmf_field), allocatable :: tracers_b4adj_target_grid(:)
70  type(esmf_field), public :: ps_target_grid
71  type(esmf_field) :: ps_b4adj_target_grid
72  type(esmf_field) :: pres_target_grid
73  type(esmf_field) :: pres_b4adj_target_grid
74  type(esmf_field), public :: temp_target_grid
75  type(esmf_field) :: temp_b4adj_target_grid
76  type(esmf_field) :: terrain_interp_to_target_grid
77  type(esmf_field), public :: u_s_target_grid
78  type(esmf_field), public :: v_s_target_grid
79  type(esmf_field) :: wind_target_grid
80  type(esmf_field) :: wind_b4adj_target_grid
81  type(esmf_field) :: wind_s_target_grid
82  type(esmf_field), public :: u_w_target_grid
83  type(esmf_field), public :: v_w_target_grid
84  type(esmf_field) :: wind_w_target_grid
85  type(esmf_field), public :: zh_target_grid
86 
87 ! Fields associated with thompson microphysics climatological tracers.
88 
89  type(esmf_field) :: qnifa_climo_b4adj_target_grid
91  type(esmf_field), public :: qnifa_climo_target_grid
94  type(esmf_field) :: qnwfa_climo_b4adj_target_grid
96  type(esmf_field), public :: qnwfa_climo_target_grid
99  type(esmf_field) :: thomp_pres_climo_b4adj_target_grid
101 
102  public :: atmosphere_driver
103 
104  contains
105 
110  subroutine atmosphere_driver(localpet)
111 
112  use mpi
113 
114  implicit none
115 
116  integer, intent(in) :: localpet
117 
118  integer :: isrctermprocessing
119  integer :: rc, n
120 
121  type(esmf_regridmethod_flag) :: method
122  type(esmf_routehandle) :: regrid_bl
123 
124  real(esmf_kind_r8), parameter :: p0=101325.0
125  real(esmf_kind_r8), parameter :: rd = 287.058
126  real(esmf_kind_r8), parameter :: grav = 9.81
127  real(esmf_kind_r8), parameter :: lapse = -6.5e-03
128 
129  real(esmf_kind_r8), parameter :: exponent = rd*lapse/grav
130  real(esmf_kind_r8), parameter :: one_over_exponent = 1.0 / exponent
131 
132  real(esmf_kind_r8), pointer :: psptr(:,:), tempptr(:,:,:)
133 
134 !-----------------------------------------------------------------------------------
135 ! Read atmospheric fields on the input grid.
136 !-----------------------------------------------------------------------------------
137 
138  call read_input_atm_data(localpet)
139 
140 !-----------------------------------------------------------------------------------
141 ! Read vertical coordinate info for target grid.
142 !-----------------------------------------------------------------------------------
143 
144  call read_vcoord_info
145 
146 !-----------------------------------------------------------------------------------
147 ! Create target grid field objects to hold data before vertical adjustment.
148 !-----------------------------------------------------------------------------------
149 
150  call create_atm_b4adj_esmf_fields
151 
152 !-----------------------------------------------------------------------------------
153 ! Horizontally interpolate. If specified, use weights from file.
154 !-----------------------------------------------------------------------------------
155 
156  isrctermprocessing = 1
157 
158  if (trim(atm_weight_file) /= "NULL") then
159 
160  print*,"- CALL FieldSMMStore FOR ATMOSPHERIC FIELDS."
161 
162  call esmf_fieldsmmstore(temp_input_grid, &
163  temp_b4adj_target_grid, &
164  atm_weight_file, &
165  routehandle=regrid_bl, &
166  srctermprocessing=isrctermprocessing, rc=rc)
167  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
168  call error_handler("IN FieldSMMStore", rc)
169 
170  else
171 
172  print*,"- CALL FieldRegridStore FOR ATMOSPHERIC FIELDS."
173 
174  method=esmf_regridmethod_bilinear
175 
176  call esmf_fieldregridstore(temp_input_grid, &
177  temp_b4adj_target_grid, &
178  polemethod=esmf_polemethod_allavg, &
179  srctermprocessing=isrctermprocessing, &
180  routehandle=regrid_bl, &
181  regridmethod=method, rc=rc)
182  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
183  call error_handler("IN FieldRegridStore", rc)
184 
185  endif
186 
187  print*,"- CALL Field_Regrid FOR TEMPERATURE."
188  call esmf_fieldregrid(temp_input_grid, &
189  temp_b4adj_target_grid, &
190  routehandle=regrid_bl, &
191  termorderflag=esmf_termorder_srcseq, &
192  rc=rc)
193  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
194  call error_handler("IN FieldRegrid", rc)
195 
196  print*,"- CALL Field_Regrid FOR PRESSURE."
197  call esmf_fieldregrid(pres_input_grid, &
198  pres_b4adj_target_grid, &
199  routehandle=regrid_bl, &
200  termorderflag=esmf_termorder_srcseq, &
201  rc=rc)
202  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
203  call error_handler("IN FieldRegrid", rc)
204 
205  do n = 1, num_tracers
206  print*,"- CALL Field_Regrid FOR TRACER ", trim(tracers(n))
207  call esmf_fieldregrid(tracers_input_grid(n), &
208  tracers_b4adj_target_grid(n), &
209  routehandle=regrid_bl, &
210  termorderflag=esmf_termorder_srcseq, rc=rc)
211  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
212  call error_handler("IN FieldRegrid", rc)
213 
214  enddo
215 
216  print*,"- CALL Field_Regrid FOR VERTICAL VELOCITY."
217  call esmf_fieldregrid(dzdt_input_grid, &
218  dzdt_b4adj_target_grid, &
219  routehandle=regrid_bl, &
220  termorderflag=esmf_termorder_srcseq, rc=rc)
221  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
222  call error_handler("IN FieldRegrid", rc)
223 
224  nullify(tempptr)
225  print*,"- CALL FieldGet FOR INPUT GRID VERTICAL VEL."
226  call esmf_fieldget(dzdt_input_grid, &
227  farrayptr=tempptr, rc=rc)
228  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
229  call error_handler("IN FieldGet", rc)
230 
231  print*, "MIN MAX W INPUT = ", minval(tempptr), maxval(tempptr)
232 
233  nullify(tempptr)
234  print*,"- CALL FieldGet FOR VERTICAL VEL B4ADJ."
235  call esmf_fieldget(dzdt_b4adj_target_grid, &
236  farrayptr=tempptr, rc=rc)
237  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
238  call error_handler("IN FieldGet", rc)
239 
240  print*, "MIN MAX W B4ADJ = ", minval(tempptr), maxval(tempptr)
241 
242  nullify(psptr)
243  print*,"- CALL FieldGet FOR INPUT SURFACE PRESSURE."
244  call esmf_fieldget(ps_input_grid, &
245  farrayptr=psptr, rc=rc)
246  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
247  call error_handler("IN FieldGet", rc)
248 
249 !------------------------------------------------------------------------------------
250 ! Assume standard lapse rate when interpolating pressure (per Phil Pegion).
251 !------------------------------------------------------------------------------------
252 
253  psptr = (psptr/p0)**exponent
254 
255  print*,"- CALL Field_Regrid FOR SURFACE PRESSURE."
256  call esmf_fieldregrid(ps_input_grid, &
257  ps_b4adj_target_grid, &
258  routehandle=regrid_bl, &
259  termorderflag=esmf_termorder_srcseq, &
260  rc=rc)
261  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
262  call error_handler("IN FieldRegrid", rc)
263 
264  nullify(psptr)
265  print*,"- CALL FieldGet FOR INPUT SURFACE PRESSURE B4ADJ."
266  call esmf_fieldget(ps_b4adj_target_grid, &
267  farrayptr=psptr, rc=rc)
268  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
269  call error_handler("IN FieldGet", rc)
270 
271  psptr = p0 * psptr**one_over_exponent
272 
273  print*,"- CALL Field_Regrid FOR TERRAIN."
274  call esmf_fieldregrid(terrain_input_grid, &
275  terrain_interp_to_target_grid, &
276  routehandle=regrid_bl, &
277  termorderflag=esmf_termorder_srcseq, &
278  rc=rc)
279  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
280  call error_handler("IN FieldRegrid", rc)
281 
282  print*,"- CALL Field_Regrid FOR 3-D WIND."
283  call esmf_fieldregrid(wind_input_grid, &
284  wind_b4adj_target_grid, &
285  routehandle=regrid_bl, &
286  termorderflag=esmf_termorder_srcseq, 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 FieldRegridRelease."
291  call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
292  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
293  call error_handler("IN FieldRegridRelease", rc)
294 
295 !-----------------------------------------------------------------------------------
296 ! Deallocate input fields.
297 !-----------------------------------------------------------------------------------
298 
299  call cleanup_input_atm_data
300 
301 !-----------------------------------------------------------------------------------
302 ! Create target grid field objects to hold data after vertical interpolation.
303 !-----------------------------------------------------------------------------------
304 
305  call create_atm_esmf_fields
306 
307 !-----------------------------------------------------------------------------------
308 ! Adjust surface pressure for terrain differences.
309 !-----------------------------------------------------------------------------------
310 
311  call newps(localpet)
312 
313 !-----------------------------------------------------------------------------------
314 ! Compute 3-d pressure based on adjusted surface pressure.
315 !-----------------------------------------------------------------------------------
316 
317  call newpr1(localpet)
318 
319 !-----------------------------------------------------------------------------------
320 ! Vertically interpolate.
321 !-----------------------------------------------------------------------------------
322 
323  call vintg
324 
325 !-----------------------------------------------------------------------------------
326 ! Compute height.
327 !-----------------------------------------------------------------------------------
328 
329  call compute_zh
330 
331 !-----------------------------------------------------------------------------------
332 ! Free up memory.
333 !-----------------------------------------------------------------------------------
334 
335  call cleanup_target_atm_b4adj_data
336 
337 !-----------------------------------------------------------------------------------
338 ! Interpolate winds to 'd' grid.
339 !-----------------------------------------------------------------------------------
340 
341  isrctermprocessing = 1
342  method=esmf_regridmethod_bilinear
343 
344  print*,"- CALL FieldRegridStore FOR 3D-WIND WEST EDGE."
345  call esmf_fieldregridstore(wind_target_grid, &
346  wind_w_target_grid, &
347  polemethod=esmf_polemethod_allavg, &
348  srctermprocessing=isrctermprocessing, &
349  routehandle=regrid_bl, &
350  extrapmethod=esmf_extrapmethod_nearest_stod, &
351  regridmethod=method, rc=rc)
352  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
353  call error_handler("IN FieldRegridStore", rc)
354 
355  print*,"- CALL Field_Regrid FOR 3-D WIND WEST EDGE."
356  call esmf_fieldregrid(wind_target_grid, &
357  wind_w_target_grid, &
358  routehandle=regrid_bl, &
359  termorderflag=esmf_termorder_srcseq, rc=rc)
360  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
361  call error_handler("IN FieldRegrid", rc)
362 
363  print*,"- CALL FieldRegridRelease."
364  call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
365  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
366  call error_handler("IN FieldRegridRelease", rc)
367 
368  isrctermprocessing = 1
369  method=esmf_regridmethod_bilinear
370 
371  print*,"- CALL FieldRegridStore FOR 3D-WIND SOUTH EDGE."
372  call esmf_fieldregridstore(wind_target_grid, &
373  wind_s_target_grid, &
374  polemethod=esmf_polemethod_allavg, &
375  srctermprocessing=isrctermprocessing, &
376  routehandle=regrid_bl, &
377  extrapmethod=esmf_extrapmethod_nearest_stod, &
378  regridmethod=method, rc=rc)
379  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
380  call error_handler("IN FieldRegridStore", rc)
381 
382  print*,"- CALL Field_Regrid FOR 3-D WIND SOUTH EDGE."
383  call esmf_fieldregrid(wind_target_grid, &
384  wind_s_target_grid, &
385  routehandle=regrid_bl, &
386  termorderflag=esmf_termorder_srcseq, rc=rc)
387  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
388  call error_handler("IN FieldRegrid", rc)
389 
390  print*,"- CALL FieldRegridRelease."
391  call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
392  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
393  call error_handler("IN FieldRegridRelease", rc)
394 
395 !-----------------------------------------------------------------------------------
396 ! Convert from 3-d to 2-d cartesian winds.
397 !-----------------------------------------------------------------------------------
398 
399  call convert_winds
400 
401 !-----------------------------------------------------------------------------------
402 ! If selected, process thompson microphysics climatological fields.
403 !-----------------------------------------------------------------------------------
404 
405  if (use_thomp_mp_climo) then
406  call read_thomp_mp_climo_data
407  call horiz_interp_thomp_mp_climo
408  call vintg_thomp_mp_climo
409  endif
410 
411 !-----------------------------------------------------------------------------------
412 ! Write target data to file.
413 !-----------------------------------------------------------------------------------
414 
415  call write_fv3_atm_header_netcdf(localpet)
416  if (regional <= 1) call write_fv3_atm_data_netcdf(localpet)
417  if (regional >= 1) call write_fv3_atm_bndy_data_netcdf(localpet)
418 
419 !-----------------------------------------------------------------------------------
420 ! Free up memory.
421 !-----------------------------------------------------------------------------------
422 
423  call cleanup_target_atm_data
424 
425  end subroutine atmosphere_driver
426 
432  subroutine create_atm_b4adj_esmf_fields
433 
434  implicit none
435 
436  integer :: rc, n
437 
438  allocate(tracers_b4adj_target_grid(num_tracers))
439 
440  do n = 1, num_tracers
441  print*,"- CALL FieldCreate FOR TARGET GRID TRACER BEFORE ADJUSTMENT ", trim(tracers(n))
442  tracers_b4adj_target_grid(n) = esmf_fieldcreate(target_grid, &
443  typekind=esmf_typekind_r8, &
444  staggerloc=esmf_staggerloc_center, &
445  ungriddedlbound=(/1/), &
446  ungriddedubound=(/lev_input/), rc=rc)
447  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
448  call error_handler("IN FieldCreate", rc)
449  enddo
450 
451  print*,"- CALL FieldCreate FOR TARGET GRID TEMPERATURE BEFORE ADJUSTMENT."
452  temp_b4adj_target_grid = esmf_fieldcreate(target_grid, &
453  typekind=esmf_typekind_r8, &
454  staggerloc=esmf_staggerloc_center, &
455  ungriddedlbound=(/1/), &
456  ungriddedubound=(/lev_input/), rc=rc)
457  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
458  call error_handler("IN FieldCreate", rc)
459 
460  print*,"- CALL FieldCreate FOR TARGET GRID PRESSURE BEFORE ADJUSTMENT."
461  pres_b4adj_target_grid = esmf_fieldcreate(target_grid, &
462  typekind=esmf_typekind_r8, &
463  staggerloc=esmf_staggerloc_center, &
464  ungriddedlbound=(/1/), &
465  ungriddedubound=(/lev_input/), rc=rc)
466  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
467  call error_handler("IN FieldCreate", rc)
468 
469  print*,"- CALL FieldCreate FOR TARGET GRID VERTICAL VELOCITY BEFORE ADJUSTMENT."
470  dzdt_b4adj_target_grid = esmf_fieldcreate(target_grid, &
471  typekind=esmf_typekind_r8, &
472  staggerloc=esmf_staggerloc_center, &
473  ungriddedlbound=(/1/), &
474  ungriddedubound=(/lev_input/), rc=rc)
475  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
476  call error_handler("IN FieldCreate", rc)
477 
478  print*,"- CALL FieldCreate FOR TARGET GRID UNSTAGGERED WINDS BEFORE ADJUSTMENT."
479  wind_b4adj_target_grid = esmf_fieldcreate(target_grid, &
480  typekind=esmf_typekind_r8, &
481  staggerloc=esmf_staggerloc_center, &
482  ungriddedlbound=(/1,1/), &
483  ungriddedubound=(/lev_input,3/), rc=rc)
484  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
485  call error_handler("IN FieldCreate", rc)
486 
487  print*,"- CALL FieldCreate FOR TARGET TERRAIN."
488  terrain_interp_to_target_grid = esmf_fieldcreate(target_grid, &
489  typekind=esmf_typekind_r8, &
490  staggerloc=esmf_staggerloc_center, rc=rc)
491  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
492  call error_handler("IN FieldCreate", rc)
493 
494  print*,"- CALL FieldCreate FOR TARGET SURFACE PRESSURE BEFORE ADJUSTMENT."
495  ps_b4adj_target_grid = esmf_fieldcreate(target_grid, &
496  typekind=esmf_typekind_r8, &
497  staggerloc=esmf_staggerloc_center, rc=rc)
498  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
499  call error_handler("IN FieldCreate", rc)
500 
501  end subroutine create_atm_b4adj_esmf_fields
502 
506  subroutine create_atm_esmf_fields
507 
508  implicit none
509 
510  integer :: rc, n
511 
512  allocate(tracers_target_grid(num_tracers))
513 
514  do n = 1, num_tracers
515  print*,"- CALL FieldCreate FOR TARGET GRID TRACERS ", trim(tracers(n))
516  tracers_target_grid(n) = esmf_fieldcreate(target_grid, &
517  typekind=esmf_typekind_r8, &
518  staggerloc=esmf_staggerloc_center, &
519  ungriddedlbound=(/1/), &
520  ungriddedubound=(/lev_target/), rc=rc)
521  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
522  call error_handler("IN FieldCreate", rc)
523  enddo
524 
525  print*,"- CALL FieldCreate FOR TARGET GRID TEMPERATURE."
526  temp_target_grid = esmf_fieldcreate(target_grid, &
527  typekind=esmf_typekind_r8, &
528  staggerloc=esmf_staggerloc_center, &
529  ungriddedlbound=(/1/), &
530  ungriddedubound=(/lev_target/), rc=rc)
531  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
532  call error_handler("IN FieldCreate", rc)
533 
534  print*,"- CALL FieldCreate FOR TARGET GRID PRESSURE."
535  pres_target_grid = esmf_fieldcreate(target_grid, &
536  typekind=esmf_typekind_r8, &
537  staggerloc=esmf_staggerloc_center, &
538  ungriddedlbound=(/1/), &
539  ungriddedubound=(/lev_target/), rc=rc)
540  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
541  call error_handler("IN FieldCreate", rc)
542 
543  print*,"- CALL FieldCreate FOR TARGET GRID VERTICAL VELOCITY."
544  dzdt_target_grid = esmf_fieldcreate(target_grid, &
545  typekind=esmf_typekind_r8, &
546  staggerloc=esmf_staggerloc_center, &
547  ungriddedlbound=(/1/), &
548  ungriddedubound=(/lev_target/), rc=rc)
549  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
550  call error_handler("IN FieldCreate", rc)
551 
552  print*,"- CALL FieldCreate FOR TARGET GRID DELP."
553  delp_target_grid = esmf_fieldcreate(target_grid, &
554  typekind=esmf_typekind_r8, &
555  staggerloc=esmf_staggerloc_center, &
556  ungriddedlbound=(/1/), &
557  ungriddedubound=(/lev_target/), rc=rc)
558  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
559  call error_handler("IN FieldCreate", rc)
560 
561  print*,"- CALL FieldCreate FOR TARGET HEIGHT."
562  zh_target_grid = esmf_fieldcreate(target_grid, &
563  typekind=esmf_typekind_r8, &
564  staggerloc=esmf_staggerloc_center, &
565  ungriddedlbound=(/1/), &
566  ungriddedubound=(/levp1_target/), rc=rc)
567  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
568  call error_handler("IN FieldCreate", rc)
569 
570  print*,"- CALL FieldCreate FOR TARGET UNSTAGGERED 3D-WIND."
571  wind_target_grid = esmf_fieldcreate(target_grid, &
572  typekind=esmf_typekind_r8, &
573  staggerloc=esmf_staggerloc_center, &
574  ungriddedlbound=(/1,1/), &
575  ungriddedubound=(/lev_target,3/), rc=rc)
576  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
577  call error_handler("IN FieldCreate", rc)
578 
579  print*,"- CALL FieldCreate FOR TARGET U_S."
580  u_s_target_grid = esmf_fieldcreate(target_grid, &
581  typekind=esmf_typekind_r8, &
582  staggerloc=esmf_staggerloc_edge2, &
583  ungriddedlbound=(/1/), &
584  ungriddedubound=(/lev_target/), rc=rc)
585  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
586  call error_handler("IN FieldCreate", rc)
587 
588  print*,"- CALL FieldCreate FOR TARGET V_S."
589  v_s_target_grid = esmf_fieldcreate(target_grid, &
590  typekind=esmf_typekind_r8, &
591  staggerloc=esmf_staggerloc_edge2, &
592  ungriddedlbound=(/1/), &
593  ungriddedubound=(/lev_target/), rc=rc)
594  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
595  call error_handler("IN FieldCreate", rc)
596 
597  print*,"- CALL FieldCreate FOR TARGET 3D-WIND_S."
598  wind_s_target_grid = esmf_fieldcreate(target_grid, &
599  typekind=esmf_typekind_r8, &
600  staggerloc=esmf_staggerloc_edge2, &
601  ungriddedlbound=(/1,1/), &
602  ungriddedubound=(/lev_target,3/), rc=rc)
603  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
604  call error_handler("IN FieldCreate", rc)
605 
606  print*,"- CALL FieldCreate FOR TARGET U_W."
607  u_w_target_grid = esmf_fieldcreate(target_grid, &
608  typekind=esmf_typekind_r8, &
609  staggerloc=esmf_staggerloc_edge1, &
610  ungriddedlbound=(/1/), &
611  ungriddedubound=(/lev_target/), rc=rc)
612  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
613  call error_handler("IN FieldCreate", rc)
614 
615  print*,"- CALL FieldCreate FOR TARGET V_W."
616  v_w_target_grid = esmf_fieldcreate(target_grid, &
617  typekind=esmf_typekind_r8, &
618  staggerloc=esmf_staggerloc_edge1, &
619  ungriddedlbound=(/1/), &
620  ungriddedubound=(/lev_target/), rc=rc)
621  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
622  call error_handler("IN FieldCreate", rc)
623 
624  print*,"- CALL FieldCreate FOR TARGET 3D-WIND_W."
625  wind_w_target_grid = esmf_fieldcreate(target_grid, &
626  typekind=esmf_typekind_r8, &
627  staggerloc=esmf_staggerloc_edge1, &
628  ungriddedlbound=(/1,1/), &
629  ungriddedubound=(/lev_target,3/), rc=rc)
630  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
631  call error_handler("IN FieldCreate", rc)
632 
633  print*,"- CALL FieldCreate FOR TARGET SURFACE PRESSURE."
634  ps_target_grid = esmf_fieldcreate(target_grid, &
635  typekind=esmf_typekind_r8, &
636  staggerloc=esmf_staggerloc_center, rc=rc)
637  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
638  call error_handler("IN FieldCreate", rc)
639 
640  end subroutine create_atm_esmf_fields
641 
645  subroutine convert_winds
646 
647  implicit none
648 
649  integer :: clb(4), cub(4)
650  integer :: i, j, k, rc
651 
652  real(esmf_kind_r8), pointer :: latptr(:,:)
653  real(esmf_kind_r8), pointer :: lonptr(:,:)
654  real(esmf_kind_r8), pointer :: uptr(:,:,:)
655  real(esmf_kind_r8), pointer :: vptr(:,:,:)
656  real(esmf_kind_r8), pointer :: windptr(:,:,:,:)
657  real(esmf_kind_r8) :: latrad, lonrad
658 
659 !-----------------------------------------------------------------------------------
660 ! Convert from 3-d cartesian to 2-cartesian winds
661 !-----------------------------------------------------------------------------------
662 
663  print*,'- CONVERT WINDS.'
664 
665  print*,"- CALL FieldGet FOR 3-D WIND_S."
666  call esmf_fieldget(wind_s_target_grid, &
667  computationallbound=clb, &
668  computationalubound=cub, &
669  farrayptr=windptr, rc=rc)
670  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
671  call error_handler("IN FieldGet", rc)
672 
673  print*,"- CALL FieldGet FOR U_S."
674  call esmf_fieldget(u_s_target_grid, &
675  farrayptr=uptr, rc=rc)
676  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
677  call error_handler("IN FieldGet", rc)
678 
679  print*,"- CALL FieldGet FOR V_S."
680  call esmf_fieldget(v_s_target_grid, &
681  farrayptr=vptr, rc=rc)
682  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
683  call error_handler("IN FieldGet", rc)
684 
685  print*,"- CALL FieldGet FOR LATITUDE_S."
686  call esmf_fieldget(latitude_s_target_grid, &
687  farrayptr=latptr, rc=rc)
688  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
689  call error_handler("IN FieldGet", rc)
690 
691  print*,"- CALL FieldGet FOR LONGITUDE_S."
692  call esmf_fieldget(longitude_s_target_grid, &
693  farrayptr=lonptr, rc=rc)
694  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
695  call error_handler("IN FieldGet", rc)
696 
697  do i = clb(1), cub(1)
698  do j = clb(2), cub(2)
699  latrad = latptr(i,j) * acos(-1.) / 180.0
700  lonrad = lonptr(i,j) * acos(-1.) / 180.0
701  do k = clb(3), cub(3)
702  uptr(i,j,k) = windptr(i,j,k,1) * cos(lonrad) + windptr(i,j,k,2) * sin(lonrad)
703  vptr(i,j,k) = -windptr(i,j,k,1) * sin(latrad) * sin(lonrad) + &
704  windptr(i,j,k,2) * sin(latrad) * cos(lonrad) + &
705  windptr(i,j,k,3) * cos(latrad)
706  enddo
707  enddo
708  enddo
709 
710 
711  print*,"- CALL FieldGet FOR 3-D WIND_W."
712  call esmf_fieldget(wind_w_target_grid, &
713  computationallbound=clb, &
714  computationalubound=cub, &
715  farrayptr=windptr, rc=rc)
716  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
717  call error_handler("IN FieldGet", rc)
718 
719  print*,"- CALL FieldGet FOR U_W."
720  call esmf_fieldget(u_w_target_grid, &
721  farrayptr=uptr, rc=rc)
722  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
723  call error_handler("IN FieldGet", rc)
724 
725  print*,"- CALL FieldGet FOR V_W."
726  call esmf_fieldget(v_w_target_grid, &
727  farrayptr=vptr, rc=rc)
728  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
729  call error_handler("IN FieldGet", rc)
730 
731  print*,"- CALL FieldGet FOR LATITUDE_W."
732  call esmf_fieldget(latitude_w_target_grid, &
733  farrayptr=latptr, rc=rc)
734  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
735  call error_handler("IN FieldGet", rc)
736 
737  print*,"- CALL FieldGet FOR LONGITUDE_W."
738  call esmf_fieldget(longitude_w_target_grid, &
739  farrayptr=lonptr, rc=rc)
740  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
741  call error_handler("IN FieldGet", rc)
742 
743  do i = clb(1), cub(1)
744  do j = clb(2), cub(2)
745  latrad = latptr(i,j) * acos(-1.) / 180.0
746  lonrad = lonptr(i,j) * acos(-1.) / 180.0
747  do k = clb(3), cub(3)
748  uptr(i,j,k) = windptr(i,j,k,1) * cos(lonrad) + windptr(i,j,k,2) * sin(lonrad)
749  vptr(i,j,k) = -windptr(i,j,k,1) * sin(latrad) * sin(lonrad) + &
750  windptr(i,j,k,2) * sin(latrad) * cos(lonrad) + &
751  windptr(i,j,k,3) * cos(latrad)
752  enddo
753  enddo
754  enddo
755 
756  end subroutine convert_winds
757 
789  subroutine newpr1(localpet)
790  implicit none
791 
792  integer, intent(in) :: localpet
793 
794  integer :: idsl, idvc, rc
795  integer :: i, j, k, clb(3), cub(3)
796 
797  real(esmf_kind_r8), parameter :: rd=287.05
798  real(esmf_kind_r8), parameter :: cp=1004.6
799  real(esmf_kind_r8), parameter :: rocp=rd/cp
800  real(esmf_kind_r8), parameter :: rocp1=rocp+1
801  real(esmf_kind_r8), parameter :: rocpr=1/rocp
802 
803  real(esmf_kind_r8), pointer :: delp_ptr(:,:,:)
804  real(esmf_kind_r8), pointer :: pptr(:,:,:) ! adjusted 3-d p.
805  real(esmf_kind_r8), pointer :: psptr(:,:) ! adjusted surface p.
806  real(esmf_kind_r8) :: ak, bk
807  real(esmf_kind_r8), allocatable :: pi(:,:,:)
808 
809  print*,"COMPUTE 3-D PRESSURE FROM ADJUSTED SURFACE PRESSURE."
810 
811  idvc = 2 ! hard wire for now.
812  idsl = 2 ! hard wire for now.
813 
814  print*,"- CALL FieldGet FOR 3-D PRES."
815  call esmf_fieldget(pres_target_grid, &
816  computationallbound=clb, &
817  computationalubound=cub, &
818  farrayptr=pptr, rc=rc)
819  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
820  call error_handler("IN FieldGet", rc)
821 
822  print*,"- CALL FieldGet FOR DELP."
823  call esmf_fieldget(delp_target_grid, &
824  computationallbound=clb, &
825  computationalubound=cub, &
826  farrayptr=delp_ptr, 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 SURFACE PRESSURE AFTER ADJUSTMENT"
831  call esmf_fieldget(ps_target_grid, &
832  farrayptr=psptr, rc=rc)
833  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
834  call error_handler("IN FieldGet", rc)
835 
836  allocate(pi(clb(1):cub(1),clb(2):cub(2),1:levp1_target))
837 
838  if(idvc.eq.2) then
839  do k=1,levp1_target
840  ak = vcoord_target(k,1)
841  bk = vcoord_target(k,2)
842  do i= clb(1), cub(1)
843  do j= clb(2), cub(2)
844  pi(i,j,k) = ak + bk*psptr(i,j)
845  enddo
846  enddo
847  enddo
848  do k=1,lev_target
849  do i= clb(1), cub(1)
850  do j= clb(2), cub(2)
851  delp_ptr(i,j,k) = pi(i,j,k) - pi(i,j,k+1)
852  enddo
853  enddo
854  enddo
855  else
856  call error_handler("PROGRAM ONLY WORKS WITH IDVC 2", 1)
857  endif
858 
859  if(idsl.eq.2) then
860  do k=1,lev_target
861  do i= clb(1), cub(1)
862  do j= clb(2), cub(2)
863  pptr(i,j,k) = (pi(i,j,k)+pi(i,j,k+1))/2.0
864  enddo
865  enddo
866  enddo
867  else
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)**rocp1-pi(i,j,k+1)**rocp1)/ &
872  (rocp1*(pi(i,j,k)-pi(i,j,k+1))))**rocpr
873  enddo
874  enddo
875  enddo
876  endif
877 
878  deallocate(pi)
879 
880  if (localpet == 0) then
881  print*,'new pres ',pptr(clb(1),clb(2),:)
882  print*,'delp ',delp_ptr(clb(1),clb(2),:)
883  endif
884 
885  end subroutine newpr1
886 
900  subroutine newps(localpet)
901 
902  implicit none
903 
904  integer, intent(in) :: localpet
905  integer :: i, j, k, ii
906  integer :: clb(3), cub(3), ls, rc
907 
908  real(esmf_kind_r8), pointer :: pptr(:,:,:)
909  real(esmf_kind_r8), pointer :: psptr(:,:)
910  real(esmf_kind_r8), pointer :: psnewptr(:,:) ! adjusted surface p.
911  real(esmf_kind_r8), pointer :: tptr(:,:,:)
912  real(esmf_kind_r8), pointer :: qptr(:,:,:)
913  real(esmf_kind_r8), pointer :: zsptr(:,:)
914  real(esmf_kind_r8), pointer :: zsnewptr(:,:)
915  real(esmf_kind_r8), allocatable :: zu(:,:)
916  real(esmf_kind_r8), parameter :: beta=-6.5e-3
917  real(esmf_kind_r8), parameter :: epsilon=1.e-9
918  real(esmf_kind_r8), parameter :: g=9.80665
919  real(esmf_kind_r8), parameter :: rd=287.05
920  real(esmf_kind_r8), parameter :: rv=461.50
921  real(esmf_kind_r8), parameter :: gor=g/rd
922  real(esmf_kind_r8), parameter :: fv=rv/rd-1.
923  real(esmf_kind_r8) :: ftv, fgam, apu, fz0
924  real(esmf_kind_r8) :: atvu, atv, fz1, fp0
925  real(esmf_kind_r8) :: apd, azd, agam, azu
926  real(esmf_kind_r8) :: atvd, fp1, gamma, pu
927  real(esmf_kind_r8) :: tvu, pd, tvd
928  real(esmf_kind_r8) :: at, aq, ap, az
929 
930  ftv(at,aq)=at*(1+fv*aq)
931  fgam(apu,atvu,apd,atvd)=-gor*log(atvd/atvu)/log(apd/apu)
932  fz0(ap,atv,azd,apd)=azd+atv/gor*log(apd/ap)
933  fz1(ap,atv,azd,apd,agam)=azd-atv/agam*((apd/ap)**(-agam/gor)-1)
934  fp0(az,azu,apu,atvu)=apu*exp(-gor/atvu*(az-azu))
935  fp1(az,azu,apu,atvu,agam)=apu*(1+agam/atvu*(az-azu))**(-gor/agam)
936 
937  print*,"- ADJUST SURFACE PRESSURE FOR NEW TERRAIN."
938 
939  print*,"- CALL FieldGet FOR 3-D PRES."
940  call esmf_fieldget(pres_b4adj_target_grid, &
941  computationallbound=clb, &
942  computationalubound=cub, &
943  farrayptr=pptr, rc=rc)
944  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
945  call error_handler("IN FieldGet", rc)
946 
947  if(localpet==0) then
948  print*,'old pres ',pptr(clb(1),clb(2),:)
949  endif
950 
951  print*,"- CALL FieldGet FOR TEMPERATURE"
952  call esmf_fieldget(temp_b4adj_target_grid, &
953  farrayptr=tptr, rc=rc)
954  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
955  call error_handler("IN FieldGet", rc)
956 
957 ! Find specific humidity in the array of tracer fields.
958 
959  do ii = 1, num_tracers
960  if (trim(tracers(ii)) == "sphum") exit
961  enddo
962 
963  print*,"- CALL FieldGet FOR SPECIFIC HUMIDITY"
964  call esmf_fieldget(tracers_b4adj_target_grid(ii), &
965  farrayptr=qptr, rc=rc)
966  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
967  call error_handler("IN FieldGet", rc)
968 
969  print*,"- CALL FieldGet FOR SURFACE PRESSURE BEFORE ADJUSTMENT"
970  call esmf_fieldget(ps_b4adj_target_grid, &
971  farrayptr=psptr, rc=rc)
972  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
973  call error_handler("IN FieldGet", rc)
974 
975  print*,"- CALL FieldGet FOR SURFACE PRESSURE AFTER ADJUSTMENT"
976  call esmf_fieldget(ps_target_grid, &
977  farrayptr=psnewptr, 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 OLD TERRAIN"
982  call esmf_fieldget(terrain_interp_to_target_grid, &
983  farrayptr=zsptr, rc=rc)
984  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
985  call error_handler("IN FieldGet", rc)
986 
987  print*,"- CALL FieldGet FOR NEW TERRAIN"
988  call esmf_fieldget(terrain_target_grid, &
989  farrayptr=zsnewptr, rc=rc)
990  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
991  call error_handler("IN FieldGet", rc)
992 
993  allocate(zu(clb(1):cub(1),clb(2):cub(2)))
994 
995 !-----------------------------------------------------------------------------------
996 ! Note, this routine was adapted from the spectral GFS which labeled the lowest
997 ! model layer as '1'.
998 !-----------------------------------------------------------------------------------
999 
1000 !-----------------------------------------------------------------------------------
1001 ! Compute surface pressure below the original ground.
1002 !-----------------------------------------------------------------------------------
1003 
1004  ls=0
1005  k=1
1006  gamma=beta
1007  do i=clb(1), cub(1)
1008  do j=clb(2), cub(2)
1009  pu=pptr(i,j,k)
1010  tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1011  zu(i,j)=fz1(pu,tvu,zsptr(i,j),psptr(i,j),gamma)
1012  if(zsnewptr(i,j).le.zu(i,j)) then
1013  pu=pptr(i,j,k)
1014  tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1015  if(abs(gamma).gt.epsilon) then
1016  psnewptr(i,j)=fp1(zsnewptr(i,j),zu(i,j),pu,tvu,gamma)
1017  else
1018  psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1019  endif
1020  else
1021  psnewptr(i,j)=0
1022  ls=ls+1
1023  endif
1024  enddo
1025  enddo
1026 
1027 !-----------------------------------------------------------------------------------
1028 ! Compute surface pressure above the original ground.
1029 !-----------------------------------------------------------------------------------
1030 
1031  do k=2,cub(3)
1032  if(ls.gt.0) then
1033  do i=clb(1),cub(1)
1034  do j=clb(2),cub(2)
1035  if(psnewptr(i,j).eq.0) then
1036  pu=pptr(i,j,k)
1037  tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1038  pd=pptr(i,j,k-1)
1039  tvd=ftv(tptr(i,j,k-1),qptr(i,j,k-1))
1040  gamma=fgam(pu,tvu,pd,tvd)
1041  if(abs(gamma).gt.epsilon) then
1042  zu(i,j)=fz1(pu,tvu,zu(i,j),pd,gamma)
1043  else
1044  zu(i,j)=fz0(pu,tvu,zu(i,j),pd)
1045  endif
1046  if(zsnewptr(i,j).le.zu(i,j)) then
1047  if(abs(gamma).gt.epsilon) then
1048  psnewptr(i,j)=fp1(zsnewptr(i,j),zu(i,j),pu,tvu,gamma)
1049  else
1050  psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1051  endif
1052  ls=ls-1
1053  endif
1054  endif
1055  enddo
1056  enddo
1057  endif
1058  enddo
1059 
1060 !-----------------------------------------------------------------------------------
1061 ! Compute surface pressure over the top.
1062 !-----------------------------------------------------------------------------------
1063 
1064 
1065  if(ls.gt.0) then
1066  k=cub(3)
1067  gamma=0
1068  do i=clb(1),cub(1)
1069  do j=clb(2),cub(2)
1070  if(psnewptr(i,j).eq.0) then
1071  pu=pptr(i,j,k)
1072  tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1073  psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1074  endif
1075  enddo
1076  enddo
1077  endif
1078 
1079  deallocate(zu)
1080 
1081  if (localpet == 0) then
1082 ! do i=clb(1),cub(1)
1083 ! do j=clb(2),cub(2)
1084  do i=clb(1),clb(1)
1085  do j=clb(2),clb(2)
1086  print*,'sfcp adjust ',(zsnewptr(i,j)-zsptr(i,j)), psptr(i,j),psnewptr(i,j)
1087  enddo
1088  enddo
1089  endif
1090 
1091  end subroutine newps
1092 
1097  subroutine read_vcoord_info
1098  implicit none
1099 
1100  integer :: istat, n, k
1101 
1102  print*
1103  print*,"OPEN VERTICAL COORD FILE: ", trim(vcoord_file_target_grid)
1104  open(14, file=trim(vcoord_file_target_grid), form='formatted', iostat=istat)
1105  if (istat /= 0) then
1106  call error_handler("OPENING VERTICAL COORD FILE", istat)
1107  endif
1108 
1109  read(14, *, iostat=istat) nvcoord_target, lev_target
1110  if (istat /= 0) then
1111  call error_handler("READING VERTICAL COORD FILE", istat)
1112  endif
1113 
1114  levp1_target = lev_target + 1
1115 
1116  allocate(vcoord_target(levp1_target, nvcoord_target))
1117  read(14, *, iostat=istat) ((vcoord_target(n,k), k=1,nvcoord_target), n=1,levp1_target)
1118  if (istat /= 0) then
1119  call error_handler("READING VERTICAL COORD FILE", istat)
1120  endif
1121 
1122  print*
1123 
1124  close(14)
1125 
1126  end subroutine read_vcoord_info
1127 
1132  subroutine horiz_interp_thomp_mp_climo
1133 
1134  implicit none
1135 
1136  integer :: isrctermprocessing, rc
1137 
1138  type(esmf_regridmethod_flag) :: method
1139  type(esmf_routehandle) :: regrid_bl
1140 
1141  isrctermprocessing=1
1142 
1143  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNIFA BEFORE ADJUSTMENT."
1144  qnifa_climo_b4adj_target_grid = esmf_fieldcreate(target_grid, &
1145  typekind=esmf_typekind_r8, &
1146  staggerloc=esmf_staggerloc_center, &
1147  ungriddedlbound=(/1/), &
1148  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1149  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1150  call error_handler("IN FieldCreate", rc)
1151 
1152  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNWFA BEFORE ADJUSTMENT."
1153  qnwfa_climo_b4adj_target_grid = esmf_fieldcreate(target_grid, &
1154  typekind=esmf_typekind_r8, &
1155  staggerloc=esmf_staggerloc_center, &
1156  ungriddedlbound=(/1/), &
1157  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1158  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1159  call error_handler("IN FieldCreate", rc)
1160 
1161  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO PRESSURE BEFORE ADJUSTMENT."
1162  thomp_pres_climo_b4adj_target_grid = esmf_fieldcreate(target_grid, &
1163  typekind=esmf_typekind_r8, &
1164  staggerloc=esmf_staggerloc_center, &
1165  ungriddedlbound=(/1/), &
1166  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1167  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1168  call error_handler("IN FieldCreate", rc)
1169 
1170  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNIFA."
1171  qnifa_climo_target_grid = esmf_fieldcreate(target_grid, &
1172  typekind=esmf_typekind_r8, &
1173  staggerloc=esmf_staggerloc_center, &
1174  ungriddedlbound=(/1/), &
1175  ungriddedubound=(/lev_target/), rc=rc)
1176  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1177  call error_handler("IN FieldCreate", rc)
1178 
1179  print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNWFA."
1180  qnwfa_climo_target_grid = esmf_fieldcreate(target_grid, &
1181  typekind=esmf_typekind_r8, &
1182  staggerloc=esmf_staggerloc_center, &
1183  ungriddedlbound=(/1/), &
1184  ungriddedubound=(/lev_target/), rc=rc)
1185  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1186  call error_handler("IN FieldCreate", rc)
1187 
1188  print*,"- CALL FieldRegridStore FOR THOMPSON CLIMO FIELDS."
1189 
1190  method=esmf_regridmethod_bilinear
1191 
1192  call esmf_fieldregridstore(qnifa_climo_input_grid, &
1193  qnifa_climo_b4adj_target_grid, &
1194  polemethod=esmf_polemethod_allavg, &
1195  srctermprocessing=isrctermprocessing, &
1196  routehandle=regrid_bl, &
1197  regridmethod=method, rc=rc)
1198  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1199  call error_handler("IN FieldRegridStore", rc)
1200 
1201  print*,"- CALL Field_Regrid FOR THOMP CLIMO QNIFA."
1202  call esmf_fieldregrid(qnifa_climo_input_grid, &
1203  qnifa_climo_b4adj_target_grid, &
1204  routehandle=regrid_bl, &
1205  termorderflag=esmf_termorder_srcseq, rc=rc)
1206  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1207  call error_handler("IN FieldRegrid", rc)
1208 
1209  print*,"- CALL Field_Regrid FOR THOMP CLIMO QNWFA."
1210  call esmf_fieldregrid(qnwfa_climo_input_grid, &
1211  qnwfa_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 PRESSURE."
1218  call esmf_fieldregrid(thomp_pres_climo_input_grid, &
1219  thomp_pres_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 FieldRegridRelease."
1226  call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
1227  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1228  call error_handler("IN FieldRegridRelease", rc)
1229 
1230 !-----------------------------------------------------------------------------------
1231 ! Free up input data memory.
1232 !-----------------------------------------------------------------------------------
1233 
1234  call cleanup_thomp_mp_climo_input_data
1235 
1236  end subroutine horiz_interp_thomp_mp_climo
1237 
1244  SUBROUTINE vintg_thomp_mp_climo
1245 
1246  implicit none
1247 
1248  INTEGER :: clb(3), cub(3), rc
1249  INTEGER :: im, km1, km2, nt
1250  INTEGER :: i, j, k
1251 
1252  REAL(ESMF_KIND_R8), ALLOCATABLE :: z1(:,:,:), z2(:,:,:)
1253  REAL(ESMF_KIND_R8), ALLOCATABLE :: c1(:,:,:,:),c2(:,:,:,:)
1254 
1255  REAL(ESMF_KIND_R8), POINTER :: qnifa1ptr(:,:,:) ! input
1256  REAL(ESMF_KIND_R8), POINTER :: qnifa2ptr(:,:,:) ! target
1257  REAL(ESMF_KIND_R8), POINTER :: qnwfa1ptr(:,:,:) ! input
1258  REAL(ESMF_KIND_R8), POINTER :: qnwfa2ptr(:,:,:) ! target
1259  REAL(ESMF_KIND_R8), POINTER :: p1ptr(:,:,:) ! input pressure
1260  REAL(ESMF_KIND_R8), POINTER :: p2ptr(:,:,:) ! target pressure
1261 
1262  print*,"- VERTICALY INTERPOLATE THOMP MP CLIMO TRACERS."
1263 
1264  print*,"- CALL FieldGet FOR 3-D THOMP PRES."
1265  call esmf_fieldget(thomp_pres_climo_b4adj_target_grid, &
1266  computationallbound=clb, &
1267  computationalubound=cub, &
1268  farrayptr=p1ptr, rc=rc)
1269  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1270  call error_handler("IN FieldGet", rc)
1271 
1272 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1273 ! The '1'/'2' arrays hold fields before/after interpolation.
1274 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1275 
1276  nt= 2 ! number of thomp tracers
1277 
1278  ALLOCATE(z1(clb(1):cub(1),clb(2):cub(2),lev_thomp_mp_climo))
1279  ALLOCATE(z2(clb(1):cub(1),clb(2):cub(2),lev_target))
1280  ALLOCATE(c1(clb(1):cub(1),clb(2):cub(2),lev_thomp_mp_climo,nt))
1281  ALLOCATE(c2(clb(1):cub(1),clb(2):cub(2),lev_target,nt))
1282 
1283  z1 = -log(p1ptr)
1284 
1285  print*,"- CALL FieldGet FOR 3-D ADJUSTED PRESS"
1286  call esmf_fieldget(pres_target_grid, &
1287  farrayptr=p2ptr, rc=rc)
1288  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1289  call error_handler("IN FieldGet", rc)
1290 
1291  z2 = -log(p2ptr)
1292 
1293 !print*,'pres check 1 ', p1ptr(clb(1),clb(2),:)
1294 !print*,'pres check 2 ', p2ptr(clb(1),clb(2),:)
1295 
1296  print*,"- CALL FieldGet FOR qnifa before vertical adjustment."
1297  call esmf_fieldget(qnifa_climo_b4adj_target_grid, &
1298  farrayptr=qnifa1ptr, rc=rc)
1299  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1300  call error_handler("IN FieldGet", rc)
1301 
1302  c1(:,:,:,1) = qnifa1ptr(:,:,:)
1303 
1304  print*,"- CALL FieldGet FOR qnwfa before vertical adjustment."
1305  call esmf_fieldget(qnwfa_climo_b4adj_target_grid, &
1306  farrayptr=qnwfa1ptr, 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(:,:,:,2) = qnwfa1ptr(:,:,:)
1311 
1312 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1313 ! PERFORM LAGRANGIAN ONE-DIMENSIONAL INTERPOLATION
1314 ! THAT IS 4TH-ORDER IN INTERIOR, 2ND-ORDER IN OUTSIDE INTERVALS
1315 ! AND 1ST-ORDER FOR EXTRAPOLATION.
1316 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1317 
1318  im = (cub(1)-clb(1)+1) * (cub(2)-clb(2)+1)
1319  km1= lev_thomp_mp_climo
1320  km2= lev_target
1321 
1322  CALL terp3(im,1,1,1,1,nt,(im*km1),(im*km2), &
1323  km1,im,im,z1,c1,km2,im,im,z2,c2)
1324 
1325  print*,"- CALL FieldGet FOR ADJUSTED climo qnifa."
1326  call esmf_fieldget(qnifa_climo_target_grid, &
1327  farrayptr=qnifa2ptr, rc=rc)
1328  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1329  call error_handler("IN FieldGet", rc)
1330 
1331  print*,"- CALL FieldGet FOR ADJUSTED climo qnwfa."
1332  call esmf_fieldget(qnwfa_climo_target_grid, &
1333  farrayptr=qnwfa2ptr, rc=rc)
1334  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1335  call error_handler("IN FieldGet", rc)
1336 
1337  DO k=1,lev_target
1338  DO i=clb(1),cub(1)
1339  DO j=clb(2),cub(2)
1340  qnifa2ptr(i,j,k) = c2(i,j,k,1)
1341  qnwfa2ptr(i,j,k) = c2(i,j,k,2)
1342  ENDDO
1343  ENDDO
1344  ENDDO
1345 
1346  DEALLOCATE (z1, z2, c1, c2)
1347 
1348  call esmf_fielddestroy(qnifa_climo_b4adj_target_grid, rc=rc)
1349  call esmf_fielddestroy(qnwfa_climo_b4adj_target_grid, rc=rc)
1350  call esmf_fielddestroy(thomp_pres_climo_b4adj_target_grid, rc=rc)
1351 
1352  END SUBROUTINE vintg_thomp_mp_climo
1353 
1367  SUBROUTINE vintg
1368  use mpi
1369 
1370  IMPLICIT NONE
1371 
1372  REAL(ESMF_KIND_R8), PARAMETER :: dltdz=-6.5e-3*287.05/9.80665
1373  REAL(ESMF_KIND_R8), PARAMETER :: dlpvdrt=-2.5e6/461.50
1374  REAL(ESMF_KIND_R8), PARAMETER :: one = 1.0_esmf_kind_r8
1375 
1376  INTEGER :: i, j, k, clb(3), cub(3), rc
1377  INTEGER :: im, km1, km2, nt, ii
1378 
1379  REAL(ESMF_KIND_R8) :: dz
1380  REAL(ESMF_KIND_R8), ALLOCATABLE :: z1(:,:,:), z2(:,:,:)
1381  REAL(ESMF_KIND_R8), ALLOCATABLE :: c1(:,:,:,:),c2(:,:,:,:)
1382 
1383  REAL(ESMF_KIND_R8), POINTER :: p1ptr(:,:,:) ! input pressure
1384  REAL(ESMF_KIND_R8), POINTER :: p2ptr(:,:,:) ! output pressure
1385  REAL(ESMF_KIND_R8), POINTER :: dzdt1ptr(:,:,:) ! input vvel
1386  REAL(ESMF_KIND_R8), POINTER :: dzdt2ptr(:,:,:) ! output vvel
1387  REAL(ESMF_KIND_R8), POINTER :: t1ptr(:,:,:) ! input temperature
1388  REAL(ESMF_KIND_R8), POINTER :: t2ptr(:,:,:) ! output temperature
1389  REAL(ESMF_KIND_R8), POINTER :: q1ptr(:,:,:) ! input tracer
1390  REAL(ESMF_KIND_R8), POINTER :: q2ptr(:,:,:) ! output tracer
1391  REAL(ESMF_KIND_R8), POINTER :: wind1ptr(:,:,:,:) ! input wind (x,y,z components)
1392  REAL(ESMF_KIND_R8), POINTER :: wind2ptr(:,:,:,:) ! input wind (x,y,z components)
1393 
1394 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1395 ! COMPUTE LOG PRESSURE INTERPOLATING COORDINATE
1396 ! AND COPY INPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS
1397 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1398 
1399  print*,"- VERTICALY INTERPOLATE FIELDS."
1400 
1401  print*,"- CALL FieldGet FOR 3-D PRES."
1402  call esmf_fieldget(pres_b4adj_target_grid, &
1403  computationallbound=clb, &
1404  computationalubound=cub, &
1405  farrayptr=p1ptr, rc=rc)
1406  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1407  call error_handler("IN FieldGet", rc)
1408 
1409 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1410 ! The '1'/'2' arrays hold fields before/after interpolation.
1411 ! Note the 'z' component of the horizontal wind will be treated as a
1412 ! tracer. So add one extra third dimension to these 3-d arrays.
1413 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1414 
1415  ALLOCATE(z1(clb(1):cub(1),clb(2):cub(2),lev_input))
1416  ALLOCATE(z2(clb(1):cub(1),clb(2):cub(2),lev_target))
1417  ALLOCATE(c1(clb(1):cub(1),clb(2):cub(2),lev_input,num_tracers+5))
1418  ALLOCATE(c2(clb(1):cub(1),clb(2):cub(2),lev_target,num_tracers+5))
1419 
1420  z1 = -log(p1ptr)
1421 
1422  print*,"- CALL FieldGet FOR 3-D ADJUSTED PRESS"
1423  call esmf_fieldget(pres_target_grid, &
1424  farrayptr=p2ptr, rc=rc)
1425  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1426  call error_handler("IN FieldGet", rc)
1427 
1428  z2 = -log(p2ptr)
1429 
1430  print*,"- CALL FieldGet FOR 3-D WIND."
1431  call esmf_fieldget(wind_b4adj_target_grid, &
1432  farrayptr=wind1ptr, rc=rc)
1433  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1434  call error_handler("IN FieldGet", rc)
1435 
1436  c1(:,:,:,1) = wind1ptr(:,:,:,1)
1437  c1(:,:,:,2) = wind1ptr(:,:,:,2)
1438  c1(:,:,:,3) = wind1ptr(:,:,:,3)
1439 
1440  print*,"- CALL FieldGet FOR VERTICAL VELOCITY."
1441  call esmf_fieldget(dzdt_b4adj_target_grid, &
1442  farrayptr=dzdt1ptr, rc=rc)
1443  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1444  call error_handler("IN FieldGet", rc)
1445 
1446  c1(:,:,:,4) = dzdt1ptr(:,:,:)
1447  print*,"MIN MAX W TARGETB4 IN VINTG = ", minval(dzdt1ptr(:,:,:)), maxval(dzdt1ptr(:,:,:))
1448 
1449  print*,"- CALL FieldGet FOR 3-D TEMP."
1450  call esmf_fieldget(temp_b4adj_target_grid, &
1451  farrayptr=t1ptr, rc=rc)
1452  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1453  call error_handler("IN FieldGet", rc)
1454 
1455  c1(:,:,:,5) = t1ptr(:,:,:)
1456 
1457  DO i = 1, num_tracers
1458 
1459  print*,"- CALL FieldGet FOR 3-D TRACERS ", trim(tracers(i))
1460  call esmf_fieldget(tracers_b4adj_target_grid(i), &
1461  farrayptr=q1ptr, rc=rc)
1462  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1463  call error_handler("IN FieldGet", rc)
1464 
1465  c1(:,:,:,5+i) = q1ptr(:,:,:)
1466 
1467  ENDDO
1468 
1469 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1470 ! PERFORM LAGRANGIAN ONE-DIMENSIONAL INTERPOLATION
1471 ! THAT IS 4TH-ORDER IN INTERIOR, 2ND-ORDER IN OUTSIDE INTERVALS
1472 ! AND 1ST-ORDER FOR EXTRAPOLATION.
1473 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1474 
1475  im = (cub(1)-clb(1)+1) * (cub(2)-clb(2)+1)
1476  km1= lev_input
1477  km2= lev_target
1478  nt= num_tracers + 1 ! treat 'z' wind as tracer.
1479 
1480  CALL terp3(im,1,1,1,1,4+nt,(im*km1),(im*km2), &
1481  km1,im,im,z1,c1,km2,im,im,z2,c2)
1482 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1483 ! COPY OUTPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS
1484 ! EXCEPT BELOW THE INPUT DOMAIN, LET TEMPERATURE INCREASE WITH A FIXED
1485 ! LAPSE RATE AND LET THE RELATIVE HUMIDITY REMAIN CONSTANT.
1486 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1487 
1488  print*,"- CALL FieldGet FOR 3-D ADJUSTED TEMP."
1489  call esmf_fieldget(temp_target_grid, &
1490  farrayptr=t2ptr, rc=rc)
1491  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1492  call error_handler("IN FieldGet", rc)
1493 
1494  print*,"- CALL FieldGet FOR ADJUSTED VERTICAL VELOCITY."
1495  call esmf_fieldget(dzdt_target_grid, &
1496  farrayptr=dzdt2ptr, rc=rc)
1497  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1498  call error_handler("IN FieldGet", rc)
1499 
1500  print*,"- CALL FieldGet FOR 3-D ADJUSTED WIND."
1501  call esmf_fieldget(wind_target_grid, &
1502  farrayptr=wind2ptr, rc=rc)
1503  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1504  call error_handler("IN FieldGet", rc)
1505 
1506  DO k=1,lev_target
1507  DO i=clb(1),cub(1)
1508  DO j=clb(2),cub(2)
1509  wind2ptr(i,j,k,1)=c2(i,j,k,1)
1510  wind2ptr(i,j,k,2)=c2(i,j,k,2)
1511  wind2ptr(i,j,k,3)=c2(i,j,k,3)
1512  dzdt2ptr(i,j,k)=c2(i,j,k,4)
1513  dz=z2(i,j,k)-z1(i,j,1)
1514  IF(dz.GE.0) THEN
1515  t2ptr(i,j,k)=c2(i,j,k,5)
1516  ELSE
1517  t2ptr(i,j,k)=c1(i,j,1,5)*exp(dltdz*dz)
1518  ENDIF
1519  ENDDO
1520  ENDDO
1521  ENDDO
1522 
1523  DO ii = 1, num_tracers
1524 
1525  print*,"- CALL FieldGet FOR 3-D TRACER ", trim(tracers(ii))
1526  call esmf_fieldget(tracers_target_grid(ii), &
1527  farrayptr=q2ptr, rc=rc)
1528  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1529  call error_handler("IN FieldGet", rc)
1530 
1531  IF (trim(tracers(ii)) == "sphum") THEN ! specific humidity
1532 
1533  DO k=1,lev_target
1534  DO i=clb(1),cub(1)
1535  DO j=clb(2),cub(2)
1536  dz=z2(i,j,k)-z1(i,j,1)
1537  IF(dz.GE.0) THEN
1538  q2ptr(i,j,k) = c2(i,j,k,5+ii)
1539  ELSE
1540  q2ptr(i,j,k) = c1(i,j,1,5+ii)*exp(dlpvdrt*(one/t2ptr(i,j,k)-one/t1ptr(i,j,1))-dz)
1541  ENDIF
1542  ENDDO
1543  ENDDO
1544  ENDDO
1545 
1546  ELSE ! all other tracers
1547 
1548  DO k=1,lev_target
1549  DO i=clb(1),cub(1)
1550  DO j=clb(2),cub(2)
1551  q2ptr(i,j,k) = c2(i,j,k,5+ii)
1552  ENDDO
1553  ENDDO
1554  ENDDO
1555 
1556  ENDIF
1557 
1558  ENDDO
1559 
1560  DEALLOCATE (z1, z2, c1, c2)
1561 
1562  END SUBROUTINE vintg
1563 
1600  SUBROUTINE terp3(IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2, &
1601  km1,kxz1,kxq1,z1,q1,km2,kxz2,kxq2,z2,q2)
1602  IMPLICIT NONE
1603  INTEGER im,ixz1,ixq1,ixz2,ixq2,nm,nxq1,nxq2
1604  INTEGER km1,kxz1,kxq1,km2,kxz2,kxq2
1605  INTEGER i,k1,k2,n
1606  INTEGER k1s(im,km2)
1607  REAL(ESMF_KIND_R8), PARAMETER :: one = 1.0_esmf_kind_r8
1608  REAL(ESMF_KIND_R8) :: z1(1+(im-1)*ixz1+(km1-1)*kxz1)
1609  REAL(ESMF_KIND_R8) :: q1(1+(im-1)*ixq1+(km1-1)*kxq1+(nm-1)*nxq1)
1610  REAL(ESMF_KIND_R8) :: z2(1+(im-1)*ixz2+(km2-1)*kxz2)
1611  REAL(ESMF_KIND_R8) :: q2(1+(im-1)*ixq2+(km2-1)*kxq2+(nm-1)*nxq2)
1612 ! REAL(ESMF_KIND_R8) :: J2(1+(IM-1)*IXQ2+(KM2-1)*KXQ2+(NM-1)*NXQ2)
1613  REAL(ESMF_KIND_R8) :: ffa(im),ffb(im),ffc(im),ffd(im)
1614  REAL(ESMF_KIND_R8) :: gga(im),ggb(im),ggc(im),ggd(im)
1615  REAL(ESMF_KIND_R8) :: z1a,z1b,z1c,z1d,q1a,q1b,q1c,q1d,z2s,q2s
1616 ! REAL(ESMF_KIND_R8) :: J2S
1617 
1618 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1619 ! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT.
1620  CALL rsearch(im,km1,ixz1,kxz1,z1,km2,ixz2,kxz2,z2,1,im,k1s)
1621 
1622 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1623 ! GENERALLY INTERPOLATE CUBICALLY WITH MONOTONIC CONSTRAINT
1624 ! FROM TWO NEAREST INPUT POINTS ON EITHER SIDE OF THE OUTPUT POINT,
1625 ! BUT WITHIN THE TWO EDGE INTERVALS INTERPOLATE LINEARLY.
1626 ! KEEP THE OUTPUT FIELDS CONSTANT OUTSIDE THE INPUT DOMAIN.
1627 
1628 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(IM,IXZ1,IXQ1,IXZ2), &
1629 !$OMP& SHARED(IXQ2,NM,NXQ1,NXQ2,KM1,KXZ1,KXQ1,Z1,Q1,KM2,KXZ2), &
1630 !$OMP& SHARED(KXQ2,Z2,Q2,K1S)
1631  DO k2=1,km2
1632  DO i=1,im
1633  k1=k1s(i,k2)
1634  IF(k1.EQ.1.OR.k1.EQ.km1-1) THEN
1635  z2s=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
1636  z1a=z1(1+(i-1)*ixz1+(k1-1)*kxz1)
1637  z1b=z1(1+(i-1)*ixz1+(k1+0)*kxz1)
1638  ffa(i)=(z2s-z1b)/(z1a-z1b)
1639  ffb(i)=(z2s-z1a)/(z1b-z1a)
1640  gga(i)=one/(z1a-z1b)
1641  ggb(i)=one/(z1b-z1a)
1642  ELSEIF(k1.GT.1.AND.k1.LT.km1-1) THEN
1643  z2s=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
1644  z1a=z1(1+(i-1)*ixz1+(k1-2)*kxz1)
1645  z1b=z1(1+(i-1)*ixz1+(k1-1)*kxz1)
1646  z1c=z1(1+(i-1)*ixz1+(k1+0)*kxz1)
1647  z1d=z1(1+(i-1)*ixz1+(k1+1)*kxz1)
1648  ffa(i)=(z2s-z1b)/(z1a-z1b)* &
1649  (z2s-z1c)/(z1a-z1c)* &
1650  (z2s-z1d)/(z1a-z1d)
1651  ffb(i)=(z2s-z1a)/(z1b-z1a)* &
1652  (z2s-z1c)/(z1b-z1c)* &
1653  (z2s-z1d)/(z1b-z1d)
1654  ffc(i)=(z2s-z1a)/(z1c-z1a)* &
1655  (z2s-z1b)/(z1c-z1b)* &
1656  (z2s-z1d)/(z1c-z1d)
1657  ffd(i)=(z2s-z1a)/(z1d-z1a)* &
1658  (z2s-z1b)/(z1d-z1b)* &
1659  (z2s-z1c)/(z1d-z1c)
1660  gga(i)= one/(z1a-z1b)* &
1661  (z2s-z1c)/(z1a-z1c)* &
1662  (z2s-z1d)/(z1a-z1d)+ &
1663  (z2s-z1b)/(z1a-z1b)* &
1664  one/(z1a-z1c)* &
1665  (z2s-z1d)/(z1a-z1d)+ &
1666  (z2s-z1b)/(z1a-z1b)* &
1667  (z2s-z1c)/(z1a-z1c)* &
1668  one/(z1a-z1d)
1669  ggb(i)= one/(z1b-z1a)* &
1670  (z2s-z1c)/(z1b-z1c)* &
1671  (z2s-z1d)/(z1b-z1d)+ &
1672  (z2s-z1a)/(z1b-z1a)* &
1673  one/(z1b-z1c)* &
1674  (z2s-z1d)/(z1b-z1d)+ &
1675  (z2s-z1a)/(z1b-z1a)* &
1676  (z2s-z1c)/(z1b-z1c)* &
1677  one/(z1b-z1d)
1678  ggc(i)= one/(z1c-z1a)* &
1679  (z2s-z1b)/(z1c-z1b)* &
1680  (z2s-z1d)/(z1c-z1d)+ &
1681  (z2s-z1a)/(z1c-z1a)* &
1682  one/(z1c-z1b)* &
1683  (z2s-z1d)/(z1c-z1d)+ &
1684  (z2s-z1a)/(z1c-z1a)* &
1685  (z2s-z1b)/(z1c-z1b)* &
1686  one/(z1c-z1d)
1687  ggd(i)= one/(z1d-z1a)* &
1688  (z2s-z1b)/(z1d-z1b)* &
1689  (z2s-z1c)/(z1d-z1c)+ &
1690  (z2s-z1a)/(z1d-z1a)* &
1691  one/(z1d-z1b)* &
1692  (z2s-z1c)/(z1d-z1c)+ &
1693  (z2s-z1a)/(z1d-z1a)* &
1694  (z2s-z1b)/(z1d-z1b)* &
1695  one/(z1d-z1c)
1696  ENDIF
1697  ENDDO
1698 
1699 ! INTERPOLATE.
1700  DO n=1,nm
1701  DO i=1,im
1702  k1=k1s(i,k2)
1703  IF(k1.EQ.0) THEN
1704  q2s=q1(1+(i-1)*ixq1+(n-1)*nxq1)
1705 ! J2S=0
1706  ELSEIF(k1.EQ.km1) THEN
1707  q2s=q1(1+(i-1)*ixq1+(km1-1)*kxq1+(n-1)*nxq1)
1708 ! J2S=0
1709  ELSEIF(k1.EQ.1.OR.k1.EQ.km1-1) THEN
1710  q1a=q1(1+(i-1)*ixq1+(k1-1)*kxq1+(n-1)*nxq1)
1711  q1b=q1(1+(i-1)*ixq1+(k1+0)*kxq1+(n-1)*nxq1)
1712  q2s=ffa(i)*q1a+ffb(i)*q1b
1713 ! J2S=GGA(I)*Q1A+GGB(I)*Q1B
1714  ELSE
1715  q1a=q1(1+(i-1)*ixq1+(k1-2)*kxq1+(n-1)*nxq1)
1716  q1b=q1(1+(i-1)*ixq1+(k1-1)*kxq1+(n-1)*nxq1)
1717  q1c=q1(1+(i-1)*ixq1+(k1+0)*kxq1+(n-1)*nxq1)
1718  q1d=q1(1+(i-1)*ixq1+(k1+1)*kxq1+(n-1)*nxq1)
1719  q2s=ffa(i)*q1a+ffb(i)*q1b+ffc(i)*q1c+ffd(i)*q1d
1720 ! J2S=GGA(I)*Q1A+GGB(I)*Q1B+GGC(I)*Q1C+GGD(I)*Q1D
1721  IF(q2s.LT.min(q1b,q1c)) THEN
1722  q2s=min(q1b,q1c)
1723 ! J2S=0
1724  ELSEIF(q2s.GT.max(q1b,q1c)) THEN
1725  q2s=max(q1b,q1c)
1726 ! J2S=0
1727  ENDIF
1728  ENDIF
1729  q2(1+(i-1)*ixq2+(k2-1)*kxq2+(n-1)*nxq2)=q2s
1730 ! J2(1+(I-1)*IXQ2+(K2-1)*KXQ2+(N-1)*NXQ2)=J2S
1731  ENDDO
1732  ENDDO
1733  ENDDO
1734 !$OMP END PARALLEL DO
1735 
1736  END SUBROUTINE terp3
1737 
1794  SUBROUTINE rsearch(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,IXL2,KXL2,L2)
1795  IMPLICIT NONE
1796 
1797  INTEGER,INTENT(IN) :: im,km1,ixz1,kxz1,km2,ixz2,kxz2,ixl2,kxl2
1798  INTEGER,INTENT(OUT) :: l2(1+(im-1)*ixl2+(km2-1)*kxl2)
1799 
1800  REAL(ESMF_KIND_R8),INTENT(IN) :: z1(1+(im-1)*ixz1+(km1-1)*kxz1)
1801  REAL(ESMF_KIND_R8),INTENT(IN) :: z2(1+(im-1)*ixz2+(km2-1)*kxz2)
1802 
1803  INTEGER :: i,k2,l
1804 
1805  REAL(ESMF_KIND_R8) :: z
1806 
1807 
1808 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1809 ! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT.
1810  DO i=1,im
1811  IF (z1(1+(i-1)*ixz1).LE.z1(1+(i-1)*ixz1+(km1-1)*kxz1)) THEN
1812 ! INPUT COORDINATE IS MONOTONICALLY ASCENDING.
1813  DO k2=1,km2
1814  z=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
1815  l=0
1816  DO
1817  IF(z.LT.z1(1+(i-1)*ixz1+l*kxz1)) EXIT
1818  l=l+1
1819  IF(l.EQ.km1) EXIT
1820  ENDDO
1821  l2(1+(i-1)*ixl2+(k2-1)*kxl2)=l
1822  ENDDO
1823  ELSE
1824 ! INPUT COORDINATE IS MONOTONICALLY DESCENDING.
1825  DO k2=1,km2
1826  z=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
1827  l=0
1828  DO
1829  IF(z.GT.z1(1+(i-1)*ixz1+l*kxz1)) EXIT
1830  l=l+1
1831  IF(l.EQ.km1) EXIT
1832  ENDDO
1833  l2(1+(i-1)*ixl2+(k2-1)*kxl2)=l
1834  ENDDO
1835  ENDIF
1836  ENDDO
1837 
1838  END SUBROUTINE rsearch
1839 
1842  subroutine compute_zh
1843 
1844  implicit none
1845 
1846  integer :: i,ii, j,k, rc, clb(2), cub(2)
1847 
1848  real(esmf_kind_r8), allocatable :: pe0(:), pn0(:)
1849  real(esmf_kind_r8), pointer :: psptr(:,:)
1850  real(esmf_kind_r8), pointer :: zhsfcptr(:,:)
1851  real(esmf_kind_r8), pointer :: zhptr(:,:,:)
1852  real(esmf_kind_r8), pointer :: tptr(:,:,:)
1853  real(esmf_kind_r8), pointer :: qptr(:,:,:)
1854  real(esmf_kind_r8) :: ak, bk, zvir, grd
1855  real(esmf_kind_r8), parameter :: grav = 9.80665
1856  real(esmf_kind_r8), parameter :: rdgas = 287.05
1857  real(esmf_kind_r8), parameter :: rvgas = 461.50
1858 
1859  print*,"- COMPUTE HEIGHT"
1860 
1861  print*,"- CALL FieldGet FOR SURFACE PRESSURE"
1862  call esmf_fieldget(ps_target_grid, &
1863  computationallbound=clb, &
1864  computationalubound=cub, &
1865  farrayptr=psptr, 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  print*,"- CALL FieldGet FOR TERRAIN HEIGHT"
1870  call esmf_fieldget(terrain_target_grid, &
1871  farrayptr=zhsfcptr, rc=rc)
1872  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1873  call error_handler("IN FieldGet", rc)
1874 
1875  print*,"- CALL FieldGet FOR HEIGHT"
1876  call esmf_fieldget(zh_target_grid, &
1877  farrayptr=zhptr, rc=rc)
1878  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1879  call error_handler("IN FieldGet", rc)
1880 
1881  print*,"- CALL FieldGet FOR TEMPERATURE"
1882  call esmf_fieldget(temp_target_grid, &
1883  farrayptr=tptr, rc=rc)
1884  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1885  call error_handler("IN FieldGet", rc)
1886 
1887  do ii = 1, num_tracers
1888  if (trim(tracers(ii)) == "sphum") exit
1889  enddo
1890 
1891  print*,"- CALL FieldGet FOR SPECIFIC HUMIDITY"
1892  call esmf_fieldget(tracers_target_grid(ii), &
1893  farrayptr=qptr, rc=rc)
1894  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1895  call error_handler("IN FieldGet", rc)
1896 
1897  grd = grav/rdgas
1898  zvir = rvgas/rdgas - 1.0_esmf_kind_r8
1899 
1900  allocate(pe0(levp1_target))
1901  allocate(pn0(levp1_target))
1902 
1903  do j = clb(2), cub(2)
1904  do i = clb(1), cub(1)
1905 
1906  do k = 1, levp1_target
1907  ak = vcoord_target(k,1)
1908  ak = max(ak, 1.e-9)
1909  bk = vcoord_target(k,2)
1910 
1911  pe0(k) = ak + bk*psptr(i,j)
1912  pn0(k) = log(pe0(k))
1913  enddo
1914 
1915  zhptr(i,j,1) = zhsfcptr(i,j)
1916 
1917  do k = 2, levp1_target
1918  zhptr(i,j,k) = zhptr(i,j,k-1)+tptr(i,j,k-1)*(1.+zvir*qptr(i,j,k-1))* &
1919  (pn0(k-1)-pn0(k))/grd
1920  enddo
1921 
1922  enddo
1923  enddo
1924 
1925  deallocate(pe0, pn0)
1926 
1927  end subroutine compute_zh
1928 
1932  subroutine cleanup_target_atm_b4adj_data
1933 
1934  implicit none
1935 
1936  integer :: i, rc
1937 
1938  print*,"- DESTROY TARGET GRID ATMOSPHERIC BEFORE ADJUSTMENT FIELDS."
1939 
1940  call esmf_fielddestroy(wind_b4adj_target_grid, rc=rc)
1941  call esmf_fielddestroy(dzdt_b4adj_target_grid, rc=rc)
1942  call esmf_fielddestroy(ps_b4adj_target_grid, rc=rc)
1943  call esmf_fielddestroy(pres_b4adj_target_grid, rc=rc)
1944  call esmf_fielddestroy(temp_b4adj_target_grid, rc=rc)
1945  call esmf_fielddestroy(terrain_interp_to_target_grid, rc=rc)
1946 
1947  do i = 1, num_tracers
1948  call esmf_fielddestroy(tracers_b4adj_target_grid(i), rc=rc)
1949  enddo
1950 
1951  deallocate(tracers_b4adj_target_grid)
1952 
1953  end subroutine cleanup_target_atm_b4adj_data
1954 
1957  subroutine cleanup_target_atm_data
1958 
1959  implicit none
1960 
1961  integer :: i, rc
1962 
1963  print*,"- DESTROY TARGET GRID ATMOSPHERIC FIELDS."
1964 
1965  call esmf_fielddestroy(delp_target_grid, rc=rc)
1966  call esmf_fielddestroy(dzdt_target_grid, rc=rc)
1967  call esmf_fielddestroy(ps_target_grid, rc=rc)
1968  call esmf_fielddestroy(pres_target_grid, rc=rc)
1969  call esmf_fielddestroy(temp_target_grid, rc=rc)
1970  call esmf_fielddestroy(u_s_target_grid, rc=rc)
1971  call esmf_fielddestroy(v_s_target_grid, rc=rc)
1972  call esmf_fielddestroy(wind_target_grid, rc=rc)
1973  call esmf_fielddestroy(wind_s_target_grid, rc=rc)
1974  call esmf_fielddestroy(wind_w_target_grid, rc=rc)
1975  call esmf_fielddestroy(u_w_target_grid, rc=rc)
1976  call esmf_fielddestroy(v_w_target_grid, rc=rc)
1977  call esmf_fielddestroy(zh_target_grid, rc=rc)
1978 
1979  do i = 1, num_tracers
1980  call esmf_fielddestroy(tracers_target_grid(i), rc=rc)
1981  enddo
1982 
1983  deallocate(tracers_target_grid)
1984 
1985  if (esmf_fieldiscreated(qnifa_climo_target_grid)) then
1986  call esmf_fielddestroy(qnifa_climo_target_grid, rc=rc)
1987  endif
1988 
1989  if (esmf_fieldiscreated(qnwfa_climo_target_grid)) then
1990  call esmf_fielddestroy(qnwfa_climo_target_grid, rc=rc)
1991  endif
1992 
1993  end subroutine cleanup_target_atm_data
1994 
1995  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:1932
atmosphere::vintg
subroutine vintg
Vertically interpolate upper-air fields.
Definition: atmosphere.F90:1367
atmosphere::vintg_thomp_mp_climo
subroutine vintg_thomp_mp_climo
Vertically interpolate atmospheric fields to target FV3 grid.
Definition: atmosphere.F90:1244
input_data::cleanup_input_atm_data
subroutine, public cleanup_input_atm_data
Free up memory associated with atm data.
Definition: input_data.F90:6472
atmosphere::newps
subroutine newps(localpet)
Computes adjusted surface pressure given a new terrain height.
Definition: atmosphere.F90:900
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:432
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:1097
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:1600
atmosphere::atmosphere_driver
subroutine, public atmosphere_driver(localpet)
Driver routine to process for atmospheric fields.
Definition: atmosphere.F90:110
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:506
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:144
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:1794
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:1132
program_setup
This module contains code to read the setup namelist file, handle the varmap file for GRIB2 data...
Definition: program_setup.F90:9
atmosphere::compute_zh
subroutine compute_zh
Compute vertical level height.
Definition: atmosphere.F90:1842
atmosphere::newpr1
subroutine newpr1(localpet)
Computes 3-D pressure given an adjusted surface pressure.
Definition: atmosphere.F90:789
atmosphere::convert_winds
subroutine convert_winds
Convert 3-d component winds to u and v.
Definition: atmosphere.F90:645
atmosphere::cleanup_target_atm_data
subroutine cleanup_target_atm_data
Cleanup target grid atmospheric field objects.
Definition: atmosphere.F90:1957