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