FV3 Bundle
fv_dynamics_tlm.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU General Public License *
3 !* This file is a part of fvGFS. *
4 !* *
5 !* fvGFS is free software; you can redistribute it and/or modify it *
6 !* and are expected to follow the terms of the GNU General Public *
7 !* License as published by the Free Software Foundation; either *
8 !* version 2 of the License, or (at your option) any later version. *
9 !* *
10 !* fvGFS is distributed in the hope that it will be useful, but *
11 !* WITHOUT ANY WARRANTY; without even the implied warranty of *
12 !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
13 !* General Public License for more details. *
14 !* *
15 !* For the full text of the GNU General Public License, *
16 !* write to: Free Software Foundation, Inc., *
17 !* 675 Mass Ave, Cambridge, MA 02139, USA. *
18 !* or see: http://www.gnu.org/licenses/gpl.html *
19 !***********************************************************************
20 module fv_dynamics_tlm_mod
21  use constants_mod, only: grav, pi=>pi_8, radius, hlv, rdgas, omega, rvgas, cp_vapor
22  use dyn_core_tlm_mod, only: dyn_core, del2_cubed, init_ijk_mem
23  use dyn_core_tlm_mod, only: dyn_core_tlm, del2_cubed_tlm
24  use fv_mapz_tlm_mod, only: compute_total_energy, lagrangian_to_eulerian, moist_cv, moist_cp
25  use fv_mapz_tlm_mod, only: compute_total_energy_tlm, lagrangian_to_eulerian_tlm
26  use fv_tracer2d_tlm_mod, only: tracer_2d, tracer_2d_1l, tracer_2d_nested
27  use fv_tracer2d_tlm_mod, only: tracer_2d_tlm, tracer_2d_1l_tlm, tracer_2d_nested_tlm
28  use fv_grid_utils_tlm_mod, only: cubed_to_latlon, c2l_ord2, g_sum
29  use fv_grid_utils_tlm_mod, only: c2l_ord2_tlm, g_sum_tlm
30  use fv_fill_nlm_mod, only: fill2d
31  use fv_mp_nlm_mod, only: is_master
32  use fv_mp_nlm_mod, only: group_halo_update_type
33  use fv_mp_tlm_mod, only: start_group_halo_update, complete_group_halo_update
34  use fv_mp_tlm_mod, only: start_group_halo_update_tlm
35  use fv_timing_nlm_mod, only: timing_on, timing_off
36  use diag_manager_mod, only: send_data
37  use fv_diagnostics_nlm_mod, only: fv_time, prt_mxm, range_check, prt_minmax
38  use mpp_domains_mod, only: mpp_update_domains, dgrid_ne, cgrid_ne, domain2d
39  use fv_mp_tlm_mod, only: mpp_update_domains_tlm
40  use mpp_mod, only: mpp_pe
41  use field_manager_mod, only: model_atmos
42  use tracer_manager_mod, only: get_tracer_index
43  use fv_sg_nlm_mod, only: neg_adj3
44  use fv_nesting_tlm_mod, only: setup_nested_grid_bcs
45  use fv_nesting_tlm_mod, only: setup_nested_grid_bcs_tlm
46  use boundary_tlm_mod, only: nested_grid_bc_apply_intt
47  use boundary_tlm_mod, only: nested_grid_bc_apply_intt_tlm
48  use fv_arrays_nlm_mod, only: fv_grid_type, fv_flags_type, fv_atmos_type, fv_nest_type, fv_diag_type, fv_grid_bounds_type
49  use fv_nwp_nudge_nlm_mod, only: do_adiabatic_init
50 !#ifdef MAPL_MODE
51 ! use fv_control_nlm_mod, only: dyn_timer, comm_timer
52 !#endif
53  use fv_arrays_tlmadm_mod, only: fv_flags_pert_type, fpp
54 
55 implicit none
56 
57 !#ifdef MAPL_MODE
58 ! ! Include the MPI library definitons:
59 ! include 'mpif.h'
60 !#endif
61 
62  logical :: rf_initialized = .false.
63  logical :: pt_initialized = .false.
64  logical :: bad_range = .false.
65  real, allocatable :: rf(:)
66  integer :: kmax=1
67  real :: agrav
68  logical, public, save :: idealtest=.false.
69 #ifdef HIWPP
70  real, allocatable:: u00(:,:,:), v00(:,:,:)
71 #endif
72 private
73 public :: fv_dynamics, fv_dynamics_tlm
74 
75 CONTAINS
76 ! Differentiation of fv_dynamics in forward (tangent) mode:
77 ! variations of useful results: q u v w delp delz pt
78 ! with respect to varying inputs: q u v w delp delz pt
79 ! RW status of diff variables: peln:(loc) q:in-out u:in-out v:in-out
80 ! w:in-out delp:in-out ua:(loc) uc:(loc) mfx:(loc)
81 ! delz:in-out mfy:(loc) omga:(loc) va:(loc) vc:(loc)
82 ! pkz:(loc) pe:(loc) pk:(loc) ps:(loc) pt:in-out
83 ! cx:(loc) cy:(loc)
84 !-----------------------------------------------------------------------
85 ! fv_dynamics :: FV dynamical core driver
86 !-----------------------------------------------------------------------
87  SUBROUTINE fv_dynamics_tlm(npx, npy, npz, nq_tot, ng, bdt, consv_te, &
88 & fill, reproduce_sum, kappa, cp_air, zvir, ptop, ks, ncnst, n_split, &
89 & q_split, u, u_tl, v, v_tl, w, w_tl, delz, delz_tl, hydrostatic, pt, &
90 & pt_tl, delp, delp_tl, q, q_tl, ps, ps_tl, pe, pe_tl, pk, pk_tl, peln&
91 & , peln_tl, pkz, pkz_tl, phis, q_con, omga, omga_tl, ua, ua_tl, va, &
92 & va_tl, uc, uc_tl, vc, vc_tl, ak, bk, mfx, mfx_tl, mfy, mfy_tl, cx, &
93 & cx_tl, cy, cy_tl, ze0, hybrid_z, gridstruct, flagstruct, flagstructp&
94 & , neststruct, idiag, bd, parent_grid, domain, time_total)
95  IMPLICIT NONE
96 ! Large time-step
97  REAL, INTENT(IN) :: bdt
98  REAL, INTENT(IN) :: consv_te
99  REAL, INTENT(IN) :: kappa, cp_air
100  REAL, INTENT(IN) :: zvir, ptop
101  REAL, INTENT(IN), OPTIONAL :: time_total
102  INTEGER, INTENT(IN) :: npx
103  INTEGER, INTENT(IN) :: npy
104  INTEGER, INTENT(IN) :: npz
105 ! transported tracers
106  INTEGER, INTENT(IN) :: nq_tot
107  INTEGER, INTENT(IN) :: ng
108  INTEGER, INTENT(IN) :: ks
109  INTEGER, INTENT(IN) :: ncnst
110 ! small-step horizontal dynamics
111  INTEGER, INTENT(IN) :: n_split
112 ! tracer
113  INTEGER, INTENT(IN) :: q_split
114  LOGICAL, INTENT(IN) :: fill
115  LOGICAL, INTENT(IN) :: reproduce_sum
116  LOGICAL, INTENT(IN) :: hydrostatic
117 ! Using hybrid_z for remapping
118  LOGICAL, INTENT(IN) :: hybrid_z
119  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
120 ! D grid zonal wind (m/s)
121  REAL, DIMENSION(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz), INTENT(INOUT) &
122 & :: u
123  REAL, DIMENSION(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz), INTENT(INOUT) &
124 & :: u_tl
125 ! D grid meridional wind (m/s)
126  REAL, DIMENSION(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz), INTENT(INOUT) &
127 & :: v
128  REAL, DIMENSION(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz), INTENT(INOUT) &
129 & :: v_tl
130 ! W (m/s)
131  REAL, INTENT(INOUT) :: w(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
132  REAL, INTENT(INOUT) :: w_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
133 ! temperature (K)
134  REAL, INTENT(INOUT) :: pt(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
135  REAL, INTENT(INOUT) :: pt_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
136 ! pressure thickness (pascal)
137  REAL, INTENT(INOUT) :: delp(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
138  REAL, INTENT(INOUT) :: delp_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
139 ! specific humidity and constituents
140  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, ncnst)
141  REAL, INTENT(INOUT) :: q_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz, ncnst&
142 & )
143 ! delta-height (m); non-hydrostatic only
144  REAL, INTENT(INOUT) :: delz(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
145  REAL, INTENT(INOUT) :: delz_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
146 ! height at edges (m); non-hydrostatic
147  REAL, INTENT(INOUT) :: ze0(bd%is:bd%is, bd%js:bd%js, 1)
148 ! ze0 no longer used
149 !-----------------------------------------------------------------------
150 ! Auxilliary pressure arrays:
151 ! The 5 vars below can be re-computed from delp and ptop.
152 !-----------------------------------------------------------------------
153 ! dyn_aux:
154 ! Surface pressure (pascal)
155  REAL, INTENT(INOUT) :: ps(bd%isd:bd%ied, bd%jsd:bd%jed)
156  REAL, INTENT(INOUT) :: ps_tl(bd%isd:bd%ied, bd%jsd:bd%jed)
157 ! edge pressure (pascal)
158  REAL, INTENT(INOUT) :: pe(bd%is-1:bd%ie+1, npz+1, bd%js-1:bd%je+1)
159  REAL, INTENT(INOUT) :: pe_tl(bd%is-1:bd%ie+1, npz+1, bd%js-1:bd%je+1&
160 & )
161 ! pe**kappa
162  REAL, INTENT(INOUT) :: pk(bd%is:bd%ie, bd%js:bd%je, npz+1)
163  REAL, INTENT(INOUT) :: pk_tl(bd%is:bd%ie, bd%js:bd%je, npz+1)
164 ! ln(pe)
165  REAL, INTENT(INOUT) :: peln(bd%is:bd%ie, npz+1, bd%js:bd%je)
166  REAL, INTENT(INOUT) :: peln_tl(bd%is:bd%ie, npz+1, bd%js:bd%je)
167 ! finite-volume mean pk
168  REAL, INTENT(INOUT) :: pkz(bd%is:bd%ie, bd%js:bd%je, npz)
169  REAL, INTENT(INOUT) :: pkz_tl(bd%is:bd%ie, bd%js:bd%je, npz)
170  REAL, INTENT(INOUT) :: q_con(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
171 !-----------------------------------------------------------------------
172 ! Others:
173 !-----------------------------------------------------------------------
174 ! Surface geopotential (g*Z_surf)
175  REAL, INTENT(INOUT) :: phis(bd%isd:bd%ied, bd%jsd:bd%jed)
176 ! Vertical pressure velocity (pa/s)
177  REAL, INTENT(INOUT) :: omga(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
178  REAL, INTENT(INOUT) :: omga_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
179 ! (uc,vc) mostly used as the C grid winds
180  REAL, INTENT(INOUT) :: uc(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
181  REAL, INTENT(INOUT) :: uc_tl(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
182  REAL, INTENT(INOUT) :: vc(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
183  REAL, INTENT(INOUT) :: vc_tl(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
184  REAL, DIMENSION(bd%isd:bd%ied, bd%jsd:bd%jed, npz), INTENT(INOUT) ::&
185 & ua, va
186  REAL, DIMENSION(bd%isd:bd%ied, bd%jsd:bd%jed, npz), INTENT(INOUT) ::&
187 & ua_tl, va_tl
188  REAL, DIMENSION(npz+1), INTENT(IN) :: ak, bk
189 ! Accumulated Mass flux arrays: the "Flux Capacitor"
190  REAL, INTENT(INOUT) :: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz)
191  REAL, INTENT(INOUT) :: mfx_tl(bd%is:bd%ie+1, bd%js:bd%je, npz)
192  REAL, INTENT(INOUT) :: mfy(bd%is:bd%ie, bd%js:bd%je+1, npz)
193  REAL, INTENT(INOUT) :: mfy_tl(bd%is:bd%ie, bd%js:bd%je+1, npz)
194 ! Accumulated Courant number arrays
195  REAL, INTENT(INOUT) :: cx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
196  REAL, INTENT(INOUT) :: cx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
197  REAL, INTENT(INOUT) :: cy(bd%isd:bd%ied, bd%js:bd%je+1, npz)
198  REAL, INTENT(INOUT) :: cy_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
199  TYPE(fv_grid_type), INTENT(INOUT), TARGET :: gridstruct
200  TYPE(fv_flags_type), INTENT(INOUT) :: flagstruct
201  TYPE(fv_flags_pert_type), INTENT(INOUT) :: flagstructp
202  TYPE(fv_nest_type), INTENT(INOUT) :: neststruct
203  TYPE(domain2d), INTENT(INOUT) :: domain
204  TYPE(fv_atmos_type), INTENT(INOUT) :: parent_grid
205  TYPE(fv_diag_type), INTENT(IN) :: idiag
206 ! Local Arrays
207  REAL :: ws(bd%is:bd%ie, bd%js:bd%je)
208  REAL :: ws_tl(bd%is:bd%ie, bd%js:bd%je)
209  REAL :: te_2d(bd%is:bd%ie, bd%js:bd%je)
210  REAL :: te_2d_tl(bd%is:bd%ie, bd%js:bd%je)
211  REAL :: teq(bd%is:bd%ie, bd%js:bd%je)
212  REAL :: teq_tl(bd%is:bd%ie, bd%js:bd%je)
213  REAL :: ps2(bd%isd:bd%ied, bd%jsd:bd%jed)
214  REAL :: ps2_tl(bd%isd:bd%ied, bd%jsd:bd%jed)
215  REAL :: m_fac(bd%is:bd%ie, bd%js:bd%je)
216  REAL :: m_fac_tl(bd%is:bd%ie, bd%js:bd%je)
217  REAL :: pfull(npz)
218  REAL, DIMENSION(bd%is:bd%ie) :: cvm
219  REAL :: dp1(bd%isd:bd%ied, bd%jsd:bd%jed, npz), dtdt_m(bd%is:bd%ie, &
220 & bd%js:bd%je, npz), cappa(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
221  REAL :: dp1_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
222  REAL(kind=8) :: psx(bd%isd:bd%ied, bd%jsd:bd%jed)
223  REAL(kind=8) :: psx_tl(bd%isd:bd%ied, bd%jsd:bd%jed)
224  REAL(kind=8) :: dpx(bd%is:bd%ie, bd%js:bd%je)
225  REAL(kind=8) :: dpx_tl(bd%is:bd%ie, bd%js:bd%je)
226  REAL :: akap, rdg, ph1, ph2, mdt, gam, amdt, u0
227  REAL :: amdt_tl, u0_tl
228  INTEGER :: kord_tracer(ncnst), kord_mt, kord_wz, kord_tm
229  INTEGER :: kord_tracer_pert(ncnst), kord_mt_pert, kord_wz_pert, &
230 & kord_tm_pert
231  INTEGER :: i, j, k, n, iq, n_map, nq, nwat, k_split
232 ! GFDL physics
233  INTEGER :: sphum
234  INTEGER, SAVE :: liq_wat=-999
235  INTEGER, SAVE :: ice_wat=-999
236  INTEGER, SAVE :: rainwat=-999
237  INTEGER, SAVE :: snowwat=-999
238  INTEGER, SAVE :: graupel=-999
239  INTEGER, SAVE :: cld_amt=-999
240  INTEGER, SAVE :: theta_d=-999
241  LOGICAL :: used, last_step, do_omega
242  INTEGER, PARAMETER :: max_packs=12
243  TYPE(group_halo_update_type), SAVE :: i_pack(max_packs)
244  INTEGER :: is, ie, js, je
245  INTEGER :: isd, ied, jsd, jed
246  REAL :: dt2
247  REAL(kind=8) :: t1, t2
248  INTEGER :: status
249  REAL :: rf(npz)
250  REAL :: gz(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
251  REAL :: gz_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
252  REAL :: pkc(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
253  REAL :: pkc_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
254  REAL :: ptc(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
255  REAL :: ptc_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
256  REAL :: crx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
257  REAL :: crx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
258  REAL :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
259  REAL :: xfx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
260  REAL :: cry(bd%isd:bd%ied, bd%js:bd%je+1, npz)
261  REAL :: cry_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
262  REAL :: yfx(bd%isd:bd%ied, bd%js:bd%je+1, npz)
263  REAL :: yfx_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
264  REAL :: divgd(bd%isd:bd%ied+1, bd%jsd:bd%jed+1, npz)
265  REAL :: divgd_tl(bd%isd:bd%ied+1, bd%jsd:bd%jed+1, npz)
266  REAL :: delpc(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
267  REAL :: delpc_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
268  REAL :: ut(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
269  REAL :: ut_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
270  REAL :: vt(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
271  REAL :: vt_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
272  REAL :: zh(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
273  REAL :: zh_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
274  REAL :: pk3(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
275  REAL :: pk3_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
276  REAL :: du(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
277  REAL :: du_tl(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
278  REAL :: dv(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
279  REAL :: dv_tl(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
280  INTRINSIC any
281  INTRINSIC log
282  INTRINSIC exp
283  INTRINSIC abs
284  INTRINSIC real
285  INTRINSIC cos
286  REAL :: abs0
287  REAL :: abs1
288  REAL :: arg1
289  REAL :: arg1_tl
290  REAL :: arg2
291  REAL :: arg2_tl
292  REAL :: result1
293  REAL :: result1_tl
294  gz = 0.0
295  pkc = 0.0
296  ptc = 0.0
297  crx = 0.0
298  xfx = 0.0
299  cry = 0.0
300  yfx = 0.0
301  divgd = 0.0
302  delpc = 0.0
303  ut = 0.0
304  vt = 0.0
305  zh = 0.0
306  pk3 = 0.0
307  du = 0.0
308  dv = 0.0
309  is = bd%is
310  ie = bd%ie
311  js = bd%js
312  je = bd%je
313  isd = bd%isd
314  ied = bd%ied
315  jsd = bd%jsd
316  jed = bd%jed
317 ! dyn_timer = 0
318 ! comm_timer = 0
319 ! cv_air = cp_air - rdgas
320  agrav = 1./grav
321  dt2 = 0.5*bdt
322  k_split = flagstruct%k_split
323  nwat = flagstruct%nwat
324  nq = nq_tot - flagstruct%dnats
325  rdg = -(rdgas*agrav)
326 !allocate ( dp1(isd:ied, jsd:jed, 1:npz) )
327 ! Begin Dynamics timer for GEOS history processing
328 !t1 = MPI_Wtime(status)
329  t1 = 0.0
330  t2 = 0.0
331 !allocate ( cappa(isd:isd,jsd:jsd,1) )
332  cappa = 0.
333 !We call this BEFORE converting pt to virtual potential temperature,
334 !since we interpolate on (regular) temperature rather than theta.
335  IF (gridstruct%nested .OR. any(neststruct%child_grids)) THEN
336  CALL timing_on('NEST_BCs')
337  CALL setup_nested_grid_bcs_tlm(npx, npy, npz, zvir, ncnst, u, u_tl&
338 & , v, v_tl, w, pt, delp, delz, q, uc, &
339 & uc_tl, vc, vc_tl, pkz, neststruct%nested&
340 & , flagstruct%inline_q, flagstruct%make_nh&
341 & , ng, gridstruct, flagstruct, neststruct&
342 & , neststruct%nest_timestep, neststruct%&
343 & tracer_nest_timestep, domain, bd, nwat)
344  IF (gridstruct%nested) CALL nested_grid_bc_apply_intt_tlm(pt, &
345 & pt_tl, 0, 0, &
346 & npx, npy, npz&
347 & , bd, 1., 1., &
348 & neststruct%&
349 & pt_bc, bctype=&
350 & neststruct%&
351 & nestbctype)
352 !Correct halo values have now been set up for BCs; we can go ahead and apply them too...
353  CALL timing_off('NEST_BCs')
354  ELSE
355  uc_tl = 0.0
356  vc_tl = 0.0
357  END IF
358  IF (flagstruct%no_dycore) THEN
359  IF (nwat .EQ. 2 .AND. (.NOT.hydrostatic)) sphum = get_tracer_index&
360 & (model_atmos, 'sphum')
361  END IF
362 !goto 911
363  IF (fpp%fpp_mapl_mode) THEN
364  SELECT CASE (nwat)
365  CASE (0)
366  sphum = 1
367 ! to cause trouble if (mis)used
368  cld_amt = -1
369  CASE (1)
370  sphum = 1
371 ! to cause trouble if (mis)used
372  liq_wat = -1
373 ! to cause trouble if (mis)used
374  ice_wat = -1
375 ! to cause trouble if (mis)used
376  rainwat = -1
377 ! to cause trouble if (mis)used
378  snowwat = -1
379 ! to cause trouble if (mis)used
380  graupel = -1
381 ! to cause trouble if (mis)used
382  cld_amt = -1
383 ! to cause trouble if (mis)used
384  theta_d = -1
385  CASE (3)
386  sphum = 1
387  liq_wat = 2
388  ice_wat = 3
389 ! to cause trouble if (mis)used
390  rainwat = -1
391 ! to cause trouble if (mis)used
392  snowwat = -1
393 ! to cause trouble if (mis)used
394  graupel = -1
395 ! to cause trouble if (mis)used
396  cld_amt = -1
397 ! to cause trouble if (mis)used
398  theta_d = -1
399  END SELECT
400  ELSE
401  IF (nwat .EQ. 0) THEN
402  sphum = 1
403 ! to cause trouble if (mis)used
404  cld_amt = -1
405  ELSE
406  sphum = get_tracer_index(model_atmos, 'sphum')
407  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
408  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
409  rainwat = get_tracer_index(model_atmos, 'rainwat')
410  snowwat = get_tracer_index(model_atmos, 'snowwat')
411  graupel = get_tracer_index(model_atmos, 'graupel')
412  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
413  END IF
414  theta_d = get_tracer_index(model_atmos, 'theta_d')
415  END IF
416  akap = kappa
417 !$OMP parallel do default(none) shared(npz,ak,bk,flagstruct,pfull) &
418 !$OMP private(ph1, ph2)
419  DO k=1,npz
420  ph1 = ak(k) + bk(k)*flagstruct%p_ref
421  ph2 = ak(k+1) + bk(k+1)*flagstruct%p_ref
422  pfull(k) = (ph2-ph1)/log(ph2/ph1)
423  END DO
424  IF (hydrostatic) THEN
425  dp1_tl = 0.0
426 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,nwat,q,q_con,sphum,liq_wat, &
427 !$OMP rainwat,ice_wat,snowwat,graupel) private(cvm)
428  DO k=1,npz
429  DO j=js,je
430  DO i=is,ie
431  dp1_tl(i, j, k) = zvir*q_tl(i, j, k, sphum)
432  dp1(i, j, k) = zvir*q(i, j, k, sphum)
433  END DO
434  END DO
435  END DO
436  ELSE
437  dp1_tl = 0.0
438 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,q,q_con,sphum,liq_wat, &
439 !$OMP rainwat,ice_wat,snowwat,graupel,pkz,flagstruct, &
440 !$OMP cappa,kappa,rdg,delp,pt,delz,nwat) &
441 !$OMP private(cvm)
442  DO k=1,npz
443  IF (flagstruct%moist_phys) THEN
444  DO j=js,je
445  DO i=is,ie
446  dp1_tl(i, j, k) = zvir*q_tl(i, j, k, sphum)
447  dp1(i, j, k) = zvir*q(i, j, k, sphum)
448  arg1_tl = (rdg*((delp_tl(i, j, k)*pt(i, j, k)+delp(i, j, k&
449 & )*pt_tl(i, j, k))*(1.+dp1(i, j, k))+delp(i, j, k)*pt(i, &
450 & j, k)*dp1_tl(i, j, k))*delz(i, j, k)-rdg*delp(i, j, k)*&
451 & pt(i, j, k)*(1.+dp1(i, j, k))*delz_tl(i, j, k))/delz(i, &
452 & j, k)**2
453  arg1 = rdg*delp(i, j, k)*pt(i, j, k)*(1.+dp1(i, j, k))/&
454 & delz(i, j, k)
455  arg2_tl = kappa*arg1_tl/arg1
456  arg2 = kappa*log(arg1)
457  pkz_tl(i, j, k) = arg2_tl*exp(arg2)
458  pkz(i, j, k) = exp(arg2)
459  END DO
460  END DO
461  ELSE
462 ! Using dry pressure for the definition of the virtual potential temperature
463 ! pkz(i,j,k) = exp( kappa*log(rdg*delp(i,j,k)*pt(i,j,k)* &
464 ! (1.-q(i,j,k,sphum))/delz(i,j,k)) )
465  DO j=js,je
466  DO i=is,ie
467  dp1_tl(i, j, k) = 0.0
468  dp1(i, j, k) = 0.
469  arg1_tl = (rdg*(delp_tl(i, j, k)*pt(i, j, k)+delp(i, j, k)&
470 & *pt_tl(i, j, k))*delz(i, j, k)-rdg*delp(i, j, k)*pt(i, j&
471 & , k)*delz_tl(i, j, k))/delz(i, j, k)**2
472  arg1 = rdg*delp(i, j, k)*pt(i, j, k)/delz(i, j, k)
473  arg2_tl = kappa*arg1_tl/arg1
474  arg2 = kappa*log(arg1)
475  pkz_tl(i, j, k) = arg2_tl*exp(arg2)
476  pkz(i, j, k) = exp(arg2)
477  END DO
478  END DO
479  END IF
480  END DO
481  END IF
482  IF (flagstruct%fv_debug) THEN
483  CALL prt_mxm('PS', ps, is, ie, js, je, ng, 1, 0.01, gridstruct%&
484 & area_64, domain)
485  CALL prt_mxm('T_dyn_b', pt, is, ie, js, je, ng, npz, 1., &
486 & gridstruct%area_64, domain)
487  IF (.NOT.hydrostatic) CALL prt_mxm('delz', delz, is, ie, js, je, &
488 & ng, npz, 1., gridstruct%area_64, &
489 & domain)
490  CALL prt_mxm('delp_b ', delp, is, ie, js, je, ng, npz, 0.01, &
491 & gridstruct%area_64, domain)
492  CALL prt_mxm('pk_b', pk, is, ie, js, je, 0, npz + 1, 1., &
493 & gridstruct%area_64, domain)
494  CALL prt_mxm('pkz_b', pkz, is, ie, js, je, 0, npz, 1., gridstruct%&
495 & area_64, domain)
496  END IF
497 !---------------------
498 ! Compute Total Energy
499 !---------------------
500  IF (consv_te .GT. 0. .AND. (.NOT.do_adiabatic_init)) THEN
501  CALL compute_total_energy_tlm(is, ie, js, je, isd, ied, jsd, jed, &
502 & npz, u, u_tl, v, v_tl, w, w_tl, delz, &
503 & delz_tl, pt, pt_tl, delp, delp_tl, q, q_tl&
504 & , dp1, dp1_tl, pe, pe_tl, peln, peln_tl, &
505 & phis, gridstruct%rsin2, gridstruct%cosa_s&
506 & , zvir, cp_air, rdgas, hlv, te_2d, &
507 & te_2d_tl, ua, va, teq, teq_tl, flagstruct%&
508 & moist_phys, nwat, sphum, liq_wat, rainwat&
509 & , ice_wat, snowwat, graupel, hydrostatic, &
510 & idiag%id_te)
511  IF (idiag%id_te .GT. 0) used = send_data(idiag%id_te, teq, fv_time&
512 & )
513 ! te_den=1.E-9*g_sum(teq, is, ie, js, je, ng, area, 0)/(grav*4.*pi*radius**2)
514 ! if(is_master()) write(*,*) 'Total Energy Density (Giga J/m**2)=',te_den
515  ELSE
516  teq_tl = 0.0
517  te_2d_tl = 0.0
518  END IF
519  IF ((flagstruct%consv_am .OR. idiag%id_amdt .GT. 0) .AND. (.NOT.&
520 & do_adiabatic_init)) THEN
521  m_fac_tl = 0.0
522  ps2_tl = 0.0
523  va_tl = 0.0
524  ua_tl = 0.0
525  CALL compute_aam_tlm(npz, is, ie, js, je, isd, ied, jsd, jed, &
526 & gridstruct, bd, ptop, ua, ua_tl, va, va_tl, u, u_tl&
527 & , v, v_tl, delp, delp_tl, teq, teq_tl, ps2, ps2_tl&
528 & , m_fac, m_fac_tl)
529  ELSE
530  ua_tl = 0.0
531  va_tl = 0.0
532  ps2_tl = 0.0
533  m_fac_tl = 0.0
534  END IF
535  IF (flagstruct%tau .GT. 0.) THEN
536  IF (gridstruct%grid_type .LT. 4) THEN
537  IF (bdt .GE. 0.) THEN
538  abs0 = bdt
539  ELSE
540  abs0 = -bdt
541  END IF
542  CALL rayleigh_super_tlm(abs0, npx, npy, npz, ks, pfull, phis, &
543 & flagstruct%tau, u, u_tl, v, v_tl, w, w_tl, pt&
544 & , pt_tl, ua, ua_tl, va, va_tl, delz, &
545 & gridstruct%agrid, cp_air, rdgas, ptop, &
546 & hydrostatic, .NOT.neststruct%nested, &
547 & flagstruct%rf_cutoff, rf, gridstruct, domain, &
548 & bd)
549  ELSE
550  IF (bdt .GE. 0.) THEN
551  abs1 = bdt
552  ELSE
553  abs1 = -bdt
554  END IF
555  CALL rayleigh_friction_tlm(abs1, npx, npy, npz, ks, pfull, &
556 & flagstruct%tau, u, u_tl, v, v_tl, w, w_tl, &
557 & pt, pt_tl, ua, ua_tl, va, va_tl, delz, &
558 & delz_tl, cp_air, rdgas, ptop, hydrostatic, &
559 & .true., flagstruct%rf_cutoff, rf, &
560 & gridstruct, domain, bd)
561  END IF
562  END IF
563 ! Convert pt to virtual potential temperature on the first timestep
564  IF (flagstruct%adiabatic) THEN
565 !$OMP parallel do default(none) shared(theta_d,is,ie,js,je,npz,pt,pkz,q)
566  DO k=1,npz
567  DO j=js,je
568  DO i=is,ie
569  pt_tl(i, j, k) = (pt_tl(i, j, k)*pkz(i, j, k)-pt(i, j, k)*&
570 & pkz_tl(i, j, k))/pkz(i, j, k)**2
571  pt(i, j, k) = pt(i, j, k)/pkz(i, j, k)
572  END DO
573  END DO
574  IF (theta_d .GT. 0) THEN
575  DO j=js,je
576  DO i=is,ie
577  q_tl(i, j, k, theta_d) = pt_tl(i, j, k)
578  q(i, j, k, theta_d) = pt(i, j, k)
579  END DO
580  END DO
581  END IF
582  END DO
583  ELSE
584 !$OMP parallel do default(none) shared(is,ie,js,je,npz,pt,dp1,pkz,q_con)
585  DO k=1,npz
586  DO j=js,je
587  DO i=is,ie
588  pt_tl(i, j, k) = ((pt_tl(i, j, k)*(1.+dp1(i, j, k))+pt(i, j&
589 & , k)*dp1_tl(i, j, k))*pkz(i, j, k)-pt(i, j, k)*(1.+dp1(i, &
590 & j, k))*pkz_tl(i, j, k))/pkz(i, j, k)**2
591  pt(i, j, k) = pt(i, j, k)*(1.+dp1(i, j, k))/pkz(i, j, k)
592  END DO
593  END DO
594  END DO
595  END IF
596  last_step = .false.
597  mdt = bdt/REAL(k_split)
598  IF (idiag%id_mdt .GT. 0 .AND. (.NOT.do_adiabatic_init)) THEN
599 !allocate ( dtdt_m(is:ie,js:je,npz) )
600 !$OMP parallel do default(none) shared(is,ie,js,je,npz,dtdt_m)
601  DO k=1,npz
602  DO j=js,je
603  DO i=is,ie
604  dtdt_m(i, j, k) = 0.
605  END DO
606  END DO
607  END DO
608  END IF
609 !DryMassRoundoffControl
610 !allocate(psx(isd:ied,jsd:jed),dpx(is:ie,js:je))
611  IF (fpp%fpp_overload_r4) THEN
612  psx_tl = 0.0_8
613  DO j=js,je
614  DO i=is,ie
615  psx_tl(i, j) = pe_tl(i, npz+1, j)
616  psx(i, j) = pe(i, npz+1, j)
617  dpx(i, j) = 0.0
618  END DO
619  END DO
620  ELSE
621  psx_tl = 0.0_8
622  END IF
623  CALL timing_on('FV_DYN_LOOP')
624  omga_tl = 0.0
625  ps_tl = 0.0
626  pk3_tl = 0.0
627  xfx_tl = 0.0
628  ws_tl = 0.0
629  gz_tl = 0.0
630  du_tl = 0.0
631  dv_tl = 0.0
632  ptc_tl = 0.0
633  ut_tl = 0.0
634  divgd_tl = 0.0
635  pkc_tl = 0.0
636  delpc_tl = 0.0
637  yfx_tl = 0.0
638  vt_tl = 0.0
639  zh_tl = 0.0
640  dpx_tl = 0.0_8
641  crx_tl = 0.0
642  cry_tl = 0.0
643 ! first level of time-split
644  DO n_map=1,k_split
645  CALL timing_on('COMM_TOTAL')
646  CALL start_group_halo_update_tlm(i_pack(1), delp, delp_tl, domain&
647 & , complete=.true.)
648  CALL start_group_halo_update_tlm(i_pack(1), pt, pt_tl, domain, &
649 & complete=.true.)
650  CALL start_group_halo_update_tlm(i_pack(8), u, u_tl, v, v_tl, &
651 & domain, gridtype=dgrid_ne)
652  CALL timing_off('COMM_TOTAL')
653 !$OMP parallel do default(none) shared(isd,ied,jsd,jed,npz,dp1,delp)
654  DO k=1,npz
655  DO j=jsd,jed
656  DO i=isd,ied
657  dp1_tl(i, j, k) = delp_tl(i, j, k)
658  dp1(i, j, k) = delp(i, j, k)
659  END DO
660  END DO
661  END DO
662  IF (n_map .EQ. k_split) last_step = .true.
663  CALL timing_on('DYN_CORE')
664  CALL dyn_core_tlm(npx, npy, npz, ng, sphum, nq, mdt, n_split, zvir&
665 & , cp_air, akap, cappa, grav, hydrostatic, u, u_tl, v, &
666 & v_tl, w, w_tl, delz, delz_tl, pt, pt_tl, q, q_tl, delp&
667 & , delp_tl, pe, pe_tl, pk, pk_tl, phis, ws, ws_tl, omga&
668 & , omga_tl, ptop, pfull, ua, ua_tl, va, va_tl, uc, &
669 & uc_tl, vc, vc_tl, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl&
670 & , cy, cy_tl, pkz, pkz_tl, peln, peln_tl, q_con, ak, bk&
671 & , dpx, dpx_tl, ks, gridstruct, flagstruct, flagstructp&
672 & , neststruct, idiag, bd, domain, n_map .EQ. 1, i_pack&
673 & , last_step, gz, gz_tl, pkc, pkc_tl, ptc, ptc_tl, crx&
674 & , crx_tl, xfx, xfx_tl, cry, cry_tl, yfx, yfx_tl, divgd&
675 & , divgd_tl, delpc, delpc_tl, ut, ut_tl, vt, vt_tl, zh&
676 & , zh_tl, pk3, pk3_tl, du, du_tl, dv, dv_tl, time_total&
677 & )
678  CALL timing_off('DYN_CORE')
679 !DryMassRoundoffControl
680  IF (last_step) THEN
681  IF (fpp%fpp_overload_r4) THEN
682  DO j=js,je
683  DO i=is,ie
684  psx_tl(i, j) = psx_tl(i, j) + dpx_tl(i, j)
685  psx(i, j) = psx(i, j) + dpx(i, j)
686  END DO
687  END DO
688  CALL timing_on('COMM_TOTAL')
689  CALL mpp_update_domains_tlm(psx, psx_tl, domain)
690  CALL timing_off('COMM_TOTAL')
691  DO j=js-1,je+1
692  DO i=is-1,ie+1
693  pe_tl(i, npz+1, j) = psx_tl(i, j)
694  pe(i, npz+1, j) = psx(i, j)
695  END DO
696  END DO
697  END IF
698  END IF
699 !deallocate(psx,dpx)
700  IF (.NOT.flagstruct%inline_q .AND. nq .NE. 0) THEN
701 !--------------------------------------------------------
702 ! Perform large-time-step scalar transport using the accumulated CFL and
703 ! mass fluxes
704  CALL timing_on('tracer_2d')
705 !!! CLEANUP: merge these two calls?
706  IF (gridstruct%nested) THEN
707  CALL tracer_2d_nested_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, &
708 & mfy, mfy_tl, cx, cx_tl, cy, cy_tl, &
709 & gridstruct, bd, domain, npx, npy, npz, nq&
710 & , flagstruct%hord_tr, q_split, mdt, idiag%&
711 & id_divg, i_pack(10), flagstruct%nord_tr, &
712 & flagstruct%trdm2, k_split, neststruct, &
713 & parent_grid, flagstructp%hord_tr_pert, &
714 & flagstructp%nord_tr_pert, flagstructp%&
715 & trdm2_pert, flagstructp%split_damp_tr)
716  ELSE IF (flagstruct%z_tracer) THEN
717  CALL tracer_2d_1l_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, &
718 & mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, &
719 & domain, npx, npy, npz, nq, flagstruct%hord_tr&
720 & , q_split, mdt, idiag%id_divg, i_pack(10), &
721 & flagstruct%nord_tr, flagstruct%trdm2, &
722 & flagstructp%hord_tr_pert, flagstructp%&
723 & nord_tr_pert, flagstructp%trdm2_pert, &
724 & flagstructp%split_damp_tr)
725  ELSE
726  CALL tracer_2d_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, &
727 & mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, &
728 & domain, npx, npy, npz, nq, flagstruct%hord_tr, &
729 & q_split, mdt, idiag%id_divg, i_pack(10), &
730 & flagstruct%nord_tr, flagstruct%trdm2, flagstructp&
731 & %hord_tr_pert, flagstructp%nord_tr_pert, &
732 & flagstructp%trdm2_pert, flagstructp%split_damp_tr&
733 & )
734  END IF
735  CALL timing_off('tracer_2d')
736  IF (flagstruct%hord_tr .LT. 8 .AND. flagstruct%moist_phys) THEN
737  CALL timing_on('Fill2D')
738  IF (liq_wat .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
739 & :ied, jsd:jed, 1, liq_wat), delp, &
740 & gridstruct%area, domain, neststruct%&
741 & nested, npx, npy)
742  IF (rainwat .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
743 & :ied, jsd:jed, 1, rainwat), delp, &
744 & gridstruct%area, domain, neststruct%&
745 & nested, npx, npy)
746  IF (ice_wat .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
747 & :ied, jsd:jed, 1, ice_wat), delp, &
748 & gridstruct%area, domain, neststruct%&
749 & nested, npx, npy)
750  IF (snowwat .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
751 & :ied, jsd:jed, 1, snowwat), delp, &
752 & gridstruct%area, domain, neststruct%&
753 & nested, npx, npy)
754  IF (graupel .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
755 & :ied, jsd:jed, 1, graupel), delp, &
756 & gridstruct%area, domain, neststruct%&
757 & nested, npx, npy)
758  CALL timing_off('Fill2D')
759  END IF
760  IF (last_step .AND. idiag%id_divg .GT. 0) THEN
761  used = send_data(idiag%id_divg, dp1, fv_time)
762  IF (flagstruct%fv_debug) CALL prt_mxm('divg', dp1, is, ie, js&
763 & , je, 0, npz, 1., gridstruct%&
764 & area_64, domain)
765  END IF
766  END IF
767  IF (npz .GT. 4) THEN
768 !------------------------------------------------------------------------
769 ! Perform vertical remapping from Lagrangian control-volume to
770 ! the Eulerian coordinate as specified by the routine set_eta.
771 ! Note that this finite-volume dycore is otherwise independent of the vertical
772 ! Eulerian coordinate.
773 !------------------------------------------------------------------------
774  DO iq=1,nq
775  kord_tracer(iq) = flagstruct%kord_tr
776 ! monotonic
777  IF (iq .EQ. cld_amt) kord_tracer(iq) = 9
778  kord_tracer_pert(iq) = flagstructp%kord_tr_pert
779 ! linear
780  IF (iq .EQ. cld_amt) kord_tracer_pert(iq) = 17
781  END DO
782  do_omega = hydrostatic .AND. last_step
783  CALL timing_on('Remapping')
784  kord_mt = flagstruct%kord_mt
785  kord_wz = flagstruct%kord_wz
786  kord_tm = flagstruct%kord_tm
787  kord_mt_pert = flagstructp%kord_mt_pert
788  kord_wz_pert = flagstructp%kord_wz_pert
789  kord_tm_pert = flagstructp%kord_tm_pert
790  IF (n_map .EQ. k_split) THEN
791  kord_mt = kord_mt_pert
792  kord_wz = kord_wz_pert
793  kord_tm = kord_tm_pert
794  kord_tracer = kord_tracer_pert
795  END IF
796  CALL lagrangian_to_eulerian_tlm(last_step, consv_te, ps, ps_tl, &
797 & pe, pe_tl, delp, delp_tl, pkz, pkz_tl&
798 & , pk, pk_tl, mdt, bdt, npz, is, ie, js&
799 & , je, isd, ied, jsd, jed, nq, nwat, &
800 & sphum, q_con, u, u_tl, v, v_tl, w, &
801 & w_tl, delz, delz_tl, pt, pt_tl, q, &
802 & q_tl, phis, zvir, cp_air, akap, cappa&
803 & , kord_mt, kord_wz, kord_tracer, &
804 & kord_tm, peln, peln_tl, te_2d, &
805 & te_2d_tl, ng, ua, ua_tl, va, omga, &
806 & omga_tl, dp1, dp1_tl, ws, ws_tl, fill&
807 & , reproduce_sum, idiag%id_mdt .GT. 0, &
808 & dtdt_m, ptop, ak, bk, pfull, &
809 & flagstruct, gridstruct, domain, &
810 & flagstruct%do_sat_adj, hydrostatic, &
811 & hybrid_z, do_omega, flagstruct%&
812 & adiabatic, do_adiabatic_init, mfx, mfy&
813 & , flagstruct%remap_option, &
814 & kord_mt_pert, kord_wz_pert, &
815 & kord_tracer_pert, kord_tm_pert)
816  CALL timing_off('Remapping')
817  IF (last_step) THEN
818  IF (.NOT.hydrostatic) THEN
819 !$OMP parallel do default(none) shared(is,ie,js,je,npz,omga,delp,delz,w)
820  DO k=1,npz
821  DO j=js,je
822  DO i=is,ie
823  omga_tl(i, j, k) = (delp_tl(i, j, k)*delz(i, j, k)-&
824 & delp(i, j, k)*delz_tl(i, j, k))*w(i, j, k)/delz(i, j&
825 & , k)**2 + delp(i, j, k)*w_tl(i, j, k)/delz(i, j, k)
826  omga(i, j, k) = delp(i, j, k)/delz(i, j, k)*w(i, j, k)
827  END DO
828  END DO
829  END DO
830  END IF
831 !--------------------------
832 ! Filter omega for physics:
833 !--------------------------
834  IF (flagstruct%nf_omega .GT. 0) CALL del2_cubed_tlm(omga, &
835 & omga_tl, 0.18*&
836 & gridstruct%&
837 & da_min, &
838 & gridstruct, &
839 & domain, npx, npy&
840 & , npz, &
841 & flagstruct%&
842 & nf_omega, bd)
843  END IF
844  END IF
845  END DO
846 ! n_map loop
847  CALL timing_off('FV_DYN_LOOP')
848  IF (idiag%id_mdt .GT. 0 .AND. (.NOT.do_adiabatic_init)) THEN
849 ! Output temperature tendency due to inline moist physics:
850 !$OMP parallel do default(none) shared(is,ie,js,je,npz,dtdt_m,bdt)
851  DO k=1,npz
852  DO j=js,je
853  DO i=is,ie
854  dtdt_m(i, j, k) = dtdt_m(i, j, k)/bdt*86400.
855  END DO
856  END DO
857  END DO
858 ! call prt_mxm('Fast DTDT (deg/Day)', dtdt_m, is, ie, js, je, 0, npz, 1., gridstruct%area_64, domain)
859  used = send_data(idiag%id_mdt, dtdt_m, fv_time)
860 !deallocate ( dtdt_m )
861  END IF
862  IF (nwat .EQ. 6) THEN
863  IF (cld_amt .GT. 0) THEN
864  CALL neg_adj3(is, ie, js, je, ng, npz, flagstruct%hydrostatic, &
865 & peln, delz, pt, delp, q(isd:ied, jsd:jed, 1, sphum), q(&
866 & isd:ied, jsd:jed, 1, liq_wat), q(isd:ied, jsd:jed, 1, &
867 & rainwat), q(isd:ied, jsd:jed, 1, ice_wat), q(isd:ied, &
868 & jsd:jed, 1, snowwat), q(isd:ied, jsd:jed, 1, graupel), q&
869 & (isd:ied, jsd:jed, 1, cld_amt), flagstruct%&
870 & check_negative)
871  ELSE
872  CALL neg_adj3(is, ie, js, je, ng, npz, flagstruct%hydrostatic, &
873 & peln, delz, pt, delp, q(isd:ied, jsd:jed, 1, sphum), q(&
874 & isd:ied, jsd:jed, 1, liq_wat), q(isd:ied, jsd:jed, 1, &
875 & rainwat), q(isd:ied, jsd:jed, 1, ice_wat), q(isd:ied, &
876 & jsd:jed, 1, snowwat), q(isd:ied, jsd:jed, 1, graupel), &
877 & check_negative=flagstruct%check_negative)
878  END IF
879  IF (flagstruct%fv_debug) THEN
880  CALL prt_mxm('T_dyn_a3', pt, is, ie, js, je, ng, npz, 1., &
881 & gridstruct%area_64, domain)
882  CALL prt_mxm('SPHUM_dyn', q(isd:ied, jsd:jed, 1, sphum), is, ie&
883 & , js, je, ng, npz, 1., gridstruct%area_64, domain)
884  CALL prt_mxm('liq_wat_dyn', q(isd:ied, jsd:jed, 1, liq_wat), is&
885 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
886  CALL prt_mxm('rainwat_dyn', q(isd:ied, jsd:jed, 1, rainwat), is&
887 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
888  CALL prt_mxm('ice_wat_dyn', q(isd:ied, jsd:jed, 1, ice_wat), is&
889 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
890  CALL prt_mxm('snowwat_dyn', q(isd:ied, jsd:jed, 1, snowwat), is&
891 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
892  CALL prt_mxm('graupel_dyn', q(isd:ied, jsd:jed, 1, graupel), is&
893 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
894  END IF
895  END IF
896  IF (((flagstruct%consv_am .OR. idiag%id_amdt .GT. 0) .OR. idiag%&
897 & id_aam .GT. 0) .AND. (.NOT.do_adiabatic_init)) THEN
898  CALL compute_aam_tlm(npz, is, ie, js, je, isd, ied, jsd, jed, &
899 & gridstruct, bd, ptop, ua, ua_tl, va, va_tl, u, u_tl&
900 & , v, v_tl, delp, delp_tl, te_2d, te_2d_tl, ps, &
901 & ps_tl, m_fac, m_fac_tl)
902  IF (idiag%id_aam .GT. 0) THEN
903  used = send_data(idiag%id_aam, te_2d, fv_time)
904  IF (prt_minmax) gam = g_sum(domain, te_2d, is, ie, js, je, ng, &
905 & gridstruct%area_64, 0)
906 !if( is_master() ) write(6,*) 'Total AAM =', gam
907  END IF
908  END IF
909  IF ((flagstruct%consv_am .OR. idiag%id_amdt .GT. 0) .AND. (.NOT.&
910 & do_adiabatic_init)) THEN
911 !$OMP parallel do default(none) shared(is,ie,js,je,te_2d,teq,dt2,ps2,ps,idiag)
912  DO j=js,je
913  DO i=is,ie
914 ! Note: the mountain torque computation contains also numerical error
915 ! The numerical error is mostly from the zonal gradient of the terrain (zxg)
916  te_2d_tl(i, j) = te_2d_tl(i, j) - teq_tl(i, j) + dt2*idiag%zxg&
917 & (i, j)*(ps2_tl(i, j)+ps_tl(i, j))
918  te_2d(i, j) = te_2d(i, j) - teq(i, j) + dt2*(ps2(i, j)+ps(i, j&
919 & ))*idiag%zxg(i, j)
920  END DO
921  END DO
922  IF (idiag%id_amdt .GT. 0) used = send_data(idiag%id_amdt, te_2d/&
923 & bdt, fv_time)
924  IF (flagstruct%consv_am .OR. prt_minmax) THEN
925  amdt_tl = g_sum_tlm(domain, te_2d, te_2d_tl, is, ie, js, je, ng&
926 & , gridstruct%area_64, 0, reproduce=.true., g_sum=amdt)
927  result1_tl = g_sum_tlm(domain, m_fac, m_fac_tl, is, ie, js, je, &
928 & ng, gridstruct%area_64, 0, reproduce=.true., g_sum=result1)
929  u0_tl = -((radius*amdt_tl*result1-radius*amdt*result1_tl)/&
930 & result1**2)
931  u0 = -(radius*amdt/result1)
932  IF (is_master() .AND. prt_minmax) WRITE(6, *) &
933 & 'Dynamic AM tendency (Hadleys)='&
934 & , amdt/(bdt*1.e18), &
935 & 'del-u (per day)=', u0*86400./&
936 & bdt
937  ELSE
938  u0_tl = 0.0
939  END IF
940 ! consv_am
941  IF (flagstruct%consv_am) THEN
942 !$OMP parallel do default(none) shared(is,ie,js,je,m_fac,u0,gridstruct)
943  DO j=js,je
944  DO i=is,ie
945  m_fac(i, j) = u0*cos(gridstruct%agrid(i, j, 2))
946  END DO
947  END DO
948 !$OMP parallel do default(none) shared(is,ie,js,je,npz,hydrostatic,pt,m_fac,ua,cp_air, &
949 !$OMP u,u0,gridstruct,v )
950  DO k=1,npz
951  DO j=js,je+1
952  DO i=is,ie
953  u_tl(i, j, k) = u_tl(i, j, k) + gridstruct%l2c_u(i, j)*&
954 & u0_tl
955  u(i, j, k) = u(i, j, k) + u0*gridstruct%l2c_u(i, j)
956  END DO
957  END DO
958  DO j=js,je
959  DO i=is,ie+1
960  v_tl(i, j, k) = v_tl(i, j, k) + gridstruct%l2c_v(i, j)*&
961 & u0_tl
962  v(i, j, k) = v(i, j, k) + u0*gridstruct%l2c_v(i, j)
963  END DO
964  END DO
965  END DO
966  END IF
967  END IF
968 !911 call cubed_to_latlon(u, v, ua, va, gridstruct, &
969 ! npx, npy, npz, 1, gridstruct%grid_type, domain, gridstruct%nested, flagstruct%c2l_ord, bd)
970 !deallocate(dp1)
971 !deallocate(cappa)
972  IF (flagstruct%fv_debug) THEN
973  CALL prt_mxm('UA', ua, is, ie, js, je, ng, npz, 1., gridstruct%&
974 & area_64, domain)
975  CALL prt_mxm('VA', va, is, ie, js, je, ng, npz, 1., gridstruct%&
976 & area_64, domain)
977  CALL prt_mxm('TA', pt, is, ie, js, je, ng, npz, 1., gridstruct%&
978 & area_64, domain)
979  IF (.NOT.hydrostatic) CALL prt_mxm('W ', w, is, ie, js, je, ng, &
980 & npz, 1., gridstruct%area_64, domain)
981  END IF
982  IF (flagstruct%range_warn) THEN
983  CALL range_check('UA_dyn', ua, is, ie, js, je, ng, npz, gridstruct&
984 & %agrid, -280., 280., bad_range)
985  CALL range_check('VA_dyn', ua, is, ie, js, je, ng, npz, gridstruct&
986 & %agrid, -280., 280., bad_range)
987  CALL range_check('TA_dyn', pt, is, ie, js, je, ng, npz, gridstruct&
988 & %agrid, 150., 335., bad_range)
989  IF (.NOT.hydrostatic) CALL range_check('W_dyn', w, is, ie, js, je&
990 & , ng, npz, gridstruct%agrid, -50.&
991 & , 100., bad_range)
992  END IF
993 ! IF (fpp%fpp_mapl_mode) dyn_timer = dyn_timer + (t2-t1)
994 !t2 = MPI_Wtime(status)
995  END SUBROUTINE fv_dynamics_tlm
996 !-----------------------------------------------------------------------
997 ! fv_dynamics :: FV dynamical core driver
998 !-----------------------------------------------------------------------
999  SUBROUTINE fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill&
1000 & , reproduce_sum, kappa, cp_air, zvir, ptop, ks, ncnst, n_split, &
1001 & q_split, u, v, w, delz, hydrostatic, pt, delp, q, ps, pe, pk, peln, &
1002 & pkz, phis, q_con, omga, ua, va, uc, vc, ak, bk, mfx, mfy, cx, cy, &
1003 & ze0, hybrid_z, gridstruct, flagstruct, flagstructp, neststruct, &
1004 & idiag, bd, parent_grid, domain, time_total)
1005  IMPLICIT NONE
1006 ! Large time-step
1007  REAL, INTENT(IN) :: bdt
1008  REAL, INTENT(IN) :: consv_te
1009  REAL, INTENT(IN) :: kappa, cp_air
1010  REAL, INTENT(IN) :: zvir, ptop
1011  REAL, INTENT(IN), OPTIONAL :: time_total
1012  INTEGER, INTENT(IN) :: npx
1013  INTEGER, INTENT(IN) :: npy
1014  INTEGER, INTENT(IN) :: npz
1015 ! transported tracers
1016  INTEGER, INTENT(IN) :: nq_tot
1017  INTEGER, INTENT(IN) :: ng
1018  INTEGER, INTENT(IN) :: ks
1019  INTEGER, INTENT(IN) :: ncnst
1020 ! small-step horizontal dynamics
1021  INTEGER, INTENT(IN) :: n_split
1022 ! tracer
1023  INTEGER, INTENT(IN) :: q_split
1024  LOGICAL, INTENT(IN) :: fill
1025  LOGICAL, INTENT(IN) :: reproduce_sum
1026  LOGICAL, INTENT(IN) :: hydrostatic
1027 ! Using hybrid_z for remapping
1028  LOGICAL, INTENT(IN) :: hybrid_z
1029  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
1030 ! D grid zonal wind (m/s)
1031  REAL, DIMENSION(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz), INTENT(INOUT) &
1032 & :: u
1033 ! D grid meridional wind (m/s)
1034  REAL, DIMENSION(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz), INTENT(INOUT) &
1035 & :: v
1036 ! W (m/s)
1037  REAL, INTENT(INOUT) :: w(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1038 ! temperature (K)
1039  REAL, INTENT(INOUT) :: pt(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1040 ! pressure thickness (pascal)
1041  REAL, INTENT(INOUT) :: delp(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1042 ! specific humidity and constituents
1043  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, ncnst)
1044 ! delta-height (m); non-hydrostatic only
1045  REAL, INTENT(INOUT) :: delz(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1046 ! height at edges (m); non-hydrostatic
1047  REAL, INTENT(INOUT) :: ze0(bd%is:bd%is, bd%js:bd%js, 1)
1048 ! ze0 no longer used
1049 !-----------------------------------------------------------------------
1050 ! Auxilliary pressure arrays:
1051 ! The 5 vars below can be re-computed from delp and ptop.
1052 !-----------------------------------------------------------------------
1053 ! dyn_aux:
1054 ! Surface pressure (pascal)
1055  REAL, INTENT(INOUT) :: ps(bd%isd:bd%ied, bd%jsd:bd%jed)
1056 ! edge pressure (pascal)
1057  REAL, INTENT(INOUT) :: pe(bd%is-1:bd%ie+1, npz+1, bd%js-1:bd%je+1)
1058 ! pe**kappa
1059  REAL, INTENT(INOUT) :: pk(bd%is:bd%ie, bd%js:bd%je, npz+1)
1060 ! ln(pe)
1061  REAL, INTENT(INOUT) :: peln(bd%is:bd%ie, npz+1, bd%js:bd%je)
1062 ! finite-volume mean pk
1063  REAL, INTENT(INOUT) :: pkz(bd%is:bd%ie, bd%js:bd%je, npz)
1064  REAL, INTENT(INOUT) :: q_con(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1065 !-----------------------------------------------------------------------
1066 ! Others:
1067 !-----------------------------------------------------------------------
1068 ! Surface geopotential (g*Z_surf)
1069  REAL, INTENT(INOUT) :: phis(bd%isd:bd%ied, bd%jsd:bd%jed)
1070 ! Vertical pressure velocity (pa/s)
1071  REAL, INTENT(INOUT) :: omga(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1072 ! (uc,vc) mostly used as the C grid winds
1073  REAL, INTENT(INOUT) :: uc(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
1074  REAL, INTENT(INOUT) :: vc(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
1075  REAL, DIMENSION(bd%isd:bd%ied, bd%jsd:bd%jed, npz), INTENT(INOUT) ::&
1076 & ua, va
1077  REAL, DIMENSION(npz+1), INTENT(IN) :: ak, bk
1078 ! Accumulated Mass flux arrays: the "Flux Capacitor"
1079  REAL, INTENT(INOUT) :: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz)
1080  REAL, INTENT(INOUT) :: mfy(bd%is:bd%ie, bd%js:bd%je+1, npz)
1081 ! Accumulated Courant number arrays
1082  REAL, INTENT(INOUT) :: cx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1083  REAL, INTENT(INOUT) :: cy(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1084  TYPE(fv_grid_type), INTENT(INOUT), TARGET :: gridstruct
1085  TYPE(fv_flags_type), INTENT(INOUT) :: flagstruct
1086  TYPE(fv_flags_pert_type), INTENT(INOUT) :: flagstructp
1087  TYPE(fv_nest_type), INTENT(INOUT) :: neststruct
1088  TYPE(domain2d), INTENT(INOUT) :: domain
1089  TYPE(fv_atmos_type), INTENT(INOUT) :: parent_grid
1090  TYPE(fv_diag_type), INTENT(IN) :: idiag
1091 ! Local Arrays
1092  REAL :: ws(bd%is:bd%ie, bd%js:bd%je)
1093  REAL :: te_2d(bd%is:bd%ie, bd%js:bd%je)
1094  REAL :: teq(bd%is:bd%ie, bd%js:bd%je)
1095  REAL :: ps2(bd%isd:bd%ied, bd%jsd:bd%jed)
1096  REAL :: m_fac(bd%is:bd%ie, bd%js:bd%je)
1097  REAL :: pfull(npz)
1098  REAL, DIMENSION(bd%is:bd%ie) :: cvm
1099  REAL :: dp1(bd%isd:bd%ied, bd%jsd:bd%jed, npz), dtdt_m(bd%is:bd%ie, &
1100 & bd%js:bd%je, npz), cappa(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1101  REAL(kind=8) :: psx(bd%isd:bd%ied, bd%jsd:bd%jed)
1102  REAL(kind=8) :: dpx(bd%is:bd%ie, bd%js:bd%je)
1103  REAL :: akap, rdg, ph1, ph2, mdt, gam, amdt, u0
1104  INTEGER :: kord_tracer(ncnst), kord_mt, kord_wz, kord_tm
1105  INTEGER :: kord_tracer_pert(ncnst), kord_mt_pert, kord_wz_pert, &
1106 & kord_tm_pert
1107  INTEGER :: i, j, k, n, iq, n_map, nq, nwat, k_split
1108 ! GFDL physics
1109  INTEGER :: sphum
1110  INTEGER, SAVE :: liq_wat=-999
1111  INTEGER, SAVE :: ice_wat=-999
1112  INTEGER, SAVE :: rainwat=-999
1113  INTEGER, SAVE :: snowwat=-999
1114  INTEGER, SAVE :: graupel=-999
1115  INTEGER, SAVE :: cld_amt=-999
1116  INTEGER, SAVE :: theta_d=-999
1117  LOGICAL :: used, last_step, do_omega
1118  INTEGER, PARAMETER :: max_packs=12
1119  TYPE(group_halo_update_type), SAVE :: i_pack(max_packs)
1120  INTEGER :: is, ie, js, je
1121  INTEGER :: isd, ied, jsd, jed
1122  REAL :: dt2
1123  REAL(kind=8) :: t1, t2
1124  INTEGER :: status
1125  REAL :: rf(npz)
1126  REAL :: gz(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
1127  REAL :: pkc(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
1128  REAL :: ptc(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1129  REAL :: crx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1130  REAL :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1131  REAL :: cry(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1132  REAL :: yfx(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1133  REAL :: divgd(bd%isd:bd%ied+1, bd%jsd:bd%jed+1, npz)
1134  REAL :: delpc(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1135  REAL :: ut(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1136  REAL :: vt(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1137  REAL :: zh(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
1138  REAL :: pk3(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1)
1139  REAL :: du(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
1140  REAL :: dv(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
1141  INTRINSIC any
1142  INTRINSIC log
1143  INTRINSIC exp
1144  INTRINSIC abs
1145  INTRINSIC real
1146  INTRINSIC cos
1147  REAL :: abs0
1148  REAL :: abs1
1149  REAL :: arg1
1150  REAL :: arg2
1151  REAL :: result1
1152  gz = 0.0
1153  pkc = 0.0
1154  ptc = 0.0
1155  crx = 0.0
1156  xfx = 0.0
1157  cry = 0.0
1158  yfx = 0.0
1159  divgd = 0.0
1160  delpc = 0.0
1161  ut = 0.0
1162  vt = 0.0
1163  zh = 0.0
1164  pk3 = 0.0
1165  du = 0.0
1166  dv = 0.0
1167  is = bd%is
1168  ie = bd%ie
1169  js = bd%js
1170  je = bd%je
1171  isd = bd%isd
1172  ied = bd%ied
1173  jsd = bd%jsd
1174  jed = bd%jed
1175 ! dyn_timer = 0
1176 ! comm_timer = 0
1177 ! cv_air = cp_air - rdgas
1178  agrav = 1./grav
1179  dt2 = 0.5*bdt
1180  k_split = flagstruct%k_split
1181  nwat = flagstruct%nwat
1182  nq = nq_tot - flagstruct%dnats
1183  rdg = -(rdgas*agrav)
1184 !allocate ( dp1(isd:ied, jsd:jed, 1:npz) )
1185 ! Begin Dynamics timer for GEOS history processing
1186 !t1 = MPI_Wtime(status)
1187  t1 = 0.0
1188  t2 = 0.0
1189 !allocate ( cappa(isd:isd,jsd:jsd,1) )
1190  cappa = 0.
1191 !We call this BEFORE converting pt to virtual potential temperature,
1192 !since we interpolate on (regular) temperature rather than theta.
1193  IF (gridstruct%nested .OR. any(neststruct%child_grids)) THEN
1194  CALL timing_on('NEST_BCs')
1195  CALL setup_nested_grid_bcs(npx, npy, npz, zvir, ncnst, u, v, w, pt&
1196 & , delp, delz, q, uc, vc, pkz, neststruct%&
1197 & nested, flagstruct%inline_q, flagstruct%&
1198 & make_nh, ng, gridstruct, flagstruct, &
1199 & neststruct, neststruct%nest_timestep, &
1200 & neststruct%tracer_nest_timestep, domain, bd, &
1201 & nwat)
1202  IF (gridstruct%nested) CALL nested_grid_bc_apply_intt(pt, 0, 0, &
1203 & npx, npy, npz, bd&
1204 & , 1., 1., &
1205 & neststruct%pt_bc, &
1206 & neststruct%&
1207 & nestbctype)
1208 !Correct halo values have now been set up for BCs; we can go ahead and apply them too...
1209  CALL timing_off('NEST_BCs')
1210  END IF
1211  IF (flagstruct%no_dycore) THEN
1212  IF (nwat .EQ. 2 .AND. (.NOT.hydrostatic)) sphum = get_tracer_index&
1213 & (model_atmos, 'sphum')
1214  END IF
1215 !goto 911
1216  IF (fpp%fpp_mapl_mode) THEN
1217  SELECT CASE (nwat)
1218  CASE (0)
1219  sphum = 1
1220 ! to cause trouble if (mis)used
1221  cld_amt = -1
1222  CASE (1)
1223  sphum = 1
1224 ! to cause trouble if (mis)used
1225  liq_wat = -1
1226 ! to cause trouble if (mis)used
1227  ice_wat = -1
1228 ! to cause trouble if (mis)used
1229  rainwat = -1
1230 ! to cause trouble if (mis)used
1231  snowwat = -1
1232 ! to cause trouble if (mis)used
1233  graupel = -1
1234 ! to cause trouble if (mis)used
1235  cld_amt = -1
1236 ! to cause trouble if (mis)used
1237  theta_d = -1
1238  CASE (3)
1239  sphum = 1
1240  liq_wat = 2
1241  ice_wat = 3
1242 ! to cause trouble if (mis)used
1243  rainwat = -1
1244 ! to cause trouble if (mis)used
1245  snowwat = -1
1246 ! to cause trouble if (mis)used
1247  graupel = -1
1248 ! to cause trouble if (mis)used
1249  cld_amt = -1
1250 ! to cause trouble if (mis)used
1251  theta_d = -1
1252  END SELECT
1253  ELSE
1254  IF (nwat .EQ. 0) THEN
1255  sphum = 1
1256 ! to cause trouble if (mis)used
1257  cld_amt = -1
1258  ELSE
1259  sphum = get_tracer_index(model_atmos, 'sphum')
1260  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
1261  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
1262  rainwat = get_tracer_index(model_atmos, 'rainwat')
1263  snowwat = get_tracer_index(model_atmos, 'snowwat')
1264  graupel = get_tracer_index(model_atmos, 'graupel')
1265  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
1266  END IF
1267  theta_d = get_tracer_index(model_atmos, 'theta_d')
1268  END IF
1269  akap = kappa
1270 !$OMP parallel do default(none) shared(npz,ak,bk,flagstruct,pfull) &
1271 !$OMP private(ph1, ph2)
1272  DO k=1,npz
1273  ph1 = ak(k) + bk(k)*flagstruct%p_ref
1274  ph2 = ak(k+1) + bk(k+1)*flagstruct%p_ref
1275  pfull(k) = (ph2-ph1)/log(ph2/ph1)
1276  END DO
1277  IF (hydrostatic) THEN
1278 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,nwat,q,q_con,sphum,liq_wat, &
1279 !$OMP rainwat,ice_wat,snowwat,graupel) private(cvm)
1280  DO k=1,npz
1281  DO j=js,je
1282  DO i=is,ie
1283  dp1(i, j, k) = zvir*q(i, j, k, sphum)
1284  END DO
1285  END DO
1286  END DO
1287  ELSE
1288 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,q,q_con,sphum,liq_wat, &
1289 !$OMP rainwat,ice_wat,snowwat,graupel,pkz,flagstruct, &
1290 !$OMP cappa,kappa,rdg,delp,pt,delz,nwat) &
1291 !$OMP private(cvm)
1292  DO k=1,npz
1293  IF (flagstruct%moist_phys) THEN
1294  DO j=js,je
1295  DO i=is,ie
1296  dp1(i, j, k) = zvir*q(i, j, k, sphum)
1297  arg1 = rdg*delp(i, j, k)*pt(i, j, k)*(1.+dp1(i, j, k))/&
1298 & delz(i, j, k)
1299  arg2 = kappa*log(arg1)
1300  pkz(i, j, k) = exp(arg2)
1301  END DO
1302  END DO
1303  ELSE
1304 ! Using dry pressure for the definition of the virtual potential temperature
1305 ! pkz(i,j,k) = exp( kappa*log(rdg*delp(i,j,k)*pt(i,j,k)* &
1306 ! (1.-q(i,j,k,sphum))/delz(i,j,k)) )
1307  DO j=js,je
1308  DO i=is,ie
1309  dp1(i, j, k) = 0.
1310  arg1 = rdg*delp(i, j, k)*pt(i, j, k)/delz(i, j, k)
1311  arg2 = kappa*log(arg1)
1312  pkz(i, j, k) = exp(arg2)
1313  END DO
1314  END DO
1315  END IF
1316  END DO
1317  END IF
1318  IF (flagstruct%fv_debug) THEN
1319  CALL prt_mxm('PS', ps, is, ie, js, je, ng, 1, 0.01, gridstruct%&
1320 & area_64, domain)
1321  CALL prt_mxm('T_dyn_b', pt, is, ie, js, je, ng, npz, 1., &
1322 & gridstruct%area_64, domain)
1323  IF (.NOT.hydrostatic) CALL prt_mxm('delz', delz, is, ie, js, je, &
1324 & ng, npz, 1., gridstruct%area_64, &
1325 & domain)
1326  CALL prt_mxm('delp_b ', delp, is, ie, js, je, ng, npz, 0.01, &
1327 & gridstruct%area_64, domain)
1328  CALL prt_mxm('pk_b', pk, is, ie, js, je, 0, npz + 1, 1., &
1329 & gridstruct%area_64, domain)
1330  CALL prt_mxm('pkz_b', pkz, is, ie, js, je, 0, npz, 1., gridstruct%&
1331 & area_64, domain)
1332  END IF
1333 !---------------------
1334 ! Compute Total Energy
1335 !---------------------
1336  IF (consv_te .GT. 0. .AND. (.NOT.do_adiabatic_init)) THEN
1337  CALL compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, npz&
1338 & , u, v, w, delz, pt, delp, q, dp1, pe, peln, &
1339 & phis, gridstruct%rsin2, gridstruct%cosa_s, &
1340 & zvir, cp_air, rdgas, hlv, te_2d, ua, va, teq, &
1341 & flagstruct%moist_phys, nwat, sphum, liq_wat, &
1342 & rainwat, ice_wat, snowwat, graupel, &
1343 & hydrostatic, idiag%id_te)
1344  IF (idiag%id_te .GT. 0) used = send_data(idiag%id_te, teq, fv_time&
1345 & )
1346 ! te_den=1.E-9*g_sum(teq, is, ie, js, je, ng, area, 0)/(grav*4.*pi*radius**2)
1347 ! if(is_master()) write(*,*) 'Total Energy Density (Giga J/m**2)=',te_den
1348  END IF
1349  IF ((flagstruct%consv_am .OR. idiag%id_amdt .GT. 0) .AND. (.NOT.&
1350 & do_adiabatic_init)) CALL compute_aam(npz, is, ie, js, je, isd, &
1351 & ied, jsd, jed, gridstruct, bd, &
1352 & ptop, ua, va, u, v, delp, teq, &
1353 & ps2, m_fac)
1354  IF (flagstruct%tau .GT. 0.) THEN
1355  IF (gridstruct%grid_type .LT. 4) THEN
1356  IF (bdt .GE. 0.) THEN
1357  abs0 = bdt
1358  ELSE
1359  abs0 = -bdt
1360  END IF
1361  CALL rayleigh_super(abs0, npx, npy, npz, ks, pfull, phis, &
1362 & flagstruct%tau, u, v, w, pt, ua, va, delz, &
1363 & gridstruct%agrid, cp_air, rdgas, ptop, hydrostatic&
1364 & , .NOT.neststruct%nested, flagstruct%rf_cutoff, rf&
1365 & , gridstruct, domain, bd)
1366  ELSE
1367  IF (bdt .GE. 0.) THEN
1368  abs1 = bdt
1369  ELSE
1370  abs1 = -bdt
1371  END IF
1372  CALL rayleigh_friction(abs1, npx, npy, npz, ks, pfull, &
1373 & flagstruct%tau, u, v, w, pt, ua, va, delz, &
1374 & cp_air, rdgas, ptop, hydrostatic, .true., &
1375 & flagstruct%rf_cutoff, rf, gridstruct, domain, &
1376 & bd)
1377  END IF
1378  END IF
1379 ! Convert pt to virtual potential temperature on the first timestep
1380  IF (flagstruct%adiabatic) THEN
1381 !$OMP parallel do default(none) shared(theta_d,is,ie,js,je,npz,pt,pkz,q)
1382  DO k=1,npz
1383  DO j=js,je
1384  DO i=is,ie
1385  pt(i, j, k) = pt(i, j, k)/pkz(i, j, k)
1386  END DO
1387  END DO
1388  IF (theta_d .GT. 0) THEN
1389  DO j=js,je
1390  DO i=is,ie
1391  q(i, j, k, theta_d) = pt(i, j, k)
1392  END DO
1393  END DO
1394  END IF
1395  END DO
1396  ELSE
1397 !$OMP parallel do default(none) shared(is,ie,js,je,npz,pt,dp1,pkz,q_con)
1398  DO k=1,npz
1399  DO j=js,je
1400  DO i=is,ie
1401  pt(i, j, k) = pt(i, j, k)*(1.+dp1(i, j, k))/pkz(i, j, k)
1402  END DO
1403  END DO
1404  END DO
1405  END IF
1406  last_step = .false.
1407  mdt = bdt/REAL(k_split)
1408  IF (idiag%id_mdt .GT. 0 .AND. (.NOT.do_adiabatic_init)) THEN
1409 !allocate ( dtdt_m(is:ie,js:je,npz) )
1410 !$OMP parallel do default(none) shared(is,ie,js,je,npz,dtdt_m)
1411  DO k=1,npz
1412  DO j=js,je
1413  DO i=is,ie
1414  dtdt_m(i, j, k) = 0.
1415  END DO
1416  END DO
1417  END DO
1418  END IF
1419 !DryMassRoundoffControl
1420 !allocate(psx(isd:ied,jsd:jed),dpx(is:ie,js:je))
1421  IF (fpp%fpp_overload_r4) THEN
1422  DO j=js,je
1423  DO i=is,ie
1424  psx(i, j) = pe(i, npz+1, j)
1425  dpx(i, j) = 0.0
1426  END DO
1427  END DO
1428  END IF
1429  CALL timing_on('FV_DYN_LOOP')
1430 ! first level of time-split
1431  DO n_map=1,k_split
1432  CALL timing_on('COMM_TOTAL')
1433  CALL start_group_halo_update(i_pack(1), delp, domain, complete=&
1434 & .true.)
1435  CALL start_group_halo_update(i_pack(1), pt, domain, complete=&
1436 & .true.)
1437  CALL start_group_halo_update(i_pack(8), u, v, domain, gridtype=&
1438 & dgrid_ne)
1439  CALL timing_off('COMM_TOTAL')
1440 !$OMP parallel do default(none) shared(isd,ied,jsd,jed,npz,dp1,delp)
1441  DO k=1,npz
1442  DO j=jsd,jed
1443  DO i=isd,ied
1444  dp1(i, j, k) = delp(i, j, k)
1445  END DO
1446  END DO
1447  END DO
1448  IF (n_map .EQ. k_split) last_step = .true.
1449  CALL timing_on('DYN_CORE')
1450  CALL dyn_core(npx, npy, npz, ng, sphum, nq, mdt, n_split, zvir, &
1451 & cp_air, akap, cappa, grav, hydrostatic, u, v, w, delz, pt&
1452 & , q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, uc&
1453 & , vc, mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, dpx, ks&
1454 & , gridstruct, flagstruct, flagstructp, neststruct, idiag, &
1455 & bd, domain, n_map .EQ. 1, i_pack, last_step, gz, pkc, ptc&
1456 & , crx, xfx, cry, yfx, divgd, delpc, ut, vt, zh, pk3, du, &
1457 & dv, time_total)
1458  CALL timing_off('DYN_CORE')
1459 !DryMassRoundoffControl
1460  IF (last_step) THEN
1461  IF (fpp%fpp_overload_r4) THEN
1462  DO j=js,je
1463  DO i=is,ie
1464  psx(i, j) = psx(i, j) + dpx(i, j)
1465  END DO
1466  END DO
1467  CALL timing_on('COMM_TOTAL')
1468  CALL mpp_update_domains(psx, domain)
1469  CALL timing_off('COMM_TOTAL')
1470  DO j=js-1,je+1
1471  DO i=is-1,ie+1
1472  pe(i, npz+1, j) = psx(i, j)
1473  END DO
1474  END DO
1475  END IF
1476  END IF
1477 !deallocate(psx,dpx)
1478  IF (.NOT.flagstruct%inline_q .AND. nq .NE. 0) THEN
1479 !--------------------------------------------------------
1480 ! Perform large-time-step scalar transport using the accumulated CFL and
1481 ! mass fluxes
1482  CALL timing_on('tracer_2d')
1483 !!! CLEANUP: merge these two calls?
1484  IF (gridstruct%nested) THEN
1485  CALL tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd&
1486 & , domain, npx, npy, npz, nq, flagstruct%&
1487 & hord_tr, q_split, mdt, idiag%id_divg, i_pack(&
1488 & 10), flagstruct%nord_tr, flagstruct%trdm2, &
1489 & k_split, neststruct, parent_grid, flagstructp%&
1490 & hord_tr_pert, flagstructp%nord_tr_pert, &
1491 & flagstructp%trdm2_pert, flagstructp%&
1492 & split_damp_tr)
1493  ELSE IF (flagstruct%z_tracer) THEN
1494  CALL tracer_2d_1l(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, &
1495 & domain, npx, npy, npz, nq, flagstruct%hord_tr, &
1496 & q_split, mdt, idiag%id_divg, i_pack(10), &
1497 & flagstruct%nord_tr, flagstruct%trdm2, flagstructp%&
1498 & hord_tr_pert, flagstructp%nord_tr_pert, &
1499 & flagstructp%trdm2_pert, flagstructp%split_damp_tr)
1500  ELSE
1501  CALL tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, &
1502 & domain, npx, npy, npz, nq, flagstruct%hord_tr, &
1503 & q_split, mdt, idiag%id_divg, i_pack(10), flagstruct%&
1504 & nord_tr, flagstruct%trdm2, flagstructp%hord_tr_pert, &
1505 & flagstructp%nord_tr_pert, flagstructp%trdm2_pert, &
1506 & flagstructp%split_damp_tr)
1507  END IF
1508  CALL timing_off('tracer_2d')
1509  IF (flagstruct%hord_tr .LT. 8 .AND. flagstruct%moist_phys) THEN
1510  CALL timing_on('Fill2D')
1511  IF (liq_wat .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
1512 & :ied, jsd:jed, 1, liq_wat), delp, &
1513 & gridstruct%area, domain, neststruct%&
1514 & nested, npx, npy)
1515  IF (rainwat .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
1516 & :ied, jsd:jed, 1, rainwat), delp, &
1517 & gridstruct%area, domain, neststruct%&
1518 & nested, npx, npy)
1519  IF (ice_wat .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
1520 & :ied, jsd:jed, 1, ice_wat), delp, &
1521 & gridstruct%area, domain, neststruct%&
1522 & nested, npx, npy)
1523  IF (snowwat .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
1524 & :ied, jsd:jed, 1, snowwat), delp, &
1525 & gridstruct%area, domain, neststruct%&
1526 & nested, npx, npy)
1527  IF (graupel .GT. 0) CALL fill2d(is, ie, js, je, ng, npz, q(isd&
1528 & :ied, jsd:jed, 1, graupel), delp, &
1529 & gridstruct%area, domain, neststruct%&
1530 & nested, npx, npy)
1531  CALL timing_off('Fill2D')
1532  END IF
1533  IF (last_step .AND. idiag%id_divg .GT. 0) THEN
1534  used = send_data(idiag%id_divg, dp1, fv_time)
1535  IF (flagstruct%fv_debug) CALL prt_mxm('divg', dp1, is, ie, js&
1536 & , je, 0, npz, 1., gridstruct%&
1537 & area_64, domain)
1538  END IF
1539  END IF
1540  IF (npz .GT. 4) THEN
1541 !------------------------------------------------------------------------
1542 ! Perform vertical remapping from Lagrangian control-volume to
1543 ! the Eulerian coordinate as specified by the routine set_eta.
1544 ! Note that this finite-volume dycore is otherwise independent of the vertical
1545 ! Eulerian coordinate.
1546 !------------------------------------------------------------------------
1547  DO iq=1,nq
1548  kord_tracer(iq) = flagstruct%kord_tr
1549 ! monotonic
1550  IF (iq .EQ. cld_amt) kord_tracer(iq) = 9
1551  kord_tracer_pert(iq) = flagstructp%kord_tr_pert
1552 ! linear
1553  IF (iq .EQ. cld_amt) kord_tracer_pert(iq) = 17
1554  END DO
1555  do_omega = hydrostatic .AND. last_step
1556  CALL timing_on('Remapping')
1557  kord_mt = flagstruct%kord_mt
1558  kord_wz = flagstruct%kord_wz
1559  kord_tm = flagstruct%kord_tm
1560  kord_mt_pert = flagstructp%kord_mt_pert
1561  kord_wz_pert = flagstructp%kord_wz_pert
1562  kord_tm_pert = flagstructp%kord_tm_pert
1563  IF (n_map .EQ. k_split) THEN
1564  kord_mt = kord_mt_pert
1565  kord_wz = kord_wz_pert
1566  kord_tm = kord_tm_pert
1567  kord_tracer = kord_tracer_pert
1568  END IF
1569  CALL lagrangian_to_eulerian(last_step, consv_te, ps, pe, delp, &
1570 & pkz, pk, mdt, bdt, npz, is, ie, js, je, &
1571 & isd, ied, jsd, jed, nq, nwat, sphum, q_con&
1572 & , u, v, w, delz, pt, q, phis, zvir, cp_air&
1573 & , akap, cappa, kord_mt, kord_wz, &
1574 & kord_tracer, kord_tm, peln, te_2d, ng, ua&
1575 & , va, omga, dp1, ws, fill, reproduce_sum, &
1576 & idiag%id_mdt .GT. 0, dtdt_m, ptop, ak, bk&
1577 & , pfull, flagstruct, gridstruct, domain, &
1578 & flagstruct%do_sat_adj, hydrostatic, &
1579 & hybrid_z, do_omega, flagstruct%adiabatic, &
1580 & do_adiabatic_init, mfx, mfy, flagstruct%&
1581 & remap_option, kord_mt_pert, kord_wz_pert, &
1582 & kord_tracer_pert, kord_tm_pert)
1583  CALL timing_off('Remapping')
1584  IF (last_step) THEN
1585  IF (.NOT.hydrostatic) THEN
1586 !$OMP parallel do default(none) shared(is,ie,js,je,npz,omga,delp,delz,w)
1587  DO k=1,npz
1588  DO j=js,je
1589  DO i=is,ie
1590  omga(i, j, k) = delp(i, j, k)/delz(i, j, k)*w(i, j, k)
1591  END DO
1592  END DO
1593  END DO
1594  END IF
1595 !--------------------------
1596 ! Filter omega for physics:
1597 !--------------------------
1598  IF (flagstruct%nf_omega .GT. 0) CALL del2_cubed(omga, 0.18*&
1599 & gridstruct%da_min, &
1600 & gridstruct, domain, &
1601 & npx, npy, npz, &
1602 & flagstruct%nf_omega&
1603 & , bd)
1604  END IF
1605  END IF
1606  END DO
1607 ! n_map loop
1608  CALL timing_off('FV_DYN_LOOP')
1609  IF (idiag%id_mdt .GT. 0 .AND. (.NOT.do_adiabatic_init)) THEN
1610 ! Output temperature tendency due to inline moist physics:
1611 !$OMP parallel do default(none) shared(is,ie,js,je,npz,dtdt_m,bdt)
1612  DO k=1,npz
1613  DO j=js,je
1614  DO i=is,ie
1615  dtdt_m(i, j, k) = dtdt_m(i, j, k)/bdt*86400.
1616  END DO
1617  END DO
1618  END DO
1619 ! call prt_mxm('Fast DTDT (deg/Day)', dtdt_m, is, ie, js, je, 0, npz, 1., gridstruct%area_64, domain)
1620  used = send_data(idiag%id_mdt, dtdt_m, fv_time)
1621 !deallocate ( dtdt_m )
1622  END IF
1623  IF (nwat .EQ. 6) THEN
1624  IF (cld_amt .GT. 0) THEN
1625  CALL neg_adj3(is, ie, js, je, ng, npz, flagstruct%hydrostatic, &
1626 & peln, delz, pt, delp, q(isd:ied, jsd:jed, 1, sphum), q(&
1627 & isd:ied, jsd:jed, 1, liq_wat), q(isd:ied, jsd:jed, 1, &
1628 & rainwat), q(isd:ied, jsd:jed, 1, ice_wat), q(isd:ied, &
1629 & jsd:jed, 1, snowwat), q(isd:ied, jsd:jed, 1, graupel), q&
1630 & (isd:ied, jsd:jed, 1, cld_amt), flagstruct%&
1631 & check_negative)
1632  ELSE
1633  CALL neg_adj3(is, ie, js, je, ng, npz, flagstruct%hydrostatic, &
1634 & peln, delz, pt, delp, q(isd:ied, jsd:jed, 1, sphum), q(&
1635 & isd:ied, jsd:jed, 1, liq_wat), q(isd:ied, jsd:jed, 1, &
1636 & rainwat), q(isd:ied, jsd:jed, 1, ice_wat), q(isd:ied, &
1637 & jsd:jed, 1, snowwat), q(isd:ied, jsd:jed, 1, graupel), &
1638 & check_negative=flagstruct%check_negative)
1639  END IF
1640  IF (flagstruct%fv_debug) THEN
1641  CALL prt_mxm('T_dyn_a3', pt, is, ie, js, je, ng, npz, 1., &
1642 & gridstruct%area_64, domain)
1643  CALL prt_mxm('SPHUM_dyn', q(isd:ied, jsd:jed, 1, sphum), is, ie&
1644 & , js, je, ng, npz, 1., gridstruct%area_64, domain)
1645  CALL prt_mxm('liq_wat_dyn', q(isd:ied, jsd:jed, 1, liq_wat), is&
1646 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
1647  CALL prt_mxm('rainwat_dyn', q(isd:ied, jsd:jed, 1, rainwat), is&
1648 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
1649  CALL prt_mxm('ice_wat_dyn', q(isd:ied, jsd:jed, 1, ice_wat), is&
1650 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
1651  CALL prt_mxm('snowwat_dyn', q(isd:ied, jsd:jed, 1, snowwat), is&
1652 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
1653  CALL prt_mxm('graupel_dyn', q(isd:ied, jsd:jed, 1, graupel), is&
1654 & , ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
1655  END IF
1656  END IF
1657  IF (((flagstruct%consv_am .OR. idiag%id_amdt .GT. 0) .OR. idiag%&
1658 & id_aam .GT. 0) .AND. (.NOT.do_adiabatic_init)) THEN
1659  CALL compute_aam(npz, is, ie, js, je, isd, ied, jsd, jed, &
1660 & gridstruct, bd, ptop, ua, va, u, v, delp, te_2d, ps, &
1661 & m_fac)
1662  IF (idiag%id_aam .GT. 0) THEN
1663  used = send_data(idiag%id_aam, te_2d, fv_time)
1664  IF (prt_minmax) gam = g_sum(domain, te_2d, is, ie, js, je, ng, &
1665 & gridstruct%area_64, 0)
1666 !if( is_master() ) write(6,*) 'Total AAM =', gam
1667  END IF
1668  END IF
1669  IF ((flagstruct%consv_am .OR. idiag%id_amdt .GT. 0) .AND. (.NOT.&
1670 & do_adiabatic_init)) THEN
1671 !$OMP parallel do default(none) shared(is,ie,js,je,te_2d,teq,dt2,ps2,ps,idiag)
1672  DO j=js,je
1673  DO i=is,ie
1674 ! Note: the mountain torque computation contains also numerical error
1675 ! The numerical error is mostly from the zonal gradient of the terrain (zxg)
1676  te_2d(i, j) = te_2d(i, j) - teq(i, j) + dt2*(ps2(i, j)+ps(i, j&
1677 & ))*idiag%zxg(i, j)
1678  END DO
1679  END DO
1680  IF (idiag%id_amdt .GT. 0) used = send_data(idiag%id_amdt, te_2d/&
1681 & bdt, fv_time)
1682  IF (flagstruct%consv_am .OR. prt_minmax) THEN
1683  amdt = g_sum(domain, te_2d, is, ie, js, je, ng, gridstruct%&
1684 & area_64, 0, .true.)
1685  result1 = g_sum(domain, m_fac, is, ie, js, je, ng, gridstruct%&
1686 & area_64, 0, .true.)
1687  u0 = -(radius*amdt/result1)
1688  IF (is_master() .AND. prt_minmax) WRITE(6, *) &
1689 & 'Dynamic AM tendency (Hadleys)='&
1690 & , amdt/(bdt*1.e18), &
1691 & 'del-u (per day)=', u0*86400./&
1692 & bdt
1693  END IF
1694 ! consv_am
1695  IF (flagstruct%consv_am) THEN
1696 !$OMP parallel do default(none) shared(is,ie,js,je,m_fac,u0,gridstruct)
1697  DO j=js,je
1698  DO i=is,ie
1699  m_fac(i, j) = u0*cos(gridstruct%agrid(i, j, 2))
1700  END DO
1701  END DO
1702 !$OMP parallel do default(none) shared(is,ie,js,je,npz,hydrostatic,pt,m_fac,ua,cp_air, &
1703 !$OMP u,u0,gridstruct,v )
1704  DO k=1,npz
1705  DO j=js,je+1
1706  DO i=is,ie
1707  u(i, j, k) = u(i, j, k) + u0*gridstruct%l2c_u(i, j)
1708  END DO
1709  END DO
1710  DO j=js,je
1711  DO i=is,ie+1
1712  v(i, j, k) = v(i, j, k) + u0*gridstruct%l2c_v(i, j)
1713  END DO
1714  END DO
1715  END DO
1716  END IF
1717  END IF
1718 !911 call cubed_to_latlon(u, v, ua, va, gridstruct, &
1719 ! npx, npy, npz, 1, gridstruct%grid_type, domain, gridstruct%nested, flagstruct%c2l_ord, bd)
1720 !deallocate(dp1)
1721 !deallocate(cappa)
1722  IF (flagstruct%fv_debug) THEN
1723  CALL prt_mxm('UA', ua, is, ie, js, je, ng, npz, 1., gridstruct%&
1724 & area_64, domain)
1725  CALL prt_mxm('VA', va, is, ie, js, je, ng, npz, 1., gridstruct%&
1726 & area_64, domain)
1727  CALL prt_mxm('TA', pt, is, ie, js, je, ng, npz, 1., gridstruct%&
1728 & area_64, domain)
1729  IF (.NOT.hydrostatic) CALL prt_mxm('W ', w, is, ie, js, je, ng, &
1730 & npz, 1., gridstruct%area_64, domain)
1731  END IF
1732  IF (flagstruct%range_warn) THEN
1733  CALL range_check('UA_dyn', ua, is, ie, js, je, ng, npz, gridstruct&
1734 & %agrid, -280., 280., bad_range)
1735  CALL range_check('VA_dyn', ua, is, ie, js, je, ng, npz, gridstruct&
1736 & %agrid, -280., 280., bad_range)
1737  CALL range_check('TA_dyn', pt, is, ie, js, je, ng, npz, gridstruct&
1738 & %agrid, 150., 335., bad_range)
1739  IF (.NOT.hydrostatic) CALL range_check('W_dyn', w, is, ie, js, je&
1740 & , ng, npz, gridstruct%agrid, -50.&
1741 & , 100., bad_range)
1742  END IF
1743 ! IF (fpp%fpp_mapl_mode) dyn_timer = dyn_timer + (t2-t1)
1744 !t2 = MPI_Wtime(status)
1745  END SUBROUTINE fv_dynamics
1746 ! Differentiation of rayleigh_super in forward (tangent) mode:
1747 ! variations of useful results: u v w ua va pt
1748 ! with respect to varying inputs: u v w ua va pt
1749  SUBROUTINE rayleigh_super_tlm(dt, npx, npy, npz, ks, pm, phis, tau, u&
1750 & , u_tl, v, v_tl, w, w_tl, pt, pt_tl, ua, ua_tl, va, va_tl, delz, &
1751 & agrid, cp, rg, ptop, hydrostatic, conserve, rf_cutoff, rf, &
1752 & gridstruct, domain, bd)
1753  IMPLICIT NONE
1754 !deallocate ( u2f )
1755  REAL, INTENT(IN) :: dt
1756 ! time scale (days)
1757  REAL, INTENT(IN) :: tau
1758  REAL, INTENT(IN) :: cp, rg, ptop, rf_cutoff
1759  INTEGER, INTENT(IN) :: npx, npy, npz, ks
1760  REAL, DIMENSION(npz), INTENT(IN) :: pm
1761  LOGICAL, INTENT(IN) :: hydrostatic
1762  LOGICAL, INTENT(IN) :: conserve
1763  TYPE(FV_GRID_BOUNDS_TYPE), INTENT(IN) :: bd
1764 ! D grid zonal wind (m/s)
1765  REAL, INTENT(INOUT) :: u(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
1766  REAL, INTENT(INOUT) :: u_tl(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
1767 ! D grid meridional wind (m/s)
1768  REAL, INTENT(INOUT) :: v(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
1769  REAL, INTENT(INOUT) :: v_tl(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
1770 ! cell center vertical wind (m/s)
1771  REAL, INTENT(INOUT) :: w(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1772  REAL, INTENT(INOUT) :: w_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1773 ! temp
1774  REAL, INTENT(INOUT) :: pt(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1775  REAL, INTENT(INOUT) :: pt_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1776 !
1777  REAL, INTENT(INOUT) :: ua(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1778  REAL, INTENT(INOUT) :: ua_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1779 !
1780  REAL, INTENT(INOUT) :: va(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1781  REAL, INTENT(INOUT) :: va_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1782 ! delta-height (m); non-hydrostatic only
1783  REAL, INTENT(INOUT) :: delz(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1784  REAL, INTENT(INOUT) :: rf(npz)
1785  REAL, INTENT(IN) :: agrid(bd%isd:bd%ied, bd%jsd:bd%jed, 2)
1786 ! Surface geopotential (g*Z_surf)
1787  REAL, INTENT(IN) :: phis(bd%isd:bd%ied, bd%jsd:bd%jed)
1788  TYPE(FV_GRID_TYPE), INTENT(IN) :: gridstruct
1789  TYPE(DOMAIN2D), INTENT(INOUT) :: domain
1790 !
1791  REAL :: u2f(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1792 ! scaling velocity
1793  REAL, PARAMETER :: u0=60.
1794  REAL, PARAMETER :: sday=86400.
1795  REAL :: rcv, tau0
1796  INTEGER :: i, j, k
1797  INTEGER :: is, ie, js, je
1798  INTEGER :: isd, ied, jsd, jed
1799  INTRINSIC log
1800  INTRINSIC sin
1801  REAL*8 :: arg1
1802  is = bd%is
1803  ie = bd%ie
1804  js = bd%js
1805  je = bd%je
1806  isd = bd%isd
1807  ied = bd%ied
1808  jsd = bd%jsd
1809  jed = bd%jed
1810  rcv = 1./(cp-rg)
1811  IF (.NOT.rf_initialized) THEN
1812  tau0 = tau*sday
1813 !allocate( rf(npz) )
1814  rf(:) = 0.
1815  k = ks + 2
1816 !if( is_master() ) write(6,*) k, 0.01*pm(k)
1817  IF (is_master()) WRITE(6, *) &
1818 & 'Rayleigh friction E-folding time (days):'
1819  DO k=1,npz
1820  IF (pm(k) .LT. rf_cutoff) THEN
1821  arg1 = 0.5*pi*log(rf_cutoff/pm(k))/log(rf_cutoff/ptop)
1822  rf(k) = dt/tau0*sin(arg1)**2
1823 !if( is_master() ) write(6,*) k, 0.01*pm(k), dt/(rf(k)*sday)
1824  kmax = k
1825  ELSE
1826  EXIT
1827  END IF
1828  END DO
1829  rf_initialized = .true.
1830  END IF
1831  CALL c2l_ord2_tlm(u, u_tl, v, v_tl, ua, ua_tl, va, va_tl, gridstruct&
1832 & , npz, gridstruct%grid_type, bd, gridstruct%nested)
1833 !allocate( u2f(isd:ied,jsd:jed,kmax) )
1834  u2f = 0.0
1835 !$OMP parallel do default(none) shared(is,ie,js,je,kmax,pm,rf_cutoff,hydrostatic,ua,va,agrid, &
1836 !$OMP u2f,rf,w)
1837  DO k=1,kmax
1838  IF (pm(k) .LT. rf_cutoff) THEN
1839  u2f(:, :, k) = 1./(1.+rf(k))
1840  ELSE
1841  u2f(:, :, k) = 1.
1842  END IF
1843  END DO
1844  CALL timing_on('COMM_TOTAL')
1845  CALL mpp_update_domains(u2f, domain)
1846  CALL timing_off('COMM_TOTAL')
1847 !$OMP parallel do default(none) shared(is,ie,js,je,kmax,pm,rf_cutoff,w,rf,u,v, &
1848 !$OMP conserve,hydrostatic,pt,ua,va,u2f,cp,rg,ptop,rcv)
1849  DO k=1,kmax
1850  IF (pm(k) .LT. rf_cutoff) THEN
1851 ! Add heat so as to conserve TE
1852  IF (conserve) THEN
1853  IF (hydrostatic) THEN
1854  DO j=js,je
1855  DO i=is,ie
1856  pt_tl(i, j, k) = pt_tl(i, j, k) + 0.5*(1.-u2f(i, j, k)**&
1857 & 2)*(2*ua(i, j, k)*ua_tl(i, j, k)+2*va(i, j, k)*va_tl(i&
1858 & , j, k))/(cp-rg*ptop/pm(k))
1859  pt(i, j, k) = pt(i, j, k) + 0.5*(ua(i, j, k)**2+va(i, j&
1860 & , k)**2)*(1.-u2f(i, j, k)**2)/(cp-rg*ptop/pm(k))
1861  END DO
1862  END DO
1863  ELSE
1864  DO j=js,je
1865  DO i=is,ie
1866  pt_tl(i, j, k) = pt_tl(i, j, k) + 0.5*(1.-u2f(i, j, k)**&
1867 & 2)*rcv*(2*ua(i, j, k)*ua_tl(i, j, k)+2*va(i, j, k)*&
1868 & va_tl(i, j, k)+2*w(i, j, k)*w_tl(i, j, k))
1869  pt(i, j, k) = pt(i, j, k) + 0.5*(ua(i, j, k)**2+va(i, j&
1870 & , k)**2+w(i, j, k)**2)*(1.-u2f(i, j, k)**2)*rcv
1871  END DO
1872  END DO
1873  END IF
1874  END IF
1875  DO j=js,je+1
1876  DO i=is,ie
1877  u_tl(i, j, k) = 0.5*(u2f(i, j-1, k)+u2f(i, j, k))*u_tl(i, j&
1878 & , k)
1879  u(i, j, k) = 0.5*(u2f(i, j-1, k)+u2f(i, j, k))*u(i, j, k)
1880  END DO
1881  END DO
1882  DO j=js,je
1883  DO i=is,ie+1
1884  v_tl(i, j, k) = 0.5*(u2f(i-1, j, k)+u2f(i, j, k))*v_tl(i, j&
1885 & , k)
1886  v(i, j, k) = 0.5*(u2f(i-1, j, k)+u2f(i, j, k))*v(i, j, k)
1887  END DO
1888  END DO
1889  IF (.NOT.hydrostatic) THEN
1890  DO j=js,je
1891  DO i=is,ie
1892  w_tl(i, j, k) = u2f(i, j, k)*w_tl(i, j, k)
1893  w(i, j, k) = u2f(i, j, k)*w(i, j, k)
1894  END DO
1895  END DO
1896  END IF
1897  END IF
1898  END DO
1899  END SUBROUTINE rayleigh_super_tlm
1900  SUBROUTINE rayleigh_super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, &
1901 & w, pt, ua, va, delz, agrid, cp, rg, ptop, hydrostatic, conserve, &
1902 & rf_cutoff, rf, gridstruct, domain, bd)
1903  IMPLICIT NONE
1904 !deallocate ( u2f )
1905  REAL, INTENT(IN) :: dt
1906 ! time scale (days)
1907  REAL, INTENT(IN) :: tau
1908  REAL, INTENT(IN) :: cp, rg, ptop, rf_cutoff
1909  INTEGER, INTENT(IN) :: npx, npy, npz, ks
1910  REAL, DIMENSION(npz), INTENT(IN) :: pm
1911  LOGICAL, INTENT(IN) :: hydrostatic
1912  LOGICAL, INTENT(IN) :: conserve
1913  TYPE(FV_GRID_BOUNDS_TYPE), INTENT(IN) :: bd
1914 ! D grid zonal wind (m/s)
1915  REAL, INTENT(INOUT) :: u(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
1916 ! D grid meridional wind (m/s)
1917  REAL, INTENT(INOUT) :: v(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
1918 ! cell center vertical wind (m/s)
1919  REAL, INTENT(INOUT) :: w(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1920 ! temp
1921  REAL, INTENT(INOUT) :: pt(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1922 !
1923  REAL, INTENT(INOUT) :: ua(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1924 !
1925  REAL, INTENT(INOUT) :: va(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1926 ! delta-height (m); non-hydrostatic only
1927  REAL, INTENT(INOUT) :: delz(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1928  REAL, INTENT(INOUT) :: rf(npz)
1929  REAL, INTENT(IN) :: agrid(bd%isd:bd%ied, bd%jsd:bd%jed, 2)
1930 ! Surface geopotential (g*Z_surf)
1931  REAL, INTENT(IN) :: phis(bd%isd:bd%ied, bd%jsd:bd%jed)
1932  TYPE(FV_GRID_TYPE), INTENT(IN) :: gridstruct
1933  TYPE(DOMAIN2D), INTENT(INOUT) :: domain
1934 !
1935  REAL :: u2f(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1936 ! scaling velocity
1937  REAL, PARAMETER :: u0=60.
1938  REAL, PARAMETER :: sday=86400.
1939  REAL :: rcv, tau0
1940  INTEGER :: i, j, k
1941  INTEGER :: is, ie, js, je
1942  INTEGER :: isd, ied, jsd, jed
1943  INTRINSIC log
1944  INTRINSIC sin
1945  REAL*8 :: arg1
1946  is = bd%is
1947  ie = bd%ie
1948  js = bd%js
1949  je = bd%je
1950  isd = bd%isd
1951  ied = bd%ied
1952  jsd = bd%jsd
1953  jed = bd%jed
1954  rcv = 1./(cp-rg)
1955  IF (.NOT.rf_initialized) THEN
1956  tau0 = tau*sday
1957 !allocate( rf(npz) )
1958  rf(:) = 0.
1959  k = ks + 2
1960 !if( is_master() ) write(6,*) k, 0.01*pm(k)
1961  IF (is_master()) WRITE(6, *) &
1962 & 'Rayleigh friction E-folding time (days):'
1963  DO k=1,npz
1964  IF (pm(k) .LT. rf_cutoff) THEN
1965  arg1 = 0.5*pi*log(rf_cutoff/pm(k))/log(rf_cutoff/ptop)
1966  rf(k) = dt/tau0*sin(arg1)**2
1967 !if( is_master() ) write(6,*) k, 0.01*pm(k), dt/(rf(k)*sday)
1968  kmax = k
1969  ELSE
1970  GOTO 100
1971  END IF
1972  END DO
1973  100 rf_initialized = .true.
1974  END IF
1975  CALL c2l_ord2(u, v, ua, va, gridstruct, npz, gridstruct%grid_type, &
1976 & bd, gridstruct%nested)
1977 !allocate( u2f(isd:ied,jsd:jed,kmax) )
1978  u2f = 0.0
1979 !$OMP parallel do default(none) shared(is,ie,js,je,kmax,pm,rf_cutoff,hydrostatic,ua,va,agrid, &
1980 !$OMP u2f,rf,w)
1981  DO k=1,kmax
1982  IF (pm(k) .LT. rf_cutoff) THEN
1983  u2f(:, :, k) = 1./(1.+rf(k))
1984  ELSE
1985  u2f(:, :, k) = 1.
1986  END IF
1987  END DO
1988  CALL timing_on('COMM_TOTAL')
1989  CALL mpp_update_domains(u2f, domain)
1990  CALL timing_off('COMM_TOTAL')
1991 !$OMP parallel do default(none) shared(is,ie,js,je,kmax,pm,rf_cutoff,w,rf,u,v, &
1992 !$OMP conserve,hydrostatic,pt,ua,va,u2f,cp,rg,ptop,rcv)
1993  DO k=1,kmax
1994  IF (pm(k) .LT. rf_cutoff) THEN
1995 ! Add heat so as to conserve TE
1996  IF (conserve) THEN
1997  IF (hydrostatic) THEN
1998  DO j=js,je
1999  DO i=is,ie
2000  pt(i, j, k) = pt(i, j, k) + 0.5*(ua(i, j, k)**2+va(i, j&
2001 & , k)**2)*(1.-u2f(i, j, k)**2)/(cp-rg*ptop/pm(k))
2002  END DO
2003  END DO
2004  ELSE
2005  DO j=js,je
2006  DO i=is,ie
2007  pt(i, j, k) = pt(i, j, k) + 0.5*(ua(i, j, k)**2+va(i, j&
2008 & , k)**2+w(i, j, k)**2)*(1.-u2f(i, j, k)**2)*rcv
2009  END DO
2010  END DO
2011  END IF
2012  END IF
2013  DO j=js,je+1
2014  DO i=is,ie
2015  u(i, j, k) = 0.5*(u2f(i, j-1, k)+u2f(i, j, k))*u(i, j, k)
2016  END DO
2017  END DO
2018  DO j=js,je
2019  DO i=is,ie+1
2020  v(i, j, k) = 0.5*(u2f(i-1, j, k)+u2f(i, j, k))*v(i, j, k)
2021  END DO
2022  END DO
2023  IF (.NOT.hydrostatic) THEN
2024  DO j=js,je
2025  DO i=is,ie
2026  w(i, j, k) = u2f(i, j, k)*w(i, j, k)
2027  END DO
2028  END DO
2029  END IF
2030  END IF
2031  END DO
2032  END SUBROUTINE rayleigh_super
2033 ! Differentiation of rayleigh_friction in forward (tangent) mode:
2034 ! variations of useful results: u v w ua delz va pt
2035 ! with respect to varying inputs: u v w ua delz va pt
2036  SUBROUTINE rayleigh_friction_tlm(dt, npx, npy, npz, ks, pm, tau, u, &
2037 & u_tl, v, v_tl, w, w_tl, pt, pt_tl, ua, ua_tl, va, va_tl, delz, &
2038 & delz_tl, cp, rg, ptop, hydrostatic, conserve, rf_cutoff, rf, &
2039 & gridstruct, domain, bd)
2040  IMPLICIT NONE
2041 !deallocate ( u2f )
2042  REAL, INTENT(IN) :: dt
2043 ! time scale (days)
2044  REAL, INTENT(IN) :: tau
2045  REAL, INTENT(IN) :: cp, rg, ptop, rf_cutoff
2046  INTEGER, INTENT(IN) :: npx, npy, npz, ks
2047  REAL, DIMENSION(npz), INTENT(IN) :: pm
2048  LOGICAL, INTENT(IN) :: hydrostatic
2049  LOGICAL, INTENT(IN) :: conserve
2050  TYPE(FV_GRID_BOUNDS_TYPE), INTENT(IN) :: bd
2051 ! D grid zonal wind (m/s)
2052  REAL, INTENT(INOUT) :: u(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
2053  REAL, INTENT(INOUT) :: u_tl(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
2054 ! D grid meridional wind (m/s)
2055  REAL, INTENT(INOUT) :: v(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
2056  REAL, INTENT(INOUT) :: v_tl(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
2057 ! cell center vertical wind (m/s)
2058  REAL, INTENT(INOUT) :: w(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2059  REAL, INTENT(INOUT) :: w_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2060 ! temp
2061  REAL, INTENT(INOUT) :: pt(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2062  REAL, INTENT(INOUT) :: pt_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2063 !
2064  REAL, INTENT(INOUT) :: ua(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2065  REAL, INTENT(INOUT) :: ua_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2066 !
2067  REAL, INTENT(INOUT) :: va(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2068  REAL, INTENT(INOUT) :: va_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2069 ! delta-height (m); non-hydrostatic only
2070  REAL, INTENT(INOUT) :: delz(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2071  REAL, INTENT(INOUT) :: delz_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2072  REAL, INTENT(INOUT) :: rf(npz)
2073  TYPE(FV_GRID_TYPE), INTENT(IN) :: gridstruct
2074  TYPE(DOMAIN2D), INTENT(INOUT) :: domain
2075 ! local:
2076  REAL :: u2f(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2077  REAL :: u2f_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2078  REAL, PARAMETER :: sday=86400.
2079 ! scaling velocity **2
2080  REAL, PARAMETER :: u000=4900.
2081  REAL :: rcv
2082  INTEGER :: i, j, k
2083  INTEGER :: is, ie, js, je
2084  INTEGER :: isd, ied, jsd, jed
2085  INTRINSIC log
2086  INTRINSIC sin
2087  INTRINSIC sqrt
2088  REAL*8 :: arg1
2089  REAL :: arg10
2090  REAL :: arg10_tl
2091  REAL :: result1
2092  REAL :: result1_tl
2093  is = bd%is
2094  ie = bd%ie
2095  js = bd%js
2096  je = bd%je
2097  isd = bd%isd
2098  ied = bd%ied
2099  jsd = bd%jsd
2100  jed = bd%jed
2101  rcv = 1./(cp-rg)
2102  IF (.NOT.rf_initialized) THEN
2103 !allocate( rf(npz) )
2104  rf = 0.0
2105  IF (is_master()) WRITE(6, *) &
2106 & 'Rayleigh friction E-folding time (days):'
2107  DO k=1,npz
2108  IF (pm(k) .LT. rf_cutoff) THEN
2109  arg1 = 0.5*pi*log(rf_cutoff/pm(k))/log(rf_cutoff/ptop)
2110  rf(k) = dt/(tau*sday)*sin(arg1)**2
2111 !if( is_master() ) write(6,*) k, 0.01*pm(k), dt/(rf(k)*sday)
2112  kmax = k
2113  ELSE
2114  EXIT
2115  END IF
2116  END DO
2117  rf_initialized = .true.
2118  END IF
2119 !allocate( u2f(isd:ied,jsd:jed,kmax) )
2120  CALL c2l_ord2_tlm(u, u_tl, v, v_tl, ua, ua_tl, va, va_tl, gridstruct&
2121 & , npz, gridstruct%grid_type, bd, gridstruct%nested)
2122  u2f_tl = 0.0
2123 !$OMP parallel do default(none) shared(is,ie,js,je,kmax,u2f,hydrostatic,ua,va,w)
2124  DO k=1,kmax
2125  IF (hydrostatic) THEN
2126  DO j=js,je
2127  DO i=is,ie
2128  u2f_tl(i, j, k) = 2*ua(i, j, k)*ua_tl(i, j, k) + 2*va(i, j, &
2129 & k)*va_tl(i, j, k)
2130  u2f(i, j, k) = ua(i, j, k)**2 + va(i, j, k)**2
2131  END DO
2132  END DO
2133  ELSE
2134  DO j=js,je
2135  DO i=is,ie
2136  u2f_tl(i, j, k) = 2*ua(i, j, k)*ua_tl(i, j, k) + 2*va(i, j, &
2137 & k)*va_tl(i, j, k) + 2*w(i, j, k)*w_tl(i, j, k)
2138  u2f(i, j, k) = ua(i, j, k)**2 + va(i, j, k)**2 + w(i, j, k)&
2139 & **2
2140  END DO
2141  END DO
2142  END IF
2143  END DO
2144  CALL timing_on('COMM_TOTAL')
2145  CALL mpp_update_domains_tlm(u2f, u2f_tl, domain)
2146  CALL timing_off('COMM_TOTAL')
2147 !$OMP parallel do default(none) shared(is,ie,js,je,kmax,conserve,hydrostatic,pt,u2f,cp,rg, &
2148 !$OMP ptop,pm,rf,delz,rcv,u,v,w)
2149  DO k=1,kmax
2150  IF (conserve) THEN
2151  IF (hydrostatic) THEN
2152  DO j=js,je
2153  DO i=is,ie
2154  arg10_tl = u2f_tl(i, j, k)/u000
2155  arg10 = u2f(i, j, k)/u000
2156  IF (arg10 .EQ. 0.0) THEN
2157  result1_tl = 0.0
2158  ELSE
2159  result1_tl = arg10_tl/(2.0*sqrt(arg10))
2160  END IF
2161  result1 = sqrt(arg10)
2162  pt_tl(i, j, k) = pt_tl(i, j, k) + 0.5*u2f_tl(i, j, k)*(1.-&
2163 & 1./(1.+rf(k)*result1)**2)/(cp-rg*ptop/pm(k)) + 0.5*u2f(i&
2164 & , j, k)*2*rf(k)*result1_tl/((cp-rg*ptop/pm(k))*(1.+rf(k)&
2165 & *result1)**3)
2166  pt(i, j, k) = pt(i, j, k) + 0.5*u2f(i, j, k)/(cp-rg*ptop/&
2167 & pm(k))*(1.-1./(1.+rf(k)*result1)**2)
2168  END DO
2169  END DO
2170  ELSE
2171  DO j=js,je
2172  DO i=is,ie
2173  delz_tl(i, j, k) = (delz_tl(i, j, k)*pt(i, j, k)-delz(i, j&
2174 & , k)*pt_tl(i, j, k))/pt(i, j, k)**2
2175  delz(i, j, k) = delz(i, j, k)/pt(i, j, k)
2176  arg10_tl = u2f_tl(i, j, k)/u000
2177  arg10 = u2f(i, j, k)/u000
2178  IF (arg10 .EQ. 0.0) THEN
2179  result1_tl = 0.0
2180  ELSE
2181  result1_tl = arg10_tl/(2.0*sqrt(arg10))
2182  END IF
2183  result1 = sqrt(arg10)
2184  pt_tl(i, j, k) = pt_tl(i, j, k) + 0.5*rcv*(u2f_tl(i, j, k)&
2185 & *(1.-1./(1.+rf(k)*result1)**2)+u2f(i, j, k)*2*rf(k)*&
2186 & result1_tl/(1.+rf(k)*result1)**3)
2187  pt(i, j, k) = pt(i, j, k) + 0.5*u2f(i, j, k)*rcv*(1.-1./(&
2188 & 1.+rf(k)*result1)**2)
2189  delz_tl(i, j, k) = delz_tl(i, j, k)*pt(i, j, k) + delz(i, &
2190 & j, k)*pt_tl(i, j, k)
2191  delz(i, j, k) = delz(i, j, k)*pt(i, j, k)
2192  END DO
2193  END DO
2194  END IF
2195  END IF
2196  DO j=js-1,je+1
2197  DO i=is-1,ie+1
2198  arg10_tl = u2f_tl(i, j, k)/u000
2199  arg10 = u2f(i, j, k)/u000
2200  IF (arg10 .EQ. 0.0) THEN
2201  result1_tl = 0.0
2202  ELSE
2203  result1_tl = arg10_tl/(2.0*sqrt(arg10))
2204  END IF
2205  result1 = sqrt(arg10)
2206  u2f_tl(i, j, k) = rf(k)*result1_tl
2207  u2f(i, j, k) = rf(k)*result1
2208  END DO
2209  END DO
2210  DO j=js,je+1
2211  DO i=is,ie
2212  u_tl(i, j, k) = (u_tl(i, j, k)*(1.+0.5*(u2f(i, j-1, k)+u2f(i, &
2213 & j, k)))-u(i, j, k)*0.5*(u2f_tl(i, j-1, k)+u2f_tl(i, j, k)))/&
2214 & (1.+0.5*(u2f(i, j-1, k)+u2f(i, j, k)))**2
2215  u(i, j, k) = u(i, j, k)/(1.+0.5*(u2f(i, j-1, k)+u2f(i, j, k)))
2216  END DO
2217  END DO
2218  DO j=js,je
2219  DO i=is,ie+1
2220  v_tl(i, j, k) = (v_tl(i, j, k)*(1.+0.5*(u2f(i-1, j, k)+u2f(i, &
2221 & j, k)))-v(i, j, k)*0.5*(u2f_tl(i-1, j, k)+u2f_tl(i, j, k)))/&
2222 & (1.+0.5*(u2f(i-1, j, k)+u2f(i, j, k)))**2
2223  v(i, j, k) = v(i, j, k)/(1.+0.5*(u2f(i-1, j, k)+u2f(i, j, k)))
2224  END DO
2225  END DO
2226  IF (.NOT.hydrostatic) THEN
2227  DO j=js,je
2228  DO i=is,ie
2229  w_tl(i, j, k) = (w_tl(i, j, k)*(1.+u2f(i, j, k))-w(i, j, k)*&
2230 & u2f_tl(i, j, k))/(1.+u2f(i, j, k))**2
2231  w(i, j, k) = w(i, j, k)/(1.+u2f(i, j, k))
2232  END DO
2233  END DO
2234  END IF
2235  END DO
2236  END SUBROUTINE rayleigh_friction_tlm
2237  SUBROUTINE rayleigh_friction(dt, npx, npy, npz, ks, pm, tau, u, v, w, &
2238 & pt, ua, va, delz, cp, rg, ptop, hydrostatic, conserve, rf_cutoff, rf&
2239 & , gridstruct, domain, bd)
2240  IMPLICIT NONE
2241 !deallocate ( u2f )
2242  REAL, INTENT(IN) :: dt
2243 ! time scale (days)
2244  REAL, INTENT(IN) :: tau
2245  REAL, INTENT(IN) :: cp, rg, ptop, rf_cutoff
2246  INTEGER, INTENT(IN) :: npx, npy, npz, ks
2247  REAL, DIMENSION(npz), INTENT(IN) :: pm
2248  LOGICAL, INTENT(IN) :: hydrostatic
2249  LOGICAL, INTENT(IN) :: conserve
2250  TYPE(FV_GRID_BOUNDS_TYPE), INTENT(IN) :: bd
2251 ! D grid zonal wind (m/s)
2252  REAL, INTENT(INOUT) :: u(bd%isd:bd%ied, bd%jsd:bd%jed+1, npz)
2253 ! D grid meridional wind (m/s)
2254  REAL, INTENT(INOUT) :: v(bd%isd:bd%ied+1, bd%jsd:bd%jed, npz)
2255 ! cell center vertical wind (m/s)
2256  REAL, INTENT(INOUT) :: w(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2257 ! temp
2258  REAL, INTENT(INOUT) :: pt(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2259 !
2260  REAL, INTENT(INOUT) :: ua(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2261 !
2262  REAL, INTENT(INOUT) :: va(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2263 ! delta-height (m); non-hydrostatic only
2264  REAL, INTENT(INOUT) :: delz(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2265  REAL, INTENT(INOUT) :: rf(npz)
2266  TYPE(FV_GRID_TYPE), INTENT(IN) :: gridstruct
2267  TYPE(DOMAIN2D), INTENT(INOUT) :: domain
2268 ! local:
2269  REAL :: u2f(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
2270  REAL, PARAMETER :: sday=86400.
2271 ! scaling velocity **2
2272  REAL, PARAMETER :: u000=4900.
2273  REAL :: rcv
2274  INTEGER :: i, j, k
2275  INTEGER :: is, ie, js, je
2276  INTEGER :: isd, ied, jsd, jed
2277  INTRINSIC log
2278  INTRINSIC sin
2279  INTRINSIC sqrt
2280  REAL*8 :: arg1
2281  REAL :: arg10
2282  REAL :: result1
2283  is = bd%is
2284  ie = bd%ie
2285  js = bd%js
2286  je = bd%je
2287  isd = bd%isd
2288  ied = bd%ied
2289  jsd = bd%jsd
2290  jed = bd%jed
2291  rcv = 1./(cp-rg)
2292  IF (.NOT.rf_initialized) THEN
2293 !allocate( rf(npz) )
2294  rf = 0.0
2295  IF (is_master()) WRITE(6, *) &
2296 & 'Rayleigh friction E-folding time (days):'
2297  DO k=1,npz
2298  IF (pm(k) .LT. rf_cutoff) THEN
2299  arg1 = 0.5*pi*log(rf_cutoff/pm(k))/log(rf_cutoff/ptop)
2300  rf(k) = dt/(tau*sday)*sin(arg1)**2
2301 !if( is_master() ) write(6,*) k, 0.01*pm(k), dt/(rf(k)*sday)
2302  kmax = k
2303  ELSE
2304  GOTO 100
2305  END IF
2306  END DO
2307  100 rf_initialized = .true.
2308  END IF
2309 !allocate( u2f(isd:ied,jsd:jed,kmax) )
2310  CALL c2l_ord2(u, v, ua, va, gridstruct, npz, gridstruct%grid_type, &
2311 & bd, gridstruct%nested)
2312 !$OMP parallel do default(none) shared(is,ie,js,je,kmax,u2f,hydrostatic,ua,va,w)
2313  DO k=1,kmax
2314  IF (hydrostatic) THEN
2315  DO j=js,je
2316  DO i=is,ie
2317  u2f(i, j, k) = ua(i, j, k)**2 + va(i, j, k)**2
2318  END DO
2319  END DO
2320  ELSE
2321  DO j=js,je
2322  DO i=is,ie
2323  u2f(i, j, k) = ua(i, j, k)**2 + va(i, j, k)**2 + w(i, j, k)&
2324 & **2
2325  END DO
2326  END DO
2327  END IF
2328  END DO
2329  CALL timing_on('COMM_TOTAL')
2330  CALL mpp_update_domains(u2f, domain)
2331  CALL timing_off('COMM_TOTAL')
2332 !$OMP parallel do default(none) shared(is,ie,js,je,kmax,conserve,hydrostatic,pt,u2f,cp,rg, &
2333 !$OMP ptop,pm,rf,delz,rcv,u,v,w)
2334  DO k=1,kmax
2335  IF (conserve) THEN
2336  IF (hydrostatic) THEN
2337  DO j=js,je
2338  DO i=is,ie
2339  arg10 = u2f(i, j, k)/u000
2340  result1 = sqrt(arg10)
2341  pt(i, j, k) = pt(i, j, k) + 0.5*u2f(i, j, k)/(cp-rg*ptop/&
2342 & pm(k))*(1.-1./(1.+rf(k)*result1)**2)
2343  END DO
2344  END DO
2345  ELSE
2346  DO j=js,je
2347  DO i=is,ie
2348  delz(i, j, k) = delz(i, j, k)/pt(i, j, k)
2349  arg10 = u2f(i, j, k)/u000
2350  result1 = sqrt(arg10)
2351  pt(i, j, k) = pt(i, j, k) + 0.5*u2f(i, j, k)*rcv*(1.-1./(&
2352 & 1.+rf(k)*result1)**2)
2353  delz(i, j, k) = delz(i, j, k)*pt(i, j, k)
2354  END DO
2355  END DO
2356  END IF
2357  END IF
2358  DO j=js-1,je+1
2359  DO i=is-1,ie+1
2360  arg10 = u2f(i, j, k)/u000
2361  result1 = sqrt(arg10)
2362  u2f(i, j, k) = rf(k)*result1
2363  END DO
2364  END DO
2365  DO j=js,je+1
2366  DO i=is,ie
2367  u(i, j, k) = u(i, j, k)/(1.+0.5*(u2f(i, j-1, k)+u2f(i, j, k)))
2368  END DO
2369  END DO
2370  DO j=js,je
2371  DO i=is,ie+1
2372  v(i, j, k) = v(i, j, k)/(1.+0.5*(u2f(i-1, j, k)+u2f(i, j, k)))
2373  END DO
2374  END DO
2375  IF (.NOT.hydrostatic) THEN
2376  DO j=js,je
2377  DO i=is,ie
2378  w(i, j, k) = w(i, j, k)/(1.+u2f(i, j, k))
2379  END DO
2380  END DO
2381  END IF
2382  END DO
2383  END SUBROUTINE rayleigh_friction
2384 ! Differentiation of compute_aam in forward (tangent) mode:
2385 ! variations of useful results: ua aam va m_fac ps
2386 ! with respect to varying inputs: u v delp ua aam va m_fac ps
2387  SUBROUTINE compute_aam_tlm(npz, is, ie, js, je, isd, ied, jsd, jed, &
2388 & gridstruct, bd, ptop, ua, ua_tl, va, va_tl, u, u_tl, v, v_tl, delp, &
2389 & delp_tl, aam, aam_tl, ps, ps_tl, m_fac, m_fac_tl)
2390  IMPLICIT NONE
2391 ! Compute vertically (mass) integrated Atmospheric Angular Momentum
2392  INTEGER, INTENT(IN) :: npz
2393  INTEGER, INTENT(IN) :: is, ie, js, je
2394  INTEGER, INTENT(IN) :: isd, ied, jsd, jed
2395  REAL, INTENT(IN) :: ptop
2396 ! D grid zonal wind (m/s)
2397  REAL, INTENT(INOUT) :: u(isd:ied, jsd:jed+1, npz)
2398  REAL, INTENT(INOUT) :: u_tl(isd:ied, jsd:jed+1, npz)
2399 ! D grid meridional wind (m/s)
2400  REAL, INTENT(INOUT) :: v(isd:ied+1, jsd:jed, npz)
2401  REAL, INTENT(INOUT) :: v_tl(isd:ied+1, jsd:jed, npz)
2402  REAL, INTENT(INOUT) :: delp(isd:ied, jsd:jed, npz)
2403  REAL, INTENT(INOUT) :: delp_tl(isd:ied, jsd:jed, npz)
2404  REAL, DIMENSION(isd:ied, jsd:jed, npz), INTENT(INOUT) :: ua, va
2405  REAL, DIMENSION(isd:ied, jsd:jed, npz), INTENT(INOUT) :: ua_tl, &
2406 & va_tl
2407  REAL, INTENT(OUT) :: aam(is:ie, js:je)
2408  REAL, INTENT(OUT) :: aam_tl(is:ie, js:je)
2409  REAL, INTENT(OUT) :: m_fac(is:ie, js:je)
2410  REAL, INTENT(OUT) :: m_fac_tl(is:ie, js:je)
2411  REAL, INTENT(OUT) :: ps(isd:ied, jsd:jed)
2412  REAL, INTENT(OUT) :: ps_tl(isd:ied, jsd:jed)
2413  TYPE(FV_GRID_BOUNDS_TYPE), INTENT(IN) :: bd
2414  TYPE(FV_GRID_TYPE), INTENT(IN) :: gridstruct
2415 ! local:
2416  REAL, DIMENSION(is:ie) :: r1, r2, dm
2417  REAL, DIMENSION(is:ie) :: dm_tl
2418  INTEGER :: i, j, k
2419  INTRINSIC cos
2420  CALL c2l_ord2_tlm(u, u_tl, v, v_tl, ua, ua_tl, va, va_tl, gridstruct&
2421 & , npz, gridstruct%grid_type, bd, gridstruct%nested)
2422  dm_tl = 0.0
2423 !$OMP parallel do default(none) shared(is,ie,js,je,npz,gridstruct,aam,m_fac,ps,ptop,delp,agrav,ua) &
2424 !$OMP private(r1, r2, dm)
2425  DO j=js,je
2426  DO i=is,ie
2427  r1(i) = radius*cos(gridstruct%agrid(i, j, 2))
2428  r2(i) = r1(i)*r1(i)
2429  aam_tl(i, j) = 0.0
2430  aam(i, j) = 0.
2431  m_fac_tl(i, j) = 0.0
2432  m_fac(i, j) = 0.
2433  ps_tl(i, j) = 0.0
2434  ps(i, j) = ptop
2435  END DO
2436  DO k=1,npz
2437  DO i=is,ie
2438  dm_tl(i) = delp_tl(i, j, k)
2439  dm(i) = delp(i, j, k)
2440  ps_tl(i, j) = ps_tl(i, j) + dm_tl(i)
2441  ps(i, j) = ps(i, j) + dm(i)
2442  dm_tl(i) = agrav*dm_tl(i)
2443  dm(i) = dm(i)*agrav
2444  aam_tl(i, j) = aam_tl(i, j) + r1(i)*ua_tl(i, j, k)*dm(i) + (r2&
2445 & (i)*omega+r1(i)*ua(i, j, k))*dm_tl(i)
2446  aam(i, j) = aam(i, j) + (r2(i)*omega+r1(i)*ua(i, j, k))*dm(i)
2447  m_fac_tl(i, j) = m_fac_tl(i, j) + r2(i)*dm_tl(i)
2448  m_fac(i, j) = m_fac(i, j) + dm(i)*r2(i)
2449  END DO
2450  END DO
2451  END DO
2452  END SUBROUTINE compute_aam_tlm
2453  SUBROUTINE compute_aam(npz, is, ie, js, je, isd, ied, jsd, jed, &
2454 & gridstruct, bd, ptop, ua, va, u, v, delp, aam, ps, m_fac)
2455  IMPLICIT NONE
2456 ! Compute vertically (mass) integrated Atmospheric Angular Momentum
2457  INTEGER, INTENT(IN) :: npz
2458  INTEGER, INTENT(IN) :: is, ie, js, je
2459  INTEGER, INTENT(IN) :: isd, ied, jsd, jed
2460  REAL, INTENT(IN) :: ptop
2461 ! D grid zonal wind (m/s)
2462  REAL, INTENT(INOUT) :: u(isd:ied, jsd:jed+1, npz)
2463 ! D grid meridional wind (m/s)
2464  REAL, INTENT(INOUT) :: v(isd:ied+1, jsd:jed, npz)
2465  REAL, INTENT(INOUT) :: delp(isd:ied, jsd:jed, npz)
2466  REAL, DIMENSION(isd:ied, jsd:jed, npz), INTENT(INOUT) :: ua, va
2467  REAL, INTENT(OUT) :: aam(is:ie, js:je)
2468  REAL, INTENT(OUT) :: m_fac(is:ie, js:je)
2469  REAL, INTENT(OUT) :: ps(isd:ied, jsd:jed)
2470  TYPE(FV_GRID_BOUNDS_TYPE), INTENT(IN) :: bd
2471  TYPE(FV_GRID_TYPE), INTENT(IN) :: gridstruct
2472 ! local:
2473  REAL, DIMENSION(is:ie) :: r1, r2, dm
2474  INTEGER :: i, j, k
2475  INTRINSIC cos
2476  CALL c2l_ord2(u, v, ua, va, gridstruct, npz, gridstruct%grid_type, &
2477 & bd, gridstruct%nested)
2478 !$OMP parallel do default(none) shared(is,ie,js,je,npz,gridstruct,aam,m_fac,ps,ptop,delp,agrav,ua) &
2479 !$OMP private(r1, r2, dm)
2480  DO j=js,je
2481  DO i=is,ie
2482  r1(i) = radius*cos(gridstruct%agrid(i, j, 2))
2483  r2(i) = r1(i)*r1(i)
2484  aam(i, j) = 0.
2485  m_fac(i, j) = 0.
2486  ps(i, j) = ptop
2487  END DO
2488  DO k=1,npz
2489  DO i=is,ie
2490  dm(i) = delp(i, j, k)
2491  ps(i, j) = ps(i, j) + dm(i)
2492  dm(i) = dm(i)*agrav
2493  aam(i, j) = aam(i, j) + (r2(i)*omega+r1(i)*ua(i, j, k))*dm(i)
2494  m_fac(i, j) = m_fac(i, j) + dm(i)*r2(i)
2495  END DO
2496  END DO
2497  END DO
2498  END SUBROUTINE compute_aam
2499 end module fv_dynamics_tlm_mod
tracer_manager_mod::get_tracer_index
Definition: tracer_manager.F90:124
boundary_tlm_mod
Definition: boundary_tlm.F90:20
fv_mp_nlm_mod
Definition: fv_mp_nlm_mod.F90:24
constants_mod
Definition: constants.F90:24
constants_mod::radius
real, parameter, public radius
Radius of the Earth [m].
Definition: constants.F90:72
field_manager_mod::model_atmos
integer, parameter, public model_atmos
Definition: field_manager.F90:299
fv_dynamics_tlm_mod::pt_initialized
logical pt_initialized
Definition: fv_dynamics_tlm.F90:63
dyn_core_tlm_mod
Definition: dyn_core_tlm.F90:20
fv_grid_utils_tlm_mod
Definition: fv_grid_utils_tlm.F90:20
boundary_tlm_mod::nested_grid_bc_apply_intt
subroutine, public nested_grid_bc_apply_intt(var_nest, istag, jstag, npx, npy, npz, bd, step, split, bc, bctype)
Definition: boundary_tlm.F90:1459
constants_mod::omega
real, parameter, public omega
Rotation rate of the Earth [1/s].
Definition: constants.F90:75
fv_grid_utils_tlm_mod::cubed_to_latlon
subroutine, public cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, km, mode, grid_type, domain, nested, c2l_ord, bd)
Definition: fv_grid_utils_tlm.F90:129
field_manager_mod
Definition: field_manager.F90:20
fv_mapz_tlm_mod::moist_cp
subroutine, public moist_cp(is, ie, isd, ied, jsd, jed, km, j, k, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, q, qd, cpm, t1)
Definition: fv_mapz_tlm.F90:7140
fv_mp_tlm_mod
Definition: fv_mp_tlm.F90:20
fv_dynamics_tlm_mod::rf
real, dimension(:), allocatable rf
Definition: fv_dynamics_tlm.F90:65
constants_mod::pi_8
real(kind=8), parameter, public pi_8
Ratio of circle circumference to diameter [N/A].
Definition: constants.F90:73
fv_arrays_nlm_mod::fv_grid_type
Definition: fv_arrays_nlm.F90:115
fv_tracer2d_tlm_mod::tracer_2d_1l_tlm
subroutine, public tracer_2d_1l_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, dpa)
Definition: fv_tracer2d_tlm.F90:55
fv_arrays_nlm_mod::fv_diag_type
Definition: fv_arrays_nlm.F90:54
fv_dynamics_tlm_mod::idealtest
logical, save, public idealtest
Definition: fv_dynamics_tlm.F90:68
mpp_domains_mod::mpp_update_domains
Definition: mpp_domains.F90:1068
fv_nwp_nudge_nlm_mod
Definition: fv_nudge_nlm.F90:26
constants_mod::hlv
real, parameter, public hlv
Latent heat of evaporation [J/kg].
Definition: constants.F90:80
fv_dynamics_tlm_mod::compute_aam
subroutine compute_aam(npz, is, ie, js, je, isd, ied, jsd, jed, gridstruct, bd, ptop, ua, va, u, v, delp, aam, ps, m_fac)
Definition: fv_dynamics_tlm.F90:2455
fv_mapz_tlm_mod
Definition: fv_mapz_tlm.F90:23
dyn_core_tlm_mod::dyn_core
subroutine, public dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, cappa, grav, hydrostatic, u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, uc, vc, mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, dpx, ks, gridstruct, flagstruct, flagstructp, neststruct, idiag, bd, domain, init_step, i_pack, end_step, gz, pkc, ptc, crx, xfx, cry, yfx, divgd, delpc, ut, vt, zh, pk3, du, dv, time_total)
Definition: dyn_core_tlm.F90:1485
fv_grid_utils_tlm_mod::c2l_ord2_tlm
subroutine, public c2l_ord2_tlm(u, u_tl, v, v_tl, ua, ua_tl, va, va_tl, gridstruct, km, grid_type, bd, do_halo)
Definition: fv_grid_utils_tlm.F90:351
fv_tracer2d_tlm_mod::tracer_2d_tlm
subroutine, public tracer_2d_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, dpa)
Definition: fv_tracer2d_tlm.F90:761
fv_grid_utils_tlm_mod::c2l_ord2
subroutine, public c2l_ord2(u, v, ua, va, gridstruct, km, grid_type, bd, do_halo)
Definition: fv_grid_utils_tlm.F90:453
dyn_core_tlm_mod::del2_cubed
subroutine, public del2_cubed(q, cd, gridstruct, domain, npx, npy, km, nmax, bd)
Definition: dyn_core_tlm.F90:5036
fv_sg_nlm_mod
Definition: fv_sg_nlm.F90:20
constants_mod::rdgas
real, parameter, public rdgas
Gas constant for dry air [J/kg/deg].
Definition: constants.F90:77
dyn_core_tlm_mod::del2_cubed_tlm
subroutine, public del2_cubed_tlm(q, q_tl, cd, gridstruct, domain, npx, npy, km, nmax, bd)
Definition: dyn_core_tlm.F90:4911
fv_arrays_nlm_mod::fv_atmos_type
Definition: fv_arrays_nlm.F90:611
constants_mod::cp_vapor
real, parameter, public cp_vapor
Specific heat capacity of water vapor at constant pressure [J/kg/deg].
Definition: constants.F90:89
fv_mp_tlm_mod::mpp_update_domains_tlm
Definition: fv_mp_tlm.F90:109
fv_dynamics_tlm_mod::rayleigh_friction_tlm
subroutine rayleigh_friction_tlm(dt, npx, npy, npz, ks, pm, tau, u, u_tl, v, v_tl, w, w_tl, pt, pt_tl, ua, ua_tl, va, va_tl, delz, delz_tl, cp, rg, ptop, hydrostatic, conserve, rf_cutoff, rf, gridstruct, domain, bd)
Definition: fv_dynamics_tlm.F90:2040
diag_manager_mod::send_data
Definition: diag_manager.F90:357
mpp_mod
Definition: mpp.F90:39
fv_nesting_tlm_mod
Definition: fv_nesting_tlm.F90:20
fv_timing_nlm_mod
Definition: fv_timing_nlm.F90:20
fv_tracer2d_tlm_mod::tracer_2d_nested
subroutine, public tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, k_split, neststruct, parent_grid, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr)
Definition: fv_tracer2d_tlm.F90:1873
fv_mp_tlm_mod::complete_group_halo_update
subroutine, public complete_group_halo_update(group, group_tl, domain)
Definition: fv_mp_tlm.F90:815
fv_dynamics_tlm_mod::rayleigh_super_tlm
subroutine rayleigh_super_tlm(dt, npx, npy, npz, ks, pm, phis, tau, u, u_tl, v, v_tl, w, w_tl, pt, pt_tl, ua, ua_tl, va, va_tl, delz, agrid, cp, rg, ptop, hydrostatic, conserve, rf_cutoff, rf, gridstruct, domain, bd)
Definition: fv_dynamics_tlm.F90:1753
fv_mapz_tlm_mod::lagrangian_to_eulerian
subroutine, public lagrangian_to_eulerian(last_step, consv, ps, pe, delp, pkz, pk, mdt, pdt, km, is, ie, js, je, isd, ied, jsd, jed, nq, nwat, sphum, q_con, u, v, w, delz, pt, q, hs, r_vir, cp, akap, cappa, kord_mt, kord_wz, kord_tr, kord_tm, peln, te0_2d, ng, ua, va, omga, te, ws, fill, reproduce_sum, out_dt, dtdt, ptop, ak, bk, pfull, flagstruct, gridstruct, domain, do_sat_adj, hydrostatic, hybrid_z, do_omega, adiabatic, do_adiabatic_init, mfx, mfy, remap_option, kord_mt_pert, kord_wz_pert, kord_tr_pert, kord_tm_pert)
Definition: fv_mapz_tlm.F90:1369
fv_tracer2d_tlm_mod
Definition: fv_tracer2d_tlm.F90:20
mpp_domains_mod
Definition: mpp_domains.F90:128
fv_mp_tlm_mod::start_group_halo_update_tlm
Definition: fv_mp_tlm.F90:64
fv_grid_utils_tlm_mod::g_sum
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
Definition: fv_grid_utils_tlm.F90:596
mpp_domains_mod::domain2d
Definition: mpp_domains.F90:378
dyn_core_tlm_mod::dyn_core_tlm
subroutine, public dyn_core_tlm(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, cappa, grav, hydrostatic, u, u_tl, v, v_tl, w, w_tl, delz, delz_tl, pt, pt_tl, q, q_tl, delp, delp_tl, pe, pe_tl, pk, pk_tl, phis, ws, ws_tl, omga, omga_tl, ptop, pfull, ua, ua_tl, va, va_tl, uc, uc_tl, vc, vc_tl, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl, cy, cy_tl, pkz, pkz_tl, peln, peln_tl, q_con, ak, bk, dpx, dpx_tl, ks, gridstruct, flagstruct, flagstructp, neststruct, idiag, bd, domain, init_step, i_pack, end_step, gz, gz_tl, pkc, pkc_tl, ptc, ptc_tl, crx, crx_tl, xfx, xfx_tl, cry, cry_tl, yfx, yfx_tl, divgd, divgd_tl, delpc, delpc_tl, ut, ut_tl, vt, vt_tl, zh, zh_tl, pk3, pk3_tl, du, du_tl, dv, dv_tl, time_total)
Definition: dyn_core_tlm.F90:104
constants_mod::rvgas
real, parameter, public rvgas
Gas constant for water vapor [J/kg/deg].
Definition: constants.F90:78
fv_mapz_tlm_mod::lagrangian_to_eulerian_tlm
subroutine, public lagrangian_to_eulerian_tlm(last_step, consv, ps, ps_tl, pe, pe_tl, delp, delp_tl, pkz, pkz_tl, pk, pk_tl, mdt, pdt, km, is, ie, js, je, isd, ied, jsd, jed, nq, nwat, sphum, q_con, u, u_tl, v, v_tl, w, w_tl, delz, delz_tl, pt, pt_tl, q, q_tl, hs, r_vir, cp, akap, cappa, kord_mt, kord_wz, kord_tr, kord_tm, peln, peln_tl, te0_2d, te0_2d_tl, ng, ua, ua_tl, va, omga, omga_tl, te, te_tl, ws, ws_tl, fill, reproduce_sum, out_dt, dtdt, ptop, ak, bk, pfull, flagstruct, gridstruct, domain, do_sat_adj, hydrostatic, hybrid_z, do_omega, adiabatic, do_adiabatic_init, mfx, mfy, remap_option, kord_mt_pert, kord_wz_pert, kord_tr_pert, kord_tm_pert)
Definition: fv_mapz_tlm.F90:79
fv_nesting_tlm_mod::setup_nested_grid_bcs_tlm
subroutine, public setup_nested_grid_bcs_tlm(npx, npy, npz, zvir, ncnst, u, u_tl, v, v_tl, w, pt, delp, delz, q, uc, uc_tl, vc, vc_tl, pkz, nested, inline_q, make_nh, ng, gridstruct, flagstruct, neststruct, nest_timestep, tracer_nest_timestep, domain, bd, nwat)
Definition: fv_nesting_tlm.F90:79
fv_dynamics_tlm_mod
Definition: fv_dynamics_tlm.F90:20
fv_dynamics_tlm_mod::bad_range
logical bad_range
Definition: fv_dynamics_tlm.F90:64
fv_timing_nlm_mod::timing_on
subroutine timing_on(blk_name)
Definition: fv_timing_nlm.F90:101
fv_mapz_tlm_mod::compute_total_energy
subroutine, public compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, u, v, w, delz, pt, delp, q, qc, pe, peln, hs, rsin2_l, cosa_s_l, r_vir, cp, rg, hlv, te_2d, ua, va, teq, moist_phys, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hydrostatic, id_te)
Definition: fv_mapz_tlm.F90:2493
fv_arrays_nlm_mod::fv_flags_type
Definition: fv_arrays_nlm.F90:236
fv_mapz_tlm_mod::compute_total_energy_tlm
subroutine, public compute_total_energy_tlm(is, ie, js, je, isd, ied, jsd, jed, km, u, u_tl, v, v_tl, w, w_tl, delz, delz_tl, pt, pt_tl, delp, delp_tl, q, q_tl, qc, qc_tl, pe, pe_tl, peln, peln_tl, hs, rsin2_l, cosa_s_l, r_vir, cp, rg, hlv, te_2d, te_2d_tl, ua, va, teq, teq_tl, moist_phys, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hydrostatic, id_te)
Definition: fv_mapz_tlm.F90:2297
fv_arrays_tlmadm_mod
Definition: fv_arrays_tlmadm.F90:20
fv_arrays_nlm_mod
Definition: fv_arrays_nlm.F90:20
fv_dynamics_tlm_mod::kmax
integer kmax
Definition: fv_dynamics_tlm.F90:66
fv_dynamics_tlm_mod::agrav
real agrav
Definition: fv_dynamics_tlm.F90:67
fv_arrays_tlmadm_mod::fpp
type(fpp_type) fpp
Definition: fv_arrays_tlmadm.F90:35
diag_manager_mod
Definition: diag_manager.F90:20
fv_nwp_nudge_nlm_mod::do_adiabatic_init
logical, public do_adiabatic_init
Definition: fv_nudge_nlm.F90:57
dyn_core_tlm_mod::init_ijk_mem
subroutine, public init_ijk_mem(i1, i2, j1, j2, km, array, var)
Definition: dyn_core_tlm.F90:5132
fv_dynamics_tlm_mod::compute_aam_tlm
subroutine compute_aam_tlm(npz, is, ie, js, je, isd, ied, jsd, jed, gridstruct, bd, ptop, ua, ua_tl, va, va_tl, u, u_tl, v, v_tl, delp, delp_tl, aam, aam_tl, ps, ps_tl, m_fac, m_fac_tl)
Definition: fv_dynamics_tlm.F90:2390
fv_dynamics_tlm_mod::rayleigh_super
subroutine rayleigh_super(dt, npx, npy, npz, ks, pm, phis, tau, u, v, w, pt, ua, va, delz, agrid, cp, rg, ptop, hydrostatic, conserve, rf_cutoff, rf, gridstruct, domain, bd)
Definition: fv_dynamics_tlm.F90:1903
fv_tracer2d_tlm_mod::tracer_2d
subroutine, public tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, dpa)
Definition: fv_tracer2d_tlm.F90:1151
fv_mp_tlm_mod::start_group_halo_update
Definition: fv_mp_tlm.F90:56
fv_dynamics_tlm_mod::rf_initialized
logical rf_initialized
Definition: fv_dynamics_tlm.F90:62
constants_mod::grav
real, parameter, public grav
Acceleration due to gravity [m/s^2].
Definition: constants.F90:76
tracer_manager_mod
Definition: tracer_manager.F90:20
fv_arrays_nlm_mod::fv_nest_type
Definition: fv_arrays_nlm.F90:534
fv_mapz_tlm_mod::moist_cv
subroutine, public moist_cv(is, ie, isd, ied, jsd, jed, km, j, k, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, q, qd, cvm, t1)
Definition: fv_mapz_tlm.F90:7050
fv_sg_nlm_mod::neg_adj3
subroutine, public neg_adj3(is, ie, js, je, ng, kbot, hydrostatic, peln, delz, pt, dp, qv, ql, qr, qi, qs, qg, qa, check_negative)
Definition: fv_sg_nlm.F90:1132
fv_dynamics_tlm_mod::fv_dynamics
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, flagstructp, neststruct, idiag, bd, parent_grid, domain, time_total)
Definition: fv_dynamics_tlm.F90:1005
fv_fill_nlm_mod
Definition: fv_fill_nlm.F90:20
boundary_tlm_mod::nested_grid_bc_apply_intt_tlm
subroutine, public nested_grid_bc_apply_intt_tlm(var_nest, var_nest_tl, istag, jstag, npx, npy, npz, bd, step, split, bc, bctype)
Definition: boundary_tlm.F90:1356
fv_diagnostics_nlm_mod::prt_mxm
subroutine, public prt_mxm(qname, q, is, ie, js, je, n_g, km, fac, area, domain)
Definition: fv_diagnostics_nlm.F90:2809
fv_tracer2d_tlm_mod::tracer_2d_1l
subroutine, public tracer_2d_1l(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, dpa)
Definition: fv_tracer2d_tlm.F90:455
fv_arrays_tlmadm_mod::fv_flags_pert_type
Definition: fv_arrays_tlmadm.F90:37
fv_fill_nlm_mod::fill2d
subroutine, public fill2d(is, ie, js, je, ng, km, q, delp, area, domain, nested, npx, npy)
Definition: fv_fill_nlm.F90:183
fv_diagnostics_nlm_mod::prt_minmax
logical, public prt_minmax
Definition: fv_diagnostics_nlm.F90:65
fv_diagnostics_nlm_mod
Definition: fv_diagnostics_nlm.F90:20
fv_diagnostics_nlm_mod::range_check
subroutine, public range_check(qname, q, is, ie, js, je, n_g, km, pos, q_low, q_hi, bad_range)
Definition: fv_diagnostics_nlm.F90:2714
fv_diagnostics_nlm_mod::fv_time
type(time_type), public fv_time
Definition: fv_diagnostics_nlm.F90:61
fv_arrays_nlm_mod::fv_grid_bounds_type
Definition: fv_arrays_nlm.F90:601
fv_nesting_tlm_mod::setup_nested_grid_bcs
subroutine, public setup_nested_grid_bcs(npx, npy, npz, zvir, ncnst, u, v, w, pt, delp, delz, q, uc, vc, pkz, nested, inline_q, make_nh, ng, gridstruct, flagstruct, neststruct, nest_timestep, tracer_nest_timestep, domain, bd, nwat)
Definition: fv_nesting_tlm.F90:342
fv_grid_utils_tlm_mod::g_sum_tlm
real function, public g_sum_tlm(domain, p, p_tl, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce, g_sum)
Definition: fv_grid_utils_tlm.F90:532
fv_tracer2d_tlm_mod::tracer_2d_nested_tlm
subroutine, public tracer_2d_nested_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, k_split, neststruct, parent_grid, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr)
Definition: fv_tracer2d_tlm.F90:1455
crtm_parameters::pi
real(fp), parameter, public pi
Definition: CRTM_Parameters.f90:86
fv_dynamics_tlm_mod::rayleigh_friction
subroutine rayleigh_friction(dt, npx, npy, npz, ks, pm, tau, u, v, w, pt, ua, va, delz, cp, rg, ptop, hydrostatic, conserve, rf_cutoff, rf, gridstruct, domain, bd)
Definition: fv_dynamics_tlm.F90:2240
fv_timing_nlm_mod::timing_off
subroutine timing_off(blk_name)
Definition: fv_timing_nlm.F90:175
fv_dynamics_tlm_mod::fv_dynamics_tlm
subroutine, public fv_dynamics_tlm(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, reproduce_sum, kappa, cp_air, zvir, ptop, ks, ncnst, n_split, q_split, u, u_tl, v, v_tl, w, w_tl, delz, delz_tl, hydrostatic, pt, pt_tl, delp, delp_tl, q, q_tl, ps, ps_tl, pe, pe_tl, pk, pk_tl, peln, peln_tl, pkz, pkz_tl, phis, q_con, omga, omga_tl, ua, ua_tl, va, va_tl, uc, uc_tl, vc, vc_tl, ak, bk, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl, cy, cy_tl, ze0, hybrid_z, gridstruct, flagstruct, flagstructp, neststruct, idiag, bd, parent_grid, domain, time_total)
Definition: fv_dynamics_tlm.F90:95