FV3DYCORE  Version 1.1.0
atmosphere.F90
Go to the documentation of this file.
1 
2 !***********************************************************************
3 !* GNU Lesser General Public License
4 !*
5 !* This file is part of the FV3 dynamical core.
6 !*
7 !* The FV3 dynamical core is free software: you can redistribute it
8 !* and/or modify it under the terms of the
9 !* GNU Lesser General Public License as published by the
10 !* Free Software Foundation, either version 3 of the License, or
11 !* (at your option) any later version.
12 !*
13 !* The FV3 dynamical core is distributed in the hope that it will be
14 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
15 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16 !* See the GNU General Public License for more details.
17 !*
18 !* You should have received a copy of the GNU Lesser General Public
19 !* License along with the FV3 dynamical core.
20 !* If not, see <http://www.gnu.org/licenses/>.
21 !***********************************************************************
22 
25 
27 
28 ! <table>
29 ! <tr>
30 ! <th>Module Name</th>
31 ! <th>Functions Included</th>
32 ! </tr>
33 ! <tr>
34 ! <td>block_control_mod</td>
35 ! <td>block_control_type</td>
36 ! </tr>
37 ! <tr>
38 ! <td>constants_mod</td>
39 ! <td>cp_air, rdgas, grav, rvgas, kappa, pstd_mks</td>
40 ! </tr>
41 ! <tr>
42 ! <td>field_manager_mod</td>
43 ! <td>MODEL_ATMOS</td>
44 ! </tr>
45 ! <tr>
46 ! <td>fms_mod</td>
47 ! <td>file_exist, open_namelist_file,close_file, error_mesg, FATAL,
48 ! check_nml_error, stdlog,write_version_number,set_domain,
49 ! mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_SUBCOMPONENT,
50 ! clock_flag_default, nullify_domain</td>
51 ! </tr>
52 ! <tr>
53 ! <td>fv_arrays_mod</td>
54 ! <td>fv_atmos_type, R_GRID</td>
55 ! </tr>
56 ! <tr>
57 ! <td>fv_control_mod</td>
58 ! <td>fv_init, fv_end, ngrids</td>
59 ! </tr>
60 ! <tr>
61 ! <td>fv_diagnostics_mod</td>
62 ! <td>fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height</td>
63 ! </tr>
64 ! <tr>
65 ! <td>fv_dynamics_mod</td>
66 ! <td>fv_dynamics</td>
67 ! </tr>
68 ! <tr>
69 ! <td>fv_eta_mod</td>
70 ! <td>get_eta_level</td>
71 ! </tr>
72 ! <tr>
73 ! <td>fv_fill_mod</td>
74 ! <td>fill_gfs</td>
75 ! </tr>
76 ! <tr>
77 ! <td>fv_mp_mod</td>
78 ! <td>switch_current_Atm</td>
79 ! </tr>
80 ! <tr>
81 ! <td>fv_nesting_mod</td>
82 ! <td>twoway_nesting</td>
83 ! </tr>
84 ! <tr>
85 ! <td>fv_nggps_diags_mod</td>
86 ! <td>fv_nggps_diag_init, fv_nggps_diag</td>
87 ! </tr>
88 ! <tr>
89 ! <td>fv_nwp_nudge_mod</td>
90 ! <td>fv_nwp_nudge_init, fv_nwp_nudge_end, do_adiabatic_init</td>
91 ! </tr>
92 ! <tr>
93 ! <td>fv_restart_mod</td>
94 ! <td>fv_restart, fv_write_restart</td>
95 ! </tr>
96 ! <tr>
97 ! <td>fv_sg_mod</td>
98 ! <td>fv_subgrid_z</td>
99 ! </tr>
100 ! <tr>
101 ! <td>fv_timing_mod</td>
102 ! <td>timing_on, timing_off</td>
103 ! </tr>
104 ! <tr>
105 ! <td>fv_update_phys_mod</td>
106 ! <td>fv_update_phys</td>
107 ! </tr>
108 ! <tr>
109 ! <td>IPD_typedefs_mod</td>
110 ! <td>IPD_data_type, kind_phys => IPD_kind_phys</td>
111 ! </tr>
112 ! <tr>
113 ! <td>mpp_mod</td>
114 ! <td>mpp_error, stdout, FATAL, NOTE, input_nml_file, mpp_root_pe,
115 ! mpp_npes, mpp_pe, mpp_chksum,mpp_get_current_pelist,
116 ! mpp_set_current_pelist</td>
117 ! </tr>
118 ! <tr>
119 ! <td>mpp_domains_mod</td>
120 ! <td>mpp_get_data_domain, mpp_get_compute_domain, domain2d, mpp_update_domains</td>
121 ! </tr>
122 ! <tr>
123 ! <td>mpp_parameter_mod</td>
124 ! <td>EUPDATE, WUPDATE, SUPDATE, NUPDATE</td>
125 ! </tr>
126 ! <tr>
127 ! <td>time_manager_mod</td>
128 ! <td>time_type, get_time, set_time, operator(+), operator(-)</td>
129 ! </tr>
130 ! <tr>
131 ! <td>tracer_manager_mod</td>
132 ! <td>get_tracer_index, get_number_tracers, NO_TRACER</td>
133 ! </tr>
134 ! <tr>
135 ! <td>xgrid_mod</td>
136 ! <td>grid_box_type</td>
137 ! </tr>
138 ! </table>
139 
140 #include <fms_platform.h>
141 
142 !-----------------
143 ! FMS modules:
144 !-----------------
145 use block_control_mod, only: block_control_type
146 use constants_mod, only: cp_air, rdgas, grav, rvgas, kappa, pstd_mks
147 use time_manager_mod, only: time_type, get_time, set_time, operator(+), &
148  operator(-), operator(/), time_type_to_real
149 use fms_mod, only: file_exist, open_namelist_file, &
150  close_file, error_mesg, fatal, &
151  check_nml_error, stdlog, &
152  write_version_number, &
153  set_domain, &
154  read_data, &
155  mpp_clock_id, mpp_clock_begin, &
156  mpp_clock_end, clock_subcomponent, &
157  clock_flag_default, nullify_domain
158 use mpp_mod, only: mpp_error, stdout, fatal, note, &
159  input_nml_file, mpp_root_pe, &
160  mpp_npes, mpp_pe, mpp_chksum, &
161  mpp_get_current_pelist, &
162  mpp_set_current_pelist
163 use mpp_parameter_mod, only: eupdate, wupdate, supdate, nupdate
164 use mpp_domains_mod, only: domain2d, mpp_update_domains
165 use xgrid_mod, only: grid_box_type
166 use field_manager_mod, only: model_atmos
167 use tracer_manager_mod, only: get_tracer_index, get_number_tracers, &
168  no_tracer
170 use ipd_typedefs, only: ipd_data_type, kind_phys => ipd_kind_phys
172 
173 !-----------------
174 ! FV core modules:
175 !-----------------
177 use fv_control_mod, only: fv_init, fv_end, ngrids
178 use fv_eta_mod, only: get_eta_level
179 use fv_fill_mod, only: fill_gfs
180 use fv_dynamics_mod, only: fv_dynamics
186 use fv_mp_mod, only: switch_current_atm
187 use fv_sg_mod, only: fv_subgrid_z
190 #ifdef MULTI_GASES
191 use multi_gases_mod, only: virq, virq_max, num_gas, ri, cpi
192 #endif
195 
196 use mpp_domains_mod, only: mpp_get_data_domain, mpp_get_compute_domain
197 
198 implicit none
199 private
200 
201 !--- driver routines
204 
205 !--- utility routines
211 ! experimental APIs not ready for use
212 ! atmosphere_tracer_postinit, &
213  atmosphere_diss_est, & ! dissipation estimate for SKEB
217 
218 !--- physics/radiation data exchange routines
220 
221 !-----------------------------------------------------------------------
222 ! version number of this module
223 ! Include variable "version" to be written to log file.
224 #include<file_version.h>
225 character(len=20) :: mod_name = 'fvGFS/atmosphere_mod'
226 
227 !---- private data ----
228  type(time_type) :: time_step_atmos
229  public atm, mytile
230 
231  !These are convenience variables for local use only, and are set to values in Atm%
232  real :: dt_atmos
233  real :: zvir
234  integer :: npx, npy, npz, ncnst, pnats
235  integer :: isc, iec, jsc, jec
236  integer :: isd, ied, jsd, jed
237  integer :: nq ! number of transported tracers
238  integer :: sec, seconds, days
240  logical :: cold_start = .false. ! used in initial condition
241 
242  integer, dimension(:), allocatable :: id_tracerdt_dyn
243  integer :: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel ! condensate species tracer indices
244 #ifdef CCPP
245  integer :: cld_amt
246 #endif
247 
248  integer :: mytile = 1
249  integer :: p_split = 1
250  integer, allocatable :: pelist(:)
251  logical, allocatable :: grids_on_this_pe(:)
252  type(fv_atmos_type), allocatable, target :: atm(:)
253 
254  integer :: id_udt_dyn, id_vdt_dyn
255 
256  real, parameter:: w0_big = 60. ! to prevent negative w-tracer diffusion
257 
258 !---dynamics tendencies for use in fv_subgrid_z and during fv_update_phys
259  real, allocatable, dimension(:,:,:) :: u_dt, v_dt, t_dt
260  real, allocatable :: pref(:,:), dum1d(:)
261 
262  logical :: first_diag = .true.
263 
264 contains
265 
266 
270  subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
271 #ifdef CCPP
272 #ifdef STATIC
273 ! For static builds, the ccpp_physics_{init,run,finalize} calls
274 ! are not pointing to code in the CCPP framework, but to auto-generated
275 ! ccpp_suite_cap and ccpp_group_*_cap modules behind a ccpp_static_api
276  use ccpp_api, only: ccpp_init
277  use ccpp_static_api, only: ccpp_physics_init
278 #else
279  use iso_c_binding, only: c_loc
280  use ccpp_api, only: ccpp_init, &
281  ccpp_physics_init, &
282  ccpp_field_add, &
283  ccpp_error
284 #endif
285  use ccpp_data, only: ccpp_suite, &
286  cdata => cdata_tile, &
287  ccpp_interstitial
288 #ifdef OPENMP
289  use omp_lib
290 #endif
291 #ifndef STATIC
292 ! Begin include auto-generated list of modules for ccpp
293 #include "ccpp_modules_fast_physics.inc"
294 ! End include auto-generated list of modules for ccpp
295 #endif
296 #endif
297  type(time_type), intent(in) :: time_init, time, time_step
298  type(grid_box_type), intent(inout) :: Grid_box
299  real(kind=kind_phys), pointer, dimension(:,:), intent(inout) :: area
300 !--- local variables ---
301  integer :: i, n
302 ! integer :: itrac
303  logical :: do_atmos_nudge
304  character(len=32) :: tracer_name, tracer_units
305  real :: ps1, ps2
306 #ifdef CCPP
307  integer :: nthreads
308  integer :: ierr
309 #endif
310 
311  current_time_in_seconds = time_type_to_real( time - time_init )
312  if (mpp_pe() == 0) write(*,"('atmosphere_init: current_time_seconds = ',f9.1)")current_time_in_seconds
313 
314  call timing_on('ATMOS_INIT')
315  allocate(pelist(mpp_npes()))
316  call mpp_get_current_pelist(pelist)
317 
318  zvir = rvgas/rdgas - 1.
319 
320 !---- compute physics/atmos time step in seconds ----
321 
322  time_step_atmos = time_step
323  call get_time (time_step_atmos, sec)
324  dt_atmos = real(sec)
325 
326 !----- initialize FV dynamical core -----
327  !NOTE do we still need the second file_exist call?
328  cold_start = (.not.file_exist('INPUT/fv_core.res.nc') .and. .not.file_exist('INPUT/fv_core.res.tile1.nc'))
329 
330  call fv_init( atm, dt_atmos, grids_on_this_pe, p_split ) ! allocates Atm components
331 
332  do n=1,ngrids
333  if (grids_on_this_pe(n)) mytile = n
334  enddo
335 
336  atm(mytile)%Time_init = time_init
337 
338  a_step = 0
339  if(atm(mytile)%flagstruct%warm_start) then
341  endif
342 
343 !----- write version and namelist to log file -----
344  call write_version_number ( 'fvGFS/ATMOSPHERE_MOD', version )
345 
346 !-----------------------------------
347 
348  npx = atm(mytile)%npx
349  npy = atm(mytile)%npy
350  npz = atm(mytile)%npz
351  ncnst = atm(mytile)%ncnst
352  pnats = atm(mytile)%flagstruct%pnats
353 
354  isc = atm(mytile)%bd%isc
355  iec = atm(mytile)%bd%iec
356  jsc = atm(mytile)%bd%jsc
357  jec = atm(mytile)%bd%jec
358 
359  isd = isc - atm(mytile)%bd%ng
360  ied = iec + atm(mytile)%bd%ng
361  jsd = jsc - atm(mytile)%bd%ng
362  jed = jec + atm(mytile)%bd%ng
363 
364  nq = ncnst-pnats
365  sphum = get_tracer_index(model_atmos, 'sphum' )
366  liq_wat = get_tracer_index(model_atmos, 'liq_wat' )
367  ice_wat = get_tracer_index(model_atmos, 'ice_wat' )
368  rainwat = get_tracer_index(model_atmos, 'rainwat' )
369  snowwat = get_tracer_index(model_atmos, 'snowwat' )
370  graupel = get_tracer_index(model_atmos, 'graupel' )
371 #ifdef CCPP
372  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
373 #endif
374 
375  if (max(sphum,liq_wat,ice_wat,rainwat,snowwat,graupel) > atm(mytile)%flagstruct%nwat) then
376  call mpp_error (fatal,' atmosphere_init: condensate species are not first in the list of &
377  &tracers defined in the field_table')
378  endif
379 
380  ! Allocate grid variables to be used to calculate gradient in 2nd order flux exchange
381  ! This data is only needed for the COARSEST grid.
382  call switch_current_atm(atm(mytile))
383 
384  allocate(grid_box%dx ( isc:iec , jsc:jec+1))
385  allocate(grid_box%dy ( isc:iec+1, jsc:jec ))
386  allocate(grid_box%area ( isc:iec , jsc:jec ))
387  allocate(grid_box%edge_w( jsc:jec+1))
388  allocate(grid_box%edge_e( jsc:jec+1))
389  allocate(grid_box%edge_s( isc:iec+1 ))
390  allocate(grid_box%edge_n( isc:iec+1 ))
391  allocate(grid_box%en1 (3, isc:iec , jsc:jec+1))
392  allocate(grid_box%en2 (3, isc:iec+1, jsc:jec ))
393  allocate(grid_box%vlon (3, isc:iec , jsc:jec ))
394  allocate(grid_box%vlat (3, isc:iec , jsc:jec ))
395  grid_box%dx ( isc:iec , jsc:jec+1) = atm(mytile)%gridstruct%dx ( isc:iec, jsc:jec+1)
396  grid_box%dy ( isc:iec+1, jsc:jec ) = atm(mytile)%gridstruct%dy ( isc:iec+1, jsc:jec )
397  grid_box%area ( isc:iec , jsc:jec ) = atm(mytile)%gridstruct%area ( isc:iec , jsc:jec )
398  grid_box%edge_w( jsc:jec+1) = atm(mytile)%gridstruct%edge_w( jsc:jec+1)
399  grid_box%edge_e( jsc:jec+1) = atm(mytile)%gridstruct%edge_e( jsc:jec+1)
400  grid_box%edge_s( isc:iec+1 ) = atm(mytile)%gridstruct%edge_s( isc:iec+1)
401  grid_box%edge_n( isc:iec+1 ) = atm(mytile)%gridstruct%edge_n( isc:iec+1)
402  grid_box%en1 (:, isc:iec , jsc:jec+1) = atm(mytile)%gridstruct%en1 (:, isc:iec , jsc:jec+1)
403  grid_box%en2 (:, isc:iec+1, jsc:jec ) = atm(mytile)%gridstruct%en2 (:, isc:iec+1, jsc:jec )
404  do i = 1,3
405  grid_box%vlon(i, isc:iec , jsc:jec ) = atm(mytile)%gridstruct%vlon (isc:iec , jsc:jec, i )
406  grid_box%vlat(i, isc:iec , jsc:jec ) = atm(mytile)%gridstruct%vlat (isc:iec , jsc:jec, i )
407  enddo
408  allocate (area(isc:iec , jsc:jec ))
409  area(isc:iec,jsc:jec) = atm(mytile)%gridstruct%area_64(isc:iec,jsc:jec)
410 
411 !----- allocate and zero out the dynamics (and accumulated) tendencies
412  allocate( u_dt(isd:ied,jsd:jed,npz), &
413  v_dt(isd:ied,jsd:jed,npz), &
414  t_dt(isc:iec,jsc:jec,npz) )
415 !--- allocate pref
416  allocate(pref(npz+1,2), dum1d(npz+1))
417 
418  call set_domain ( atm(mytile)%domain )
419  call fv_restart(atm(mytile)%domain, atm, dt_atmos, seconds, days, cold_start, atm(mytile)%gridstruct%grid_type, grids_on_this_pe)
420 
421  fv_time = time
422 
423 !----- initialize atmos_axes and fv_dynamics diagnostics
424  !I've had trouble getting this to work with multiple grids at a time; worth revisiting?
425  call fv_diag_init(atm(mytile:mytile), atm(mytile)%atmos_axes, time, npx, npy, npz, atm(mytile)%flagstruct%p_ref)
426 
427 !---------- reference profile -----------
428  ps1 = 101325.
429  ps2 = 81060.
430  pref(npz+1,1) = ps1
431  pref(npz+1,2) = ps2
432  call get_eta_level ( npz, ps1, pref(1,1), dum1d, atm(mytile)%ak, atm(mytile)%bk )
433  call get_eta_level ( npz, ps2, pref(1,2), dum1d, atm(mytile)%ak, atm(mytile)%bk )
434 
435 ! --- initialize clocks for dynamics, physics_down and physics_up
436  id_dynam = mpp_clock_id('FV dy-core', flags = clock_flag_default, grain=clock_subcomponent )
437  id_subgridz = mpp_clock_id('FV subgrid_z',flags = clock_flag_default, grain=clock_subcomponent )
438  id_fv_diag = mpp_clock_id('FV Diag', flags = clock_flag_default, grain=clock_subcomponent )
439 
440  call timing_off('ATMOS_INIT')
441 
442 #ifdef CCPP
443  ! Do CCPP fast physics initialization before call to adiabatic_init (since this calls fv_dynamics)
444 
445  ! Initialize the cdata structure
446  call ccpp_init(trim(ccpp_suite), cdata, ierr)
447  if (ierr/=0) then
448  cdata%errmsg = ' atmosphere_dynamics: error in ccpp_init: ' // trim(cdata%errmsg)
449  call mpp_error (fatal, cdata%errmsg)
450  end if
451 
452  ! For fast physics running over the entire domain, block and thread
453  ! number are not used; set to safe values
454  cdata%blk_no = 1
455  cdata%thrd_no = 1
456 
457  ! Create shared data type for fast and slow physics, one for each thread
458 #ifdef OPENMP
459  nthreads = omp_get_max_threads()
460 #else
461  nthreads = 1
462 #endif
463  ! Create interstitial data type for fast physics; for multi-gases physics,
464  ! pass q(:,:,:,1:num_gas) as qvi, otherwise pass q(:,:,:,1:1) as 4D array
465  call ccpp_interstitial%create(atm(mytile)%bd%is, atm(mytile)%bd%ie, atm(mytile)%bd%isd, atm(mytile)%bd%ied, &
466  atm(mytile)%bd%js, atm(mytile)%bd%je, atm(mytile)%bd%jsd, atm(mytile)%bd%jed, &
467  atm(mytile)%npz, atm(mytile)%ng, &
468  dt_atmos, p_split, atm(mytile)%flagstruct%k_split, &
469  zvir, atm(mytile)%flagstruct%p_ref, atm(mytile)%ak, atm(mytile)%bk, &
470  cld_amt>0, kappa, atm(mytile)%flagstruct%hydrostatic, &
471  atm(mytile)%flagstruct%do_sat_adj, &
472  atm(mytile)%delp, atm(mytile)%delz, atm(mytile)%gridstruct%area_64, &
473  atm(mytile)%peln, atm(mytile)%phis, atm(mytile)%pkz, atm(mytile)%pt, &
474 #ifdef MULTI_GASES
475  atm(mytile)%q(:,:,:,1:max(1,num_gas)), &
476 #else
477  atm(mytile)%q(:,:,:,1:1), &
478 #endif
479  atm(mytile)%q(:,:,:,sphum), atm(mytile)%q(:,:,:,liq_wat), &
480  atm(mytile)%q(:,:,:,ice_wat), atm(mytile)%q(:,:,:,rainwat), &
481  atm(mytile)%q(:,:,:,snowwat), atm(mytile)%q(:,:,:,graupel), &
482  atm(mytile)%q(:,:,:,cld_amt), atm(mytile)%q_con, nthreads, &
483  atm(mytile)%flagstruct%nwat, &
484 #ifdef MULTI_GASES
485  ngas=num_gas, rilist=ri, cpilist=cpi, &
486 #endif
487  mpirank=mpp_pe(), mpiroot=mpp_root_pe())
488 
489 #ifndef STATIC
490 ! Populate cdata structure with fields required to run fast physics (auto-generated).
491 #include "ccpp_fields_fast_physics.inc"
492 #endif
493 
494  if (atm(mytile)%flagstruct%do_sat_adj) then
495  ! Initialize fast physics
496 #ifdef STATIC
497  call ccpp_physics_init(cdata, suite_name=trim(ccpp_suite), group_name="fast_physics", ierr=ierr)
498 #else
499  call ccpp_physics_init(cdata, group_name="fast_physics", ierr=ierr)
500 #endif
501  if (ierr/=0) then
502  cdata%errmsg = ' atmosphere_dynamics: error in ccpp_physics_init for group fast_physics: ' // trim(cdata%errmsg)
503  call mpp_error (fatal, cdata%errmsg)
504  end if
505  end if
506 #endif
507 
508 ! --- initiate the start for a restarted regional forecast
509  if ( atm(mytile)%gridstruct%regional .and. atm(mytile)%flagstruct%warm_start ) then
510 
511  call start_regional_restart(atm(1), &
512  isc, iec, jsc, jec, &
513  isd, ied, jsd, jed )
514  endif
515 
516  if ( atm(mytile)%flagstruct%na_init>0 ) then
517  call nullify_domain ( )
518  if ( .not. atm(mytile)%flagstruct%hydrostatic ) then
519  call prt_maxmin('Before adi: W', atm(mytile)%w, isc, iec, jsc, jec, atm(mytile)%ng, npz, 1.)
520  endif
521  call adiabatic_init(zvir,atm(mytile)%flagstruct%nudge_dz,time)
522  if ( .not. atm(mytile)%flagstruct%hydrostatic ) then
523  call prt_maxmin('After adi: W', atm(mytile)%w, isc, iec, jsc, jec, atm(mytile)%ng, npz, 1.)
524 ! Not nested?
525  call prt_height('na_ini Z500', isc,iec, jsc,jec, 3, npz, 500.e2, atm(mytile)%phis, atm(mytile)%delz, &
526  atm(mytile)%peln, atm(mytile)%gridstruct%area_64(isc:iec,jsc:jec), atm(mytile)%gridstruct%agrid_64(isc:iec,jsc:jec,2))
527  endif
528  else
529  call mpp_error(note,'No adiabatic initialization correction in use')
530  endif
531 
532 #ifdef DEBUG
533  call nullify_domain()
534  call fv_diag(atm(mytile:mytile), zvir, time, -1)
535 #endif
536 
537  n = mytile
538  call switch_current_atm(atm(n))
539 
540  end subroutine atmosphere_init
541 
542 
545  subroutine p_adi(km, ng, ifirst, ilast, jfirst, jlast, ptop, &
546  delp, pt, ps, pe, peln, pk, pkz, hydrostatic)
547 ! Given (ptop, delp) computes (ps, pk, pe, peln, pkz)
548 ! Input:
549  integer, intent(in):: km, ng
550  integer, intent(in):: ifirst, ilast
551  integer, intent(in):: jfirst, jlast
552  logical, intent(in):: hydrostatic
553  real, intent(in):: ptop
554  real, intent(in):: pt(ifirst-ng:ilast+ng,jfirst-ng:jlast+ng, km)
555  real, intent(in):: delp(ifirst-ng:ilast+ng,jfirst-ng:jlast+ng, km)
556 ! Output:
557  real, intent(out) :: ps(ifirst-ng:ilast+ng, jfirst-ng:jlast+ng)
558  real, intent(out) :: pk(ifirst:ilast, jfirst:jlast, km+1)
559  real, intent(out) :: pe(ifirst-1:ilast+1,km+1,jfirst-1:jlast+1)
560  real, intent(out) :: peln(ifirst:ilast, km+1, jfirst:jlast)
561  real, intent(out) :: pkz(ifirst:ilast, jfirst:jlast, km)
562 ! Local
563  real pek
564  integer i, j, k
565 
566  pek = ptop ** kappa
567 !$OMP parallel do default (none) &
568 !$OMP shared (ifirst,ilast,jfirst,jlast,km,ptop,pek,pe,pk, &
569 !$OMP ps,delp,peln,hydrostatic,pkz) &
570 !$OMP private (j, i, k)
571  do j=jfirst,jlast
572  do i=ifirst,ilast
573  pe(i,1,j) = ptop
574  pk(i,j,1) = pek
575  enddo
576 
577  do k=2,km+1
578  do i=ifirst,ilast
579  pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
580  peln(i,k,j) = log(pe(i,k,j))
581  pk(i,j,k) = exp( kappa*peln(i,k,j) )
582  enddo
583  enddo
584 
585  do i=ifirst,ilast
586  ps(i,j) = pe(i,km+1,j)
587  enddo
588 
589  if ( hydrostatic ) then
590  do k=1,km
591  do i=ifirst,ilast
592  pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
593  enddo
594  enddo
595  endif
596  enddo
597 
598  end subroutine p_adi
599 
600 
603  subroutine atmosphere_dynamics ( Time )
604  type(time_type),intent(in) :: Time
605  integer :: n, psc, atmos_time_step
606  integer :: k, w_diff, nt_dyn, n_split_loc, seconds, days
607 
608  type(time_type) :: atmos_time
609 
610 !---- Call FV dynamics -----
611 
612  call mpp_clock_begin (id_dynam)
613 
614  n = mytile
615 
616  call get_time (time, seconds, days)
617 ! if (seconds < 10800 .and. days == 0) then
618 ! n_split_loc = (Atm(n)%flagstruct%n_split * 3) / 2
619  if (seconds < nint(3600*atm(n)%flagstruct%fhouri) .and. atm(n)%flagstruct%fac_n_spl > 1.0) then
620  n_split_loc = nint(atm(n)%flagstruct%n_split * atm(n)%flagstruct%fac_n_spl)
621  else
622  n_split_loc = atm(n)%flagstruct%n_split
623  endif
624 
625 ! write(0,*)' before calling init n_split_loc=',n_split_loc,' seconds=',seconds,' days=',days,&
626 ! ' n_split=',Atm(mytile)%flagstruct%n_split,' mytile=',mytile
627 
628  a_step = a_step + 1
629 !
630 !*** If this is a regional run then read in the next boundary data when it is time.
631 !
632  if(atm(n)%flagstruct%regional)then
633 
634  call read_new_bc_data(atm(n), time, time_step_atmos, p_split, &
635  isd, ied, jsd, jed )
636  endif
637 
638  do psc=1,abs(p_split)
639  p_step = psc
640  call timing_on('fv_dynamics')
641 !uc/vc only need be same on coarse grid? However BCs do need to be the same
642  call fv_dynamics(npx, npy, npz, nq, atm(n)%ng, dt_atmos/real(abs(p_split)),&
643  Atm(n)%flagstruct%consv_te, Atm(n)%flagstruct%fill, &
644  Atm(n)%flagstruct%reproduce_sum, kappa, cp_air, zvir, &
645  Atm(n)%ptop, Atm(n)%ks, nq, &
646  n_split_loc, Atm(n)%flagstruct%q_split, &
647 ! Atm(n)%flagstruct%n_split, Atm(n)%flagstruct%q_split, &
648  atm(n)%u, atm(n)%v, atm(n)%w, atm(n)%delz, &
649  atm(n)%flagstruct%hydrostatic, &
650  atm(n)%pt , atm(n)%delp, atm(n)%q, atm(n)%ps, &
651  atm(n)%pe, atm(n)%pk, atm(n)%peln, &
652  atm(n)%pkz, atm(n)%phis, atm(n)%q_con, &
653  atm(n)%omga, atm(n)%ua, atm(n)%va, atm(n)%uc, &
654  atm(n)%vc, atm(n)%ak, atm(n)%bk, atm(n)%mfx, &
655  atm(n)%mfy , atm(n)%cx, atm(n)%cy, atm(n)%ze0, &
656  atm(n)%flagstruct%hybrid_z, &
657  atm(n)%gridstruct, atm(n)%flagstruct, &
658  atm(n)%neststruct, atm(n)%idiag, atm(n)%bd, &
659  atm(n)%parent_grid, atm(n)%domain,atm(n)%diss_est)
660 
661  call timing_off('fv_dynamics')
662 
663  if (ngrids > 1 .and. (psc < p_split .or. p_split < 0)) then
664  call timing_on('TWOWAY_UPDATE')
665  call twoway_nesting(atm, ngrids, grids_on_this_pe, zvir)
666  call timing_off('TWOWAY_UPDATE')
667  endif
668 
669  end do !p_split
670  call mpp_clock_end (id_dynam)
671 
672 !-----------------------------------------------------
673 !--- COMPUTE SUBGRID Z
674 !-----------------------------------------------------
675 !--- zero out tendencies
676  call mpp_clock_begin (id_subgridz)
677  u_dt(:,:,:) = 0.
678  v_dt(:,:,:) = 0.
679  t_dt(:,:,:) = 0.
680 
681  w_diff = get_tracer_index(model_atmos, 'w_diff' )
682  if ( atm(n)%flagstruct%fv_sg_adj > 0 ) then
683  nt_dyn = nq
684  if ( w_diff /= no_tracer ) then
685  nt_dyn = nq - 1
686  endif
687  call fv_subgrid_z(isd, ied, jsd, jed, isc, iec, jsc, jec, atm(n)%npz, &
688  nt_dyn, dt_atmos, atm(n)%flagstruct%fv_sg_adj, &
689  atm(n)%flagstruct%nwat, atm(n)%delp, atm(n)%pe, &
690  atm(n)%peln, atm(n)%pkz, atm(n)%pt, atm(n)%q, &
691  atm(n)%ua, atm(n)%va, atm(n)%flagstruct%hydrostatic,&
692  atm(n)%w, atm(n)%delz, u_dt, v_dt, t_dt, &
693  atm(n)%flagstruct%n_sponge)
694  endif
695 
696 #ifdef USE_Q_DT
697  if ( .not. atm(n)%flagstruct%hydrostatic .and. w_diff /= no_tracer ) then
698 !$OMP parallel do default (none) &
699 !$OMP shared (isc, iec, jsc, jec, w_diff, n, Atm, q_dt) &
700 !$OMP private (k)
701  do k=1, atm(n)%npz
702  atm(n)%q(isc:iec,jsc:jec,k,w_diff) = atm(n)%w(isc:iec,jsc:jec,k) + w0_big
703  q_dt(:,:,k,w_diff) = 0.
704  enddo
705  endif
706 #endif
707 
708  call mpp_clock_end (id_subgridz)
709 
710  end subroutine atmosphere_dynamics
711 
712 
715  subroutine atmosphere_end (Time, Grid_box)
716 #ifdef CCPP
717 #ifdef STATIC
718 ! For static builds, the ccpp_physics_{init,run,finalize} calls
719 ! are not pointing to code in the CCPP framework, but to auto-generated
720 ! ccpp_suite_cap and ccpp_group_*_cap modules behind a ccpp_static_api
721  use ccpp_static_api, only: ccpp_physics_finalize
722  use ccpp_data, only: ccpp_suite
723 #else
724  use ccpp_api, only: ccpp_physics_finalize
725 #endif
726  use ccpp_data, only: cdata => cdata_tile
727 #endif
728  type(time_type), intent(in) :: time
729  type(grid_box_type), intent(inout) :: Grid_box
730 
731 #ifdef CCPP
732  integer :: ierr
733 
734  if (atm(mytile)%flagstruct%do_sat_adj) then
735  ! Finalize fast physics
736 #ifdef STATIC
737  call ccpp_physics_finalize(cdata, suite_name=trim(ccpp_suite), group_name="fast_physics", ierr=ierr)
738 #else
739  call ccpp_physics_finalize(cdata, group_name="fast_physics", ierr=ierr)
740 #endif
741  if (ierr/=0) then
742  cdata%errmsg = ' atmosphere_dynamics: error in ccpp_physics_finalize for group fast_physics: ' // trim(cdata%errmsg)
743  call mpp_error (fatal, cdata%errmsg)
744  end if
745  end if
746 #endif
747 
748  call nullify_domain ( )
749  if (first_diag) then
750  call timing_on('FV_DIAG')
751  call fv_diag(atm(mytile:mytile), zvir, fv_time, atm(mytile)%flagstruct%print_freq)
753  first_diag = .false.
754  call timing_off('FV_DIAG')
755  endif
756 
758  deallocate (atm)
759 
760  deallocate( u_dt, v_dt, t_dt, pref, dum1d )
761 
762  end subroutine atmosphere_end
763 
764 
769  subroutine atmosphere_restart(timestamp)
770  character(len=*), intent(in) :: timestamp
771 
772  call fv_write_restart(atm, grids_on_this_pe, timestamp)
773 
774  end subroutine atmosphere_restart
775 
776 
780  subroutine atmosphere_resolution (i_size, j_size, global)
781  integer, intent(out) :: i_size, j_size
782  logical, intent(in), optional :: global
783  logical :: local
784 
785  local = .true.
786  if( PRESENT(global) ) local = .NOT.global
787 
788  if( local ) then
789  i_size = iec - isc + 1
790  j_size = jec - jsc + 1
791  else
792  i_size = npx - 1
793  j_size = npy - 1
794  end if
795 
796  end subroutine atmosphere_resolution
797 
798 
800  subroutine atmosphere_pref (p_ref)
801  real, dimension(:,:), intent(inout) :: p_ref
802 
803  p_ref = pref
804 
805  end subroutine atmosphere_pref
806 
807 
808  subroutine atmosphere_control_data (i1, i2, j1, j2, kt, p_hydro, hydro, tile_num)
809  integer, intent(out) :: i1, i2, j1, j2, kt
810  logical, intent(out), optional :: p_hydro, hydro
811  integer, intent(out), optional :: tile_num
812  i1 = atm(mytile)%bd%isc
813  i2 = atm(mytile)%bd%iec
814  j1 = atm(mytile)%bd%jsc
815  j2 = atm(mytile)%bd%jec
816  kt = atm(mytile)%npz
817 
818  if (present(tile_num)) tile_num = atm(mytile)%tile
819  if (present(p_hydro)) p_hydro = atm(mytile)%flagstruct%phys_hydrostatic
820  if (present( hydro)) hydro = atm(mytile)%flagstruct%hydrostatic
821 
822  end subroutine atmosphere_control_data
823 
824 
827  subroutine atmosphere_grid_ctr (lon, lat)
828  real(kind=kind_phys), intent(out) :: lon(:,:), lat(:,:)
829 ! Local data:
830  integer i,j
831 
832  do j=jsc,jec
833  do i=isc,iec
834  lon(i-isc+1,j-jsc+1) = atm(mytile)%gridstruct%agrid_64(i,j,1)
835  lat(i-isc+1,j-jsc+1) = atm(mytile)%gridstruct%agrid_64(i,j,2)
836  enddo
837  end do
838 
839  end subroutine atmosphere_grid_ctr
840 
841 
844  subroutine atmosphere_grid_bdry (blon, blat, global)
845  real, intent(out) :: blon(:,:), blat(:,:)
846  logical, intent(in), optional :: global
847 ! Local data:
848  integer i,j
849 
850  if( PRESENT(global) ) then
851  if (global) call mpp_error(fatal, '==> global grid is no longer available &
852  & in the Cubed Sphere')
853  endif
854 
855  do j=jsc,jec+1
856  do i=isc,iec+1
857  blon(i-isc+1,j-jsc+1) = atm(mytile)%gridstruct%grid(i,j,1)
858  blat(i-isc+1,j-jsc+1) = atm(mytile)%gridstruct%grid(i,j,2)
859  enddo
860  end do
861 
862  end subroutine atmosphere_grid_bdry
863 
864 
865  subroutine set_atmosphere_pelist ()
866  call mpp_set_current_pelist(atm(mytile)%pelist, no_sync=.true.)
867  end subroutine set_atmosphere_pelist
868 
869 
874  subroutine atmosphere_domain ( fv_domain, layout, regional, nested, pelist )
875  type(domain2d), intent(out) :: fv_domain
876  integer, intent(out) :: layout(2)
877  logical, intent(out) :: regional
878  logical, intent(out) :: nested
879  integer, pointer, intent(out) :: pelist(:)
880 ! returns the domain2d variable associated with the coupling grid
881 ! note: coupling is done using the mass/temperature grid with no halos
882 
883  fv_domain = atm(mytile)%domain_for_coupler
884  layout(1:2) = atm(mytile)%layout(1:2)
885  regional = atm(mytile)%flagstruct%regional
886  nested = ngrids > 1
887  call set_atmosphere_pelist()
888  pelist => atm(mytile)%pelist
889 
890  end subroutine atmosphere_domain
891 
892 
895  subroutine atmosphere_diag_axes ( axes )
896  integer, intent(out) :: axes (:)
897 
898 !----- returns the axis indices for the atmospheric (mass) grid -----
899  if ( size(axes(:)) < 0 .or. size(axes(:)) > 4 ) call error_mesg ( &
900  'get_atmosphere_axes in atmosphere_mod', &
901  'size of argument is incorrect', fatal )
902 
903  axes(1:size(axes(:))) = atm(mytile)%atmos_axes (1:size(axes(:)))
904 
905  end subroutine atmosphere_diag_axes
906 
907 
912  subroutine atmosphere_etalvls (ak, bk, flip)
913  real(kind=kind_phys), pointer, dimension(:), intent(inout) :: ak, bk
914  logical, intent(in) :: flip
915 
916  allocate(ak(npz+1))
917  allocate(bk(npz+1))
918 
919  if (flip) then
920  ak(1:npz+1) = atm(mytile)%ak(npz+1:1:-1)
921  bk(1:npz+1) = atm(mytile)%bk(npz+1:1:-1)
922  else
923  ak(1:npz+1) = atm(mytile)%ak(1:npz+1)
924  bk(1:npz+1) = atm(mytile)%bk(1:npz+1)
925  endif
926  end subroutine atmosphere_etalvls
927 
928 
934  subroutine atmosphere_hgt (hgt, position, relative, flip)
935  real(kind=kind_phys), pointer, dimension(:,:,:), intent(inout) :: hgt
936  character(len=5), intent(in) :: position
937  logical, intent(in) :: relative
938  logical, intent(in) :: flip
939  !--- local variables ---
940  integer:: lev, k, j, i, npx, npy
941  real(kind=kind_phys), dimension(isc:iec,jsc:jec,1:npz+1) :: z
942  real(kind=kind_phys), dimension(isc:iec,jsc:jec,1:npz) :: dz
943 
944  if ((position .ne. "layer") .and. (position .ne. "level")) then
945  call mpp_error (fatal, 'atmosphere_hgt:: incorrect position specification')
946  endif
947 
948  npx = iec-isc+1
949  npy = jec-jsc+1
950  z = 0
951  dz = 0
952 
953  if (atm(mytile)%flagstruct%hydrostatic) then
954  !--- generate dz using hydrostatic assumption
955  dz(isc:iec,jsc:jec,1:npz) = (rdgas/grav)*atm(mytile)%pt(isc:iec,jsc:jec,1:npz) &
956  * (atm(mytile)%peln(isc:iec,1:npz,jsc:jec) &
957  - atm(mytile)%peln(isc:iec,2:npz+1,jsc:jec))
958  else
959  !--- use non-hydrostatic delz directly
960  dz(isc:iec,jsc:jec,1:npz) = atm(mytile)%delz(isc:iec,jsc:jec,1:npz)
961  endif
962 
963  !--- calculate geometric heights at the interfaces (levels)
964  !--- if needed, flip the indexing during this step
965  if (flip) then
966  if (.not. relative) z(isc:iec,jsc:jec,1) = atm(mytile)%phis(isc:iec,jsc:jec)/grav
967  do k = 2,npz+1
968  z(isc:iec,jsc:jec,k) = z(isc:iec,jsc:jec,k-1) - dz(isc:iec,jsc:jec,npz+2-k)
969  enddo
970  else
971  if (.not. relative) z(isc:iec,jsc:jec,npz+1) = atm(mytile)%phis(isc:iec,jsc:jec)/grav
972  do k = npz,1,-1
973  z(isc:iec,jsc:jec,k) = z(isc:iec,jsc:jec,k+1) - dz(isc:iec,jsc:jec,k)
974  enddo
975  endif
976 
977  !--- allocate and set either the level or layer height for return
978  if (position == "level") then
979  allocate (hgt(npx,npy,npz+1))
980  hgt(1:npx,1:npy,1:npz+1) = z(isc:iec,jsc:jec,1:npz+1)
981  elseif (position == "layer") then
982  allocate (hgt(npx,npy,npz))
983  hgt(1:npx,1:npy,1:npz) = 0.5d0 * (z(isc:iec,jsc:jec,1:npz) + z(isc:iec,jsc:jec,2:npz+1))
984  endif
985 
986  end subroutine atmosphere_hgt
987 
988 
994  subroutine atmosphere_scalar_field_halo (data, halo, isize, jsize, ksize, data_p)
995  !--------------------------------------------------------------------
996  ! data - output array to return the field with halo (i,j,k)
997  ! optionally input for field already in (i,j,k) form
998  ! sized to include the halo of the field (+ 2*halo)
999  ! halo - size of the halo (must be less than 3)
1000  ! ied - horizontal resolution in i-dir with haloes
1001  ! jed - horizontal resolution in j-dir with haloes
1002  ! ksize - vertical resolution
1003  ! data_p - optional input field in packed format (ix,k)
1004  !--------------------------------------------------------------------
1005  !--- interface variables ---
1006  real(kind=kind_phys), dimension(1:isize,1:jsize,ksize), intent(inout) :: data
1009  integer, intent(in) :: halo
1010  integer, intent(in) :: isize
1011  integer, intent(in) :: jsize
1012  integer, intent(in) :: ksize
1013  real(kind=kind_phys), dimension(:,:), optional, intent(in) :: data_p
1014  !--- local variables ---
1015  integer :: i, j, k
1016  integer :: ic, jc
1017  character(len=44) :: modname = 'atmosphere_mod::atmosphere_scalar_field_halo'
1018  integer :: mpp_flags
1019 
1020  !--- perform error checking
1021  if (halo .gt. 3) call mpp_error(fatal, modname//.gt.' - halo3 requires extending the MPP domain to support')
1022  ic = isize - 2 * halo
1023  jc = jsize - 2 * halo
1024 
1025  !--- if packed data is present, unpack it into the two-dimensional data array
1026  if (present(data_p)) then
1027  if (ic*jc .ne. size(data_p,1)) call mpp_error(fatal, modname//' - incorrect sizes for incoming &
1028  &variables data and data_p')
1029  data = 0.
1030 !$OMP parallel do default (none) &
1031 !$OMP shared (data, data_p, halo, ic, jc, ksize) &
1032 !$OMP private (i, j, k)
1033  do k = 1, ksize
1034  do j = 1, jc
1035  do i = 1, ic
1036  data(i+halo, j+halo, k) = data_p(i + (j-1)*ic, k)
1037  enddo
1038  enddo
1039  enddo
1040  endif
1041 
1042  mpp_flags = eupdate + wupdate + supdate + nupdate
1043  if (halo == 1) then
1044  call mpp_update_domains(data, atm(mytile)%domain_for_coupler, flags=mpp_flags, complete=.true.)
1045  elseif (halo == 3) then
1046  call mpp_update_domains(data, atm(mytile)%domain, flags=mpp_flags, complete=.true.)
1047  else
1048  call mpp_error(fatal, modname//' - unsupported halo size')
1049  endif
1050 
1051  !--- fill the halo points when at a corner of the cubed-sphere tile
1052  !--- interior domain corners are handled correctly
1053  if ( (isc==1) .or. (jsc==1) .or. (iec==npx-1) .or. (jec==npy-1) ) then
1054  do k = 1, ksize
1055  do j=1,halo
1056  do i=1,halo
1057  if ((isc== 1) .and. (jsc== 1)) data(halo+1-j ,halo+1-i ,k) = data(halo+i ,halo+1-j ,k) !SW Corner
1058  if ((isc== 1) .and. (jec==npy-1)) data(halo+1-j ,halo+jc+i,k) = data(halo+i ,halo+jc+j,k) !NW Corner
1059  if ((iec==npx-1) .and. (jsc== 1)) data(halo+ic+j,halo+1-i ,k) = data(halo+ic-i+1,halo+1-j ,k) !SE Corner
1060  if ((iec==npx-1) .and. (jec==npy-1)) data(halo+ic+j,halo+jc+i,k) = data(halo+ic-i+1,halo+jc+j,k) !NE Corner
1061  enddo
1062  enddo
1063  enddo
1064  endif
1065 
1066  return
1067  end subroutine atmosphere_scalar_field_halo
1068 
1069 
1070  subroutine atmosphere_diss_est (npass)
1072  !--- interface variables ---
1073  integer, intent(in) :: npass
1074  !--- local variables
1075  integer:: k
1076 
1077  !horizontally smooth dissiapation estimate for SKEB
1078  ! 3 passes before taking absolute value
1079  do k = 1,min(3,npass)
1080  call del2_cubed(atm(mytile)%diss_est, 0.25*atm(mytile)%gridstruct%da_min, atm(mytile)%gridstruct, &
1081  atm(mytile)%domain, npx, npy, npz, 3, atm(mytile)%bd)
1082  enddo
1083 
1084  atm(mytile)%diss_est=abs(atm(mytile)%diss_est)
1085 
1086  do k = 4,npass
1087  call del2_cubed(atm(mytile)%diss_est, 0.25*atm(mytile)%gridstruct%da_min, atm(mytile)%gridstruct, &
1088  atm(mytile)%domain, npx, npy, npz, 3, atm(mytile)%bd)
1089  enddo
1090  ! provide back sqrt of dissipation estimate
1091  atm(mytile)%diss_est=sqrt(abs(atm(mytile)%diss_est))
1092 
1093  end subroutine atmosphere_diss_est
1094 
1095 
1100  subroutine atmosphere_nggps_diag (Time, init, ltavg,avg_max_length)
1101  type(time_type), intent(in) :: Time
1102  logical, optional, intent(in) :: init, ltavg
1103  real, optional, intent(in) :: avg_max_length
1104  if (PRESENT(init)) then
1105  if (init) then
1106  call fv_nggps_diag_init(atm(mytile:mytile), atm(mytile)%atmos_axes, time)
1107  return
1108  else
1109  call mpp_error(fatal, 'atmosphere_nggps_diag - calling with init present, but set to .false.')
1110  endif
1111  endif
1112  if (PRESENT(ltavg)) then
1113  if (ltavg) then
1114  call fv_nggps_tavg(atm(mytile:mytile), time_step_atmos,avg_max_length,zvir)
1115  return
1116  endif
1117  else
1118  call fv_nggps_diag(atm(mytile:mytile), zvir, time)
1119  endif
1120  end subroutine atmosphere_nggps_diag
1121 
1122 
1123 !--- Need to know the formulation of the mixing ratio being imported into FV3
1124 !--- in order to adjust it in a consistent manner for advection
1125 !rab subroutine atmosphere_tracer_postinit (IPD_Data, Atm_block)
1126 !rab !--- interface variables ---
1127 !rab type(IPD_data_type), intent(in) :: IPD_Data(:)
1128 !rab type(block_control_type), intent(in) :: Atm_block
1129 !rab !--- local variables ---
1130 !rab integer :: i, j, ix, k, k1, n, nwat, nb, blen
1131 !rab real(kind=kind_phys) :: qwat
1132 !rab
1133 !rab if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers')
1134 !rab
1135 !rab n = mytile
1136 !rab nwat = Atm(n)%flagstruct%nwat
1137 !rab
1138 !rab!$OMP parallel do default (none) &
1139 !rab!$OMP shared (Atm_block, Atm, IPD_Data, npz, nq, ncnst, n, nwat) &
1140 !rab!$OMP private (nb, k, k1, ix, i, j, qwat)
1141 !rab do nb = 1,Atm_block%nblks
1142 !rab do k = 1, npz
1143 !rab k1 = npz+1-k !reverse the k direction
1144 !rab do ix = 1, Atm_block%blksz(nb)
1145 !rab i = Atm_block%index(nb)%ii(ix)
1146 !rab j = Atm_block%index(nb)%jj(ix)
1147 !rab qwat = sum(Atm(n)%q(i,j,k1,1:nwat))
1148 !rab Atm(n)%q(i,j,k1,1:nq) = Atm(n)%q(i,j,k1,1:nq) + IPD_Data(nb)%Stateout%gq0(ix,k,1:nq) * (1.0 - qwat)
1149 !rab if (nq .gt. ncnst) then
1150 !rab Atm(n)%qdiag(i,j,k1,nq+1:ncnst) = Atm(n)%qdiag(i,j,k1,nq+1:ncnst) + IPD_Data(nb)%Stateout%gq0(ix,k,nq+1:ncnst)
1151 !rab endif
1152 !rab enddo
1153 !rab enddo
1154 !rab enddo
1155 !rab
1156 !rab call mpp_update_domains (Atm(n)%q, Atm(n)%domain, complete=.true.)
1157 !rab
1158 !rab return
1159 !rab end subroutine atmosphere_tracer_postinit
1160 
1161 
1162  subroutine get_bottom_mass ( t_bot, tr_bot, p_bot, z_bot, p_surf, slp )
1163 !--------------------------------------------------------------
1164 ! returns temp, sphum, pres, height at the lowest model level
1165 ! and surface pressure
1166 !--------------------------------------------------------------
1167  !--- interface variables ---
1168  real, intent(out), dimension(isc:iec,jsc:jec):: t_bot, p_bot, z_bot, p_surf
1169  real, intent(out), optional, dimension(isc:iec,jsc:jec):: slp
1170  real, intent(out), dimension(isc:iec,jsc:jec,nq):: tr_bot
1171  !--- local variables ---
1172  integer :: i, j, m, k, kr
1173  real :: rrg, sigtop, sigbot
1174  real, dimension(isc:iec,jsc:jec) :: tref
1175  real, parameter :: tlaps = 6.5e-3
1176 
1177  rrg = rdgas / grav
1178 
1179  do j=jsc,jec
1180  do i=isc,iec
1181  p_surf(i,j) = atm(mytile)%ps(i,j)
1182  t_bot(i,j) = atm(mytile)%pt(i,j,npz)
1183  p_bot(i,j) = atm(mytile)%delp(i,j,npz)/(atm(mytile)%peln(i,npz+1,j)-atm(mytile)%peln(i,npz,j))
1184  z_bot(i,j) = rrg*t_bot(i,j)*(1. - atm(mytile)%pe(i,npz,j)/p_bot(i,j))
1185 #ifdef MULTI_GASES
1186  z_bot(i,j) = z_bot(i,j)*virq(atm(mytile)%q(i,j,npz,:))
1187 #else
1188  z_bot(i,j) = z_bot(i,j)*(1.+zvir*atm(mytile)%q(i,j,npz,1))
1189 #endif
1190  enddo
1191  enddo
1192 
1193  if ( present(slp) ) then
1194  ! determine 0.8 sigma reference level
1195  sigtop = atm(mytile)%ak(1)/pstd_mks+atm(mytile)%bk(1)
1196  do k = 1, npz
1197  sigbot = atm(mytile)%ak(k+1)/pstd_mks+atm(mytile)%bk(k+1)
1198  if (sigbot+sigtop > 1.6) then
1199  kr = k
1200  exit
1201  endif
1202  sigtop = sigbot
1203  enddo
1204  do j=jsc,jec
1205  do i=isc,iec
1206  ! sea level pressure
1207  tref(i,j) = atm(mytile)%pt(i,j,kr) * (atm(mytile)%delp(i,j,kr)/ &
1208  ((atm(mytile)%peln(i,kr+1,j)-atm(mytile)%peln(i,kr,j))*atm(mytile)%ps(i,j)))**(-rrg*tlaps)
1209  slp(i,j) = atm(mytile)%ps(i,j)*(1.+tlaps*atm(mytile)%phis(i,j)/(tref(i,j)*grav))**(1./(rrg*tlaps))
1210  enddo
1211  enddo
1212  endif
1213 
1214 ! Copy tracers
1215  do m=1,nq
1216  do j=jsc,jec
1217  do i=isc,iec
1218  tr_bot(i,j,m) = atm(mytile)%q(i,j,npz,m)
1219  enddo
1220  enddo
1221  enddo
1222 
1223  end subroutine get_bottom_mass
1224 
1225 
1226  subroutine get_bottom_wind ( u_bot, v_bot )
1227 !-----------------------------------------------------------
1228 ! returns u and v on the mass grid at the lowest model level
1229 !-----------------------------------------------------------
1230  !--- interface variables ---
1231  real, intent(out), dimension(isc:iec,jsc:jec):: u_bot, v_bot
1232  !--- local variables ---
1233  integer i, j
1234 
1235  do j=jsc,jec
1236  do i=isc,iec
1237  u_bot(i,j) = atm(mytile)%u_srf(i,j)
1238  v_bot(i,j) = atm(mytile)%v_srf(i,j)
1239  enddo
1240  enddo
1241 
1242  end subroutine get_bottom_wind
1243 
1244 
1248  subroutine atmosphere_get_bottom_layer (Atm_block, DYCORE_Data)
1249 !--------------------------------------------------------------
1250 ! Returns in blocked format:
1251 ! temp, pres, winds, height, and tracers at the bottom layer
1252 ! surface pressure
1253 ! sea-level pressure
1254 !--------------------------------------------------------------
1255  !--- interface variables ---
1256  type(block_control_type), intent(in) :: Atm_block
1257  type(dycore_data_type), intent(inout) :: DYCORE_Data(:)
1258  !--- local variables ---
1259  integer :: ix, i, j, nt, k, nb
1260  real :: rrg, sigtop, sigbot, tref
1261  !--- local parameters ---
1262  real, parameter :: tlaps = 6.5e-3
1263  !--- local saved variables ---
1264  integer, save :: kr
1265  logical, save :: first_time = .true.
1266 
1267  rrg = rdgas / grav
1268 
1269  if (first_time) then
1270  print *, 'calculating slp kr value'
1271  ! determine 0.8 sigma reference level
1272  sigtop = atm(mytile)%ak(1)/pstd_mks+atm(mytile)%bk(1)
1273  do k = 1, npz
1274  sigbot = atm(mytile)%ak(k+1)/pstd_mks+atm(mytile)%bk(k+1)
1275  if (sigbot+sigtop > 1.6) then
1276  kr = k
1277  exit
1278  endif
1279  sigtop = sigbot
1280  enddo
1281  !--- allocate the DYCORE_Data container
1282  do nb = 1,atm_block%nblks
1283  call dycore_data(nb)%Coupling%create (atm_block%blksz(nb), nq)
1284  enddo
1285  first_time = .false.
1286  endif
1287 
1288 !$OMP parallel do default (none) &
1289 !$OMP shared (Atm_block, DYCORE_Data, Atm, mytile, npz, kr, rrg, zvir, nq) &
1290 !$OMP private (nb, ix, i, j, tref, nt)
1291  do nb = 1,atm_block%nblks
1292  do ix = 1,atm_block%blksz(nb)
1293  i = atm_block%index(nb)%ii(ix)
1294  j = atm_block%index(nb)%jj(ix)
1295  !--- surface pressure
1296  dycore_data(nb)%Coupling%p_srf(ix) = atm(mytile)%ps(i,j)
1297  !--- bottom layer temperature, pressure, & winds
1298  dycore_data(nb)%Coupling%t_bot(ix) = atm(mytile)%pt(i,j,npz)
1299  dycore_data(nb)%Coupling%p_bot(ix) = atm(mytile)%delp(i,j,npz)/(atm(mytile)%peln(i,npz+1,j)-atm(mytile)%peln(i,npz,j))
1300  dycore_data(nb)%Coupling%u_bot(ix) = atm(mytile)%u_srf(i,j)
1301  dycore_data(nb)%Coupling%v_bot(ix) = atm(mytile)%v_srf(i,j)
1302  !--- bottom layer height based on hydrostatic assumptions
1303  dycore_data(nb)%Coupling%z_bot(ix) = rrg*dycore_data(nb)%Coupling%t_bot(ix) * &
1304  (1. - atm(mytile)%pe(i,npz,j)/dycore_data(nb)%Coupling%p_bot(ix))
1305 #ifdef MULTI_GASES
1306  dycore_data(nb)%Coupling%z_bot(ix) = dycore_data(nb)%Coupling%z_bot(ix)*virq(atm(mytile)%q(i,j,npz,:))
1307 #else
1308  dycore_data(nb)%Coupling%z_bot(ix) = dycore_data(nb)%Coupling%z_bot(ix)*(1.+zvir*atm(mytile)%q(i,j,npz,1))
1309 #endif
1310  !--- sea level pressure
1311  tref = atm(mytile)%pt(i,j,kr) * (atm(mytile)%delp(i,j,kr)/ &
1312  ((atm(mytile)%peln(i,kr+1,j)-atm(mytile)%peln(i,kr,j))*atm(mytile)%ps(i,j)))**(-rrg*tlaps)
1313  dycore_data(nb)%Coupling%slp(ix) = atm(mytile)%ps(i,j)*(1.+tlaps*atm(mytile)%phis(i,j)/(tref*grav))**(1./(rrg*tlaps))
1314  enddo
1315 
1316  !--- bottom layer tracers
1317  do nt = 1,nq
1318  do ix = 1,atm_block%blksz(nb)
1319  i = atm_block%index(nb)%ii(ix)
1320  j = atm_block%index(nb)%jj(ix)
1321  dycore_data(nb)%Coupling%tr_bot(ix,nt) = atm(mytile)%q(i,j,npz,nt)
1322  enddo
1323  enddo
1324  enddo
1325 
1326  end subroutine atmosphere_get_bottom_layer
1327 
1328 
1329  subroutine get_stock_pe(index, value)
1330  integer, intent(in) :: index
1331  real, intent(out) :: value
1332 
1333 #ifdef USE_STOCK
1334  include 'stock.inc'
1335 #endif
1336 
1337  real wm(isc:iec,jsc:jec)
1338  integer i,j,k
1339  real, pointer :: area(:,:)
1340 
1341  area => atm(mytile)%gridstruct%area
1342 
1343  select case (index)
1344 
1345 #ifdef USE_STOCK
1346  case (istock_water)
1347 #else
1348  case (1)
1349 #endif
1350 
1351 !----------------------
1352 ! Perform vertical sum:
1353 !----------------------
1354  wm = 0.
1355  do j=jsc,jec
1356  do k=1,npz
1357  do i=isc,iec
1358 ! Warning: the following works only with AM2 physics: water vapor; cloud water, cloud ice.
1359  wm(i,j) = wm(i,j) + atm(mytile)%delp(i,j,k) * ( atm(mytile)%q(i,j,k,1) + &
1360  atm(mytile)%q(i,j,k,2) + &
1361  atm(mytile)%q(i,j,k,3) )
1362  enddo
1363  enddo
1364  enddo
1365 
1366 !----------------------
1367 ! Horizontal sum:
1368 !----------------------
1369  value = 0.
1370  do j=jsc,jec
1371  do i=isc,iec
1372  value = value + wm(i,j)*area(i,j)
1373  enddo
1374  enddo
1375  value = value/grav
1376 
1377  case default
1378  value = 0.0
1379  end select
1380 
1381  end subroutine get_stock_pe
1382 
1383 
1386  subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc)
1387  type(time_type), intent(in) :: Time
1388  type(ipd_data_type), intent(in) :: IPD_Data(:)
1389  type(iau_external_data_type), intent(in) :: IAU_Data
1390  type(block_control_type), intent(in) :: Atm_block
1391  logical, intent(in) :: flip_vc
1392  !--- local variables ---
1393  type(time_type) :: Time_prev, Time_next
1394  integer :: i, j, ix, k, k1, n, w_diff, nt_dyn, iq
1395  integer :: nb, blen, nwat, dnats, nq_adv
1396  real(kind=kind_phys):: rcp, q0, qwat(nq), qt, rdt
1397  time_prev = time
1398  time_next = time + time_step_atmos
1399  rdt = 1.d0 / dt_atmos
1400 
1401  n = mytile
1402  nwat = atm(n)%flagstruct%nwat
1403  dnats = atm(mytile)%flagstruct%dnats
1404  nq_adv = nq - dnats
1405 
1406  if( nq<3 ) call mpp_error(fatal, 'GFS phys must have 3 interactive tracers')
1407 
1408  if (iau_data%in_interval) then
1409 ! IAU increments are in units of 1/sec
1410 
1411 ! add analysis increment to u,v,t tendencies
1412 ! directly update delp with analysis increment
1413  do k = 1, npz
1414  do j = jsc,jec
1415  do i = isc,iec
1416  u_dt(i,j,k) = u_dt(i,j,k) + iau_data%ua_inc(i,j,k)
1417  v_dt(i,j,k) = v_dt(i,j,k) + iau_data%va_inc(i,j,k)
1418  t_dt(i,j,k) = t_dt(i,j,k) + iau_data%temp_inc(i,j,k)
1419  atm(n)%delp(i,j,k) = atm(n)%delp(i,j,k) + iau_data%delp_inc(i,j,k)*dt_atmos
1420  enddo
1421  enddo
1422  enddo
1423  if (.not. atm(mytile)%flagstruct%hydrostatic) then
1424  do k = 1, npz
1425  do j = jsc,jec
1426  do i = isc,iec
1427  atm(n)%delz(i,j,k) = atm(n)%delz(i,j,k) + iau_data%delz_inc(i,j,k)*dt_atmos
1428  enddo
1429  enddo
1430  enddo
1431  endif
1432 ! add analysis increment to tracers to output from physics
1433  do nb = 1,atm_block%nblks
1434  !if (nb.EQ.1) print*,'in block_update',IAU_Data%in_interval,IAU_Data%temp_inc(isc,jsc,30)
1435  blen = atm_block%blksz(nb)
1436  do k = 1, npz
1437  if(flip_vc) then
1438  k1 = npz+1-k !reverse the k direction
1439  else
1440  k1 = k
1441  endif
1442  do ix = 1, blen
1443  i = atm_block%index(nb)%ii(ix)
1444  j = atm_block%index(nb)%jj(ix)
1445  ipd_data(nb)%Stateout%gq0(ix,k,:) = ipd_data(nb)%Stateout%gq0(ix,k,:) + iau_data%tracer_inc(i,j,k1,:)*dt_atmos
1446  enddo
1447  enddo
1448  enddo
1449  endif
1450 
1451  call set_domain ( atm(mytile)%domain )
1452 
1453  call timing_on('GFS_TENDENCIES')
1454 !--- put u/v tendencies into haloed arrays u_dt and v_dt
1455 !$OMP parallel do default (none) &
1456 !$OMP shared (rdt, n, nq, dnats, npz, ncnst, nwat, mytile, u_dt, v_dt, t_dt,&
1457 !$OMP Atm, IPD_Data, Atm_block, sphum, liq_wat, rainwat, ice_wat, &
1458 #ifdef MULTI_GASES
1459 !$OMP num_gas, &
1460 #endif
1461 !$OMP snowwat, graupel, nq_adv, flip_vc) &
1462 !$OMP private (nb, blen, i, j, k, k1, ix, q0, qwat, qt)
1463  do nb = 1,atm_block%nblks
1464 
1465 !SJL: perform vertical filling to fix the negative humidity if the SAS convection scheme is used
1466 ! This call may be commented out if RAS or other positivity-preserving CPS is used.
1467  blen = atm_block%blksz(nb)
1468  call fill_gfs(blen, npz, ipd_data(nb)%Statein%prsi, ipd_data(nb)%Stateout%gq0, 1.e-9_kind_phys)
1469 
1470  do k = 1, npz
1471  if(flip_vc) then
1472  k1 = npz+1-k !reverse the k direction
1473  else
1474  k1 = k
1475  endif
1476  do ix = 1, blen
1477  i = atm_block%index(nb)%ii(ix)
1478  j = atm_block%index(nb)%jj(ix)
1479  u_dt(i,j,k1) = u_dt(i,j,k1) + (ipd_data(nb)%Stateout%gu0(ix,k) - ipd_data(nb)%Statein%ugrs(ix,k)) * rdt
1480  v_dt(i,j,k1) = v_dt(i,j,k1) + (ipd_data(nb)%Stateout%gv0(ix,k) - ipd_data(nb)%Statein%vgrs(ix,k)) * rdt
1481 ! t_dt(i,j,k1) = (IPD_Data(nb)%Stateout%gt0(ix,k) - IPD_Data(nb)%Statein%tgrs(ix,k)) * rdt
1482  t_dt(i,j,k1) = t_dt(i,j,k1) + (ipd_data(nb)%Stateout%gt0(ix,k) - ipd_data(nb)%Statein%tgrs(ix,k)) * rdt
1483 ! SJL notes:
1484 ! ---- DO not touch the code below; dry mass conservation may change due to 64bit <-> 32bit conversion
1485 ! GFS total air mass = dry_mass + water_vapor (condensate excluded)
1486 ! GFS mixing ratios = tracer_mass / (dry_mass + vapor_mass)
1487 ! FV3 total air mass = dry_mass + [water_vapor + condensate ]
1488 ! FV3 mixing ratios = tracer_mass / (dry_mass+vapor_mass+cond_mass)
1489  if(flip_vc) then
1490  q0 = ipd_data(nb)%Statein%prsi(ix,k) - ipd_data(nb)%Statein%prsi(ix,k+1)
1491  else
1492  q0 = ipd_data(nb)%Statein%prsi(ix,k+1) - ipd_data(nb)%Statein%prsi(ix,k)
1493  endif
1494  qwat(1:nq_adv) = q0*ipd_data(nb)%Stateout%gq0(ix,k,1:nq_adv)
1495 ! **********************************************************************************************************
1496 ! Dry mass: the following way of updating delp is key to mass conservation with hybrid 32-64 bit computation
1497 ! **********************************************************************************************************
1498 ! The following example is for 2 water species.
1499 ! q0 = Atm(n)%delp(i,j,k1)*(1.-(Atm(n)%q(i,j,k1,1)+Atm(n)%q(i,j,k1,2))) + q1 + q2
1500 #ifdef MULTI_GASES
1501  q0 = atm(n)%delp(i,j,k1)*(1.-sum(atm(n)%q(i,j,k1,1:max(nwat,num_gas)))) + sum(qwat(1:max(nwat,num_gas)))
1502 #else
1503  qt = sum(qwat(1:nwat))
1504  q0 = atm(n)%delp(i,j,k1)*(1.-sum(atm(n)%q(i,j,k1,1:nwat))) + qt
1505 #endif
1506  atm(n)%delp(i,j,k1) = q0
1507  atm(n)%q(i,j,k1,1:nq_adv) = qwat(1:nq_adv) / q0
1508 ! if (dnats .gt. 0) Atm(n)%q(i,j,k1,nq_adv+1:nq) = IPD_Data(nb)%Stateout%gq0(ix,k,nq_adv+1:nq)
1509  enddo
1510  enddo
1511 
1512  !--- diagnostic tracers are assumed to be updated in-place
1513  !--- SHOULD THESE DIAGNOSTIC TRACERS BE MASS ADJUSTED???
1514  !--- See Note in statein...
1515  do iq = nq+1, ncnst
1516  do k = 1, npz
1517  if(flip_vc) then
1518  k1 = npz+1-k !reverse the k direction
1519  else
1520  k1 = k
1521  endif
1522  do ix = 1, blen
1523  i = atm_block%index(nb)%ii(ix)
1524  j = atm_block%index(nb)%jj(ix)
1525  atm(mytile)%qdiag(i,j,k1,iq) = ipd_data(nb)%Stateout%gq0(ix,k,iq)
1526  enddo
1527  enddo
1528  enddo
1529 
1530  enddo ! nb-loop
1531 
1532  call timing_off('GFS_TENDENCIES')
1533 
1534  w_diff = get_tracer_index(model_atmos, 'w_diff' )
1535  nt_dyn = ncnst-pnats !nothing more than nq
1536  if ( w_diff /= no_tracer ) then
1537  nt_dyn = nt_dyn - 1
1538  endif
1539 
1540 !--- adjust w and heat tendency for non-hydrostatic case
1541 #ifdef USE_Q_DT
1542  if ( .not.atm(n)%flagstruct%hydrostatic .and. w_diff /= no_tracer ) then
1543  rcp = 1. / cp_air
1544 !$OMP parallel do default (none) &
1545 !$OMP shared (jsc, jec, isc, iec, n, w_diff, Atm, q_dt, t_dt, rcp, dt_atmos) &
1546 !$OMP private (i, j, k)
1547  do k=1, atm(n)%npz
1548  do j=jsc, jec
1549  do i=isc, iec
1550  atm(n)%q(i,j,k,w_diff) = q_dt(i,j,k,w_diff) ! w tendency due to phys
1551 ! Heating due to loss of KE (vertical diffusion of w)
1552  t_dt(i,j,k) = t_dt(i,j,k) - q_dt(i,j,k,w_diff)*rcp*&
1553  (atm(n)%w(i,j,k)+0.5*dt_atmos*q_dt(i,j,k,w_diff))
1554  atm(n)%w(i,j,k) = atm(n)%w(i,j,k) + dt_atmos*atm(n)%q(i,j,k,w_diff)
1555  enddo
1556  enddo
1557  enddo
1558  endif
1559 #endif
1560 
1561  call mpp_clock_begin (id_dynam)
1562  call timing_on('FV_UPDATE_PHYS')
1563  call fv_update_phys( dt_atmos, isc, iec, jsc, jec, isd, ied, jsd, jed, atm(n)%ng, nt_dyn, &
1564  atm(n)%u, atm(n)%v, atm(n)%w, atm(n)%delp, atm(n)%pt, &
1565  atm(n)%q, atm(n)%qdiag, &
1566  atm(n)%ua, atm(n)%va, atm(n)%ps, atm(n)%pe, atm(n)%peln, &
1567  atm(n)%pk, atm(n)%pkz, atm(n)%ak, atm(n)%bk, atm(n)%phis, &
1568  atm(n)%u_srf, atm(n)%v_srf, atm(n)%ts, atm(n)%delz, &
1569  atm(n)%flagstruct%hydrostatic, u_dt, v_dt, t_dt, &
1570  .true., time_next, atm(n)%flagstruct%nudge, atm(n)%gridstruct, &
1571  atm(n)%gridstruct%agrid(:,:,1), atm(n)%gridstruct%agrid(:,:,2), &
1572  atm(n)%npx, atm(n)%npy, atm(n)%npz, atm(n)%flagstruct, &
1573  atm(n)%neststruct, atm(n)%bd, atm(n)%domain, atm(n)%ptop)
1574  call timing_off('FV_UPDATE_PHYS')
1575  call mpp_clock_end (id_dynam)
1576 
1577 !--- nesting update after updating atmospheric variables with
1578 !--- physics tendencies
1579  if (ngrids > 1 .and. p_split > 0) then
1580  call timing_on('TWOWAY_UPDATE')
1582  call timing_off('TWOWAY_UPDATE')
1583  endif
1584  call nullify_domain()
1585 
1586  !---- diagnostics for FV dynamics -----
1587  if (atm(mytile)%flagstruct%print_freq /= -99) then
1588  call mpp_clock_begin(id_fv_diag)
1589 
1590  fv_time = time_next
1591  call get_time (fv_time, seconds, days)
1592 
1593  call nullify_domain()
1594  call timing_on('FV_DIAG')
1595  call fv_diag(atm(mytile:mytile), zvir, fv_time, atm(mytile)%flagstruct%print_freq)
1596  first_diag = .false.
1597  call timing_off('FV_DIAG')
1598 
1599  call mpp_clock_end(id_fv_diag)
1600  endif
1601 
1602  end subroutine atmosphere_state_update
1603 
1604 
1608  subroutine adiabatic_init(zvir,nudge_dz,time)
1609  type(time_type),intent(in) :: time
1610  real, allocatable, dimension(:,:,:):: u0, v0, t0, dz0, dp0
1611  real, intent(in) :: zvir
1612  logical, intent(inout):: nudge_dz
1613 ! real, parameter:: wt = 1. ! was 2.
1614  real, parameter:: wt = 2.
1615 !***********
1616 ! Haloe Data
1617 !***********
1618  real, parameter:: q1_h2o = 2.2e-6
1619  real, parameter:: q7_h2o = 3.8e-6
1620  real, parameter:: q100_h2o = 3.8e-6
1621  real, parameter:: q1000_h2o = 3.1e-6
1622  real, parameter:: q2000_h2o = 2.8e-6
1623  real, parameter:: q3000_h2o = 3.0e-6
1624  real:: xt, p00, q00
1625  integer:: isc, iec, jsc, jec, npz
1626  integer:: m, n, i,j,k, ngc, n_split_loc, days
1627 
1628  character(len=80) :: errstr
1629 
1630  xt = 1./(1.+wt)
1631 
1632  write(errstr,'(A, I4, A)') 'Performing adiabatic init', atm(mytile)%flagstruct%na_init, ' times'
1633  call mpp_error(note, errstr)
1634  sphum = get_tracer_index(model_atmos, 'sphum' )
1635 
1636  npz = atm(mytile)%npz
1637 
1638  isc = atm(mytile)%bd%isc
1639  iec = atm(mytile)%bd%iec
1640  jsc = atm(mytile)%bd%jsc
1641  jec = atm(mytile)%bd%jec
1642 
1643  ngc = atm(mytile)%ng
1644  isd = isc - ngc
1645  ied = iec + ngc
1646  jsd = jsc - ngc
1647  jed = jec + ngc
1648 
1649  call timing_on('adiabatic_init')
1650  do_adiabatic_init = .true.
1651 
1652  allocate ( u0(isc:iec, jsc:jec+1, npz) )
1653  allocate ( v0(isc:iec+1,jsc:jec, npz) )
1654  allocate (dp0(isc:iec,jsc:jec, npz) )
1655 
1656  if ( atm(mytile)%flagstruct%hydrostatic ) nudge_dz = .false.
1657 
1658  if ( nudge_dz ) then
1659  allocate (dz0(isc:iec,jsc:jec, npz) )
1660  else
1661  allocate ( t0(isc:iec,jsc:jec, npz) )
1662  endif
1663 
1664 !$omp parallel do default (none) &
1665 !$omp shared (nudge_dz, npz, jsc, jec, isc, iec, n, sphum, u0, v0, t0, dz0, dp0, Atm, zvir, mytile) &
1666 !$omp private (k, j, i)
1667  do k=1,npz
1668  do j=jsc,jec+1
1669  do i=isc,iec
1670  u0(i,j,k) = atm(mytile)%u(i,j,k)
1671  enddo
1672  enddo
1673  do j=jsc,jec
1674  do i=isc,iec+1
1675  v0(i,j,k) = atm(mytile)%v(i,j,k)
1676  enddo
1677  enddo
1678  if ( nudge_dz ) then
1679  do j=jsc,jec
1680  do i=isc,iec
1681  dp0(i,j,k) = atm(mytile)%delp(i,j,k)
1682  dz0(i,j,k) = atm(mytile)%delz(i,j,k)
1683  enddo
1684  enddo
1685  else
1686  do j=jsc,jec
1687  do i=isc,iec
1688 #ifdef MULTI_GASES
1689  t0(i,j,k) = atm(mytile)%pt(i,j,k)*virq(atm(mytile)%q(i,j,k,:)) ! virt T
1690 #else
1691  t0(i,j,k) = atm(mytile)%pt(i,j,k)*(1.+zvir*atm(mytile)%q(i,j,k,sphum)) ! virt T
1692 #endif
1693  dp0(i,j,k) = atm(mytile)%delp(i,j,k)
1694  enddo
1695  enddo
1696  endif
1697  enddo
1698 
1699  call get_time (time, seconds, days)
1700  if (seconds < nint(3600*atm(mytile)%flagstruct%fhouri) .and. atm(mytile)%flagstruct%fac_n_spl > 1.0) then
1701  n_split_loc = nint(atm(mytile)%flagstruct%n_split * atm(mytile)%flagstruct%fac_n_spl)
1702  else
1703  n_split_loc = atm(mytile)%flagstruct%n_split
1704  endif
1705 
1706 ! write(0,*)' before calling init n_split_loc=',n_split_loc,' seconds=',seconds,' days=',days,&
1707 ! ' n_split=',Atm(mytile)%flagstruct%n_split,' mytile=',mytile
1708 
1709  do m=1,atm(mytile)%flagstruct%na_init
1710 ! Forward call
1711  call fv_dynamics(atm(mytile)%npx, atm(mytile)%npy, npz, nq, atm(mytile)%ng, dt_atmos, 0., &
1712  atm(mytile)%flagstruct%fill, atm(mytile)%flagstruct%reproduce_sum, kappa, cp_air, zvir, &
1713 ! Atm(mytile)%ptop, Atm(mytile)%ks, nq, Atm(mytile)%flagstruct%n_split, &
1714  atm(mytile)%ptop, atm(mytile)%ks, nq, n_split_loc, &
1715  atm(mytile)%flagstruct%q_split, atm(mytile)%u, atm(mytile)%v, atm(mytile)%w, &
1716  atm(mytile)%delz, atm(mytile)%flagstruct%hydrostatic, &
1717  atm(mytile)%pt, atm(mytile)%delp, atm(mytile)%q, atm(mytile)%ps, &
1718  atm(mytile)%pe, atm(mytile)%pk, atm(mytile)%peln, atm(mytile)%pkz, atm(mytile)%phis, &
1719  atm(mytile)%q_con, atm(mytile)%omga, atm(mytile)%ua, atm(mytile)%va, atm(mytile)%uc, atm(mytile)%vc, &
1720  atm(mytile)%ak, atm(mytile)%bk, atm(mytile)%mfx, atm(mytile)%mfy, &
1721  atm(mytile)%cx, atm(mytile)%cy, atm(mytile)%ze0, atm(mytile)%flagstruct%hybrid_z, &
1722  atm(mytile)%gridstruct, atm(mytile)%flagstruct, &
1723  atm(mytile)%neststruct, atm(mytile)%idiag, atm(mytile)%bd, atm(mytile)%parent_grid, &
1724  atm(mytile)%domain,atm(mytile)%diss_est)
1725 ! Backward
1726  call fv_dynamics(atm(mytile)%npx, atm(mytile)%npy, npz, nq, atm(mytile)%ng, -dt_atmos, 0., &
1727  atm(mytile)%flagstruct%fill, atm(mytile)%flagstruct%reproduce_sum, kappa, cp_air, zvir, &
1728 ! Atm(mytile)%ptop, Atm(mytile)%ks, nq, Atm(mytile)%flagstruct%n_split, &
1729  atm(mytile)%ptop, atm(mytile)%ks, nq, n_split_loc, &
1730  atm(mytile)%flagstruct%q_split, atm(mytile)%u, atm(mytile)%v, atm(mytile)%w, &
1731  atm(mytile)%delz, atm(mytile)%flagstruct%hydrostatic, &
1732  atm(mytile)%pt, atm(mytile)%delp, atm(mytile)%q, atm(mytile)%ps, &
1733  atm(mytile)%pe, atm(mytile)%pk, atm(mytile)%peln, atm(mytile)%pkz, atm(mytile)%phis, &
1734  atm(mytile)%q_con, atm(mytile)%omga, atm(mytile)%ua, atm(mytile)%va, atm(mytile)%uc, atm(mytile)%vc, &
1735  atm(mytile)%ak, atm(mytile)%bk, atm(mytile)%mfx, atm(mytile)%mfy, &
1736  atm(mytile)%cx, atm(mytile)%cy, atm(mytile)%ze0, atm(mytile)%flagstruct%hybrid_z, &
1737  atm(mytile)%gridstruct, atm(mytile)%flagstruct, &
1738  atm(mytile)%neststruct, atm(mytile)%idiag, atm(mytile)%bd, atm(mytile)%parent_grid, &
1739  atm(mytile)%domain,atm(mytile)%diss_est)
1740 !Nudging back to IC
1741 !$omp parallel do default (none) &
1742 !$omp shared (pref, npz, jsc, jec, isc, iec, n, sphum, Atm, u0, v0, t0, dp0, xt, zvir, mytile, nudge_dz, dz0) &
1743 !$omp private (i, j, k, p00, q00)
1744  do k=1,npz
1745  do j=jsc,jec+1
1746  do i=isc,iec
1747  atm(mytile)%u(i,j,k) = xt*(atm(mytile)%u(i,j,k) + wt*u0(i,j,k))
1748  enddo
1749  enddo
1750  do j=jsc,jec
1751  do i=isc,iec+1
1752  atm(mytile)%v(i,j,k) = xt*(atm(mytile)%v(i,j,k) + wt*v0(i,j,k))
1753  enddo
1754  enddo
1755  if( atm(mytile)%flagstruct%nudge_qv ) then
1756 ! SJL note: Nudging water vaport towards HALOE climatology:
1757 ! In case of better IC (IFS) this step may not be necessary
1758  p00 = atm(mytile)%pe(isc,k,jsc)
1759  if ( p00 < 30.e2 ) then
1760  if ( p00 < 1. ) then
1761  q00 = q1_h2o
1762  elseif ( p00 <= 7. .and. p00 >= 1. ) then
1763  q00 = q1_h2o + (q7_h2o-q1_h2o)*log(pref(k,1)/1.)/log(7.)
1764  elseif ( p00 < 100. .and. p00 >= 7. ) then
1765  q00 = q7_h2o + (q100_h2o-q7_h2o)*log(pref(k,1)/7.)/log(100./7.)
1766  elseif ( p00 < 1000. .and. p00 >= 100. ) then
1767  q00 = q100_h2o + (q1000_h2o-q100_h2o)*log(pref(k,1)/1.e2)/log(10.)
1768  elseif ( p00 < 2000. .and. p00 >= 1000. ) then
1769  q00 = q1000_h2o + (q2000_h2o-q1000_h2o)*log(pref(k,1)/1.e3)/log(2.)
1770  else
1771  q00 = q2000_h2o + (q3000_h2o-q2000_h2o)*log(pref(k,1)/2.e3)/log(1.5)
1772  endif
1773  do j=jsc,jec
1774  do i=isc,iec
1775  atm(mytile)%q(i,j,k,sphum) = xt*(atm(mytile)%q(i,j,k,sphum) + wt*q00)
1776  enddo
1777  enddo
1778  endif
1779  endif
1780  if ( nudge_dz ) then
1781  do j=jsc,jec
1782  do i=isc,iec
1783  atm(mytile)%delp(i,j,k) = xt*(atm(mytile)%delp(i,j,k) + wt*dp0(i,j,k))
1784  atm(mytile)%delz(i,j,k) = xt*(atm(mytile)%delz(i,j,k) + wt*dz0(i,j,k))
1785  enddo
1786  enddo
1787  else
1788  do j=jsc,jec
1789  do i=isc,iec
1790 #ifdef MULTI_GASES
1791  atm(mytile)%pt(i,j,k) = xt*(atm(mytile)%pt(i,j,k) + wt*t0(i,j,k)/virq(atm(mytile)%q(i,j,k,:)))
1792 #else
1793  atm(mytile)%pt(i,j,k) = xt*(atm(mytile)%pt(i,j,k) + wt*t0(i,j,k)/(1.+zvir*atm(mytile)%q(i,j,k,sphum)))
1794 #endif
1795  atm(mytile)%delp(i,j,k) = xt*(atm(mytile)%delp(i,j,k) + wt*dp0(i,j,k))
1796  enddo
1797  enddo
1798  endif
1799 
1800  enddo
1801 
1802 ! Backward
1803  call fv_dynamics(atm(mytile)%npx, atm(mytile)%npy, npz, nq, atm(mytile)%ng, -dt_atmos, 0., &
1804  atm(mytile)%flagstruct%fill, atm(mytile)%flagstruct%reproduce_sum, kappa, cp_air, zvir, &
1805 ! Atm(mytile)%ptop, Atm(mytile)%ks, nq, Atm(mytile)%flagstruct%n_split, &
1806  atm(mytile)%ptop, atm(mytile)%ks, nq, n_split_loc, &
1807  atm(mytile)%flagstruct%q_split, atm(mytile)%u, atm(mytile)%v, atm(mytile)%w, &
1808  atm(mytile)%delz, atm(mytile)%flagstruct%hydrostatic, &
1809  atm(mytile)%pt, atm(mytile)%delp, atm(mytile)%q, atm(mytile)%ps, &
1810  atm(mytile)%pe, atm(mytile)%pk, atm(mytile)%peln, atm(mytile)%pkz, atm(mytile)%phis, &
1811  atm(mytile)%q_con, atm(mytile)%omga, atm(mytile)%ua, atm(mytile)%va, atm(mytile)%uc, atm(mytile)%vc, &
1812  atm(mytile)%ak, atm(mytile)%bk, atm(mytile)%mfx, atm(mytile)%mfy, &
1813  atm(mytile)%cx, atm(mytile)%cy, atm(mytile)%ze0, atm(mytile)%flagstruct%hybrid_z, &
1814  atm(mytile)%gridstruct, atm(mytile)%flagstruct, &
1815  atm(mytile)%neststruct, atm(mytile)%idiag, atm(mytile)%bd, atm(mytile)%parent_grid, &
1816  atm(mytile)%domain,atm(mytile)%diss_est)
1817 ! Forward call
1818  call fv_dynamics(atm(mytile)%npx, atm(mytile)%npy, npz, nq, atm(mytile)%ng, dt_atmos, 0., &
1819  atm(mytile)%flagstruct%fill, atm(mytile)%flagstruct%reproduce_sum, kappa, cp_air, zvir, &
1820 ! Atm(mytile)%ptop, Atm(mytile)%ks, nq, Atm(mytile)%flagstruct%n_split, &
1821  atm(mytile)%ptop, atm(mytile)%ks, nq, n_split_loc, &
1822  atm(mytile)%flagstruct%q_split, atm(mytile)%u, atm(mytile)%v, atm(mytile)%w, &
1823  atm(mytile)%delz, atm(mytile)%flagstruct%hydrostatic, &
1824  atm(mytile)%pt, atm(mytile)%delp, atm(mytile)%q, atm(mytile)%ps, &
1825  atm(mytile)%pe, atm(mytile)%pk, atm(mytile)%peln, atm(mytile)%pkz, atm(mytile)%phis, &
1826  atm(mytile)%q_con, atm(mytile)%omga, atm(mytile)%ua, atm(mytile)%va, atm(mytile)%uc, atm(mytile)%vc, &
1827  atm(mytile)%ak, atm(mytile)%bk, atm(mytile)%mfx, atm(mytile)%mfy, &
1828  atm(mytile)%cx, atm(mytile)%cy, atm(mytile)%ze0, atm(mytile)%flagstruct%hybrid_z, &
1829  atm(mytile)%gridstruct, atm(mytile)%flagstruct, &
1830  atm(mytile)%neststruct, atm(mytile)%idiag, atm(mytile)%bd, atm(mytile)%parent_grid, &
1831  atm(mytile)%domain,atm(mytile)%diss_est)
1832 ! Nudging back to IC
1833 !$omp parallel do default (none) &
1834 !$omp shared (nudge_dz,npz, jsc, jec, isc, iec, n, sphum, Atm, u0, v0, t0, dz0, dp0, xt, zvir, mytile) &
1835 !$omp private (i, j, k)
1836  do k=1,npz
1837  do j=jsc,jec+1
1838  do i=isc,iec
1839  atm(mytile)%u(i,j,k) = xt*(atm(mytile)%u(i,j,k) + wt*u0(i,j,k))
1840  enddo
1841  enddo
1842  do j=jsc,jec
1843  do i=isc,iec+1
1844  atm(mytile)%v(i,j,k) = xt*(atm(mytile)%v(i,j,k) + wt*v0(i,j,k))
1845  enddo
1846  enddo
1847  if ( nudge_dz ) then
1848  do j=jsc,jec
1849  do i=isc,iec
1850  atm(mytile)%delp(i,j,k) = xt*(atm(mytile)%delp(i,j,k) + wt*dp0(i,j,k))
1851  atm(mytile)%delz(i,j,k) = xt*(atm(mytile)%delz(i,j,k) + wt*dz0(i,j,k))
1852  enddo
1853  enddo
1854  else
1855  do j=jsc,jec
1856  do i=isc,iec
1857 #ifdef MULTI_GASES
1858  atm(mytile)%pt(i,j,k) = xt*(atm(mytile)%pt(i,j,k) + wt*t0(i,j,k)/virq(atm(mytile)%q(i,j,k,:)))
1859 #else
1860  atm(mytile)%pt(i,j,k) = xt*(atm(mytile)%pt(i,j,k) + wt*t0(i,j,k)/(1.+zvir*atm(mytile)%q(i,j,k,sphum)))
1861 #endif
1862  atm(mytile)%delp(i,j,k) = xt*(atm(mytile)%delp(i,j,k) + wt*dp0(i,j,k))
1863  enddo
1864  enddo
1865  endif
1866  enddo
1867 
1868  enddo
1869 
1870  deallocate ( u0 )
1871  deallocate ( v0 )
1872  deallocate (dp0 )
1873  if ( allocated(t0) ) deallocate ( t0 )
1874  if ( allocated(dz0) ) deallocate ( dz0 )
1875 
1876  do_adiabatic_init = .false.
1877  call timing_off('adiabatic_init')
1878 
1879  end subroutine adiabatic_init
1880 
1881 
1882 
1889 #if defined(OVERLOAD_R4)
1890 #define _DBL_(X) DBLE(X)
1891 #define _RL_(X) REAL(X,KIND=4)
1892 #else
1893 #define _DBL_(X) X
1894 #define _RL_(X) X
1895 #endif
1896  subroutine atmos_phys_driver_statein (IPD_Data, Atm_block,flip_vc)
1897  type(ipd_data_type), intent(inout) :: ipd_data(:)
1898  type(block_control_type), intent(in) :: atm_block
1899  logical, intent(in) :: flip_vc
1900 !--------------------------------------
1901 ! Local GFS-phys consistent parameters:
1902 !--------------------------------------
1903  real(kind=kind_phys), parameter :: p00 = 1.e5
1904  real(kind=kind_phys), parameter :: qmin = 1.0e-10
1905  real(kind=kind_phys) :: pk0inv, ptop, pktop
1906  real(kind=kind_phys) :: rTv, dm, qgrs_rad
1907  integer :: nb, blen, npz, i, j, k, ix, k1, kz, dnats, nq_adv
1908 
1909 !!! NOTES: lmh 6nov15
1910 !!! - "Layer" means "layer mean", ie. the average value in a layer
1911 !!! - "Level" means "level interface", ie the point values at the top or bottom of a layer
1912 
1913  ptop = _dbl_(_rl_(atm(mytile)%ak(1)))
1914  pktop = (ptop/p00)**kappa
1915  pk0inv = (1.0_kind_phys/p00)**kappa
1916 
1917  npz = atm_block%npz
1918  dnats = atm(mytile)%flagstruct%dnats
1919  nq_adv = nq - dnats
1920 
1921 !---------------------------------------------------------------------
1922 ! use most up to date atmospheric properties when running serially
1923 !---------------------------------------------------------------------
1924 !$OMP parallel do default (none) &
1925 !$OMP shared (Atm_block, Atm, IPD_Data, npz, nq, ncnst, sphum, liq_wat, &
1926 !$OMP ice_wat, rainwat, snowwat, graupel, pk0inv, ptop, &
1927 !$OMP pktop, zvir, mytile, dnats, nq_adv, flip_vc) &
1928 !$OMP private (dm, nb, blen, i, j, ix, k1, kz, rTv, qgrs_rad)
1929 
1930  do nb = 1,atm_block%nblks
1931 ! gas_phase_mass <-- prsl
1932 ! log(pe) <-- prsik
1933 
1934  blen = atm_block%blksz(nb)
1935 
1936  !-- level interface geopotential height (relative to the surface)
1937  if(flip_vc) then
1938  ipd_data(nb)%Statein%phii(:,1) = 0.0_kind_phys
1939  else
1940  ipd_data(nb)%Statein%phii(:,npz+1) = 0.0_kind_phys
1941  endif
1942  ipd_data(nb)%Statein%prsik(:,:) = 1.e25_kind_phys
1943 
1944  do k = 1, npz
1945  !Indices for FV's vertical coordinate, for which 1 = top
1946  !here, k is the index for GFS's vertical coordinate, for which 1 =
1947  !bottom
1948  kz = npz+1-k
1949  if(flip_vc) then
1950  k1 = kz ! flipping the index
1951  else
1952  k1 = k
1953  endif
1954  do ix = 1, blen
1955  i = atm_block%index(nb)%ii(ix)
1956  j = atm_block%index(nb)%jj(ix)
1957 
1958  ipd_data(nb)%Statein%tgrs(ix,k) = _dbl_(_rl_(atm(mytile)%pt(i,j,k1)))
1959  ipd_data(nb)%Statein%ugrs(ix,k) = _dbl_(_rl_(atm(mytile)%ua(i,j,k1)))
1960  ipd_data(nb)%Statein%vgrs(ix,k) = _dbl_(_rl_(atm(mytile)%va(i,j,k1)))
1961  ipd_data(nb)%Statein%vvl(ix,k) = _dbl_(_rl_(atm(mytile)%omga(i,j,k1)))
1962  ipd_data(nb)%Statein%prsl(ix,k) = _dbl_(_rl_(atm(mytile)%delp(i,j,k1))) ! Total mass
1963  if (atm(mytile)%flagstruct%do_skeb)ipd_data(nb)%Statein%diss_est(ix,k) = _dbl_(_rl_(atm(mytile)%diss_est(i,j,k1)))
1964 
1965  if(flip_vc) then
1966  if (.not.atm(mytile)%flagstruct%hydrostatic .and. (.not.atm(mytile)%flagstruct%use_hydro_pressure)) &
1967  ipd_data(nb)%Statein%phii(ix,k+1) = ipd_data(nb)%Statein%phii(ix,k) - _dbl_(_rl_(atm(mytile)%delz(i,j,k1)*grav))
1968  else
1969  if (.not.atm(mytile)%flagstruct%hydrostatic .and. (.not.atm(mytile)%flagstruct%use_hydro_pressure)) &
1970  ipd_data(nb)%Statein%phii(ix,kz) = ipd_data(nb)%Statein%phii(ix,kz+1) - _dbl_(_rl_(atm(mytile)%delz(i,j,kz)*grav))
1971  endif
1972 
1973 ! Convert to tracer mass:
1974  ipd_data(nb)%Statein%qgrs(ix,k,1:nq_adv) = _dbl_(_rl_(atm(mytile)%q(i,j,k1,1:nq_adv))) &
1975  * ipd_data(nb)%Statein%prsl(ix,k)
1976  if (dnats > 0) &
1977  ipd_data(nb)%Statein%qgrs(ix,k,nq_adv+1:nq) = _dbl_(_rl_(atm(mytile)%q(i,j,k1,nq_adv+1:nq)))
1978  !--- SHOULD THESE BE CONVERTED TO MASS SINCE THE DYCORE DOES NOT TOUCH THEM IN ANY WAY???
1979  !--- See Note in state update...
1980  if ( ncnst > nq) &
1981  ipd_data(nb)%Statein%qgrs(ix,k,nq+1:ncnst) = _dbl_(_rl_(atm(mytile)%qdiag(i,j,k1,nq+1:ncnst)))
1982 ! Remove the contribution of condensates to delp (mass):
1983  if ( atm(mytile)%flagstruct%nwat == 6 ) then
1984  ipd_data(nb)%Statein%prsl(ix,k) = ipd_data(nb)%Statein%prsl(ix,k) &
1985  - ipd_data(nb)%Statein%qgrs(ix,k,liq_wat) &
1986  - ipd_data(nb)%Statein%qgrs(ix,k,ice_wat) &
1987  - ipd_data(nb)%Statein%qgrs(ix,k,rainwat) &
1988  - ipd_data(nb)%Statein%qgrs(ix,k,snowwat) &
1989  - ipd_data(nb)%Statein%qgrs(ix,k,graupel)
1990  else !variable condensate numbers
1991  ipd_data(nb)%Statein%prsl(ix,k) = ipd_data(nb)%Statein%prsl(ix,k) &
1992  - sum(ipd_data(nb)%Statein%qgrs(ix,k,2:atm(mytile)%flagstruct%nwat))
1993  endif
1994  enddo
1995  enddo
1996 
1997 ! Re-compute pressure (dry_mass + water_vapor) derived fields:
1998  if(flip_vc) then
1999  do i=1,blen
2000  ipd_data(nb)%Statein%prsi(i,npz+1) = ptop
2001  enddo
2002  do k=npz,1,-1
2003  do i=1,blen
2004  ipd_data(nb)%Statein%prsi(i,k) = ipd_data(nb)%Statein%prsi(i,k+1) + ipd_data(nb)%Statein%prsl(i,k)
2005  ipd_data(nb)%Statein%prsik(i,k) = log( ipd_data(nb)%Statein%prsi(i,k) )
2006 ! Redefine mixing ratios for GFS == tracer_mass / (dry_air_mass + water_vapor_mass)
2007  ipd_data(nb)%Statein%qgrs(i,k,1:nq_adv) = ipd_data(nb)%Statein%qgrs(i,k,1:nq_adv) &
2008  / ipd_data(nb)%Statein%prsl(i,k)
2009  enddo
2010  enddo
2011  else
2012  do i=1,blen
2013  ipd_data(nb)%Statein%prsi(i, 1) = ptop
2014  enddo
2015  do k=1,npz
2016  do i=1,blen
2017  ipd_data(nb)%Statein%prsi(i,k+1) = ipd_data(nb)%Statein%prsi(i,k) + ipd_data(nb)%Statein%prsl(i,k)
2018  ipd_data(nb)%Statein%prsik(i,k) = log( ipd_data(nb)%Statein%prsi(i,k) )
2019 ! Redefine mixing ratios for GFS == tracer_mass / (dry_air_mass + water_vapor_mass)
2020  ipd_data(nb)%Statein%qgrs(i,k,1:nq_adv) = ipd_data(nb)%Statein%qgrs(i,k,1:nq_adv) &
2021  / ipd_data(nb)%Statein%prsl(i,k)
2022  enddo
2023  enddo
2024  endif
2025 
2026  do i=1,blen
2027  ipd_data(nb)%Statein%pgr(i) = ipd_data(nb)%Statein%prsi(i,1) ! surface pressure for GFS
2028  ipd_data(nb)%Statein%prsik(i,npz+1) = log(ptop)
2029  enddo
2030 
2031  do k=1,npz
2032  do i=1,blen
2033 ! Geo-potential at interfaces:
2034 #ifdef MULTI_GASES
2035  rtv = rdgas*ipd_data(nb)%Statein%tgrs(i,k)*virq_max(ipd_data(nb)%Statein%qgrs(i,k,:),qmin)
2036 #else
2037  qgrs_rad = max(qmin,ipd_data(nb)%Statein%qgrs(i,k,sphum))
2038  rtv = rdgas*ipd_data(nb)%Statein%tgrs(i,k)*(1.+zvir*qgrs_rad)
2039 #endif
2040  if ( atm(mytile)%flagstruct%hydrostatic .or. atm(mytile)%flagstruct%use_hydro_pressure ) &
2041  ipd_data(nb)%Statein%phii(i,k+1) = ipd_data(nb)%Statein%phii(i,k) &
2042  + rtv*(ipd_data(nb)%Statein%prsik(i,k) &
2043  - ipd_data(nb)%Statein%prsik(i,k+1))
2044 ! Layer mean pressure by perfect gas law:
2045  dm = ipd_data(nb)%Statein%prsl(i,k)
2046  ipd_data(nb)%Statein%prsl(i,k) = dm*rtv/(ipd_data(nb)%Statein%phii(i,k+1) &
2047  - ipd_data(nb)%Statein%phii(i,k))
2048 
2049 !!! Ensure subgrid MONOTONICITY of Pressure: SJL 09/11/2016
2050  if ( .not.atm(mytile)%flagstruct%hydrostatic ) then
2051 ! If violated, replaces it with hydrostatic pressure
2052  ipd_data(nb)%Statein%prsl(i,k) = min(ipd_data(nb)%Statein%prsl(i,k), &
2053  ipd_data(nb)%Statein%prsi(i,k) - 0.01*dm)
2054  ipd_data(nb)%Statein%prsl(i,k) = max(ipd_data(nb)%Statein%prsl(i,k), &
2055  ipd_data(nb)%Statein%prsi(i,k+1) + 0.01*dm)
2056  endif
2057  enddo
2058  enddo
2059 
2060  do k = 1,npz
2061  do i=1,blen
2062 ! Exner function layer center: large sensitivity to non-hydro runs with moist kappa
2063  ipd_data(nb)%Statein%prslk(i,k) = exp( kappa*log(ipd_data(nb)%Statein%prsl(i,k)/p00) )
2064 !-- layer center geopotential; geometric midpoint
2065  ipd_data(nb)%Statein%phil(i,k) = 0.5_kind_phys*(ipd_data(nb)%Statein%phii(i,k) &
2066  + ipd_data(nb)%Statein%phii(i,k+1))
2067  enddo
2068  enddo
2069 
2070 ! Compute Exner function at layer "interfaces"
2071  do i=1,blen
2072 ! Top & Bottom edges computed hydrostaticlly
2073  ipd_data(nb)%Statein%prsik(i, 1) = exp( kappa*ipd_data(nb)%Statein%prsik(i,1) )*pk0inv ! bottom
2074  ipd_data(nb)%Statein%prsik(i,npz+1) = pktop ! TOA
2075  enddo
2076 
2077  if ( atm(mytile)%flagstruct%hydrostatic .or. atm(mytile)%flagstruct%use_hydro_pressure ) then
2078  do k=2,npz
2079  do i=1,blen
2080  ipd_data(nb)%Statein%prsik(i,k) = exp( kappa*ipd_data(nb)%Statein%prsik(i,k) )*pk0inv
2081  enddo
2082  enddo
2083  endif
2084  enddo
2085 
2086  end subroutine atmos_phys_driver_statein
2087 
2088 end module atmosphere_mod
subroutine, public atmosphere_end(Time, Grid_box)
The subroutine &#39;atmosphere_end&#39; is an API for the termination of the FV3 dynamical core responsible f...
Definition: atmosphere.F90:716
subroutine, public prt_height(qname, is, ie, js, je, ng, km, press, phis, delz, peln, area, lat)
logical first_diag
Definition: atmosphere.F90:262
subroutine, public fv_restart(fv_domain, Atm, dt_atmos, seconds, days, cold_start, grid_type, grids_on_this_pe)
The subroutine &#39;fv_restart&#39; initializes the model state, including prognaostic variables and several ...
Definition: fv_restart.F90:196
subroutine, public fv_update_phys(dt, is, ie, js, je, isd, ied, jsd, jed, ng, nq, u, v, w, delp, pt, q, qdiag, ua, va, ps, pe, peln, pk, pkz, ak, bk, phis, u_srf, v_srf, ts, delz, hydrostatic, u_dt, v_dt, t_dt, moist_phys, Time, nudge, gridstruct, lona, lata, npx, npy, npz, flagstruct, neststruct, bd, domain, ptop, q_dt)
subroutine, public atmosphere_dynamics(Time)
The subroutine &#39;atmosphere_dynamics&#39; is an API for the main driver of the FV3 dynamical core responsi...
Definition: atmosphere.F90:604
subroutine, public atmosphere_hgt(hgt, position, relative, flip)
The subroutine &#39;atmosphere_hgt&#39; is an API to return the height coordinate. By default, the vertical dimension assumes the standard FV3 convention of TOA (k=1) to Surface (k=npz). There are options to choose location [level (interface) or layer] and absolute vs. relative height (zero-based).
Definition: atmosphere.F90:935
logical, public do_adiabatic_init
Definition: fv_nudge.F90:138
integer id_udt_dyn
Definition: atmosphere.F90:254
subroutine get_stock_pe(index, value)
real, public current_time_in_seconds
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
subroutine timing_off(blk_name)
The subroutine &#39;timing_off&#39; stops a timer.
Definition: fv_timing.F90:180
subroutine, public fv_write_restart(Atm, grids_on_this_pe, timestamp)
The subroutine &#39;fv_write_restart&#39; writes restart files to disk.
subroutine, public fv_init(Atm, dt_atmos, grids_on_this_pe, p_split)
The subroutine &#39;fv_init&#39; initializes FV3.
Definition: fv_control.F90:373
integer, public p_step
subroutine, public atmosphere_get_bottom_layer(Atm_block, DYCORE_Data)
The subroutine &#39;atmosphere_get_bottom_layer&#39; is an API to provide the bottom layer quantities needed ...
integer id_vdt_dyn
Definition: atmosphere.F90:254
integer, public num_gas
Definition: multi_gases.F90:44
subroutine, public fv_nggps_tavg(Atm, Time_step_atmos, avg_max_length, zvir)
The module &#39;multi_gases&#39; peforms multi constitutents computations.
Definition: multi_gases.F90:25
pure real function, public virq_max(q, qmin)
subroutine, public atmosphere_diss_est(npass)
subroutine adiabatic_init(zvir, nudge_dz, time)
The subroutine &#39;adiabatic_init&#39; is an optional step during initialization to pre-condition a solution...
subroutine, public get_eta_level(npz, p_s, pf, ph, ak, bk, pscale)
The subroutine &#39;get_eta_level&#39; returns the interface and layer-mean pressures for reference...
Definition: fv_eta.F90:2588
subroutine, public atmosphere_control_data(i1, i2, j1, j2, kt, p_hydro, hydro, tile_num)
Definition: atmosphere.F90:809
subroutine, public atmos_phys_driver_statein(IPD_Data, Atm_block, flip_vc)
The subroutine &#39;atmos_phys_driver_statein&#39; is an API to populate the IPD_DataStatein container with t...
integer, public ngrids
Definition: fv_control.F90:348
logical, dimension(:), allocatable grids_on_this_pe
Definition: atmosphere.F90:251
The module &#39;fv_update_phys&#39; applies physics tendencies consistent with the FV3 discretization and def...
integer, parameter, public r_grid
Definition: fv_arrays.F90:35
The module &#39;fv_dynamics&#39; is the top-level routine for the dynamical core.
Definition: fv_dynamics.F90:26
real, dimension(:), allocatable ri
Definition: multi_gases.F90:48
The module &#39;fv_sg&#39; performs FV sub-grid mixing.
Definition: fv_sg.F90:54
The module fv_nwp_nudge contains routines for nudging to input analyses. note This module is currentl...
Definition: fv_nudge.F90:32
integer, dimension(:), allocatable id_tracerdt_dyn
Definition: atmosphere.F90:242
type(time_type) time_step_atmos
Definition: atmosphere.F90:228
real, dimension(:,:,:), allocatable u_dt
Definition: atmosphere.F90:259
subroutine get_bottom_mass(t_bot, tr_bot, p_bot, z_bot, p_surf, slp)
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
pure real function, public virq(q)
real, dimension(:), allocatable dum1d
Definition: atmosphere.F90:260
integer id_subgridz
Definition: atmosphere.F90:239
subroutine, public fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
subroutine, public fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, reproduce_sum, kappa, cp_air, zvir, ptop, ks, ncnst, n_split, q_split, u, v, w, delz, hydrostatic, pt, delp, q, ps, pe, pk, peln, pkz, phis, q_con, omga, ua, va, uc, vc, ak, bk, mfx, mfy, cx, cy, ze0, hybrid_z, gridstruct, flagstruct, neststruct, idiag, bd, parent_grid, domain, diss_est, time_total)
The module &#39;fv_nggps_diags&#39; computes output diagnostics entirely on 3D pressure levels.
incremental analysis update module
Definition: fv_iau_mod.F90:38
subroutine, public fv_nggps_diag(Atm, zvir, Time)
subroutine get_bottom_wind(u_bot, v_bot)
subroutine, public fv_nggps_diag_init(Atm, axes, Time)
subroutine, public atmosphere_diag_axes(axes)
The subroutine &#39;atmosphere_diag_axes&#39; is an API to return the axis indices for the atmospheric (mass)...
Definition: atmosphere.F90:896
subroutine, public fv_nwp_nudge_end
The subroutine &#39;fv_nwp_nudge_end&#39; terminates the nudging module.
Definition: fv_nudge.F90:2085
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
logical cold_start
Definition: atmosphere.F90:240
subroutine, public set_atmosphere_pelist()
Definition: atmosphere.F90:866
real, dimension(:,:,:), allocatable v_dt
Definition: atmosphere.F90:259
The module &#39;fv_eta&#39; contains routine to set up the reference (Eulerian) pressure coordinate.
Definition: fv_eta.F90:25
subroutine, public twoway_nesting(Atm, ngrids, grids_on_this_pe, zvir)
The subroutine&#39;twoway_nesting&#39; performs a two-way update of nested-grid data onto the parent grid...
integer, public mytile
Definition: atmosphere.F90:248
subroutine, public fv_diag(Atm, zvir, Time, print_freq)
The module &#39;atmosphere&#39; provides the interface for the Cubed-Sphere FV dynamical core.
Definition: atmosphere.F90:26
real, parameter w0_big
Definition: atmosphere.F90:256
subroutine, public atmosphere_resolution(i_size, j_size, global)
The subroutine &#39;atmospehre_resolution&#39; is an API to return the local extents of the current MPI-rank ...
Definition: atmosphere.F90:781
subroutine, public atmosphere_etalvls(ak, bk, flip)
The subroutine &#39;atmosphere_etalvls&#39; is an API to return the ak/bk pairs used to compute the eta or pr...
Definition: atmosphere.F90:913
subroutine, public atmosphere_init(Time_init, Time, Time_step, Grid_box, area)
The subroutine &#39;atmosphere_init&#39; is an API to initialize the FV3 dynamical core, including the grid s...
Definition: atmosphere.F90:271
subroutine, public del2_cubed(q, cd, gridstruct, domain, npx, npy, km, nmax, bd)
The subroutine &#39;del2-cubed&#39; filters the omega field for the physics.
Definition: dyn_core.F90:2421
type(fv_atmos_type), dimension(:), allocatable, target, public atm
Definition: atmosphere.F90:252
subroutine, public fill_gfs(im, km, pe2, q, q_min)
The subroutine &#39;fill_gfs&#39; is for mass-conservative filling of nonphysical negative values in the trac...
Definition: fv_fill.F90:163
integer, dimension(:), allocatable pelist
Definition: atmosphere.F90:250
subroutine p_adi(km, ng, ifirst, ilast, jfirst, jlast, ptop, delp, pt, ps, pe, peln, pk, pkz, hydrostatic)
The subroutine &#39;p_adi&#39; computes (ps, pk, pe, peln, pkz) given (ptop, delp).
Definition: atmosphere.F90:547
type(time_type), public fv_time
real, dimension(:), allocatable cpi
Definition: multi_gases.F90:49
subroutine, public atmosphere_pref(p_ref)
The subroutine &#39;atmosphere_pref&#39; is an API to return the reference pressure.
Definition: atmosphere.F90:801
subroutine timing_on(blk_name)
The subroutine &#39;timing_on&#39; starts a timer.
Definition: fv_timing.F90:116
@ The module &#39;fv_diagnostics&#39; contains routines to compute diagnosic fields.
subroutine, public atmosphere_restart(timestamp)
The subroutine &#39;atmosphere_restart&#39; is an API to save restart information at a given timestamp...
Definition: atmosphere.F90:770
The module &#39;dyn_core&#39; peforms the Lagrangian acoustic dynamics described by .
Definition: dyn_core.F90:28
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
character(len=20) mod_name
Definition: atmosphere.F90:225
subroutine, public start_regional_restart(Atm, isc, iec, jsc, jec, isd, ied, jsd, jed)
integer, public a_step
subroutine, public fv_subgrid_z(isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, tau, nwat, delp, pe, peln, pkz, ta, qa, ua, va, hydrostatic, w, delz, u_dt, v_dt, t_dt, q_dt, k_bot)
Definition: fv_sg.F90:680
real, dimension(:,:), allocatable pref
Definition: atmosphere.F90:260
subroutine, public atmosphere_scalar_field_halo(data, halo, isize, jsize, ksize, data_p)
The subroutine &#39;atmosphere_scalar_field_halo&#39; is an API to return halo information of the current MPI...
Definition: atmosphere.F90:995
subroutine, public fv_nwp_nudge_init(time, axes, npz, zvir, ak, bk, ts, phis, gridstruct, ks, npx, neststruct, bd)
The subroutine &#39;fv_nwp_nudge_init&#39; initializes the nudging module.
Definition: fv_nudge.F90:1257
subroutine, public atmosphere_nggps_diag(Time, init, ltavg, avg_max_length)
The subroutine &#39;atmosphere_nggps_diag&#39; is an API to trigger output of diagnostics in NCEP/EMC format...
integer id_dynam
Definition: atmosphere.F90:239
subroutine, public atmosphere_state_update(Time, IPD_Data, IAU_Data, Atm_block, flip_vc)
The subroutine &#39;atmosphere_state_update&#39; is an API to apply tendencies and compute a consistent progn...
The module &#39;FV3_control&#39; is for initialization and termination of the model, and controls namelist pa...
Definition: fv_control.F90:29
subroutine, public atmosphere_grid_bdry(blon, blat, global)
The subroutine &#39;atmosphere_grid_bdry&#39; is an API to returns the longitude and latitude finite volume e...
Definition: atmosphere.F90:845
subroutine, public atmosphere_grid_ctr(lon, lat)
The subroutine &#39;atmosphere_grid_ctr&#39; is an API that returns the longitude and latitude cell centers o...
Definition: atmosphere.F90:828
subroutine, public fv_end(Atm, grids_on_this_pe)
The subroutine &#39;fv_end&#39; terminates FV3, deallocates memory, saves restart files, and stops I/O...
Definition: fv_control.F90:596
subroutine, public read_new_bc_data(Atm, Time, Time_step_atmos, p_split, isd, ied, jsd, jed)
integer id_fv_diag
Definition: atmosphere.F90:239
subroutine, public atmosphere_domain(fv_domain, layout, regional, nested, pelist)
The subroutine &#39;atmosphere_domain&#39; is an API to return the "domain2d" variable associated with the co...
Definition: atmosphere.F90:875
The module &#39;fv_nesting&#39; is a collection of routines pertaining to grid nesting .
Definition: fv_nesting.F90:27
real, dimension(:,:,:), allocatable t_dt
Definition: atmosphere.F90:259