Line data Source code
1 : module dyn_comp
2 :
3 : ! CAM interfaces to the SE Dynamical Core
4 :
5 : use shr_kind_mod, only: r8=>shr_kind_r8, shr_kind_cl
6 : use physconst, only: pi
7 : use spmd_utils, only: iam, masterproc
8 : use constituents, only: pcnst, cnst_get_ind, cnst_name, cnst_longname, &
9 : cnst_read_iv, qmin, cnst_type, tottnam, &
10 : cnst_is_a_water_species
11 : use cam_control_mod, only: initial_run
12 : use cam_initfiles, only: initial_file_get_id, topo_file_get_id, pertlim
13 : use phys_control, only: use_gw_front, use_gw_front_igw
14 : use dyn_grid, only: ini_grid_name, timelevel, hvcoord, edgebuf, &
15 : ini_grid_hdim_name
16 :
17 : use cam_grid_support, only: cam_grid_id, cam_grid_get_gcid, &
18 : cam_grid_dimensions, &
19 : cam_grid_get_latvals, cam_grid_get_lonvals, &
20 : max_hcoordname_len
21 : use cam_map_utils, only: iMap
22 :
23 : use inic_analytic, only: analytic_ic_active, analytic_ic_set_ic
24 : use dyn_tests_utils, only: vcoord=>vc_dry_pressure
25 :
26 : use cam_history, only: outfld, hist_fld_active, fieldname_len
27 : use cam_history_support, only: max_fieldname_len
28 : use time_manager, only: get_step_size
29 :
30 : use ncdio_atm, only: infld
31 : use pio, only: file_desc_t, pio_seterrorhandling, PIO_BCAST_ERROR, &
32 : pio_inq_dimid, pio_inq_dimlen, PIO_NOERR
33 :
34 : use infnan, only: isnan
35 : use cam_logfile, only: iulog
36 : use cam_abortutils, only: endrun
37 : use shr_sys_mod, only: shr_sys_flush
38 :
39 : use parallel_mod, only: par
40 : use hybrid_mod, only: hybrid_t
41 : use dimensions_mod, only: nelemd, nlev, np, npsq, ntrac, nc, fv_nphys
42 : use dimensions_mod, only: qsize, use_cslam
43 : use element_mod, only: element_t, elem_state_t
44 : use fvm_control_volume_mod, only: fvm_struct
45 : use se_dyn_time_mod, only: nsplit
46 : use edge_mod, only: initEdgeBuffer, edgeVpack, edgeVunpack, FreeEdgeBuffer
47 : use edgetype_mod, only: EdgeBuffer_t
48 : use bndry_mod, only: bndry_exchange
49 :
50 : implicit none
51 : private
52 : save
53 :
54 : public :: &
55 : dyn_import_t, &
56 : dyn_export_t, &
57 : dyn_readnl, &
58 : dyn_register, &
59 : dyn_init, &
60 : dyn_run, &
61 : dyn_final
62 :
63 : type dyn_import_t
64 : type (element_t), pointer :: elem(:) => null()
65 : type (fvm_struct), pointer :: fvm(:) => null()
66 : end type dyn_import_t
67 :
68 : type dyn_export_t
69 : type (element_t), pointer :: elem(:) => null()
70 : type (fvm_struct), pointer :: fvm(:) => null()
71 : end type dyn_export_t
72 :
73 : ! Namelist
74 : logical, public, protected :: write_restart_unstruct
75 :
76 : ! Frontogenesis indices
77 : integer, public :: frontgf_idx = -1
78 : integer, public :: frontga_idx = -1
79 :
80 : interface read_dyn_var
81 : module procedure read_dyn_field_2d
82 : module procedure read_dyn_field_3d
83 : end interface read_dyn_var
84 :
85 : real(r8), parameter :: rad2deg = 180.0_r8 / pi
86 : real(r8), parameter :: deg2rad = pi / 180.0_r8
87 : real(r8), parameter :: rarea_sphere = 1.0_r8 / (4.0_r8*PI)
88 :
89 : !===============================================================================
90 : contains
91 : !===============================================================================
92 :
93 1536 : subroutine dyn_readnl(NLFileName)
94 : use air_composition,only: thermodynamic_active_species_num
95 : use namelist_utils, only: find_group_name
96 : use namelist_mod, only: homme_set_defaults, homme_postprocess_namelist
97 : use units, only: getunit, freeunit
98 : use spmd_utils, only: masterproc, masterprocid, mpicom, npes
99 : use spmd_utils, only: mpi_real8, mpi_integer, mpi_character, mpi_logical
100 : use dyn_grid, only: se_write_grid_file, se_grid_filename, se_write_gll_corners
101 : use dp_mapping, only: nphys_pts
102 : use native_mapping, only: native_mapping_readnl
103 :
104 : use control_mod, only: hypervis_subcycle, hypervis_subcycle_sponge
105 : use control_mod, only: hypervis_subcycle_q, statefreq, runtype
106 : use control_mod, only: nu, nu_div, nu_p, nu_q, nu_top, qsplit, rsplit
107 : use control_mod, only: vert_remap_uvTq_alg, vert_remap_tracer_alg
108 : use control_mod, only: tstep_type, rk_stage_user
109 : use control_mod, only: ftype, limiter_option, partmethod
110 : use control_mod, only: topology, variable_nsplit
111 : use control_mod, only: fine_ne, hypervis_power, hypervis_scaling
112 : use control_mod, only: max_hypervis_courant, statediag_numtrac,refined_mesh
113 : use control_mod, only: molecular_diff, pgf_formulation, dribble_in_rsplit_loop
114 : use control_mod, only: sponge_del4_nu_div_fac, sponge_del4_nu_fac, sponge_del4_lev
115 : use dimensions_mod, only: ne, npart
116 : use dimensions_mod, only: large_Courant_incr
117 : use dimensions_mod, only: fvm_supercycling, fvm_supercycling_jet
118 : use dimensions_mod, only: kmin_jet, kmax_jet
119 : use params_mod, only: SFCURVE
120 : use parallel_mod, only: initmpi
121 : use thread_mod, only: initomp, max_num_threads
122 : use thread_mod, only: horz_num_threads, vert_num_threads, tracer_num_threads
123 : ! Dummy argument
124 : character(len=*), intent(in) :: NLFileName
125 :
126 : ! Local variables
127 : integer :: unitn, ierr,k
128 :
129 : ! SE Namelist variables
130 : integer :: se_fine_ne
131 : integer :: se_ftype
132 : integer :: se_statediag_numtrac
133 : integer :: se_fv_nphys
134 : real(r8) :: se_hypervis_power
135 : real(r8) :: se_hypervis_scaling
136 : integer :: se_hypervis_subcycle
137 : integer :: se_hypervis_subcycle_sponge
138 : integer :: se_hypervis_subcycle_q
139 : integer :: se_limiter_option
140 : real(r8) :: se_max_hypervis_courant
141 : character(len=SHR_KIND_CL) :: se_mesh_file
142 : integer :: se_ne
143 : integer :: se_npes
144 : integer :: se_nsplit
145 : real(r8) :: se_nu
146 : real(r8) :: se_nu_div
147 : real(r8) :: se_nu_p
148 : real(r8) :: se_nu_top
149 : real(r8) :: se_sponge_del4_nu_fac
150 : real(r8) :: se_sponge_del4_nu_div_fac
151 : integer :: se_sponge_del4_lev
152 : integer :: se_qsplit
153 : logical :: se_refined_mesh
154 : integer :: se_rsplit
155 : integer :: se_statefreq
156 : integer :: se_tstep_type
157 : character(len=32) :: se_vert_remap_T
158 : character(len=32) :: se_vert_remap_uvTq_alg
159 : character(len=32) :: se_vert_remap_tracer_alg
160 : integer :: se_horz_num_threads
161 : integer :: se_vert_num_threads
162 : integer :: se_tracer_num_threads
163 : logical :: se_write_restart_unstruct
164 : logical :: se_large_Courant_incr
165 : integer :: se_fvm_supercycling
166 : integer :: se_fvm_supercycling_jet
167 : integer :: se_kmin_jet
168 : integer :: se_kmax_jet
169 : real(r8) :: se_molecular_diff
170 : integer :: se_pgf_formulation
171 : integer :: se_dribble_in_rsplit_loop
172 : namelist /dyn_se_inparm/ &
173 : se_fine_ne, & ! For refined meshes
174 : se_ftype, & ! forcing type
175 : se_statediag_numtrac, &
176 : se_fv_nphys, &
177 : se_hypervis_power, &
178 : se_hypervis_scaling, &
179 : se_hypervis_subcycle, &
180 : se_hypervis_subcycle_sponge, &
181 : se_hypervis_subcycle_q, &
182 : se_limiter_option, &
183 : se_max_hypervis_courant, &
184 : se_mesh_file, & ! Refined mesh definition file
185 : se_ne, &
186 : se_npes, &
187 : se_nsplit, & ! # of dyn steps per physics timestep
188 : se_nu, &
189 : se_nu_div, &
190 : se_nu_p, &
191 : se_nu_top, &
192 : se_sponge_del4_nu_fac, &
193 : se_sponge_del4_nu_div_fac, &
194 : se_sponge_del4_lev, &
195 : se_qsplit, &
196 : se_refined_mesh, &
197 : se_rsplit, &
198 : se_statefreq, & ! number of steps per printstate call
199 : se_tstep_type, &
200 : se_vert_remap_T, &
201 : se_vert_remap_uvTq_alg, &
202 : se_vert_remap_tracer_alg, &
203 : se_write_grid_file, &
204 : se_grid_filename, &
205 : se_write_gll_corners, &
206 : se_horz_num_threads, &
207 : se_vert_num_threads, &
208 : se_tracer_num_threads, &
209 : se_write_restart_unstruct, &
210 : se_large_Courant_incr, &
211 : se_fvm_supercycling, &
212 : se_fvm_supercycling_jet, &
213 : se_kmin_jet, &
214 : se_kmax_jet, &
215 : se_molecular_diff, &
216 : se_pgf_formulation, &
217 : se_dribble_in_rsplit_loop
218 : !--------------------------------------------------------------------------
219 :
220 : ! defaults for variables not set by build-namelist
221 1536 : se_fine_ne = -1
222 1536 : se_hypervis_power = 0
223 1536 : se_hypervis_scaling = 0
224 1536 : se_max_hypervis_courant = 1.0e99_r8
225 1536 : se_mesh_file = ''
226 1536 : se_npes = npes
227 1536 : se_write_restart_unstruct = .false.
228 :
229 : ! Read the namelist (dyn_se_inparm)
230 1536 : call MPI_barrier(mpicom, ierr)
231 1536 : if (masterproc) then
232 2 : write(iulog, *) "dyn_readnl: reading dyn_se_inparm namelist..."
233 2 : unitn = getunit()
234 2 : open( unitn, file=trim(NLFileName), status='old' )
235 2 : call find_group_name(unitn, 'dyn_se_inparm', status=ierr)
236 2 : if (ierr == 0) then
237 2 : read(unitn, dyn_se_inparm, iostat=ierr)
238 2 : if (ierr /= 0) then
239 0 : call endrun('dyn_readnl: ERROR reading dyn_se_inparm namelist')
240 : end if
241 : end if
242 2 : close(unitn)
243 2 : call freeunit(unitn)
244 : end if
245 :
246 : ! Broadcast namelist values to all PEs
247 1536 : call MPI_bcast(se_fine_ne, 1, mpi_integer, masterprocid, mpicom, ierr)
248 1536 : call MPI_bcast(se_ftype, 1, mpi_integer, masterprocid, mpicom, ierr)
249 1536 : call MPI_bcast(se_statediag_numtrac, 1, mpi_integer, masterprocid, mpicom, ierr)
250 1536 : call MPI_bcast(se_hypervis_power, 1, mpi_real8, masterprocid, mpicom, ierr)
251 1536 : call MPI_bcast(se_hypervis_scaling, 1, mpi_real8, masterprocid, mpicom, ierr)
252 1536 : call MPI_bcast(se_hypervis_subcycle, 1, mpi_integer, masterprocid, mpicom, ierr)
253 1536 : call MPI_bcast(se_hypervis_subcycle_sponge, 1, mpi_integer, masterprocid, mpicom, ierr)
254 1536 : call MPI_bcast(se_hypervis_subcycle_q, 1, mpi_integer, masterprocid, mpicom, ierr)
255 1536 : call MPI_bcast(se_limiter_option, 1, mpi_integer, masterprocid, mpicom, ierr)
256 1536 : call MPI_bcast(se_max_hypervis_courant, 1, mpi_real8, masterprocid, mpicom, ierr)
257 1536 : call MPI_bcast(se_mesh_file, SHR_KIND_CL, mpi_character, masterprocid, mpicom, ierr)
258 1536 : call MPI_bcast(se_ne, 1, mpi_integer, masterprocid, mpicom, ierr)
259 1536 : call MPI_bcast(se_npes, 1, mpi_integer, masterprocid, mpicom, ierr)
260 1536 : call MPI_bcast(se_nsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
261 1536 : call MPI_bcast(se_nu, 1, mpi_real8, masterprocid, mpicom, ierr)
262 1536 : call MPI_bcast(se_nu_div, 1, mpi_real8, masterprocid, mpicom, ierr)
263 1536 : call MPI_bcast(se_nu_p, 1, mpi_real8, masterprocid, mpicom, ierr)
264 1536 : call MPI_bcast(se_nu_top, 1, mpi_real8, masterprocid, mpicom, ierr)
265 1536 : call MPI_bcast(se_sponge_del4_nu_fac, 1, mpi_real8, masterprocid, mpicom, ierr)
266 1536 : call MPI_bcast(se_sponge_del4_nu_div_fac, 1, mpi_real8, masterprocid, mpicom, ierr)
267 1536 : call MPI_bcast(se_sponge_del4_lev, 1, mpi_integer, masterprocid, mpicom, ierr)
268 1536 : call MPI_bcast(se_qsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
269 1536 : call MPI_bcast(se_refined_mesh, 1, mpi_logical, masterprocid, mpicom, ierr)
270 1536 : call MPI_bcast(se_rsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
271 1536 : call MPI_bcast(se_statefreq, 1, mpi_integer, masterprocid, mpicom, ierr)
272 1536 : call MPI_bcast(se_tstep_type, 1, mpi_integer, masterprocid, mpicom, ierr)
273 1536 : call MPI_bcast(se_vert_remap_T, 32, mpi_character, masterprocid, mpicom, ierr)
274 1536 : call MPI_bcast(se_vert_remap_uvTq_alg, 32, mpi_character, masterprocid, mpicom, ierr)
275 1536 : call MPI_bcast(se_vert_remap_tracer_alg, 32, mpi_character, masterprocid, mpicom, ierr)
276 1536 : call MPI_bcast(se_fv_nphys, 1, mpi_integer, masterprocid, mpicom, ierr)
277 1536 : call MPI_bcast(se_write_grid_file, 16, mpi_character, masterprocid, mpicom, ierr)
278 1536 : call MPI_bcast(se_grid_filename, shr_kind_cl, mpi_character, masterprocid, mpicom, ierr)
279 1536 : call MPI_bcast(se_write_gll_corners, 1, mpi_logical, masterprocid, mpicom, ierr)
280 1536 : call MPI_bcast(se_horz_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
281 1536 : call MPI_bcast(se_vert_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
282 1536 : call MPI_bcast(se_tracer_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
283 1536 : call MPI_bcast(se_write_restart_unstruct, 1, mpi_logical, masterprocid, mpicom, ierr)
284 1536 : call MPI_bcast(se_large_Courant_incr, 1, mpi_logical, masterprocid, mpicom, ierr)
285 1536 : call MPI_bcast(se_fvm_supercycling, 1, mpi_integer, masterprocid, mpicom, ierr)
286 1536 : call MPI_bcast(se_fvm_supercycling_jet, 1, mpi_integer, masterprocid, mpicom, ierr)
287 1536 : call MPI_bcast(se_kmin_jet, 1, mpi_integer, masterprocid, mpicom, ierr)
288 1536 : call MPI_bcast(se_kmax_jet, 1, mpi_integer, masterprocid, mpicom, ierr)
289 1536 : call MPI_bcast(se_molecular_diff, 1, mpi_real8, masterprocid, mpicom, ierr)
290 1536 : call MPI_bcast(se_pgf_formulation, 1, mpi_integer, masterprocid, mpicom, ierr)
291 1536 : call MPI_bcast(se_dribble_in_rsplit_loop, 1, mpi_integer, masterprocid, mpicom, ierr)
292 1536 : if (se_npes <= 0) then
293 0 : call endrun('dyn_readnl: ERROR: se_npes must be > 0')
294 : end if
295 :
296 : ! Initialize the SE structure that holds the MPI decomposition information
297 1536 : par = initmpi(se_npes)
298 1536 : call initomp()
299 :
300 :
301 1536 : if (se_fvm_supercycling < 0) se_fvm_supercycling = se_rsplit
302 1536 : if (se_fvm_supercycling_jet < 0) se_fvm_supercycling_jet = se_rsplit
303 :
304 : ! Go ahead and enforce ne = 0 for refined mesh runs
305 1536 : if (se_refined_mesh) then
306 0 : se_ne = 0
307 : end if
308 :
309 : ! Set HOMME defaults
310 1536 : call homme_set_defaults()
311 : ! Set HOMME variables not in CAM's namelist but with different CAM defaults
312 1536 : partmethod = SFCURVE
313 1536 : npart = se_npes
314 : ! CAM requires forward-in-time, subcycled dynamics
315 : ! RK2 3 stage tracers, sign-preserving conservative
316 1536 : rk_stage_user = 3
317 1536 : topology = "cube"
318 : ! Finally, set the HOMME variables which have different names
319 1536 : fine_ne = se_fine_ne
320 1536 : ftype = se_ftype
321 1536 : statediag_numtrac = MIN(se_statediag_numtrac,pcnst)
322 1536 : hypervis_power = se_hypervis_power
323 1536 : hypervis_scaling = se_hypervis_scaling
324 1536 : hypervis_subcycle = se_hypervis_subcycle
325 1536 : if (hypervis_subcycle_sponge<0) then
326 0 : hypervis_subcycle_sponge = hypervis_subcycle
327 : else
328 1536 : hypervis_subcycle_sponge = se_hypervis_subcycle_sponge
329 : end if
330 1536 : hypervis_subcycle_q = se_hypervis_subcycle_q
331 1536 : limiter_option = se_limiter_option
332 1536 : max_hypervis_courant = se_max_hypervis_courant
333 1536 : refined_mesh = se_refined_mesh
334 1536 : ne = se_ne
335 1536 : nsplit = se_nsplit
336 1536 : nu = se_nu
337 1536 : nu_div = se_nu_div
338 1536 : nu_p = se_nu_p
339 1536 : nu_q = se_nu_p !for tracer-wind consistency nu_q must me equal to nu_p
340 1536 : nu_top = se_nu_top
341 1536 : sponge_del4_nu_fac = se_sponge_del4_nu_fac
342 1536 : sponge_del4_nu_div_fac = se_sponge_del4_nu_div_fac
343 1536 : sponge_del4_lev = se_sponge_del4_lev
344 1536 : qsplit = se_qsplit
345 1536 : rsplit = se_rsplit
346 1536 : statefreq = se_statefreq
347 1536 : tstep_type = se_tstep_type
348 1536 : vert_remap_uvTq_alg = set_vert_remap(se_vert_remap_T, se_vert_remap_uvTq_alg)
349 1536 : vert_remap_tracer_alg = set_vert_remap(se_vert_remap_T, se_vert_remap_tracer_alg)
350 1536 : fv_nphys = se_fv_nphys
351 1536 : large_Courant_incr = se_large_Courant_incr
352 1536 : fvm_supercycling = se_fvm_supercycling
353 1536 : fvm_supercycling_jet = se_fvm_supercycling_jet
354 1536 : kmin_jet = se_kmin_jet
355 1536 : kmax_jet = se_kmax_jet
356 1536 : variable_nsplit = .false.
357 1536 : molecular_diff = se_molecular_diff
358 1536 : pgf_formulation = se_pgf_formulation
359 1536 : dribble_in_rsplit_loop = se_dribble_in_rsplit_loop
360 1536 : if (fv_nphys > 0) then
361 : ! Use finite volume physics grid and CSLAM for tracer advection
362 1536 : nphys_pts = fv_nphys*fv_nphys
363 1536 : qsize = thermodynamic_active_species_num ! number tracers advected by GLL
364 1536 : ntrac = pcnst ! number tracers advected by CSLAM
365 1536 : use_cslam = .true.
366 : else
367 : ! Use GLL grid for physics and tracer advection
368 0 : nphys_pts = npsq
369 0 : qsize = pcnst
370 0 : ntrac = 0
371 0 : use_cslam = .false.
372 : end if
373 :
374 1536 : if (rsplit < 1) then
375 0 : call endrun('dyn_readnl: rsplit must be > 0')
376 : end if
377 :
378 : ! if restart or branch run
379 1536 : if (.not. initial_run) then
380 768 : runtype = 1
381 : end if
382 :
383 : ! HOMME wants 'none' to indicate no mesh file
384 1536 : if (len_trim(se_mesh_file) == 0) then
385 1536 : se_mesh_file = 'none'
386 1536 : if (se_refined_mesh) then
387 0 : call endrun('dyn_readnl ERROR: se_refined_mesh=.true. but no se_mesh_file')
388 : end if
389 : end if
390 1536 : call homme_postprocess_namelist(se_mesh_file, par)
391 :
392 : ! Set threading numbers to reasonable values
393 1536 : if ((se_horz_num_threads == 0) .and. (se_vert_num_threads == 0) .and. (se_tracer_num_threads == 0)) then
394 : ! The user has not set any threading values, choose defaults
395 1536 : se_horz_num_threads = 1
396 1536 : se_vert_num_threads = max_num_threads
397 1536 : se_tracer_num_threads = se_vert_num_threads
398 : end if
399 1536 : if (se_horz_num_threads < 1) then
400 0 : se_horz_num_threads = 1
401 : end if
402 1536 : if (se_vert_num_threads < 1) then
403 0 : se_vert_num_threads = 1
404 : end if
405 1536 : if (se_tracer_num_threads < 1) then
406 0 : se_tracer_num_threads = 1
407 : end if
408 1536 : horz_num_threads = se_horz_num_threads
409 1536 : vert_num_threads = se_vert_num_threads
410 1536 : tracer_num_threads = se_tracer_num_threads
411 :
412 1536 : write_restart_unstruct = se_write_restart_unstruct
413 :
414 1536 : if (se_kmin_jet<0 ) kmin_jet = 1
415 1536 : if (se_kmax_jet<0 ) kmax_jet = nlev
416 :
417 1536 : if (masterproc) then
418 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_ftype = ',ftype
419 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_statediag_numtrac = ',statediag_numtrac
420 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_hypervis_subcycle = ',se_hypervis_subcycle
421 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_hypervis_subcycle_sponge = ',se_hypervis_subcycle_sponge
422 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_hypervis_subcycle_q = ',se_hypervis_subcycle_q
423 2 : write(iulog, '(a,l4)') 'dyn_readnl: se_large_Courant_incr = ',se_large_Courant_incr
424 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_limiter_option = ',se_limiter_option
425 2 : if (.not. se_refined_mesh) then
426 2 : write(iulog, '(a,i0)')'dyn_readnl: se_ne = ',se_ne
427 : end if
428 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_npes = ',se_npes
429 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_nsplit = ',se_nsplit
430 : !
431 : ! se_nu<0 then coefficients are set automatically in module global_norms_mod
432 : !
433 2 : if (se_nu_div>0) &
434 0 : write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu = ',se_nu
435 2 : if (se_nu_div>0) &
436 0 : write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_div = ',se_nu_div
437 2 : if (se_nu_p>0) then
438 0 : write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_p = ',se_nu_p
439 0 : write(iulog, '(a)') 'Note that nu_q must be the same as nu_p for mass / tracer inconsistency'
440 : end if
441 2 : write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_top = ',se_nu_top
442 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_qsplit = ',se_qsplit
443 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_rsplit = ',se_rsplit
444 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_statefreq = ',se_statefreq
445 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_tstep_type = ',se_tstep_type
446 2 : write(iulog, '(a,a)') 'dyn_readnl: se_vert_remap_T = ',trim(se_vert_remap_T)
447 2 : write(iulog, '(a,a)') 'dyn_readnl: se_vert_remap_uvTq_alg = ',trim(se_vert_remap_uvTq_alg)
448 2 : write(iulog, '(a,a)') 'dyn_readnl: se_vert_remap_tracer_alg = ',trim(se_vert_remap_tracer_alg)
449 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_fvm_supercycling = ',fvm_supercycling
450 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_fvm_supercycling_jet = ',fvm_supercycling_jet
451 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_kmin_jet = ',kmin_jet
452 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_kmax_jet = ',kmax_jet
453 :
454 2 : write(iulog, *) 'dyn_readnl: se_sponge_del4_nu_fac = ',se_sponge_del4_nu_fac
455 2 : if (se_sponge_del4_nu_fac < 0) write(iulog, '(a)') ' (automatically set based on model top location)'
456 2 : write(iulog, *) 'dyn_readnl: se_sponge_del4_nu_div_fac = ',se_sponge_del4_nu_div_fac
457 2 : if (se_sponge_del4_nu_div_fac < 0) write(iulog, '(a)') ' (automatically set based on model top location)'
458 2 : write(iulog, *) 'dyn_readnl: se_sponge_del4_lev = ',se_sponge_del4_lev
459 2 : if (se_sponge_del4_lev < 0) write(iulog, '(a)') ' (automatically set based on model top location)'
460 :
461 2 : if (se_refined_mesh) then
462 0 : write(iulog, '(a)') 'dyn_readnl: Refined mesh simulation'
463 0 : write(iulog, '(a)') 'dyn_readnl: se_mesh_file = ',trim(se_mesh_file)
464 0 : if (hypervis_power /= 0) then
465 0 : write(iulog, '(a)') 'Using scalar viscosity (Zarzycki et al 2014 JClim)'
466 0 : write(iulog, '(a,e11.4)') 'dyn_readnl: se_hypervis_power = ',se_hypervis_power, ', (tensor hyperviscosity)'
467 0 : write(iulog, '(a,e11.4)') 'dyn_readnl: se_max_hypervis_courant = ',se_max_hypervis_courant
468 : end if
469 0 : if (hypervis_scaling /= 0) then
470 0 : write(iulog, '(a)') 'Using tensor viscosity (Guba et al., 2014)'
471 0 : write(iulog, '(a,e11.4)') 'dyn_readnl: se_hypervis_scaling = ',se_hypervis_scaling
472 : end if
473 : end if
474 :
475 2 : if (use_cslam) then
476 2 : write(iulog, '(a)') 'dyn_readnl: physics will run on FVM points; advection by CSLAM'
477 2 : write(iulog,'(a,i0)') 'dyn_readnl: se_fv_nphys = ', fv_nphys
478 : else
479 0 : write(iulog, '(a)') 'dyn_readnl: physics will run on SE GLL points'
480 : end if
481 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_horz_num_threads = ',horz_num_threads
482 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_vert_num_threads = ',vert_num_threads
483 2 : write(iulog, '(a,i0)') 'dyn_readnl: se_tracer_num_threads = ',tracer_num_threads
484 2 : if (trim(se_write_grid_file) == 'SCRIP') then
485 0 : write(iulog,'(2a)') "dyn_readnl: write SCRIP grid file = ", trim(se_grid_filename)
486 : else
487 2 : write(iulog,'(a)') "dyn_readnl: do not write grid file"
488 : end if
489 2 : write(iulog,'(a,l1)') 'dyn_readnl: write gll corners to SEMapping.nc = ', &
490 4 : se_write_gll_corners
491 2 : write(iulog,'(a,l1)') 'dyn_readnl: write restart data on unstructured grid = ', &
492 4 : se_write_restart_unstruct
493 :
494 2 : write(iulog, '(a,e9.2)') 'dyn_readnl: se_molecular_diff = ', molecular_diff
495 : end if
496 :
497 1536 : call native_mapping_readnl(NLFileName)
498 :
499 : !---------------------------------------------------------------------------
500 : contains
501 : !---------------------------------------------------------------------------
502 :
503 3072 : integer function set_vert_remap( remap_T, remap_alg )
504 :
505 : ! Convert namelist input strings to the internally used integers.
506 :
507 : character(len=*), intent(in) :: remap_T ! scheme for remapping temperature
508 : character(len=*), intent(in) :: remap_alg ! remapping algorithm
509 :
510 : ! check valid remap_T values:
511 3072 : if (remap_T /= 'thermal_energy_over_P' .and. remap_T /= 'Tv_over_logP') then
512 0 : write(iulog,*)'set_vert_remap: invalid remap_T= ',trim(remap_T)
513 0 : call endrun('set_vert_remap: invalid remap_T')
514 : end if
515 :
516 : select case (remap_alg)
517 : case ('PPM_bc_mirror')
518 0 : set_vert_remap = 1
519 : case ('PPM_bc_PCoM')
520 0 : set_vert_remap = 2
521 : case ('PPM_bc_linear_extrapolation')
522 1536 : set_vert_remap = 10
523 : case ('FV3_PPM')
524 0 : if (remap_T == 'thermal_energy_over_P') then
525 : set_vert_remap = -4
526 : else
527 0 : set_vert_remap = -40
528 : end if
529 : case ('FV3_CS')
530 1536 : if (remap_T == 'thermal_energy_over_P') then
531 : set_vert_remap = -9
532 : else
533 0 : set_vert_remap = -90
534 : end if
535 : case ('FV3_CS_2dz_filter')
536 0 : if (remap_T == 'thermal_energy_over_P') then
537 : set_vert_remap = -10
538 : else
539 0 : set_vert_remap = -100
540 : end if
541 : case ('FV3_non_monotone_CS_2dz_filter')
542 0 : if (remap_T == 'thermal_energy_over_P') then
543 : set_vert_remap = -11
544 : else
545 0 : set_vert_remap = -110
546 : end if
547 : case default
548 0 : write(iulog,*)'set_vert_remap: invalid remap_alg= ',trim(remap_alg)
549 3072 : call endrun('set_vert_remap: invalid remap_alg')
550 : end select
551 :
552 1536 : end function set_vert_remap
553 :
554 : end subroutine dyn_readnl
555 :
556 : !=========================================================================================
557 :
558 1536 : subroutine dyn_register()
559 :
560 : use physics_buffer, only: pbuf_add_field, dtype_r8
561 : use ppgrid, only: pcols, pver
562 :
563 : ! These fields are computed by the dycore and passed to the physics via the
564 : ! physics buffer.
565 :
566 1536 : if (use_gw_front .or. use_gw_front_igw) then
567 : call pbuf_add_field("FRONTGF", "global", dtype_r8, (/pcols,pver/), &
568 1536 : frontgf_idx)
569 : call pbuf_add_field("FRONTGA", "global", dtype_r8, (/pcols,pver/), &
570 1536 : frontga_idx)
571 : end if
572 :
573 1536 : end subroutine dyn_register
574 :
575 : !=========================================================================================
576 :
577 3072 : subroutine dyn_init(dyn_in, dyn_out)
578 1536 : use prim_advance_mod, only: prim_advance_init
579 : use dyn_grid, only: elem, fvm
580 : use cam_pio_utils, only: clean_iodesc_list
581 : use physconst, only: cpair, pstd
582 : use air_composition, only: thermodynamic_active_species_num, thermodynamic_active_species_idx
583 : use air_composition, only: thermodynamic_active_species_idx_dycore
584 : use air_composition, only: thermodynamic_active_species_liq_idx,thermodynamic_active_species_ice_idx
585 : use air_composition, only: thermodynamic_active_species_liq_idx_dycore,thermodynamic_active_species_ice_idx_dycore
586 : use air_composition, only: thermodynamic_active_species_liq_num, thermodynamic_active_species_ice_num
587 : use cam_history, only: addfld, add_default, horiz_only, register_vector_field
588 : use gravity_waves_sources, only: gws_init
589 :
590 : use thread_mod, only: horz_num_threads
591 : use hybrid_mod, only: get_loop_ranges, config_thread_region
592 : use dimensions_mod, only: nu_scale_top
593 : use dimensions_mod, only: ksponge_end, kmvis_ref, kmcnd_ref,rho_ref,km_sponge_factor
594 : use dimensions_mod, only: cnst_name_gll, cnst_longname_gll
595 : use dimensions_mod, only: irecons_tracer_lev, irecons_tracer, kord_tr, kord_tr_cslam
596 : use prim_driver_mod, only: prim_init2
597 : use control_mod, only: molecular_diff, nu_top
598 : use test_fvm_mapping, only: test_mapping_addfld
599 : use phys_control, only: phys_getopts
600 : use cam_thermo, only: get_molecular_diff_coef_reference
601 : use control_mod, only: vert_remap_uvTq_alg, vert_remap_tracer_alg
602 : use std_atm_profile, only: std_atm_height
603 : use dyn_tests_utils, only: vc_dycore, vc_dry_pressure, string_vc, vc_str_lgth
604 : use cam_budget, only: cam_budget_em_snapshot, cam_budget_em_register, thermo_budget_history
605 :
606 : ! Dummy arguments:
607 : type(dyn_import_t), intent(out) :: dyn_in
608 : type(dyn_export_t), intent(out) :: dyn_out
609 :
610 : ! Local variables
611 : integer :: nets, nete, ie, k, kmol_end, mfound
612 : real(r8), parameter :: Tinit = 300.0_r8
613 : real(r8) :: press(1), ptop, tref,z(1)
614 :
615 : type(hybrid_t) :: hybrid
616 :
617 : integer :: ixcldice, ixcldliq
618 : integer :: m_cnst, m
619 :
620 : ! variables for initializing energy and axial angular momentum diagnostics
621 : integer, parameter :: num_stages = 14
622 : character (len = 4), dimension(num_stages) :: stage = (/"dED","dAF","dBD","dBL","dAL","dAD","dAR","dBF","dBH","dCH","dAH","dBS","dAS","p2d"/)
623 : character (len = 70),dimension(num_stages) :: stage_txt = (/&
624 : " end of previous dynamics ",& !dED
625 : " from previous remapping or state passed to dynamics",& !dAF - state in beginning of nsplit loop
626 : " state after applying CAM forcing ",& !dBD - state after applyCAMforcing
627 : " before floating dynamics ",& !dBL
628 : " after floating dynamics ",& !dAL
629 : " before vertical remapping ",& !dAD - state before vertical remapping
630 : " after vertical remapping ",& !dAR - state at end of nsplit loop
631 : " state passed to parameterizations ",& !dBF
632 : " state before hypervis ",& !dBH
633 : " state after hypervis but before adding heating term",& !dCH
634 : " state after hypervis ",& !dAH
635 : " state before sponge layer diffusion ",& !dBS - state before sponge del2
636 : " state after sponge layer diffusion ",& !dAS - state after sponge del2
637 : " phys2dyn mapping errors (requires ftype-1) " & !p2d - for assessing phys2dyn mapping errors
638 : /)
639 :
640 : integer :: istage
641 : character (len=vc_str_lgth) :: vc_str
642 :
643 : logical :: history_budget ! output tendencies and state variables for budgets
644 : integer :: budget_hfile_num
645 :
646 : character(len=*), parameter :: sub = 'dyn_init'
647 :
648 : real(r8) :: km_sponge_factor_local(nlev+1)
649 : !----------------------------------------------------------------------------
650 1536 : vc_dycore = vc_dry_pressure
651 1536 : if (masterproc) then
652 2 : call string_vc(vc_dycore,vc_str)
653 2 : write(iulog,*) sub//': vertical coordinate dycore : ',trim(vc_str)
654 : end if
655 : ! Now allocate and set condenstate vars
656 4608 : allocate(cnst_name_gll(qsize)) ! constituent names for gll tracers
657 4608 : allocate(cnst_longname_gll(qsize)) ! long name of constituents for gll tracers
658 :
659 4608 : allocate(kord_tr(qsize))
660 10752 : kord_tr(:) = vert_remap_tracer_alg
661 1536 : if (use_cslam) then
662 4608 : allocate(kord_tr_cslam(ntrac))
663 64512 : kord_tr_cslam(:) = vert_remap_tracer_alg
664 : end if
665 :
666 :
667 10752 : do m=1,qsize
668 : !
669 : ! The "_gll" index variables below are used to keep track of condensate-loading tracers
670 : ! since they are not necessarily indexed contiguously and not necessarily in the same
671 : ! order (physics is in charge of the order)
672 : !
673 : ! if running with CSLAM then the SE (gll) condensate-loading water tracers are always
674 : ! indexed contiguously (q,cldliq,cldice,rain,snow,graupel) - see above
675 : !
676 : ! CSLAM tracers are always indexed as in physics
677 : ! of no CSLAM then SE tracers are always indexed as in physics
678 : !
679 10752 : if (use_cslam) then
680 : !
681 : ! note that in this case qsize = thermodynamic_active_species_num
682 : !
683 9216 : thermodynamic_active_species_idx_dycore(m) = m
684 9216 : kord_tr_cslam(thermodynamic_active_species_idx(m)) = vert_remap_uvTq_alg
685 9216 : kord_tr(m) = vert_remap_uvTq_alg
686 9216 : cnst_name_gll (m) = cnst_name (thermodynamic_active_species_idx(m))
687 9216 : cnst_longname_gll(m) = cnst_longname(thermodynamic_active_species_idx(m))
688 : else
689 : !
690 : ! if not running with CSLAM then the condensate-loading water tracers are not necessarily
691 : ! indexed contiguously (are indexed as in physics)
692 : !
693 0 : if (m.le.thermodynamic_active_species_num) then
694 0 : thermodynamic_active_species_idx_dycore(m) = thermodynamic_active_species_idx(m)
695 0 : kord_tr(thermodynamic_active_species_idx_dycore(m)) = vert_remap_uvTq_alg
696 : end if
697 0 : cnst_name_gll (m) = cnst_name (m)
698 0 : cnst_longname_gll(m) = cnst_longname(m)
699 : end if
700 : end do
701 :
702 4608 : do m=1,thermodynamic_active_species_liq_num
703 3072 : if (use_cslam) then
704 21504 : do mfound=1,qsize
705 21504 : if (TRIM(cnst_name(thermodynamic_active_species_liq_idx(m)))==TRIM(cnst_name_gll(mfound))) then
706 3072 : thermodynamic_active_species_liq_idx_dycore(m) = mfound
707 : end if
708 : end do
709 : else
710 0 : thermodynamic_active_species_liq_idx_dycore(m) = thermodynamic_active_species_liq_idx(m)
711 : end if
712 4608 : if (masterproc) then
713 4 : write(iulog,*) sub//": m,thermodynamic_active_species_idx_liq_dycore: ",m,thermodynamic_active_species_liq_idx_dycore(m)
714 : end if
715 : end do
716 6144 : do m=1,thermodynamic_active_species_ice_num
717 4608 : if (use_cslam) then
718 32256 : do mfound=1,qsize
719 32256 : if (TRIM(cnst_name(thermodynamic_active_species_ice_idx(m)))==TRIM(cnst_name_gll(mfound))) then
720 4608 : thermodynamic_active_species_ice_idx_dycore(m) = mfound
721 : end if
722 : end do
723 : else
724 0 : thermodynamic_active_species_ice_idx_dycore(m) = thermodynamic_active_species_ice_idx(m)
725 : end if
726 6144 : if (masterproc) then
727 6 : write(iulog,*) sub//": m,thermodynamic_active_species_idx_ice_dycore: ",m,thermodynamic_active_species_ice_idx_dycore(m)
728 : end if
729 : end do
730 :
731 : !
732 : ! Initialize the import/export objects
733 : !
734 1536 : if(iam < par%nprocs) then
735 1536 : dyn_in%elem => elem
736 1536 : dyn_in%fvm => fvm
737 :
738 1536 : dyn_out%elem => elem
739 1536 : dyn_out%fvm => fvm
740 : else
741 0 : nullify(dyn_in%elem)
742 0 : nullify(dyn_in%fvm)
743 0 : nullify(dyn_out%elem)
744 0 : nullify(dyn_out%fvm)
745 : end if
746 :
747 1536 : call set_phis(dyn_in)
748 :
749 1536 : if (initial_run) then
750 768 : call read_inidat(dyn_in)
751 768 : call clean_iodesc_list()
752 : end if
753 : !
754 : ! initialize diffusion in dycore
755 : !
756 1536 : kmol_end = 0
757 1536 : if (molecular_diff>0) then
758 : !
759 : ! molecular diffusion and thermal conductivity reference values
760 : !
761 0 : if (masterproc) write(iulog,*) sub//": initialize molecular diffusion reference profiles"
762 0 : tref = 1000._r8 !mean value at model top for solar max
763 0 : km_sponge_factor = molecular_diff
764 0 : km_sponge_factor_local = molecular_diff
765 : !
766 : ! get rho, kmvis and kmcnd at mid-levels
767 : !
768 : call get_molecular_diff_coef_reference(tref,&
769 : (hvcoord%hyam(:)+hvcoord%hybm(:))*hvcoord%ps0,km_sponge_factor,&
770 0 : kmvis_ref,kmcnd_ref,rho_ref)
771 :
772 0 : if (masterproc) then
773 0 : write(iulog,*) "Molecular viscosity and thermal conductivity reference profile"
774 0 : write(iulog,*) "k, p, z, km_sponge_factor, kmvis_ref/rho_ref, kmcnd_ref/(cp*rho_ref):"
775 : end if
776 0 : do k=1,nlev
777 : ! only apply molecular viscosity where viscosity is > 1000 m/s^2
778 0 : if (MIN(kmvis_ref(k)/rho_ref(k),kmcnd_ref(k)/(cpair*rho_ref(k)))>1000.0_r8) then
779 0 : if (masterproc) then
780 0 : press = hvcoord%hyam(k)*hvcoord%ps0+hvcoord%hybm(k)*pstd
781 0 : call std_atm_height(press,z)
782 0 : write(iulog,'(i3,5e11.4)') k,press, z,km_sponge_factor(k),kmvis_ref(k)/rho_ref(k),kmcnd_ref(k)/(cpair*rho_ref(k))
783 : end if
784 0 : kmol_end = k
785 : else
786 0 : kmvis_ref(k) = 1.0_r8
787 0 : kmcnd_ref(k) = 1.0_r8
788 : end if
789 : end do
790 : else
791 : ! -1.0E6 is an arbitrary unrealistic value. But it is used in the calculation
792 : ! of a diagnostic quantity in global_norms_mod so can't be set to huge or nan.
793 144384 : kmvis_ref(:) = -1.0E6_r8
794 144384 : kmcnd_ref(:) = -1.0E6_r8
795 144384 : rho_ref(:) = -1.0E6_r8
796 : end if
797 : !
798 144384 : irecons_tracer_lev(:) = irecons_tracer !use high-order CSLAM in all layers
799 : !
800 : ! compute scaling of traditional sponge layer damping (following cd_core.F90 in CAM-FV)
801 : !
802 1536 : nu_scale_top(:) = 0.0_r8
803 1536 : if (nu_top>0) then
804 1536 : ptop = hvcoord%hyai(1)*hvcoord%ps0
805 1536 : if (ptop>300.0_r8) then
806 : !
807 : ! for low tops the tanh formulae below makes the sponge excessively deep
808 : !
809 0 : nu_scale_top(1) = 4.0_r8
810 0 : nu_scale_top(2) = 2.0_r8
811 0 : nu_scale_top(3) = 1.0_r8
812 0 : ksponge_end = 3
813 1536 : else if (ptop>100.0_r8) then
814 : !
815 : ! CAM6 top (~225 Pa) or CAM7 low top
816 : !
817 : ! For backwards compatibility numbers below match tanh profile
818 : ! used in FV
819 : !
820 0 : nu_scale_top(1) = 4.4_r8
821 0 : nu_scale_top(2) = 1.3_r8
822 0 : nu_scale_top(3) = 3.9_r8
823 0 : ksponge_end = 3
824 1536 : else if (ptop>1e-1_r8) then
825 : !
826 : ! CAM7 FMT
827 : !
828 1536 : nu_scale_top(1) = 3.0_r8
829 1536 : nu_scale_top(2) = 1.0_r8
830 1536 : nu_scale_top(3) = 0.1_r8
831 1536 : nu_scale_top(4) = 0.05_r8
832 1536 : ksponge_end = 4
833 0 : else if (ptop>1e-4_r8) then
834 : !
835 : ! WACCM and WACCM-x
836 : !
837 0 : nu_scale_top(1) = 5.0_r8
838 0 : nu_scale_top(2) = 5.0_r8
839 0 : nu_scale_top(3) = 5.0_r8
840 0 : nu_scale_top(4) = 2.0_r8
841 0 : nu_scale_top(5) = 1.0_r8
842 0 : nu_scale_top(6) = 0.1_r8
843 0 : ksponge_end = 6
844 : end if
845 : else
846 0 : ksponge_end = 0
847 : end if
848 1536 : ksponge_end = MAX(MAX(ksponge_end,1),kmol_end)
849 1536 : if (masterproc) then
850 2 : write(iulog,*) sub//": ksponge_end = ",ksponge_end
851 2 : write(iulog,*) sub//": sponge layer Laplacian damping"
852 2 : write(iulog,*) "k, p, z, nu_scale_top, nu (actual Laplacian damping coefficient)"
853 2 : if (nu_top>0) then
854 12 : do k=1,ksponge_end+1
855 20 : press = (hvcoord%hyam(k)+hvcoord%hybm(k))*hvcoord%ps0
856 10 : call std_atm_height(press,z)
857 20 : write(iulog,'(i3,4e11.4)') k,press,z,&
858 32 : nu_scale_top(k),nu_scale_top(k)*nu_top
859 : end do
860 : end if
861 : end if
862 :
863 1536 : if (iam < par%nprocs) then
864 1536 : call prim_advance_init(par,elem)
865 : !$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete)
866 1536 : hybrid = config_thread_region(par,'horizontal')
867 1536 : call get_loop_ranges(hybrid, ibeg=nets, iend=nete)
868 1536 : call prim_init2(elem, fvm, hybrid, nets, nete, TimeLevel, hvcoord)
869 : !$OMP END PARALLEL
870 :
871 1536 : if (use_gw_front .or. use_gw_front_igw) call gws_init(elem)
872 : end if ! iam < par%nprocs
873 :
874 3072 : call addfld ('nu_kmvis', (/ 'lev' /), 'A', '', 'Molecular viscosity Laplacian coefficient' , gridname='GLL')
875 3072 : call addfld ('nu_kmcnd', (/ 'lev' /), 'A', '', 'Thermal conductivity Laplacian coefficient' , gridname='GLL')
876 3072 : call addfld ('nu_kmcnd_dp',(/ 'lev' /), 'A', '', 'Thermal conductivity like Laplacian coefficient on dp', gridname='GLL')
877 :
878 :
879 : ! Forcing from physics on the GLL grid
880 3072 : call addfld ('FU', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind forcing term on GLL grid', gridname='GLL')
881 3072 : call addfld ('FV', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind forcing term on GLL grid',gridname='GLL')
882 1536 : call register_vector_field('FU', 'FV')
883 3072 : call addfld ('FT', (/ 'lev' /), 'A', 'K/s', 'Temperature forcing term on GLL grid',gridname='GLL')
884 :
885 : ! Tracer forcing on fvm (CSLAM) grid and internal CSLAM pressure fields
886 1536 : if (use_cslam) then
887 64512 : do m = 1, ntrac
888 62976 : call addfld (trim(cnst_name(m))//'_fvm', (/ 'lev' /), 'I', 'kg/kg', &
889 188928 : trim(cnst_longname(m)), gridname='FVM')
890 :
891 62976 : call addfld ('F'//trim(cnst_name(m))//'_fvm', (/ 'lev' /), 'I', 'kg/kg/s', &
892 62976 : trim(cnst_longname(m))//' mixing ratio forcing term (q_new-q_old) on fvm grid', &
893 253440 : gridname='FVM')
894 : end do
895 :
896 3072 : call addfld ('dp_fvm' ,(/ 'lev' /), 'I', 'Pa','CSLAM Pressure level thickness', gridname='FVM')
897 1536 : call addfld ('PSDRY_fvm',horiz_only, 'I','Pa','CSLAM dry surface pressure' , gridname='FVM')
898 : end if
899 :
900 10752 : do m_cnst = 1, qsize
901 0 : call addfld ('F'//trim(cnst_name_gll(m_cnst))//'_gll', (/ 'lev' /), 'I', 'kg/kg/s', &
902 19968 : trim(cnst_longname_gll(m_cnst))//' mixing ratio forcing term (q_new-q_old) on GLL grid', gridname='GLL')
903 : end do
904 :
905 : ! Energy diagnostics and axial angular momentum diagnostics
906 1536 : call addfld ('ABS_dPSdt', horiz_only, 'A', 'Pa/s', 'Absolute surface pressure tendency',gridname='GLL')
907 :
908 1536 : if (use_cslam) then
909 : #ifdef waccm_debug
910 : call addfld ('CSLAM_gamma', (/ 'lev' /), 'A', '', 'Courant number from CSLAM', gridname='FVM')
911 : #endif
912 1536 : call addfld ('WV_PDC', horiz_only, 'A', 'kg/m2','Total column water vapor lost in physics-dynamics coupling',gridname='FVM')
913 1536 : call addfld ('WL_PDC', horiz_only, 'A', 'kg/m2','Total column cloud water lost in physics-dynamics coupling',gridname='FVM')
914 1536 : call addfld ('WI_PDC', horiz_only, 'A', 'kg/m2','Total column cloud ice lost in physics-dynamics coupling' ,gridname='FVM')
915 1536 : call addfld ('TT_PDC', horiz_only, 'A', 'kg/m2','Total column test tracer lost in physics-dynamics coupling',gridname='FVM')
916 : else
917 0 : call addfld ('WV_PDC', horiz_only, 'A', 'kg/m2','Total column water vapor lost in physics-dynamics coupling',gridname='GLL')
918 0 : call addfld ('WL_PDC', horiz_only, 'A', 'kg/m2','Total column cloud water lost in physics-dynamics coupling',gridname='GLL')
919 0 : call addfld ('WI_PDC', horiz_only, 'A', 'kg/m2','Total column cloud ice lost in physics-dynamics coupling' ,gridname='GLL')
920 0 : call addfld ('TT_PDC', horiz_only, 'A', 'kg/m2','Total column test tracer lost in physics-dynamics coupling',gridname='GLL')
921 : end if
922 :
923 1536 : if (thermo_budget_history) then
924 : ! Register stages for budgets
925 0 : do istage = 1, num_stages
926 0 : call cam_budget_em_snapshot(TRIM(ADJUSTL(stage(istage))), 'dyn', &
927 0 : longname=TRIM(ADJUSTL(stage_txt(istage))))
928 : end do
929 : !
930 : ! Register tendency (difference) budgets
931 : !
932 : call cam_budget_em_register('dEdt_floating_dyn' ,'dAL','dBL','dyn','dif', &
933 0 : longname="dE/dt floating dynamics (dAL-dBL)" )
934 : call cam_budget_em_register('dEdt_vert_remap' ,'dAR','dAD','dyn','dif', &
935 0 : longname="dE/dt vertical remapping (dAR-dAD)" )
936 : call cam_budget_em_register('dEdt_phys_tot_in_dyn','dBD','dAF','dyn','dif', &
937 0 : longname="dE/dt physics tendency in dynamics (dBD-dAF)" )
938 : call cam_budget_em_register('dEdt_del4' ,'dCH','dBH','dyn','dif', &
939 0 : longname="dE/dt del4 (dCH-dBH)" )
940 : call cam_budget_em_register('dEdt_del4_fric_heat','dAH','dCH','dyn','dif', &
941 0 : longname="dE/dt del4 frictional heating (dAH-dCH)" )
942 : call cam_budget_em_register('dEdt_del4_tot' ,'dAH','dBH','dyn','dif', &
943 0 : longname="dE/dt del4 + del4 frictional heating (dAH-dBH)" )
944 : call cam_budget_em_register('dEdt_del2_sponge' ,'dAS','dBS','dyn','dif', &
945 0 : longname="dE/dt del2 sponge (dAS-dBS)" )
946 : !
947 : ! Register derived budgets
948 : !
949 : call cam_budget_em_register('dEdt_dycore' ,'dEdt_floating_dyn','dEdt_vert_remap' ,'dyn','sum', &
950 0 : longname="dE/dt adiabatic dynamics" )
951 : call cam_budget_em_register('dEdt_del2_del4_tot' ,'dEdt_del4_tot' ,'dEdt_del2_sponge' ,'dyn','sum', &
952 0 : longname="dE/dt explicit diffusion total" )
953 : call cam_budget_em_register('dEdt_residual' ,'dEdt_floating_dyn','dEdt_del2_del4_tot','dyn','dif',&
954 0 : longname="dE/dt residual (dEdt_floating_dyn-dEdt_del2_del4_tot)" )
955 : end if
956 : !
957 : ! add dynamical core tracer tendency output
958 : !
959 1536 : if (use_cslam) then
960 64512 : do m = 1, pcnst
961 125952 : call addfld(tottnam(m),(/ 'lev' /),'A','kg/kg/s',trim(cnst_name(m))//' horz + vert', &
962 253440 : gridname='FVM')
963 : end do
964 : else
965 0 : do m = 1, pcnst
966 0 : call addfld(tottnam(m),(/ 'lev' /),'A','kg/kg/s',trim(cnst_name(m))//' horz + vert', &
967 0 : gridname='GLL')
968 : end do
969 : end if
970 1536 : call phys_getopts(history_budget_out=history_budget, history_budget_histfile_num_out=budget_hfile_num)
971 1536 : if ( history_budget ) then
972 0 : call cnst_get_ind('CLDLIQ', ixcldliq)
973 0 : call cnst_get_ind('CLDICE', ixcldice)
974 0 : call add_default(tottnam( 1), budget_hfile_num, ' ')
975 0 : call add_default(tottnam(ixcldliq), budget_hfile_num, ' ')
976 0 : call add_default(tottnam(ixcldice), budget_hfile_num, ' ')
977 : end if
978 :
979 1536 : call test_mapping_addfld
980 3072 : end subroutine dyn_init
981 :
982 : !=========================================================================================
983 :
984 14592 : subroutine dyn_run(dyn_state)
985 1536 : use air_composition, only: thermodynamic_active_species_num, dry_air_species_num
986 : use air_composition, only: thermodynamic_active_species_idx_dycore
987 : use prim_driver_mod, only: prim_run_subcycle
988 : use dimensions_mod, only: cnst_name_gll
989 : use se_dyn_time_mod, only: tstep, nsplit, timelevel_qdp, tevolve
990 : use hybrid_mod, only: config_thread_region, get_loop_ranges
991 : use control_mod, only: qsplit, rsplit, ftype_conserve
992 : use thread_mod, only: horz_num_threads
993 :
994 : type(dyn_export_t), intent(inout) :: dyn_state
995 :
996 : type(hybrid_t) :: hybrid
997 : integer :: tl_f
998 : integer :: n
999 : integer :: nets, nete
1000 : integer :: i, ie, j, k, m, nq, m_cnst
1001 : integer :: n0_qdp, nsplit_local
1002 : logical :: ldiag
1003 :
1004 : real(r8) :: ftmp(npsq,nlev,3)
1005 : real(r8) :: dtime
1006 : real(r8) :: rec2dt, pdel
1007 :
1008 14592 : real(r8), allocatable, dimension(:,:,:) :: ps_before
1009 14592 : real(r8), allocatable, dimension(:,:,:) :: abs_ps_tend
1010 29184 : real (kind=r8) :: omega_cn(2,nelemd) !min and max of vertical Courant number
1011 : !----------------------------------------------------------------------------
1012 :
1013 : #ifdef debug_coupling
1014 : return
1015 : #endif
1016 :
1017 14592 : nsplit_local = nsplit
1018 14592 : tevolve = 0._r8
1019 :
1020 14592 : if (iam >= par%nprocs) return
1021 :
1022 14592 : ldiag = hist_fld_active('ABS_dPSdt')
1023 14592 : if (ldiag) then
1024 0 : allocate(ps_before(np,np,nelemd))
1025 0 : allocate(abs_ps_tend(np,np,nelemd))
1026 :
1027 : end if
1028 :
1029 : !$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n,ie,m,i,j,k,ftmp)
1030 14592 : hybrid = config_thread_region(par,'horizontal')
1031 14592 : call get_loop_ranges(hybrid, ibeg=nets, iend=nete)
1032 :
1033 14592 : dtime = get_step_size()
1034 14592 : rec2dt = 1._r8/dtime
1035 :
1036 14592 : tl_f = TimeLevel%n0 ! timelevel which was adjusted by physics
1037 14592 : call TimeLevel_Qdp(TimeLevel, qsplit, n0_qdp)!get n0_qdp for diagnostics call
1038 :
1039 : ! output physics forcing
1040 14592 : if (hist_fld_active('FU') .or. hist_fld_active('FV') .or.hist_fld_active('FT')) then
1041 0 : do ie = nets, nete
1042 0 : do k = 1, nlev
1043 0 : do j = 1, np
1044 0 : do i = 1, np
1045 0 : ftmp(i+(j-1)*np,k,1) = dyn_state%elem(ie)%derived%FM(i,j,1,k)
1046 0 : ftmp(i+(j-1)*np,k,2) = dyn_state%elem(ie)%derived%FM(i,j,2,k)
1047 0 : ftmp(i+(j-1)*np,k,3) = dyn_state%elem(ie)%derived%FT(i,j,k)
1048 : end do
1049 : end do
1050 : end do
1051 :
1052 0 : call outfld('FU', ftmp(:,:,1), npsq, ie)
1053 0 : call outfld('FV', ftmp(:,:,2), npsq, ie)
1054 0 : call outfld('FT', ftmp(:,:,3), npsq, ie)
1055 : end do
1056 : end if
1057 :
1058 102144 : do m = 1, qsize
1059 102144 : if (hist_fld_active('F'//trim(cnst_name_gll(m))//'_gll')) then
1060 0 : do ie = nets, nete
1061 0 : call outfld('F'//trim(cnst_name_gll(m))//'_gll',&
1062 0 : RESHAPE(dyn_state%elem(ie)%derived%FQ(:,:,:,m), (/np*np,nlev/)), npsq, ie)
1063 : end do
1064 : end if
1065 : end do
1066 :
1067 : ! convert elem(ie)%derived%fq to mass tendency
1068 14592 : if (.not.use_cslam) then
1069 0 : do ie = nets, nete
1070 0 : do m = 1, qsize
1071 0 : do k = 1, nlev
1072 0 : do j = 1, np
1073 0 : do i = 1, np
1074 0 : dyn_state%elem(ie)%derived%FQ(i,j,k,m) = dyn_state%elem(ie)%derived%FQ(i,j,k,m)* &
1075 0 : rec2dt*dyn_state%elem(ie)%state%dp3d(i,j,k,tl_f)
1076 : end do
1077 : end do
1078 : end do
1079 : end do
1080 : end do
1081 : end if
1082 :
1083 14592 : if (ftype_conserve>0.and..not.use_cslam) then
1084 0 : do ie = nets, nete
1085 0 : do k=1,nlev
1086 0 : do j=1,np
1087 0 : do i = 1, np
1088 0 : pdel = dyn_state%elem(ie)%state%dp3d(i,j,k,tl_f)
1089 0 : do nq=dry_air_species_num+1,thermodynamic_active_species_num
1090 0 : m_cnst = thermodynamic_active_species_idx_dycore(nq)
1091 0 : pdel = pdel + (dyn_state%elem(ie)%state%qdp(i,j,k,m_cnst,n0_qdp)+dyn_state%elem(ie)%derived%FQ(i,j,k,m_cnst)*dtime)
1092 : end do
1093 0 : dyn_state%elem(ie)%derived%FDP(i,j,k) = pdel
1094 : end do
1095 : end do
1096 : end do
1097 : end do
1098 : end if
1099 :
1100 14592 : if (use_cslam) then
1101 117192 : do ie = nets, nete
1102 4323792 : do m = 1, ntrac
1103 395523000 : do k = 1, nlev
1104 1569061800 : do j = 1, nc
1105 5085779400 : do i = 1, nc
1106 0 : dyn_state%fvm(ie)%fc(i,j,k,m) = dyn_state%fvm(ie)%fc(i,j,k,m)* &
1107 4694565600 : rec2dt!*dyn_state%fvm(ie)%dp_fvm(i,j,k)
1108 : end do
1109 : end do
1110 : end do
1111 : end do
1112 : end do
1113 : end if
1114 :
1115 14592 : if (ldiag) then
1116 0 : abs_ps_tend(:,:,nets:nete) = 0.0_r8
1117 : endif
1118 :
1119 43776 : do n = 1, nsplit_local
1120 :
1121 29184 : if (ldiag) then
1122 0 : do ie = nets, nete
1123 0 : ps_before(:,:,ie) = dyn_state%elem(ie)%state%psdry(:,:)
1124 : end do
1125 : end if
1126 :
1127 : ! forward-in-time RK, with subcycling
1128 : call prim_run_subcycle(dyn_state%elem, dyn_state%fvm, hybrid, nets, nete, &
1129 29184 : tstep, TimeLevel, hvcoord, n, omega_cn)
1130 :
1131 43776 : if (ldiag) then
1132 0 : do ie = nets, nete
1133 0 : abs_ps_tend(:,:,ie) = abs_ps_tend(:,:,ie) + &
1134 0 : ABS(ps_before(:,:,ie)-dyn_state%elem(ie)%state%psdry(:,:)) &
1135 0 : /(tstep*qsplit*rsplit)
1136 : end do
1137 : end if
1138 :
1139 : end do
1140 :
1141 14592 : if (ldiag) then
1142 0 : do ie=nets,nete
1143 0 : abs_ps_tend(:,:,ie)=abs_ps_tend(:,:,ie)/DBLE(nsplit)
1144 0 : call outfld('ABS_dPSdt',RESHAPE(abs_ps_tend(:,:,ie),(/npsq/)),npsq,ie)
1145 : end do
1146 : end if
1147 :
1148 : !$OMP END PARALLEL
1149 :
1150 : if (ldiag) then
1151 0 : deallocate(ps_before,abs_ps_tend)
1152 : endif
1153 : ! output vars on CSLAM fvm grid
1154 14592 : call write_dyn_vars(dyn_state)
1155 :
1156 14592 : end subroutine dyn_run
1157 :
1158 : !===============================================================================
1159 :
1160 0 : subroutine dyn_final(DYN_STATE, RESTART_FILE)
1161 :
1162 : type (elem_state_t), target :: DYN_STATE
1163 : character(LEN=*) , intent(IN) :: RESTART_FILE
1164 :
1165 14592 : end subroutine dyn_final
1166 :
1167 : !===============================================================================
1168 :
1169 768 : subroutine read_inidat(dyn_in)
1170 : use air_composition, only: thermodynamic_active_species_num, dry_air_species_num
1171 : use shr_sys_mod, only: shr_sys_flush
1172 : use hycoef, only: hyai, hybi, ps0
1173 : use const_init, only: cnst_init_default
1174 :
1175 : use element_mod, only: timelevels
1176 : use fvm_mapping, only: dyn2fvm_mass_vars
1177 : use control_mod, only: runtype
1178 : use prim_driver_mod, only: prim_set_dry_mass
1179 : use air_composition, only: thermodynamic_active_species_idx
1180 : use cam_initfiles, only: scale_dry_air_mass
1181 : ! Arguments
1182 : type (dyn_import_t), target, intent(inout) :: dyn_in ! dynamics import
1183 :
1184 : ! Local variables
1185 :
1186 768 : integer(iMap), pointer :: ldof(:) ! Basic (2D) grid dof
1187 :
1188 : type(file_desc_t), pointer :: fh_ini, fh_topo
1189 :
1190 768 : type(element_t), pointer :: elem(:)
1191 :
1192 768 : real(r8), allocatable :: qtmp(:,:,:,:,:) ! (np,np,nlev,nelemd,n)
1193 768 : real(r8), allocatable :: dbuf2(:,:) ! (npsq,nelemd)
1194 768 : real(r8), allocatable :: dbuf3(:,:,:) ! (npsq,nlev,nelemd)
1195 768 : real(r8), allocatable :: phis_tmp(:,:) ! (npsp,nelemd)
1196 768 : real(r8), allocatable :: factor_array(:,:,:,:) ! (np,np,nlev,nelemd)
1197 768 : logical, allocatable :: pmask(:) ! (npsq*nelemd) unique grid vals
1198 :
1199 : character(len=max_hcoordname_len):: grid_name
1200 768 : real(r8), allocatable :: latvals(:)
1201 768 : real(r8), allocatable :: lonvals(:)
1202 768 : real(r8), pointer :: latvals_deg(:)
1203 768 : real(r8), pointer :: lonvals_deg(:)
1204 :
1205 : integer :: ie, k, t
1206 : character(len=max_fieldname_len) :: fieldname, fieldname2
1207 : logical :: found
1208 : logical :: inic_wet ! true if initial condition is based on
1209 : ! wet pressure and water species
1210 : integer :: kptr, m_cnst
1211 768 : type(EdgeBuffer_t) :: edge
1212 :
1213 : integer :: rndm_seed_sz
1214 768 : integer, allocatable :: rndm_seed(:)
1215 : integer :: dims(2)
1216 : integer :: pio_errtype
1217 : real(r8) :: pertval
1218 : integer :: i, j, indx, nq
1219 : integer :: dyn_cols
1220 : character(len=128) :: errmsg
1221 : character(len=*), parameter :: sub='READ_INIDAT'
1222 :
1223 : real(r8) :: dp_tmp, pstmp(np,np)
1224 :
1225 : ! Variables for analytic initial conditions
1226 768 : integer, allocatable :: glob_ind(:)
1227 768 : integer, allocatable :: m_ind(:)
1228 768 : real(r8), allocatable :: dbuf4(:,:,:,:)
1229 : !----------------------------------------------------------------------------
1230 :
1231 1536 : fh_ini => initial_file_get_id()
1232 768 : fh_topo => topo_file_get_id()
1233 :
1234 768 : if (iam < par%nprocs) then
1235 768 : elem => dyn_in%elem
1236 : else
1237 0 : nullify(elem)
1238 : end if
1239 :
1240 3072 : allocate(qtmp(np,np,nlev,nelemd,pcnst))
1241 432647856 : qtmp = 0._r8
1242 :
1243 : ! Set mask to indicate which columns are active
1244 768 : nullify(ldof)
1245 768 : call cam_grid_get_gcid(cam_grid_id(ini_grid_name), ldof)
1246 2304 : allocate(pmask(npsq*nelemd))
1247 87168 : pmask(:) = (ldof /= 0)
1248 :
1249 : ! lat/lon needed in radians
1250 768 : latvals_deg => cam_grid_get_latvals(cam_grid_id(ini_grid_name))
1251 768 : lonvals_deg => cam_grid_get_lonvals(cam_grid_id(ini_grid_name))
1252 2304 : allocate(latvals(np*np*nelemd))
1253 1536 : allocate(lonvals(np*np*nelemd))
1254 87168 : latvals(:) = latvals_deg(:)*deg2rad
1255 87168 : lonvals(:) = lonvals_deg(:)*deg2rad
1256 :
1257 : ! Set PIO to return error codes when reading data from IC file.
1258 768 : call pio_seterrorhandling(fh_ini, PIO_BCAST_ERROR, pio_errtype)
1259 :
1260 : ! Get the number of columns in the global GLL grid.
1261 768 : call cam_grid_dimensions(ini_grid_name, dims)
1262 768 : dyn_cols = dims(1)
1263 :
1264 : ! Set ICs. Either from analytic expressions or read from file.
1265 :
1266 768 : if (analytic_ic_active() .and. (iam < par%nprocs)) then
1267 :
1268 : ! PHIS has already been set by set_phis. Get local copy for
1269 : ! possible use in setting T and PS in the analytic IC code.
1270 0 : allocate(phis_tmp(npsq,nelemd))
1271 0 : do ie = 1, nelemd
1272 : k = 1
1273 0 : do j = 1, np
1274 0 : do i = 1, np
1275 0 : phis_tmp(k,ie) = elem(ie)%state%phis(i,j)
1276 0 : k = k + 1
1277 : end do
1278 : end do
1279 : end do
1280 :
1281 0 : inic_wet = .false.
1282 0 : allocate(glob_ind(npsq * nelemd))
1283 0 : j = 1
1284 0 : do ie = 1, nelemd
1285 0 : do i = 1, npsq
1286 : ! Create a global(ish) column index
1287 0 : glob_ind(j) = elem(ie)%GlobalId
1288 0 : j = j + 1
1289 : end do
1290 : end do
1291 :
1292 : ! First, initialize all the variables, then assign
1293 0 : allocate(dbuf4(npsq, nlev, nelemd, (qsize + 4)))
1294 0 : dbuf4 = 0.0_r8
1295 0 : allocate(m_ind(qsize))
1296 0 : do m_cnst = 1, qsize
1297 0 : m_ind(m_cnst) = m_cnst
1298 : end do
1299 :
1300 : ! Init tracers on the GLL grid. Note that analytic_ic_set_ic makes
1301 : ! use of cnst_init_default for the tracers except water vapor.
1302 :
1303 : call analytic_ic_set_ic(vcoord, latvals, lonvals, glob_ind, &
1304 : PS=dbuf4(:,1,:,(qsize+1)), U=dbuf4(:,:,:,(qsize+2)), &
1305 : V=dbuf4(:,:,:,(qsize+3)), T=dbuf4(:,:,:,(qsize+4)), &
1306 0 : Q=dbuf4(:,:,:,1:qsize), m_cnst=m_ind, mask=pmask(:), &
1307 0 : PHIS_IN=PHIS_tmp)
1308 0 : deallocate(m_ind)
1309 0 : deallocate(glob_ind)
1310 0 : deallocate(phis_tmp)
1311 0 : do ie = 1, nelemd
1312 : indx = 1
1313 0 : do j = 1, np
1314 0 : do i = 1, np
1315 : ! PS
1316 0 : elem(ie)%state%psdry(i,j) = dbuf4(indx, 1, ie, (qsize+1))
1317 : ! U
1318 0 : elem(ie)%state%v(i,j,1,:,1) = dbuf4(indx, :, ie, (qsize+2))
1319 : ! V
1320 0 : elem(ie)%state%v(i,j,2,:,1) = dbuf4(indx, :, ie, (qsize+3))
1321 : ! T
1322 0 : elem(ie)%state%T(i,j,:,1) = dbuf4(indx, :, ie, (qsize+4))
1323 0 : indx = indx + 1
1324 : end do
1325 : end do
1326 : end do
1327 :
1328 : ! Tracers to be advected on GLL grid.
1329 : ! Note that fvm tracers are initialized below.
1330 0 : do m_cnst = 1, qsize
1331 0 : do ie = 1, nelemd
1332 0 : qtmp(:,:,:,ie,m_cnst) = 0.0_r8
1333 : indx = 1
1334 0 : do j = 1, np
1335 0 : do i = 1, np
1336 : ! Set qtmp at the unique columns only
1337 0 : if (pmask(((ie - 1) * npsq) + indx)) then
1338 0 : qtmp(i,j,:,ie,m_cnst) = dbuf4(indx, :, ie, m_cnst)
1339 : end if
1340 0 : indx = indx + 1
1341 : end do
1342 : end do
1343 : end do
1344 : end do
1345 0 : deallocate(dbuf4)
1346 :
1347 :
1348 : else
1349 :
1350 : ! Read ICs from file. Assume all fields in the initial file are on the GLL grid.
1351 :
1352 2304 : allocate(dbuf2(npsq,nelemd))
1353 2304 : allocate(dbuf3(npsq,nlev,nelemd))
1354 :
1355 : ! Check that columns in IC file match grid definition.
1356 768 : call check_file_layout(fh_ini, elem, dyn_cols, 'ncdata', .true.)
1357 :
1358 : ! Read 2-D field
1359 :
1360 768 : fieldname = 'PS'
1361 768 : fieldname2 = 'PSDRY'
1362 768 : if (dyn_field_exists(fh_ini, trim(fieldname), required=.false.)) then
1363 768 : inic_wet = .true.
1364 768 : call read_dyn_var(trim(fieldname), fh_ini, ini_grid_hdim_name, dbuf2)
1365 0 : elseif (dyn_field_exists(fh_ini, trim(fieldname2), required=.false.)) then
1366 0 : inic_wet = .false.
1367 0 : call read_dyn_var(trim(fieldname2), fh_ini, ini_grid_hdim_name, dbuf2)
1368 : else
1369 0 : call endrun(trim(sub)//': PS or PSDRY must be on GLL grid')
1370 : end if
1371 : #ifndef planet_mars
1372 768 : if (iam < par%nprocs) then
1373 94872 : if (minval(dbuf2, mask=reshape(pmask, (/npsq,nelemd/))) < 10000._r8) then
1374 768 : call endrun(trim(sub)//': Problem reading ps or psdry field -- bad values')
1375 : end if
1376 : end if
1377 : #endif
1378 6168 : do ie = 1, nelemd
1379 : indx = 1
1380 27768 : do j = 1, np
1381 113400 : do i = 1, np
1382 86400 : elem(ie)%state%psdry(i,j) = dbuf2(indx,ie) ! can be either wet or dry ps
1383 108000 : indx = indx + 1
1384 : end do
1385 : end do
1386 : end do
1387 :
1388 : ! Read in 3-D fields
1389 :
1390 768 : if (dyn_field_exists(fh_ini, 'U')) then
1391 768 : call read_dyn_var('U', fh_ini, ini_grid_hdim_name, dbuf3)
1392 : else
1393 0 : call endrun(trim(sub)//': U not found')
1394 : end if
1395 6168 : do ie = 1, nelemd
1396 64805400 : elem(ie)%state%v = 0.0_r8
1397 : indx = 1
1398 27768 : do j = 1, np
1399 113400 : do i = 1, np
1400 8121600 : elem(ie)%state%v(i,j,1,:,1) = dbuf3(indx,:,ie)
1401 108000 : indx = indx + 1
1402 : end do
1403 : end do
1404 : end do
1405 :
1406 768 : if (dyn_field_exists(fh_ini, 'V')) then
1407 768 : call read_dyn_var('V', fh_ini, ini_grid_hdim_name, dbuf3)
1408 : else
1409 0 : call endrun(trim(sub)//': V not found')
1410 : end if
1411 6168 : do ie = 1, nelemd
1412 : indx = 1
1413 27768 : do j = 1, np
1414 113400 : do i = 1, np
1415 8121600 : elem(ie)%state%v(i,j,2,:,1) = dbuf3(indx,:,ie)
1416 108000 : indx = indx + 1
1417 : end do
1418 : end do
1419 : end do
1420 :
1421 768 : if (dyn_field_exists(fh_ini, 'T')) then
1422 768 : call read_dyn_var('T', fh_ini, ini_grid_hdim_name, dbuf3)
1423 : else
1424 0 : call endrun(trim(sub)//': T not found')
1425 : end if
1426 6168 : do ie=1,nelemd
1427 31660200 : elem(ie)%state%T = 0.0_r8
1428 : indx = 1
1429 27768 : do j = 1, np
1430 113400 : do i = 1, np
1431 8121600 : elem(ie)%state%T(i,j,:,1) = dbuf3(indx,:,ie)
1432 108000 : indx = indx + 1
1433 : end do
1434 : end do
1435 : end do
1436 :
1437 768 : if (pertlim .ne. 0.0_r8) then
1438 0 : if (masterproc) then
1439 0 : write(iulog,*) trim(sub), ': Adding random perturbation bounded', &
1440 0 : 'by +/- ', pertlim, ' to initial temperature field'
1441 : end if
1442 :
1443 0 : call random_seed(size=rndm_seed_sz)
1444 0 : allocate(rndm_seed(rndm_seed_sz))
1445 :
1446 0 : do ie = 1, nelemd
1447 : ! seed random number generator based on element ID
1448 : ! (possibly include a flag to allow clock-based random seeding)
1449 0 : rndm_seed = elem(ie)%GlobalId
1450 0 : call random_seed(put=rndm_seed)
1451 0 : do i = 1, np
1452 0 : do j = 1, np
1453 0 : do k = 1, nlev
1454 0 : call random_number(pertval)
1455 0 : pertval = 2.0_r8*pertlim*(0.5_r8 - pertval)
1456 0 : elem(ie)%state%T(i,j,k,1) = elem(ie)%state%T(i,j,k,1)*(1.0_r8 + pertval)
1457 : end do
1458 : end do
1459 : end do
1460 : end do
1461 :
1462 0 : deallocate(rndm_seed)
1463 : end if
1464 :
1465 : ! Cleanup
1466 768 : deallocate(dbuf2)
1467 768 : deallocate(dbuf3)
1468 :
1469 : end if ! analytic_ic_active
1470 :
1471 : ! Read in or cold-initialize all the tracer fields.
1472 : ! Data is read in on the GLL grid.
1473 : ! Both GLL and FVM tracer fields are initialized based on the
1474 : ! dimension qsize or ntrac for GLL or FVM tracers respectively.
1475 : ! Data is only read in on GLL so if FVM tracers are active,
1476 : ! interpolation is performed.
1477 : !
1478 : ! If analytic ICs are being used, we allow constituents in an initial
1479 : ! file to overwrite mixing ratios set by the default constituent initialization
1480 : ! except for the water species.
1481 :
1482 768 : if (ntrac > qsize) then
1483 768 : if (ntrac < pcnst) then
1484 0 : write(errmsg, '(a,3(i0,a))') ': ntrac (',ntrac,') > qsize (',qsize, &
1485 0 : ') but < pcnst (',pcnst,')'
1486 0 : call endrun(trim(sub)//errmsg)
1487 : end if
1488 0 : else if (qsize < pcnst) then
1489 0 : write(errmsg, '(a,2(i0,a))') ': qsize (',qsize,') < pcnst (',pcnst,')'
1490 0 : call endrun(trim(sub)//errmsg)
1491 : end if
1492 :
1493 : ! If using analytic ICs the initial file only needs the horizonal grid
1494 : ! dimension checked in the case that the file contains constituent mixing
1495 : ! ratios.
1496 3072 : do m_cnst = 1, pcnst
1497 3072 : if (cnst_read_iv(m_cnst) .and. .not. cnst_is_a_water_species(cnst_name(m_cnst))) then
1498 768 : if (dyn_field_exists(fh_ini, trim(cnst_name(m_cnst)), required=.false.)) then
1499 768 : call check_file_layout(fh_ini, elem, dyn_cols, 'ncdata', .true.)
1500 768 : exit
1501 : end if
1502 : end if
1503 : end do
1504 :
1505 2304 : allocate(dbuf3(npsq,nlev,nelemd))
1506 :
1507 32256 : do m_cnst = 1, pcnst
1508 :
1509 31488 : if (analytic_ic_active() .and. cnst_is_a_water_species(cnst_name(m_cnst))) cycle
1510 :
1511 31488 : found = .false.
1512 31488 : if (cnst_read_iv(m_cnst)) then
1513 31488 : found = dyn_field_exists(fh_ini, trim(cnst_name(m_cnst)), required=.false.)
1514 : end if
1515 :
1516 31488 : if (found) then
1517 31488 : call read_dyn_var(trim(cnst_name(m_cnst)), fh_ini, ini_grid_hdim_name, dbuf3)
1518 : else
1519 0 : call cnst_init_default(m_cnst, latvals, lonvals, dbuf3, pmask)
1520 : end if
1521 :
1522 253656 : do ie = 1, nelemd
1523 : ! Copy tracers defined on GLL grid into Eulerian array
1524 : ! Make sure tracers have at least minimum value
1525 20843088 : do k=1, nlev
1526 : indx = 1
1527 103172400 : do j = 1, np
1528 432394200 : do i = 1, np
1529 : ! Set qtmp at the unique columns only: zero non-unique columns
1530 329443200 : if (pmask(((ie - 1) * npsq) + indx)) then
1531 185319426 : qtmp(i,j, k, ie, m_cnst) = max(qmin(m_cnst),dbuf3(indx,k,ie))
1532 : else
1533 144123774 : qtmp(i,j, k, ie, m_cnst) = 0.0_r8
1534 : end if
1535 411804000 : indx = indx + 1
1536 : end do
1537 : end do
1538 : end do
1539 : end do
1540 :
1541 : end do ! pcnst
1542 :
1543 : ! Cleanup
1544 768 : deallocate(dbuf3)
1545 :
1546 : ! Put the error handling back the way it was
1547 768 : call pio_seterrorhandling(fh_ini, pio_errtype)
1548 :
1549 : ! Cleanup
1550 768 : deallocate(pmask)
1551 768 : deallocate(latvals)
1552 768 : deallocate(lonvals)
1553 :
1554 768 : if (associated(ldof)) then
1555 768 : deallocate(ldof)
1556 : nullify(ldof)
1557 : end if
1558 :
1559 : ! once we've read or initialized all the fields we do a boundary exchange to
1560 : ! update the redundent columns in the dynamics
1561 768 : if(iam < par%nprocs) then
1562 768 : call initEdgeBuffer(par, edge, elem, (3+pcnst)*nlev + 2 )
1563 : end if
1564 6168 : do ie = 1, nelemd
1565 5400 : kptr = 0
1566 5400 : call edgeVpack(edge, elem(ie)%state%psdry,1,kptr,ie)
1567 5400 : kptr = kptr + 1
1568 5400 : call edgeVpack(edge, elem(ie)%state%v(:,:,:,:,1),2*nlev,kptr,ie)
1569 5400 : kptr = kptr + (2 * nlev)
1570 5400 : call edgeVpack(edge, elem(ie)%state%T(:,:,:,1),nlev,kptr,ie)
1571 5400 : kptr = kptr + nlev
1572 432621768 : call edgeVpack(edge, qtmp(:,:,:,ie,:),nlev*pcnst,kptr,ie)
1573 : end do
1574 768 : if(iam < par%nprocs) then
1575 768 : call bndry_exchange(par,edge,location='read_inidat')
1576 : end if
1577 6168 : do ie = 1, nelemd
1578 5400 : kptr = 0
1579 5400 : call edgeVunpack(edge, elem(ie)%state%psdry,1,kptr,ie)
1580 5400 : kptr = kptr + 1
1581 5400 : call edgeVunpack(edge, elem(ie)%state%v(:,:,:,:,1),2*nlev,kptr,ie)
1582 5400 : kptr = kptr + (2 * nlev)
1583 5400 : call edgeVunpack(edge, elem(ie)%state%T(:,:,:,1),nlev,kptr,ie)
1584 5400 : kptr = kptr + nlev
1585 865237368 : call edgeVunpack(edge, qtmp(:,:,:,ie,:),nlev*pcnst,kptr,ie)
1586 : end do
1587 :
1588 768 : if (inic_wet) then
1589 : !
1590 : ! convert to dry
1591 : !
1592 : ! (this has to be done after edge-exchange since shared points between elements are only
1593 : ! initialized in one element and not the other!)
1594 : !
1595 768 : if (par%masterproc) then
1596 1 : write(iulog,*) 'Convert specific/wet mixing ratios to dry'
1597 : end if
1598 :
1599 2304 : allocate(factor_array(np,np,nlev,nelemd))
1600 : !
1601 : ! compute: factor_array = 1/(1-sum(q))
1602 : !
1603 10552368 : factor_array(:,:,:,:) = 1.0_r8
1604 6168 : do ie = 1, nelemd
1605 38568 : do k = dry_air_species_num+1, thermodynamic_active_species_num
1606 32400 : m_cnst = thermodynamic_active_species_idx(k)
1607 63315000 : factor_array(:,:,:,ie) = factor_array(:,:,:,ie) - qtmp(:,:,:,ie,m_cnst)
1608 : end do
1609 : end do
1610 10552368 : factor_array(:,:,:,:) = 1.0_r8/factor_array(:,:,:,:)
1611 :
1612 32256 : do m_cnst = 1, pcnst
1613 32256 : if (cnst_type(m_cnst) == 'wet') then
1614 67848 : do ie = 1, nelemd
1615 5592048 : do k = 1, nlev
1616 27680400 : do j = 1, np
1617 116008200 : do i = 1, np
1618 :
1619 : ! convert wet mixing ratio to dry
1620 88387200 : qtmp(i,j,k,ie,m_cnst) = qtmp(i,j,k,ie,m_cnst) * factor_array(i,j,k,ie)
1621 :
1622 : ! truncate negative values if they were not analytically specified
1623 110484000 : if (.not. analytic_ic_active()) then
1624 88387200 : qtmp(i,j,k,ie,m_cnst) = max(qmin(m_cnst), qtmp(i,j,k,ie,m_cnst))
1625 : end if
1626 : end do
1627 : end do
1628 : end do
1629 : end do
1630 : end if
1631 : end do
1632 :
1633 : ! initialize dp3d and qdp
1634 : !
1635 : ! compute: factor_array = 1/(1+sum(q))
1636 :
1637 10552368 : factor_array(:,:,:,:) = 1.0_r8
1638 6168 : do ie = 1, nelemd
1639 38568 : do k = dry_air_species_num+1, thermodynamic_active_species_num
1640 32400 : m_cnst = thermodynamic_active_species_idx(k)
1641 63315000 : factor_array(:,:,:,ie) = factor_array(:,:,:,ie) + qtmp(:,:,:,ie,m_cnst)
1642 : end do
1643 : end do
1644 10552368 : factor_array(:,:,:,:) = 1.0_r8/factor_array(:,:,:,:)
1645 6168 : do ie = 1, nelemd
1646 : ! pstmp is the wet ps
1647 113400 : pstmp = elem(ie)%state%psdry(:,:)
1648 : ! start accumulating the dry air pressure differences across each layer
1649 113400 : elem(ie)%state%psdry(:,:) = hyai(1)*ps0
1650 508368 : do k=1,nlev
1651 2516400 : do j = 1,np
1652 10546200 : do i = 1,np
1653 16070400 : dp_tmp = ((hyai(k+1) - hyai(k))*ps0) + &
1654 24105600 : ((hybi(k+1) - hybi(k))*pstmp(i,j))
1655 8035200 : if (.not. analytic_ic_active()) then
1656 :
1657 : ! if analytic_ic then the surface pressure is already dry
1658 : ! (note that it is not correct to convert to moist pressure
1659 : ! in analytic_ic and not have the #ifndef statement here
1660 : ! since the dry levels are in a different location than
1661 : ! what is obtained from algorithm below)
1662 :
1663 : ! convert dp_tmp to dry
1664 8035200 : dp_tmp = dp_tmp*factor_array(i,j,k,ie)
1665 : end if
1666 :
1667 32140800 : elem(ie)%state%dp3d(i,j,k,:) = dp_tmp
1668 :
1669 : ! compute dry surface pressure; note that at this point
1670 : !
1671 : ! dp3d .NE. (hyai(k+1) - hyai(k))*ps0 + (hybi(k+1) - hybi(k))*ps(i,j)
1672 :
1673 10044000 : elem(ie)%state%psdry(i,j) = elem(ie)%state%psdry(i,j)+elem(ie)%state%dp3d(i,j,k,1)
1674 : end do
1675 : end do
1676 : end do
1677 : end do
1678 :
1679 768 : deallocate(factor_array)
1680 :
1681 : else
1682 :
1683 : ! initial condition is based on dry surface pressure and constituents
1684 : !
1685 : ! we only need to initialize state%dp3d
1686 :
1687 0 : do ie = 1, nelemd
1688 0 : do k = 1, nlev
1689 0 : do j = 1, np
1690 0 : do i = 1, np
1691 0 : elem(ie)%state%dp3d(i,j,k,:) = (hyai(k+1) - hyai(k))*ps0 + &
1692 0 : (hybi(k+1) - hybi(k))*elem(ie)%state%psdry(i,j)
1693 : end do
1694 : end do
1695 : end do
1696 : end do
1697 : end if
1698 :
1699 : ! If scale_dry_air_mass > 0.0 then scale dry air mass to scale_dry_air_mass global average dry pressure
1700 768 : if (runtype == 0) then
1701 768 : if (scale_dry_air_mass > 0.0_r8) then
1702 768 : if (iam < par%nprocs) then
1703 768 : call prim_set_dry_mass(elem, hvcoord, scale_dry_air_mass, qtmp)
1704 : end if
1705 : end if
1706 : end if
1707 :
1708 : ! store Q values:
1709 : !
1710 : ! if CSLAM is NOT active then state%Qdp for all constituents
1711 : ! if CSLAM active then we only advect water vapor and condensate
1712 : ! loading tracers in state%qdp
1713 :
1714 768 : if (use_cslam) then
1715 6168 : do ie = 1, nelemd
1716 38568 : do nq = 1, thermodynamic_active_species_num
1717 32400 : m_cnst = thermodynamic_active_species_idx(nq)
1718 3051000 : do k = 1, nlev
1719 15098400 : do j = 1, np
1720 63277200 : do i = 1, np
1721 0 : elem(ie)%state%Qdp(i,j,k,nq,:) = &
1722 108475200 : elem(ie)%state%dp3d(i,j,k,1)*qtmp(i,j,k,ie,m_cnst)
1723 : end do
1724 : end do
1725 : end do
1726 : end do
1727 : end do
1728 : else
1729 0 : do ie = 1, nelemd
1730 0 : do m_cnst = 1, qsize
1731 0 : do k = 1, nlev
1732 0 : do j = 1, np
1733 0 : do i = 1, np
1734 0 : elem(ie)%state%Qdp(i,j,k,m_cnst,:)=&
1735 0 : elem(ie)%state%dp3d(i,j,k,1)*qtmp(i,j,k,ie,m_cnst)
1736 : end do
1737 : end do
1738 : end do
1739 : end do
1740 : end do
1741 : end if
1742 :
1743 : ! interpolate fvm tracers and fvm pressure variables
1744 :
1745 768 : if (use_cslam) then
1746 768 : if (par%masterproc) then
1747 1 : write(iulog,*) 'Initializing dp_fvm from spectral element dp'
1748 : end if
1749 :
1750 6168 : do ie = 1, nelemd
1751 :
1752 : ! note that the area over fvm cells as computed from subcell_integration is up to 1.0E-6
1753 : ! different than the areas (exact) computed by CSLAM
1754 : !
1755 : ! Map the constituents which are also to be transported by dycore
1756 0 : call dyn2fvm_mass_vars(elem(ie)%state%dp3d(:,:,:,1),elem(ie)%state%psdry(:,:),&
1757 5400 : qtmp(:,:,:,ie,1:ntrac),&
1758 0 : dyn_in%fvm(ie)%dp_fvm(1:nc,1:nc,:),dyn_in%fvm(ie)%psC(1:nc,1:nc),&
1759 0 : dyn_in%fvm(ie)%c(1:nc,1:nc,:,1:ntrac),&
1760 981472368 : ntrac,elem(ie)%metdet,dyn_in%fvm(ie)%inv_se_area_sphere(1:nc,1:nc))
1761 : end do
1762 :
1763 768 : if(par%masterproc) then
1764 1 : write(iulog,*) 'FVM tracers, FVM pressure variables and se_area_sphere initialized.'
1765 : end if
1766 :
1767 : end if ! (use_cslam)
1768 :
1769 : ! Cleanup
1770 768 : deallocate(qtmp)
1771 :
1772 6168 : do ie = 1, nelemd
1773 16968 : do t = 2, timelevels
1774 43200000 : elem(ie)%state%v(:,:,:,:,t) = elem(ie)%state%v(:,:,:,:,1)
1775 21108600 : elem(ie)%state%T(:,:,:,t) = elem(ie)%state%T(:,:,:,1)
1776 : end do
1777 : end do
1778 :
1779 768 : if(iam < par%nprocs) then
1780 768 : call FreeEdgeBuffer(edge)
1781 : end if
1782 :
1783 1536 : end subroutine read_inidat
1784 :
1785 :
1786 : !========================================================================================
1787 :
1788 1536 : subroutine set_phis(dyn_in)
1789 :
1790 : ! Set PHIS according to the following rules.
1791 : !
1792 : ! 1) If a topo file is specified use it. This option has highest precedence.
1793 : ! 2) If not using topo file, but analytic_ic option is on, use analytic phis.
1794 : ! 3) Set phis = 0.0.
1795 : !
1796 : ! If using the physics grid then the topo file will be on that grid since its
1797 : ! contents are primarily for the physics parameterizations, and the values of
1798 : ! PHIS should be consistent with the values of sub-grid variability (e.g., SGH)
1799 : ! which are computed on the physics grid. In this case phis on the physics grid
1800 : ! will be interpolated to the GLL grid.
1801 :
1802 :
1803 : ! Arguments
1804 : type (dyn_import_t), target, intent(inout) :: dyn_in ! dynamics import
1805 :
1806 : ! local variables
1807 : type(file_desc_t), pointer :: fh_topo
1808 :
1809 1536 : type(element_t), pointer :: elem(:)
1810 :
1811 1536 : real(r8), allocatable :: phis_tmp(:,:) ! (npsp,nelemd)
1812 1536 : real(r8), allocatable :: phis_phys_tmp(:,:) ! (fv_nphys**2,nelemd)
1813 :
1814 : integer :: i, ie, indx, j, kptr
1815 : integer :: ierr, pio_errtype
1816 :
1817 : character(len=max_fieldname_len) :: fieldname
1818 : character(len=max_fieldname_len) :: fieldname_gll
1819 : character(len=max_hcoordname_len):: grid_name
1820 : integer :: dims(2)
1821 : integer :: dyn_cols
1822 : integer :: ncol_did
1823 : integer :: ncol_size
1824 :
1825 1536 : integer(iMap), pointer :: ldof(:) ! Basic (2D) grid dof
1826 1536 : logical, allocatable :: pmask(:) ! (npsq*nelemd) unique columns
1827 :
1828 : ! Variables for analytic initial conditions
1829 1536 : integer, allocatable :: glob_ind(:)
1830 : logical, allocatable :: pmask_phys(:)
1831 1536 : real(r8), pointer :: latvals_deg(:)
1832 1536 : real(r8), pointer :: lonvals_deg(:)
1833 1536 : real(r8), allocatable :: latvals(:)
1834 1536 : real(r8), allocatable :: lonvals(:)
1835 : real(r8), allocatable :: latvals_phys(:)
1836 : real(r8), allocatable :: lonvals_phys(:)
1837 :
1838 : character(len=*), parameter :: sub='set_phis'
1839 : !----------------------------------------------------------------------------
1840 :
1841 3072 : fh_topo => topo_file_get_id()
1842 :
1843 1536 : if (iam < par%nprocs) then
1844 1536 : elem => dyn_in%elem
1845 : else
1846 0 : nullify(elem)
1847 : end if
1848 :
1849 4608 : allocate(phis_tmp(npsq,nelemd))
1850 185136 : phis_tmp = 0.0_r8
1851 :
1852 1536 : if (use_cslam) then
1853 6144 : allocate(phis_phys_tmp(fv_nphys**2,nelemd))
1854 109536 : phis_phys_tmp = 0.0_r8
1855 12336 : do ie=1,nelemd
1856 53245536 : elem(ie)%sub_elem_mass_flux=0.0_r8
1857 : #ifdef waccm_debug
1858 : dyn_in%fvm(ie)%CSLAM_gamma = 0.0_r8
1859 : #endif
1860 : end do
1861 : end if
1862 :
1863 : ! Set mask to indicate which columns are active in GLL grid.
1864 1536 : nullify(ldof)
1865 1536 : call cam_grid_get_gcid(cam_grid_id('GLL'), ldof)
1866 4608 : allocate(pmask(npsq*nelemd))
1867 174336 : pmask(:) = (ldof /= 0)
1868 1536 : deallocate(ldof)
1869 :
1870 1536 : if (associated(fh_topo)) then
1871 :
1872 : ! Set PIO to return error flags.
1873 1536 : call pio_seterrorhandling(fh_topo, PIO_BCAST_ERROR, pio_errtype)
1874 :
1875 : ! Set name of grid object which will be used to read data from file
1876 : ! into internal data structure via PIO.
1877 1536 : if (.not.use_cslam) then
1878 0 : grid_name = 'GLL'
1879 : else
1880 1536 : grid_name = 'physgrid_d'
1881 : end if
1882 :
1883 : ! Get number of global columns from the grid object and check that
1884 : ! it matches the file data.
1885 1536 : call cam_grid_dimensions(grid_name, dims)
1886 1536 : dyn_cols = dims(1)
1887 :
1888 : ! The dimension of the unstructured grid in the TOPO file is 'ncol'.
1889 1536 : ierr = pio_inq_dimid(fh_topo, 'ncol', ncol_did)
1890 1536 : if (ierr /= PIO_NOERR) then
1891 0 : call endrun(sub//': dimension ncol not found in bnd_topo file')
1892 : end if
1893 1536 : ierr = pio_inq_dimlen(fh_topo, ncol_did, ncol_size)
1894 1536 : if (ncol_size /= dyn_cols) then
1895 0 : if (masterproc) then
1896 0 : write(iulog,*) sub//': ncol_size=', ncol_size, ' : dyn_cols=', dyn_cols
1897 : end if
1898 0 : call endrun(sub//': ncol size in bnd_topo file does not match grid definition')
1899 : end if
1900 :
1901 1536 : fieldname = 'PHIS'
1902 1536 : fieldname_gll = 'PHIS_gll'
1903 1536 : if (use_cslam.and.dyn_field_exists(fh_topo, trim(fieldname_gll),required=.false.)) then
1904 : !
1905 : ! If physgrid it is recommended to read in PHIS on the GLL grid and then
1906 : ! map to the physgrid in d_p_coupling
1907 : !
1908 : ! This requires a topo file with PHIS_gll on it ...
1909 : !
1910 1536 : if (masterproc) then
1911 2 : write(iulog, *) "Reading in PHIS on GLL grid (mapped to physgrid in d_p_coupling)"
1912 : end if
1913 1536 : call read_dyn_var(fieldname_gll, fh_topo, 'ncol_gll', phis_tmp)
1914 0 : else if (dyn_field_exists(fh_topo, trim(fieldname))) then
1915 0 : if (.not.use_cslam) then
1916 0 : if (masterproc) then
1917 0 : write(iulog, *) "Reading in PHIS"
1918 : end if
1919 0 : call read_dyn_var(fieldname, fh_topo, 'ncol', phis_tmp)
1920 : else
1921 : !
1922 : ! For backwards compatibility we allow reading in PHIS on the physgrid
1923 : ! which is then mapped to the GLL grid and back to the physgrid in d_p_coupling
1924 : ! (the latter is to avoid noise in derived quantities such as PSL)
1925 : !
1926 0 : if (masterproc) then
1927 0 : write(iulog, *) "Reading in PHIS on physgrid"
1928 0 : write(iulog, *) "Recommended to read in PHIS on GLL grid"
1929 : end if
1930 0 : call read_phys_field_2d(fieldname, fh_topo, 'ncol', phis_phys_tmp)
1931 : call map_phis_from_physgrid_to_gll(dyn_in%fvm, elem, phis_phys_tmp, &
1932 0 : phis_tmp, pmask)
1933 0 : deallocate(phis_phys_tmp)
1934 : end if
1935 : else
1936 0 : call endrun(sub//': Could not find PHIS field on input datafile')
1937 : end if
1938 :
1939 : ! Put the error handling back the way it was
1940 1536 : call pio_seterrorhandling(fh_topo, pio_errtype)
1941 :
1942 0 : else if (analytic_ic_active() .and. (iam < par%nprocs)) then
1943 :
1944 : ! lat/lon needed in radians
1945 0 : latvals_deg => cam_grid_get_latvals(cam_grid_id('GLL'))
1946 0 : lonvals_deg => cam_grid_get_lonvals(cam_grid_id('GLL'))
1947 0 : allocate(latvals(np*np*nelemd))
1948 0 : allocate(lonvals(np*np*nelemd))
1949 0 : latvals(:) = latvals_deg(:)*deg2rad
1950 0 : lonvals(:) = lonvals_deg(:)*deg2rad
1951 :
1952 0 : allocate(glob_ind(npsq*nelemd))
1953 0 : j = 1
1954 0 : do ie = 1, nelemd
1955 0 : do i = 1, npsq
1956 : ! Create a global(ish) column index
1957 0 : glob_ind(j) = elem(ie)%GlobalId
1958 0 : j = j + 1
1959 : end do
1960 : end do
1961 : call analytic_ic_set_ic(vcoord, latvals, lonvals, glob_ind, &
1962 0 : PHIS_OUT=phis_tmp, mask=pmask(:))
1963 0 : deallocate(glob_ind)
1964 :
1965 : end if
1966 :
1967 1536 : deallocate(pmask)
1968 :
1969 : ! Set PHIS in element objects
1970 12336 : do ie = 1, nelemd
1971 226800 : elem(ie)%state%phis = 0.0_r8
1972 : indx = 1
1973 55536 : do j = 1, np
1974 226800 : do i = 1, np
1975 172800 : elem(ie)%state%phis(i,j) = phis_tmp(indx, ie)
1976 216000 : indx = indx + 1
1977 : end do
1978 : end do
1979 : end do
1980 1536 : deallocate(phis_tmp)
1981 :
1982 : ! boundary exchange to update the redundent columns in the element objects
1983 12336 : do ie = 1, nelemd
1984 10800 : kptr = 0
1985 12336 : call edgeVpack(edgebuf, elem(ie)%state%phis, 1, kptr, ie)
1986 : end do
1987 1536 : if(iam < par%nprocs) then
1988 1536 : call bndry_exchange(par, edgebuf, location=sub)
1989 : end if
1990 12336 : do ie = 1, nelemd
1991 10800 : kptr = 0
1992 12336 : call edgeVunpack(edgebuf, elem(ie)%state%phis,1,kptr,ie)
1993 : end do
1994 :
1995 3840 : end subroutine set_phis
1996 :
1997 : !========================================================================================
1998 :
1999 1536 : subroutine check_file_layout(file, elem, dyn_cols, file_desc, dyn_ok)
2000 :
2001 : ! This routine is only called when data will be read from the initial file. It is not
2002 : ! called when the initial file is only supplying vertical coordinate info.
2003 :
2004 : type(file_desc_t), pointer :: file
2005 : type(element_t), pointer :: elem(:)
2006 : integer, intent(in) :: dyn_cols
2007 : character(len=*), intent(in) :: file_desc
2008 : logical, intent(in) :: dyn_ok ! .true. iff ncol_d is okay
2009 :
2010 : integer :: ncol_did, ncol_size
2011 : integer :: ierr
2012 : integer :: ie, i, j
2013 : integer :: indx
2014 3072 : real(r8) :: dbuf2(npsq, nelemd)
2015 : logical :: found
2016 : character(len=max_fieldname_len) :: coordname
2017 :
2018 : character(len=*), parameter :: sub = 'check_file_layout'
2019 : !----------------------------------------------------------------------------
2020 :
2021 : ! Check that number of columns in IC file matches grid definition.
2022 1536 : if (trim(ini_grid_hdim_name) == 'none') then
2023 : call endrun(sub//': ERROR: no horizontal dimension in initial data file. &
2024 0 : &Cannot read data from file')
2025 : end if
2026 :
2027 1536 : ierr = pio_inq_dimid(file, trim(ini_grid_hdim_name), ncol_did)
2028 1536 : if (ierr /= PIO_NOERR) then
2029 : call endrun(sub//': ERROR: '//trim(ini_grid_hdim_name)//' dimension not found in ' &
2030 0 : //trim(file_desc)//' file')
2031 : end if
2032 :
2033 1536 : ierr = pio_inq_dimlen(file, ncol_did, ncol_size)
2034 1536 : if (ncol_size /= dyn_cols) then
2035 0 : if (masterproc) then
2036 0 : write(iulog, '(a,2(a,i0))') trim(sub), ': ncol_size=', ncol_size, &
2037 0 : ' : dyn_cols=', dyn_cols
2038 : end if
2039 0 : call endrun(sub//': ERROR: dimension '//trim(ini_grid_hdim_name)//' size not same as in ncdata file')
2040 : end if
2041 :
2042 : ! Set coordinate name associated with ini_grid_hdim_name.
2043 1536 : if (trim(ini_grid_hdim_name) == 'ncol') then
2044 0 : coordname = 'lat'
2045 : else
2046 1536 : coordname = 'lat_d'
2047 : end if
2048 :
2049 : !! Check to make sure file is in correct order
2050 1536 : call read_dyn_var(coordname, file, ini_grid_hdim_name, dbuf2)
2051 1536 : found = .true.
2052 12336 : do ie = 1, nelemd
2053 10800 : indx = 1
2054 55536 : do j = 1, np
2055 226800 : do i = 1, np
2056 172800 : if (abs(dbuf2(indx,ie)) > 1.e-12_r8) then
2057 96484 : if (abs((elem(ie)%spherep(i,j)%lat*rad2deg - dbuf2(indx,ie)) / &
2058 : dbuf2(indx,ie)) > 1.0e-10_r8) then
2059 : write(iulog, '(2a,4(i0,a),f11.5,a,f11.5)') &
2060 0 : "ncdata file latitudes not in correct column order", &
2061 0 : ' on task ', iam, ': elem(', ie, ')%spherep(', i, &
2062 0 : ', ', j, ')%lat = ', elem(ie)%spherep(i,j)%lat, &
2063 0 : ' /= ', dbuf2(indx, ie)*deg2rad
2064 0 : call shr_sys_flush(iulog)
2065 0 : found = .false.
2066 : end if
2067 : end if
2068 216000 : indx = indx + 1
2069 : end do
2070 : end do
2071 : end do
2072 1536 : if (.not. found) then
2073 0 : call endrun("ncdata file latitudes not in correct column order")
2074 : end if
2075 :
2076 1536 : if (trim(ini_grid_hdim_name) == 'ncol') then
2077 0 : coordname = 'lon'
2078 : else
2079 1536 : coordname = 'lon_d'
2080 : end if
2081 :
2082 1536 : call read_dyn_var(coordname, file, ini_grid_hdim_name, dbuf2)
2083 12336 : do ie = 1, nelemd
2084 10800 : indx = 1
2085 55536 : do j = 1, np
2086 226800 : do i = 1, np
2087 172800 : if (abs(dbuf2(indx,ie)) > 1.e-12_r8) then
2088 97200 : if (abs((elem(ie)%spherep(i,j)%lon*rad2deg - dbuf2(indx,ie)) / &
2089 : dbuf2(indx,ie)) > 1.0e-10_r8) then
2090 : write(iulog, '(2a,4(i0,a),f11.5,a,f11.5)') &
2091 0 : "ncdata file longitudes not in correct column order", &
2092 0 : ' on task ', iam, ': elem(', ie, ')%spherep(', i, &
2093 0 : ', ', j, ')%lon = ', elem(ie)%spherep(i,j)%lon, &
2094 0 : ' /= ', dbuf2(indx, ie)*deg2rad
2095 0 : call shr_sys_flush(iulog)
2096 0 : found = .false.
2097 : end if
2098 : end if
2099 216000 : indx = indx + 1
2100 : end do
2101 : end do
2102 : end do
2103 1536 : if (.not. found) then
2104 0 : call endrun("ncdata file longitudes not in correct column order")
2105 : end if
2106 1536 : end subroutine check_file_layout
2107 :
2108 : !========================================================================================
2109 :
2110 36864 : logical function dyn_field_exists(fh, fieldname, required)
2111 :
2112 : use pio, only: var_desc_t, PIO_inq_varid
2113 : use pio, only: PIO_NOERR
2114 :
2115 : type(file_desc_t), intent(in) :: fh
2116 : character(len=*), intent(in) :: fieldname
2117 : logical, optional, intent(in) :: required
2118 :
2119 : ! Local variables
2120 : logical :: found
2121 : logical :: field_required
2122 : integer :: ret
2123 : type(var_desc_t) :: varid
2124 : character(len=128) :: errormsg
2125 : !--------------------------------------------------------------------------
2126 :
2127 36864 : if (present(required)) then
2128 34560 : field_required = required
2129 : else
2130 : field_required = .true.
2131 : end if
2132 :
2133 36864 : ret = PIO_inq_varid(fh, trim(fieldname), varid)
2134 36864 : found = (ret == PIO_NOERR)
2135 36864 : if (.not. found) then
2136 0 : if (field_required) then
2137 0 : write(errormsg, *) trim(fieldname),' was not present in the input file.'
2138 0 : call endrun('DYN_FIELD_EXISTS: '//errormsg)
2139 : end if
2140 : end if
2141 :
2142 36864 : dyn_field_exists = found
2143 :
2144 36864 : end function dyn_field_exists
2145 :
2146 : !========================================================================================
2147 :
2148 5376 : subroutine read_dyn_field_2d(fieldname, fh, dimname, buffer)
2149 :
2150 : ! Dummy arguments
2151 : character(len=*), intent(in) :: fieldname
2152 : type(file_desc_t), intent(inout) :: fh
2153 : character(len=*), intent(in) :: dimname
2154 : real(r8), intent(inout) :: buffer(:, :)
2155 :
2156 : ! Local variables
2157 : logical :: found
2158 : real(r8) :: fillvalue
2159 : !----------------------------------------------------------------------------
2160 :
2161 647976 : buffer = 0.0_r8
2162 : call infld(trim(fieldname), fh, dimname, 1, npsq, 1, nelemd, buffer, &
2163 5376 : found, gridname=ini_grid_name, fillvalue=fillvalue)
2164 5376 : if(.not. found) then
2165 0 : call endrun('READ_DYN_FIELD_2D: Could not find '//trim(fieldname)//' field on input datafile')
2166 : end if
2167 :
2168 : ! This code allows use of compiler option to set uninitialized values
2169 : ! to NaN. In that case infld can return NaNs where the element GLL points
2170 : ! are not "unique columns"
2171 : ! Set NaNs or fillvalue points to zero
2172 647976 : where (isnan(buffer) .or. (buffer==fillvalue)) buffer = 0.0_r8
2173 :
2174 5376 : end subroutine read_dyn_field_2d
2175 :
2176 : !========================================================================================
2177 :
2178 33792 : subroutine read_dyn_field_3d(fieldname, fh, dimname, buffer)
2179 :
2180 : ! Dummy arguments
2181 : character(len=*), intent(in) :: fieldname
2182 : type(file_desc_t), intent(inout) :: fh
2183 : character(len=*), intent(in) :: dimname
2184 : real(r8), intent(inout) :: buffer(:,:,:)
2185 :
2186 : ! Local variables
2187 : logical :: found
2188 : real(r8) :: fillvalue
2189 : !----------------------------------------------------------------------------
2190 :
2191 375916992 : buffer = 0.0_r8
2192 : call infld(trim(fieldname), fh, dimname, 'lev', 1, npsq, 1, nlev, &
2193 33792 : 1, nelemd, buffer, found, gridname=ini_grid_name, fillvalue=fillvalue)
2194 33792 : if(.not. found) then
2195 0 : call endrun('READ_DYN_FIELD_3D: Could not find '//trim(fieldname)//' field on input datafile')
2196 : end if
2197 :
2198 : ! This code allows use of compiler option to set uninitialized values
2199 : ! to NaN. In that case infld can return NaNs where the element GLL points
2200 : ! are not "unique columns"
2201 : ! Set NaNs or fillvalue points to zero
2202 375916992 : where (isnan(buffer) .or. (buffer == fillvalue)) buffer = 0.0_r8
2203 :
2204 33792 : end subroutine read_dyn_field_3d
2205 :
2206 : !========================================================================================
2207 :
2208 0 : subroutine read_phys_field_2d(fieldname, fh, dimname, buffer)
2209 :
2210 : ! Dummy arguments
2211 : character(len=*), intent(in) :: fieldname
2212 : type(file_desc_t), intent(inout) :: fh
2213 : character(len=*), intent(in) :: dimname
2214 : real(r8), intent(inout) :: buffer(:, :)
2215 :
2216 : ! Local variables
2217 : logical :: found
2218 : !----------------------------------------------------------------------------
2219 :
2220 : call infld(trim(fieldname), fh, dimname, 1, fv_nphys**2, 1, nelemd, buffer, &
2221 0 : found, gridname='physgrid_d')
2222 0 : if(.not. found) then
2223 0 : call endrun('READ_PHYS_FIELD_2D: Could not find '//trim(fieldname)//' field on input datafile')
2224 : end if
2225 :
2226 0 : end subroutine read_phys_field_2d
2227 :
2228 : !========================================================================================
2229 :
2230 0 : subroutine map_phis_from_physgrid_to_gll(fvm,elem,phis_phys_tmp,phis_tmp,pmask)
2231 :
2232 : use hybrid_mod, only: get_loop_ranges, config_thread_region
2233 : use dimensions_mod, only: nhc_phys
2234 : use fvm_mapping, only: phys2dyn
2235 : use thread_mod, only: horz_num_threads
2236 :
2237 : type(element_t), intent(inout) :: elem(:)
2238 : type (fvm_struct), intent(in) :: fvm(:)
2239 : real(r8) , intent(in) :: phis_phys_tmp(fv_nphys**2,nelemd) !physgrid phis
2240 : real(r8) , intent(inout) :: phis_tmp(npsq,nelemd) !gll phis
2241 : logical , intent(in) :: pmask(npsq*nelemd)
2242 :
2243 : type(hybrid_t) :: hybrid
2244 : integer :: nets, nete, ie,i,j,indx
2245 0 : real(r8), allocatable :: fld_phys(:,:,:,:,:),fld_gll(:,:,:,:,:)
2246 : logical :: llimiter(1)
2247 : !----------------------------------------------------------------------------
2248 :
2249 : !!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,ie)
2250 : !hybrid = config_thread_region(par,'horizontal')
2251 0 : hybrid = config_thread_region(par,'serial')
2252 :
2253 0 : call get_loop_ranges(hybrid, ibeg=nets, iend=nete)
2254 :
2255 0 : allocate(fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,1,1,nets:nete))
2256 0 : allocate(fld_gll(np,np,1,1,nets:nete))
2257 0 : fld_phys = 0.0_r8
2258 0 : do ie = nets, nete
2259 0 : fld_phys(1:fv_nphys,1:fv_nphys,1,1,ie) = RESHAPE(phis_phys_tmp(:,ie),(/fv_nphys,fv_nphys/))
2260 : end do
2261 0 : llimiter = .true.
2262 0 : call phys2dyn(hybrid,elem,fld_phys,fld_gll,nets,nete,1,1,fvm,llimiter,halo_filled=.false.)
2263 0 : do ie = nets,nete
2264 : indx = 1
2265 0 : do j = 1, np
2266 0 : do i = 1, np
2267 0 : if (pmask(((ie - 1) * npsq) + indx)) then
2268 0 : phis_tmp(indx,ie) = fld_gll(i,j,1,1,ie)
2269 : else
2270 0 : phis_tmp(indx,ie) = 0.0_r8
2271 : end if
2272 0 : indx = indx + 1
2273 : end do
2274 : end do
2275 : end do
2276 0 : deallocate(fld_phys)
2277 0 : deallocate(fld_gll)
2278 : !!$OMP END PARALLEL
2279 0 : end subroutine map_phis_from_physgrid_to_gll
2280 :
2281 : !========================================================================================
2282 :
2283 14592 : subroutine write_dyn_vars(dyn_out)
2284 :
2285 : type (dyn_export_t), intent(inout) :: dyn_out ! Dynamics export container
2286 :
2287 : character(len=fieldname_len) :: tfname
2288 : integer :: ie, m
2289 : !----------------------------------------------------------------------------
2290 :
2291 14592 : if (use_cslam) then
2292 117192 : do ie = 1, nelemd
2293 0 : call outfld('dp_fvm', RESHAPE(dyn_out%fvm(ie)%dp_fvm(1:nc,1:nc,:), &
2294 102600 : (/nc*nc,nlev/)), nc*nc, ie)
2295 0 : call outfld('PSDRY_fvm', RESHAPE(dyn_out%fvm(ie)%psc(1:nc,1:nc), &
2296 102600 : (/nc*nc/)), nc*nc, ie)
2297 4323792 : do m = 1, ntrac
2298 4206600 : tfname = trim(cnst_name(m))//'_fvm'
2299 0 : call outfld(tfname, RESHAPE(dyn_out%fvm(ie)%c(1:nc,1:nc,:,m), &
2300 4206600 : (/nc*nc,nlev/)), nc*nc, ie)
2301 :
2302 4206600 : tfname = 'F'//trim(cnst_name(m))//'_fvm'
2303 8413200 : call outfld(tfname, RESHAPE(dyn_out%fvm(ie)%fc(1:nc,1:nc,:,m),&
2304 12722400 : (/nc*nc,nlev/)), nc*nc, ie)
2305 : end do
2306 : end do
2307 : end if
2308 :
2309 0 : end subroutine write_dyn_vars
2310 :
2311 : !=========================================================================================
2312 0 : end module dyn_comp
|
