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