Line data Source code
1 : !#define _DBG_ print *,"file: ",__FILE__," line: ",__LINE__," ithr: ",hybrid%ithr
2 : #define _DBG_
3 : module prim_driver_mod
4 : use shr_kind_mod, only: r8=>shr_kind_r8
5 : use cam_logfile, only: iulog
6 : use cam_abortutils, only: endrun
7 : use dimensions_mod, only: np, nlev, nelem, nelemd, GlobalUniqueCols, qsize, nc,nhc
8 : use hybrid_mod, only: hybrid_t, config_thread_region, PrintHybrid
9 : use derivative_mod, only: derivative_t
10 : use fvm_control_volume_mod, only: fvm_struct
11 :
12 : use element_mod, only: element_t, timelevels, allocate_element_desc
13 : use thread_mod , only: horz_num_threads, vert_num_threads, tracer_num_threads
14 : use thread_mod , only: omp_set_nested
15 : use perf_mod, only: t_startf, t_stopf
16 : use prim_init, only: gp, fvm_corners, fvm_points
17 :
18 : implicit none
19 : private
20 : public :: prim_init2, prim_run_subcycle, prim_finalize
21 : public :: prim_set_dry_mass
22 : contains
23 :
24 : !=============================================================================!
25 :
26 1536 : subroutine prim_init2(elem, fvm, hybrid, nets, nete, tl, hvcoord)
27 : use dimensions_mod, only: irecons_tracer, fvm_supercycling
28 : use dimensions_mod, only: fv_nphys, nc
29 : use parallel_mod, only: syncmp
30 : use se_dyn_time_mod, only: timelevel_t, tstep, phys_tscale, nsplit, TimeLevel_Qdp
31 : use se_dyn_time_mod, only: nsplit_baseline,rsplit_baseline
32 : use prim_state_mod, only: prim_printstate
33 : use control_mod, only: runtype, topology, rsplit, qsplit, rk_stage_user, &
34 : nu, nu_q, nu_div, hypervis_subcycle, hypervis_subcycle_q, &
35 : hypervis_subcycle_sponge, variable_nsplit
36 : use fvm_mod, only: fill_halo_fvm,ghostBufQnhc_h
37 : use thread_mod, only: omp_get_thread_num
38 : use global_norms_mod, only: print_cfl
39 : use hybvcoord_mod, only: hvcoord_t
40 : use prim_advection_mod, only: prim_advec_init2,deriv
41 : use prim_advance_mod, only: compute_omega
42 : use physconst, only: rga, cappa, cpair, tref, lapse_rate
43 : use cam_thermo, only: get_dp_ref
44 : use physconst, only: pstd
45 :
46 : type (element_t), intent(inout) :: elem(:)
47 : type (fvm_struct), intent(inout) :: fvm(:)
48 : type (hybrid_t), intent(in) :: hybrid
49 :
50 : type (TimeLevel_t), intent(inout) :: tl ! time level struct
51 : type (hvcoord_t), intent(inout) :: hvcoord ! hybrid vertical coordinate struct
52 :
53 : integer, intent(in) :: nets ! starting thread element number (private)
54 : integer, intent(in) :: nete ! ending thread element number (private)
55 :
56 :
57 : ! ==================================
58 : ! Local variables
59 : ! ==================================
60 :
61 : ! variables used to calculate CFL
62 : real (kind=r8) :: dtnu ! timestep*viscosity parameter
63 : real (kind=r8) :: dt_dyn_del2_sponge
64 : real (kind=r8) :: dt_tracer_vis ! viscosity timestep used in tracers
65 : real (kind=r8) :: dt_dyn_vis ! viscosity timestep
66 : real (kind=r8) :: dt_remap ! remapping timestep
67 :
68 : real (kind=r8) :: dp,dp0,T1,T0,pmid_ref(np,np)
69 3072 : real (kind=r8) :: ps_ref(np,np,nets:nete)
70 :
71 : integer :: i,j,k,ie,t,q
72 : integer :: n0,n0_qdp
73 :
74 :
75 12336 : do ie=nets,nete
76 43200000 : elem(ie)%derived%FM=0.0_r8
77 21103200 : elem(ie)%derived%FT=0.0_r8
78 126631536 : elem(ie)%derived%FQ=0.0_r8
79 : end do
80 :
81 : ! ==========================
82 : ! begin executable code
83 : ! ==========================
84 : !call prim_advance_init(hybrid%par,elem)
85 :
86 : ! compute most restrictive dt*nu for use by variable res viscosity:
87 : ! compute timestep seen by viscosity operator:
88 1536 : dt_dyn_vis = tstep
89 1536 : dt_dyn_del2_sponge = tstep
90 1536 : dt_tracer_vis=tstep*qsplit
91 1536 : dt_remap=dt_tracer_vis*rsplit
92 : ! compute most restrictive condition:
93 : ! note: dtnu ignores subcycling
94 1536 : dtnu=max(dt_dyn_vis*max(nu,nu_div), dt_tracer_vis*nu_q)
95 : ! compute actual viscosity timesteps with subcycling
96 1536 : dt_tracer_vis = dt_tracer_vis/hypervis_subcycle_q
97 1536 : dt_dyn_vis = dt_dyn_vis/hypervis_subcycle
98 1536 : dt_dyn_del2_sponge = dt_dyn_del2_sponge/hypervis_subcycle_sponge
99 1536 : if (variable_nsplit) then
100 0 : nsplit_baseline=nsplit
101 0 : rsplit_baseline=rsplit
102 : end if
103 : ! ==================================
104 : ! Initialize derivative structure
105 : ! ==================================
106 1536 : call Prim_Advec_Init2(fvm_corners, fvm_points)
107 1536 : if (fv_nphys>0.and.nc.ne.fv_nphys) then
108 : !
109 : ! need to fill halo for dp_coupling for fvm2phys mapping
110 : !
111 0 : call fill_halo_fvm(ghostBufQnhc_h,elem,fvm,hybrid,nets,nete,nhc,1,nlev,nlev)
112 : end if
113 : ! !$OMP BARRIER
114 : ! if (hybrid%ithr==0) then
115 : ! call syncmp(hybrid%par)
116 : ! end if
117 : ! !$OMP BARRIER
118 :
119 1536 : if (topology /= "cube") then
120 0 : call endrun('Error: only cube topology supported for primaitve equations')
121 : endif
122 :
123 : ! CAM has set tstep based on dtime before calling prim_init2(),
124 : ! so only now does HOMME learn the timstep. print them out:
125 : call print_cfl(elem,hybrid,nets,nete,dtnu,&
126 : !p top and p mid levels
127 : hvcoord%hyai(1)*hvcoord%ps0,hvcoord%hyam(:)*hvcoord%ps0+hvcoord%hybm(:)*pstd,&
128 : !dt_remap,dt_tracer_fvm,dt_tracer_se
129 : tstep*qsplit*rsplit,tstep*qsplit*fvm_supercycling,tstep*qsplit,&
130 : !dt_dyn,dt_dyn_visco,dt_tracer_visco, dt_phys
131 144384 : tstep,dt_dyn_vis,dt_dyn_del2_sponge,dt_tracer_vis,tstep*nsplit*qsplit*rsplit)
132 :
133 1536 : if (hybrid%masterthread) then
134 2 : if (phys_tscale/=0) then
135 0 : write(iulog,'(a,2f9.2)') "CAM physics timescale: ",phys_tscale
136 : endif
137 2 : write(iulog,'(a,2f9.2)') "CAM dtime (dt_phys): ",tstep*nsplit*qsplit*rsplit
138 :
139 2 : write(iulog,*) "CAM-SE uses dry-mass vertical coordinates"
140 : end if
141 :
142 1536 : n0=tl%n0
143 1536 : call TimeLevel_Qdp( tl, qsplit, n0_qdp)
144 1536 : call compute_omega(hybrid,n0,n0_qdp,elem,deriv,nets,nete,dt_remap,hvcoord)
145 : !
146 : ! pre-compute pressure-level thickness reference profile
147 : !
148 12336 : do ie=nets,nete
149 10800 : call get_dp_ref(hvcoord%hyai, hvcoord%hybi, hvcoord%ps0, elem(ie)%state%phis(:,:), &
150 23136 : elem(ie)%derived%dp_ref(:,:,:), ps_ref(:,:,ie))
151 : end do
152 : !
153 : ! pre-compute reference temperature profile (Simmons and Jiabin, 1991, QJRMS, Section 2a
154 : ! doi: https://doi.org/10.1002/qj.49711749703c)
155 : !
156 : ! Tref = T0+T1*Exner
157 : ! T1 = .0065*Tref*Cp/g ! = ~191
158 : ! T0 = Tref-T1 ! = ~97
159 : !
160 1536 : T1 = lapse_rate*Tref*cpair*rga
161 1536 : T0 = Tref-T1
162 12336 : do ie=nets,nete
163 1016736 : do k=1,nlev
164 21092400 : pmid_ref =hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*ps_ref(:,:,ie)
165 1004400 : dp0 = ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + &
166 1004400 : ( hvcoord%hybi(k+1) - hvcoord%hybi(k) )*hvcoord%ps0
167 1015200 : if (hvcoord%hybm(k)>0) then
168 11340000 : elem(ie)%derived%T_ref(:,:,k) = T0+T1*(pmid_ref/hvcoord%ps0)**cappa
169 : !
170 : ! pel@ucar.edu: resolved noise issue over Antartica
171 : !
172 11340000 : elem(ie)%derived%dp_ref(:,:,k) = elem(ie)%derived%dp_ref(:,:,k)-dp0
173 : else
174 9752400 : elem(ie)%derived%T_ref(:,:,k) = 0.0_r8
175 : end if
176 : end do
177 : end do
178 :
179 1536 : if (hybrid%masterthread) write(iulog,*) "initial state:"
180 1536 : call prim_printstate(elem, tl, hybrid,nets,nete, fvm)
181 :
182 1536 : end subroutine prim_init2
183 :
184 : !=======================================================================================================!
185 :
186 :
187 738816 : subroutine prim_run_subcycle(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord,nsubstep, single_column, omega_cn)
188 : !
189 : ! advance all variables (u,v,T,ps,Q,C) from time t to t + dt_q
190 : !
191 : ! for the RK schemes:
192 : ! input:
193 : ! tl%nm1 not used
194 : ! tl%n0 data at time t
195 : ! tl%np1 new values at t+dt_q
196 : !
197 : ! then we update timelevel pointers:
198 : ! tl%nm1 = tl%n0
199 : ! tl%n0 = tl%np1
200 : ! so that:
201 : ! tl%nm1 tracers: t dynamics: t+(qsplit-1)*dt
202 : ! tl%n0 time t + dt_q
203 : !
204 : ! for the implicit schemes:
205 : !
206 : ! input:
207 : ! tl%nm1 variables at t-1 level are stored fro BDF2 scheme
208 : ! tl%n0 data at time t
209 : ! tl%np1 new values at t+dt_q
210 : ! generally dt_q = t for BDF2, so its t+1
211 : !
212 : ! then we update timelevel pointers:
213 : ! tl%nm1 = tl%n0
214 : ! tl%n0 = tl%np1
215 : ! so that:
216 : ! tl%nm1 tracers: t dynamics: t+(qsplit-1)*dt
217 : ! tl%n0 time t + dt_q
218 : !
219 : !
220 1536 : use hybvcoord_mod, only : hvcoord_t
221 : use se_dyn_time_mod, only: TimeLevel_t, timelevel_update, timelevel_qdp, nsplit
222 : use control_mod, only: statefreq,qsplit, rsplit, variable_nsplit, dribble_in_rsplit_loop
223 : use prim_advance_mod, only: applycamforcing
224 : use prim_advance_mod, only: tot_energy_dyn,compute_omega
225 : use prim_state_mod, only: prim_printstate, adjust_nsplit
226 : use prim_advection_mod, only: vertical_remap, deriv
227 : use thread_mod, only: omp_get_thread_num
228 : use perf_mod , only: t_startf, t_stopf
229 : use fvm_mod , only: fill_halo_fvm, ghostBufQnhc_h
230 : use dimensions_mod, only: use_cslam,fv_nphys
231 : use fvm_mapping, only: cslam2gll
232 : type (element_t) , intent(inout) :: elem(:)
233 : type(fvm_struct), intent(inout) :: fvm(:)
234 : type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared)
235 : type (hvcoord_t), intent(in) :: hvcoord ! hybrid vertical coordinate struct
236 : integer, intent(in) :: nets ! starting thread element number (private)
237 : integer, intent(in) :: nete ! ending thread element number (private)
238 : real(kind=r8), intent(in) :: dt ! "timestep dependent" timestep
239 : type (TimeLevel_t), intent(inout):: tl
240 : integer, intent(in) :: nsubstep ! nsubstep = 1 .. nsplit
241 : logical, intent(in) :: single_column
242 : real (kind=r8) , intent(inout):: omega_cn(2,nets:nete) !min and max of vertical Courant number
243 :
244 : real(kind=r8) :: dt_q, dt_remap, dt_phys
245 : integer :: ie, q,k,n0_qdp,np1_qdp,r, nstep_end,region_num_threads,i,j
246 : real (kind=r8) :: dp_np1(np,np)
247 1477632 : real (kind=r8) :: dp_start(np,np,nlev+1,nets:nete),dp_end(np,np,nlev,nets:nete)
248 : logical :: compute_diagnostics
249 : ! ===================================
250 : ! Main timestepping loop
251 : ! ===================================
252 738816 : dt_q = dt*qsplit
253 738816 : nstep_end = tl%nstep + qsplit
254 738816 : dt_remap=dt_q*rsplit
255 738816 : nstep_end = tl%nstep + qsplit*rsplit ! nstep at end of this routine
256 738816 : dt_phys = nsplit*dt_remap
257 :
258 : ! compute diagnostics for STDOUT
259 738816 : compute_diagnostics=.false.
260 :
261 738816 : if (statefreq>0) then
262 0 : if (MODULO(nstep_end,statefreq)==0 .or. nstep_end==tl%nstep0) then
263 0 : compute_diagnostics=.true.
264 : endif
265 : end if
266 : !
267 : ! initialize variables for computing vertical Courant number
268 : !
269 738816 : if (variable_nsplit.or.compute_diagnostics) then
270 0 : if (nsubstep==1) then
271 0 : do ie=nets,nete
272 0 : omega_cn(1,ie) = 0.0_r8
273 0 : omega_cn(2,ie) = 0.0_r8
274 : end do
275 : end if
276 0 : do ie=nets,nete
277 0 : dp_start(:,:,1:nlev,ie) = elem(ie)%state%dp3d(:,:,:,tl%n0)
278 0 : dp_start(:,:,nlev+1,ie) = elem(ie)%state%dp3d(:,:,nlev,tl%n0)
279 : end do
280 : endif
281 :
282 :
283 738816 : call TimeLevel_Qdp( tl, qsplit, n0_qdp)
284 :
285 738816 : if (dribble_in_rsplit_loop==0) then
286 738816 : call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dAF')
287 738816 : call ApplyCAMForcing(elem,fvm,tl%n0,n0_qdp,dt_remap,dt_phys,nets,nete,nsubstep)
288 738816 : call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dBD')
289 : end if
290 2955264 : do r=1,rsplit
291 2216448 : if (r.ne.1) call TimeLevel_update(tl,"leapfrog")
292 : !
293 : ! if nsplit==1 and physics time-step is long then there will be noise in the
294 : ! pressure field; hence "dripple" in tendencies
295 : !
296 2216448 : if (dribble_in_rsplit_loop==1) then
297 0 : call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dAF')
298 0 : call ApplyCAMForcing(elem,fvm,tl%n0,n0_qdp,dt,dt_phys,nets,nete,MAX(nsubstep,r))
299 0 : call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dBD')
300 : end if
301 : !
302 : ! right after physics overwrite Qdp with CSLAM values
303 : !
304 2216448 : if (use_cslam.and.nsubstep==1.and.r==1) then
305 369408 : call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dAF')
306 369408 : call cslam2gll(elem, fvm, hybrid,nets,nete, tl%n0, n0_qdp)
307 369408 : call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dBD')
308 : end if
309 2216448 : call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dBL')
310 2216448 : if (single_column) then
311 : ! Single Column Case
312 : ! Loop over rsplit vertically lagrangian timesteps
313 0 : call prim_step_scm(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord,r)
314 : else
315 2216448 : call prim_step(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord,r,nsubstep==nsplit,dt_remap)
316 : end if
317 2955264 : call tot_energy_dyn(elem,fvm,nets,nete,tl%np1,n0_qdp,'dAL')
318 : enddo
319 :
320 :
321 : ! defer final timelevel update until after remap and diagnostics
322 738816 : call TimeLevel_Qdp( tl, qsplit, n0_qdp, np1_qdp)
323 :
324 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
325 : !
326 : ! apply vertical remap
327 : ! always for tracers
328 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329 :
330 738816 : call tot_energy_dyn(elem,fvm,nets,nete,tl%np1,np1_qdp,'dAD')
331 :
332 738816 : if (variable_nsplit.or.compute_diagnostics) then
333 : !
334 : ! initialize variables for computing vertical Courant number
335 : !
336 0 : do ie=nets,nete
337 0 : dp_end(:,:,:,ie) = elem(ie)%state%dp3d(:,:,:,tl%np1)
338 : end do
339 : end if
340 738816 : call t_startf('vertical_remap')
341 738816 : call vertical_remap(hybrid,elem,fvm,hvcoord,tl%np1,np1_qdp,nets,nete)
342 738816 : call t_stopf('vertical_remap')
343 :
344 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
345 : ! time step is complete.
346 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
347 738816 : call tot_energy_dyn(elem,fvm,nets,nete,tl%np1,np1_qdp,'dAR')
348 :
349 738816 : if (nsubstep==nsplit.and. .not. single_column) then
350 369408 : call compute_omega(hybrid,tl%np1,np1_qdp,elem,deriv,nets,nete,dt_remap,hvcoord)
351 : end if
352 :
353 : ! now we have:
354 : ! u(nm1) dynamics at t+dt_remap - 2*dt
355 : ! u(n0) dynamics at t+dt_remap - dt
356 : ! u(np1) dynamics at t+dt_remap
357 : !
358 : ! Q(1) Q at t+dt_remap
359 :
360 :
361 :
362 : ! =================================
363 : ! update dynamics time level pointers
364 : ! =================================
365 738816 : call TimeLevel_update(tl,"leapfrog")
366 : ! note: time level update for fvm tracers takes place in fvm_mod
367 :
368 : ! now we have:
369 : ! u(nm1) dynamics at t+dt_remap - dt (Robert-filtered)
370 : ! u(n0) dynamics at t+dt_remap
371 : ! u(np1) undefined
372 :
373 :
374 : !
375 : ! Compute vertical Courant numbers
376 : !
377 738816 : if (variable_nsplit.or.compute_diagnostics) then
378 0 : do ie=nets,nete
379 0 : do k=1,nlev
380 0 : do j=1,np
381 0 : do i=1,np
382 0 : if (dp_end(i,j,k,ie)<dp_start(i,j,k,ie)) then
383 0 : omega_cn(1,ie) = MIN((dp_end(i,j,k,ie)-dp_start(i,j,k,ie))/dp_start(i,j,k,ie),omega_cn(1,ie))
384 0 : omega_cn(2,ie) = MAX((dp_end(i,j,k,ie)-dp_start(i,j,k,ie))/dp_start(i,j,k,ie),omega_cn(2,ie))
385 : else
386 0 : omega_cn(1,ie) = MIN((dp_end(i,j,k,ie)-dp_start(i,j,k,ie))/dp_start(i,j,k+1,ie),omega_cn(1,ie))
387 0 : omega_cn(2,ie) = MAX((dp_end(i,j,k,ie)-dp_start(i,j,k,ie))/dp_start(i,j,k+1,ie),omega_cn(2,ie))
388 : end if
389 : end do
390 : end do
391 : end do
392 : end do
393 0 : if (nsubstep==nsplit.and.variable_nsplit) then
394 0 : call t_startf('adjust_nsplit')
395 0 : call adjust_nsplit(elem, tl, hybrid,nets,nete, fvm, omega_cn)
396 0 : call t_stopf('adjust_nsplit')
397 : end if
398 : end if
399 :
400 : ! ============================================================
401 : ! Print some diagnostic information
402 : ! ============================================================
403 0 : if (compute_diagnostics) then
404 0 : call prim_printstate(elem, tl, hybrid,nets,nete, fvm, omega_cn)
405 : end if
406 :
407 738816 : if (use_cslam.and.nsubstep==nsplit.and.nc.ne.fv_nphys) then
408 : !
409 : ! fill the fvm halo for mapping in d_p_coupling if
410 : ! physics grid resolution is different than fvm resolution
411 : !
412 0 : call fill_halo_fvm(ghostBufQnhc_h, elem,fvm,hybrid,nets,nete,nhc,1,nlev,nlev)
413 : end if
414 :
415 738816 : end subroutine prim_run_subcycle
416 :
417 :
418 2216448 : subroutine prim_step(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord, rstep, last_step,dt_remap)
419 : !
420 : ! Take qsplit dynamics steps and one tracer step
421 : ! for vertically lagrangian option, this subroutine does only the horizontal step
422 : !
423 : ! input:
424 : ! tl%nm1 not used
425 : ! tl%n0 data at time t
426 : ! tl%np1 new values at t+dt_q
427 : !
428 : ! then we update timelevel pointers:
429 : ! tl%nm1 = tl%n0
430 : ! tl%n0 = tl%np1
431 : ! so that:
432 : ! tl%nm1 tracers: t dynamics: t+(qsplit-1)*dt
433 : ! tl%n0 time t + dt_q
434 : !
435 738816 : use hybvcoord_mod, only: hvcoord_t
436 : use se_dyn_time_mod, only: TimeLevel_t, timelevel_update
437 : use control_mod, only: statefreq, qsplit, nu_p
438 : use thread_mod, only: omp_get_thread_num
439 : use prim_advance_mod, only: prim_advance_exp
440 : use prim_advection_mod, only: prim_advec_tracers_remap, prim_advec_tracers_fvm, deriv
441 : use derivative_mod, only: subcell_integration
442 : use hybrid_mod, only: set_region_num_threads, config_thread_region, get_loop_ranges
443 : use dimensions_mod, only: use_cslam,fvm_supercycling,fvm_supercycling_jet
444 : use dimensions_mod, only: kmin_jet, kmax_jet
445 : use fvm_mod, only: ghostBufQnhc_vh,ghostBufQ1_vh, ghostBufFlux_vh
446 : use fvm_mod, only: ghostBufQ1_h,ghostBufQnhcJet_h, ghostBufFluxJet_h
447 : use se_dyn_time_mod, only: timelevel_qdp
448 : use fvm_mapping, only: cslam2gll
449 : #ifdef waccm_debug
450 : use cam_history, only: outfld
451 : #endif
452 :
453 :
454 : type (element_t) , intent(inout) :: elem(:)
455 : type(fvm_struct), intent(inout) :: fvm(:)
456 : type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared)
457 : type (hvcoord_t), intent(in) :: hvcoord ! hybrid vertical coordinate struct
458 : integer, intent(in) :: nets ! starting thread element number (private)
459 : integer, intent(in) :: nete ! ending thread element number (private)
460 : real(kind=r8), intent(in) :: dt ! "timestep dependent" timestep
461 : type (TimeLevel_t), intent(inout) :: tl
462 : integer, intent(in) :: rstep ! vertical remap subcycling step
463 : logical, intent(in) :: last_step! last step before d_p_coupling
464 : real(kind=r8), intent(in) :: dt_remap
465 :
466 : type (hybrid_t):: hybridnew,hybridnew2
467 : real(kind=r8) :: st, st1, dp, dt_q
468 : integer :: ie,t,q,k,i,j,n, n_Q
469 : integer :: ithr
470 : integer :: region_num_threads
471 : integer :: kbeg,kend
472 : integer :: n0_qdp, np1_qdp
473 :
474 : real (kind=r8) :: tempdp3d(np,np), x
475 : real (kind=r8) :: tempmass(nc,nc)
476 : real (kind=r8) :: tempflux(nc,nc,4)
477 :
478 : real (kind=r8) :: dp_np1(np,np)
479 :
480 :
481 2216448 : dt_q = dt*qsplit
482 : ! ===============
483 : ! initialize mean flux accumulation variables and save some variables at n0
484 : ! for use by advection
485 : ! ===============
486 17800848 : do ie=nets,nete
487 62337600000 : elem(ie)%derived%vn0=0 ! mean horizontal mass flux
488 30451917600 : elem(ie)%derived%omega=0
489 15584400 : if (nu_p>0) then
490 30451917600 : elem(ie)%derived%dpdiss_ave=0
491 30451917600 : elem(ie)%derived%dpdiss_biharmonic=0
492 : endif
493 :
494 : ! dp at time t: use floating lagrangian levels:
495 30454134048 : elem(ie)%derived%dp(:,:,:)=elem(ie)%state%dp3d(:,:,:,tl%n0)
496 : enddo
497 :
498 : ! ===============
499 : ! Dynamical Step
500 : ! ===============
501 2216448 : n_Q = tl%n0 ! n_Q = timelevel of FV tracers at time t. need to save this
502 : ! SE tracers only carry 2 timelevels
503 :
504 2216448 : call t_startf('prim_advance_exp')
505 : ! ithr = 0 ! omp_get_thread_num()
506 : ! vybrid = hybrid_create(hybrid%par,ithr)
507 :
508 : call prim_advance_exp(elem, fvm, deriv, hvcoord, &
509 2216448 : hybrid, dt, tl, nets, nete)
510 :
511 2216448 : call t_stopf('prim_advance_exp')
512 :
513 2216448 : do n=2,qsplit
514 0 : call TimeLevel_update(tl,"leapfrog")
515 :
516 0 : call t_startf('prim_advance_exp')
517 :
518 : call prim_advance_exp(elem, fvm, deriv, hvcoord, &
519 0 : hybrid, dt, tl, nets, nete)
520 :
521 2216448 : call t_stopf('prim_advance_exp')
522 :
523 : ! defer final timelevel update until after Q update.
524 : enddo
525 : #ifdef HOMME_TEST_SUB_ELEMENT_MASS_FLUX
526 : if (use_cslam.and.rstep==1) then
527 : do ie=nets,nete
528 : do k=1,nlev
529 : tempdp3d = elem(ie)%state%dp3d(:,:,k,tl%np1) - &
530 : elem(ie)%derived%dp(:,:,k)
531 : call subcell_integration(tempdp3d, np, nc, elem(ie)%metdet,tempmass)
532 : tempflux = dt_q*elem(ie)%sub_elem_mass_flux(:,:,:,k)
533 : do i=1,nc
534 : do j=1,nc
535 : x = SUM(tempflux(i,j,:))
536 : if (ABS(tempmass(i,j)).lt.1e-11_r8 .and. 1e-11_r8.lt.ABS(x)) then
537 : write(iulog,*) __FILE__,__LINE__,"**CSLAM mass-flux ERROR***",ie,k,i,j,tempmass(i,j),x
538 : call endrun('**CSLAM mass-flux ERROR***')
539 : elseif (1e-5_r8.lt.ABS((tempmass(i,j)-x)/tempmass(i,j))) then
540 : write(iulog,*) __FILE__,__LINE__,"**CSLAM mass-flux ERROR**",ie,k,i,j,tempmass(i,j),x,&
541 : ABS((tempmass(i,j)-x)/tempmass(i,j))
542 : call endrun('**CSLAM mass-flux ERROR**')
543 : endif
544 : end do
545 : end do
546 : end do
547 : end do
548 : end if
549 : #endif
550 : ! current dynamics state variables:
551 : ! derived%dp = dp at start of timestep
552 : ! derived%vn0 = mean horiz. flux: U*dp
553 : ! rsplit>0
554 : ! state%v(:,:,:,np1) = velocity on lagrangian levels
555 : ! state%dp3d(:,:,:,np1) = dp3d
556 : !
557 :
558 :
559 : ! ===============
560 : ! Tracer Advection.
561 : ! in addition, this routine will apply the DSS to:
562 : ! derived%omega =
563 : ! Tracers are always vertically lagrangian.
564 : ! ===============
565 : ! Advect tracers if their count is > 0.
566 : ! special case in CAM: if CSLAM tracers are turned on , qsize=1 but this tracer should
567 : ! not be advected. This will be cleaned up when the physgrid is merged into CAM trunk
568 : ! Currently advecting all species
569 2216448 : if (.not.use_cslam) then
570 0 : call t_startf('prim_advec_tracers_remap')
571 0 : region_num_threads=tracer_num_threads
572 0 : call omp_set_nested(.true.)
573 : !$OMP PARALLEL NUM_THREADS(region_num_threads), DEFAULT(SHARED), PRIVATE(hybridnew)
574 0 : hybridnew = config_thread_region(hybrid,'tracer')
575 0 : call Prim_Advec_Tracers_remap(elem, deriv,hvcoord,hybridnew,dt_q,tl,nets,nete)
576 : !$OMP END PARALLEL
577 0 : call omp_set_nested(.false.)
578 0 : call t_stopf('prim_advec_tracers_remap')
579 : else
580 : !
581 : ! only run fvm transport every fvm_supercycling rstep
582 : !
583 : ! FVM transport
584 : !
585 2216448 : if ((mod(rstep,fvm_supercycling) == 0).and.(mod(rstep,fvm_supercycling_jet) == 0)) then
586 :
587 : ! call omp_set_nested(.true.)
588 : ! !$OMP PARALLEL NUM_THREADS(vert_num_threads), DEFAULT(SHARED), PRIVATE(hybridnew2,kbeg,kend)
589 : ! hybridnew2 = config_thread_region(hybrid,'vertical')
590 : ! call get_loop_ranges(hybridnew2,kbeg=kbeg,kend=kend)
591 : call Prim_Advec_Tracers_fvm(elem,fvm,hvcoord,hybrid,&
592 738816 : dt_q,tl,nets,nete,ghostBufQnhc_vh,ghostBufQ1_vh, ghostBufFlux_vh,1,nlev)
593 : ! !$OMP END PARALLEL
594 : ! call omp_set_nested(.false.)
595 : !
596 : ! to avoid accumulation of truncation error overwrite CSLAM surface pressure with SE
597 : ! surface pressure
598 : !
599 5933616 : do ie=nets,nete
600 : !
601 : ! overwrite PSDRY on CSLAM grid with SE PSDRY integrated over CSLAM control volume
602 : !
603 : ! call subcell_integration(elem(ie)%state%psdry(:,:), np, nc, elem(ie)%metdet,fvm(ie)%psc)
604 : ! fvm(ie)%psc = fvm(ie)%psc*fvm(ie)%inv_se_area_sphere
605 : !
606 : ! Update CSLAM surface pressure
607 : !
608 21518016 : do j=1,nc
609 67532400 : do i=1,nc
610 4410385200 : fvm(ie)%psc(i,j) = sum(fvm(ie)%dp_fvm(i,j,:)) + hvcoord%hyai(1)*hvcoord%ps0
611 : end do
612 : end do
613 : end do
614 738816 : call TimeLevel_Qdp( tl, qsplit, n0_qdp, np1_qdp)
615 738816 : if (.not.last_step) call cslam2gll(elem, fvm, hybrid,nets,nete, tl%np1, np1_qdp)
616 1477632 : else if ((mod(rstep,fvm_supercycling_jet) == 0)) then
617 : !
618 : ! shorter fvm time-step in jet region
619 : !
620 : call Prim_Advec_Tracers_fvm(elem,fvm,hvcoord,hybrid,&
621 0 : dt_q,tl,nets,nete,ghostBufQnhcJet_h,ghostBufQ1_h, ghostBufFluxJet_h,kmin_jet,kmax_jet)
622 : end if
623 :
624 : #ifdef waccm_debug
625 : do ie=nets,nete
626 : call outfld('CSLAM_gamma', RESHAPE(fvm(ie)%CSLAM_gamma(:,:,:,1), &
627 : (/nc*nc,nlev/)), nc*nc, ie)
628 : end do
629 : #endif
630 : endif
631 :
632 2216448 : end subroutine prim_step
633 0 : subroutine prim_step_scm(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord, rstep)
634 : !
635 : ! prim_step version for single column model (SCM)
636 : ! Here we simply want to compute the floating level tendency
637 : ! based on the prescribed large scale vertical velocity
638 : ! Take qsplit dynamics steps and one tracer step
639 : ! for vertically lagrangian option, this subroutine does only
640 : ! the horizontal step
641 : !
642 : ! input:
643 : ! tl%nm1 not used
644 : ! tl%n0 data at time t
645 : ! tl%np1 new values at t+dt_q
646 : !
647 : ! then we update timelevel pointers:
648 : ! tl%nm1 = tl%n0
649 : ! tl%n0 = tl%np1
650 : ! so that:
651 : ! tl%nm1 tracers: t dynamics: t+(qsplit-1)*dt
652 : ! tl%n0 time t + dt_q
653 : !
654 2216448 : use hybvcoord_mod, only: hvcoord_t
655 : use se_dyn_time_mod, only: TimeLevel_t, timelevel_update
656 : use control_mod, only: statefreq, qsplit, nu_p
657 : use prim_advection_mod, only: deriv
658 : use hybrid_mod, only: config_thread_region, get_loop_ranges
659 :
660 : type (element_t) , intent(inout) :: elem(:)
661 : type(fvm_struct), intent(inout) :: fvm(:)
662 : type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared)
663 : type (hvcoord_t), intent(in) :: hvcoord ! hybrid vertical coordinate struct
664 : integer, intent(in) :: nets ! starting thread element number (private)
665 : integer, intent(in) :: nete ! ending thread element number (private)
666 : real(kind=r8), intent(in) :: dt ! "timestep dependent" timestep
667 : type (TimeLevel_t), intent(inout) :: tl
668 : integer, intent(in) :: rstep ! vertical remap subcycling step
669 :
670 : integer :: ie,n
671 :
672 : ! ===============
673 : ! initialize mean flux accumulation variables and save some variables at n0
674 : ! for use by advection
675 : ! ===============
676 0 : do ie=nets,nete
677 0 : elem(ie)%derived%vn0=0 ! mean horizontal mass flux
678 0 : if (nu_p>0) then
679 0 : elem(ie)%derived%dpdiss_ave=0
680 0 : elem(ie)%derived%dpdiss_biharmonic=0
681 : endif
682 0 : elem(ie)%derived%dp(:,:,:)=elem(ie)%state%dp3d(:,:,:,tl%n0)
683 : enddo
684 :
685 : ! ===============
686 : ! Dynamical Step
687 : ! ===============
688 :
689 0 : call t_startf('set_prescribed_scm')
690 :
691 : call set_prescribed_scm(elem, fvm, deriv, hvcoord, &
692 0 : hybrid, dt, tl, nets, nete)
693 :
694 0 : call t_stopf('set_prescribed_scm')
695 :
696 0 : do n=2,qsplit
697 0 : call TimeLevel_update(tl,"leapfrog")
698 :
699 0 : call t_startf('set_prescribed_scm')
700 :
701 : call set_prescribed_scm(elem, fvm, deriv, hvcoord, &
702 0 : hybrid, dt, tl, nets, nete)
703 :
704 0 : call t_stopf('set_prescribed_scm')
705 : enddo
706 :
707 0 : end subroutine prim_step_scm
708 : !=======================================================================================================!
709 :
710 :
711 0 : subroutine prim_finalize(hybrid)
712 : type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared)
713 :
714 : ! ==========================
715 : ! end of the hybrid program
716 : ! ==========================
717 0 : end subroutine prim_finalize
718 :
719 : !=========================================================================================
720 :
721 768 : subroutine prim_set_dry_mass(elem, hvcoord,initial_global_ave_dry_ps,q)
722 : use element_mod, only: element_t
723 : use hybvcoord_mod , only: hvcoord_t
724 : use dimensions_mod, only: nelemd, nlev, np
725 : use constituents, only: cnst_type, qmin, pcnst
726 : use cam_logfile, only: iulog
727 : use spmd_utils, only: masterproc
728 :
729 : type (element_t) , intent(inout):: elem(:)
730 : type (hvcoord_t) , intent(in) :: hvcoord
731 : real (kind=r8), intent(in) :: initial_global_ave_dry_ps
732 : real (kind=r8), intent(inout):: q(np,np,nlev,nelemd,pcnst)
733 :
734 : ! local
735 : real (kind=r8) :: global_ave_ps_inic,dp_tmp, factor(np,np,nlev)
736 : integer :: ie, i, j ,k, m_cnst
737 :
738 768 : if (initial_global_ave_dry_ps == 0) return;
739 :
740 768 : call get_global_ave_surface_pressure(elem, global_ave_ps_inic)
741 :
742 6168 : do ie=1,nelemd
743 113400 : elem(ie)%state%psdry(:,:)=elem(ie)%state%psdry(:,:)*(initial_global_ave_dry_ps/global_ave_ps_inic)
744 507600 : do k=1,nlev
745 2516400 : do j = 1,np
746 10546200 : do i = 1,np
747 16070400 : dp_tmp = ((hvcoord%hyai(k+1) - hvcoord%hyai(k))*hvcoord%ps0)+&
748 24105600 : ((hvcoord%hybi(k+1) - hvcoord%hybi(k))*elem(ie)%state%psdry(i,j))
749 8035200 : factor(i,j,k) = elem(ie)%state%dp3d(i,j,k,1)/dp_tmp
750 34149600 : elem(ie)%state%dp3d(i,j,k,:) = dp_tmp
751 : end do
752 : end do
753 : end do
754 : !
755 : ! conserve initial condition mass of 'wet' tracers (following dryairm.F90 for FV dycore)
756 : ! and conserve mixing ratio (not mass) of 'dry' tracers
757 : !
758 65568 : do m_cnst=1,pcnst
759 64800 : if (cnst_type(m_cnst).ne.'dry') then
760 5583600 : do k=1,nlev
761 27680400 : do j = 1,np
762 116008200 : do i = 1,np
763 88387200 : q(i,j,k,ie,m_cnst) = q(i,j,k,ie,m_cnst)*factor(i,j,k)
764 110484000 : q(i,j,k,ie,m_cnst) = max(qmin(m_cnst),q(i,j,k,ie,m_cnst))
765 : end do
766 : end do
767 : end do
768 : end if
769 : end do
770 : end do
771 768 : if (masterproc) then
772 1 : write (iulog,*) "------ info from prim_set_dry_mass -----------------------------------------------------------"
773 1 : write (iulog,*) "Scaling dry surface pressure to global average of = ",&
774 2 : initial_global_ave_dry_ps/100.0_r8,"hPa"
775 1 : write (iulog,*) "Average dry surface pressure in initial condition = ",global_ave_ps_inic/100.0_r8,"hPa"
776 1 : write (iulog,*) "Average dry surface pressure change = ",&
777 2 : initial_global_ave_dry_ps-global_ave_ps_inic,"Pa"
778 1 : write (iulog,*) "Mixing ratios that are wet have been scaled so that total mass of tracer is conserved"
779 1 : write (iulog,*) "Mixing ratios that are dry have not been changed (mass not conserved in scaling process)"
780 1 : write (iulog,*) "------ end info from prim_set_dry_mass -------------------------------------------------------"
781 : endif
782 : end subroutine prim_set_dry_mass
783 :
784 768 : subroutine get_global_ave_surface_pressure(elem, global_ave_ps_inic)
785 : use element_mod , only : element_t
786 : use dimensions_mod , only : np
787 : use global_norms_mod , only : global_integral
788 : use hybrid_mod , only : config_thread_region, get_loop_ranges, hybrid_t
789 : use parallel_mod , only : par
790 :
791 : type (element_t) , intent(in) :: elem(:)
792 : real (kind=r8), intent(out) :: global_ave_ps_inic
793 :
794 : ! local
795 768 : real (kind=r8), allocatable :: tmp(:,:,:)
796 : type (hybrid_t) :: hybrid
797 : integer :: ie, nets, nete
798 :
799 : !JMD $OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n)
800 : !JMD hybrid = config_thread_region(par,'horizontal')
801 768 : hybrid = config_thread_region(par,'serial')
802 768 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
803 2304 : allocate(tmp(np,np,nets:nete))
804 :
805 6168 : do ie=nets,nete
806 114168 : tmp(:,:,ie)=elem(ie)%state%psdry(:,:)
807 : enddo
808 768 : global_ave_ps_inic = global_integral(elem, tmp(:,:,nets:nete),hybrid,np,nets,nete)
809 768 : deallocate(tmp)
810 768 : end subroutine get_global_ave_surface_pressure
811 :
812 0 : subroutine set_prescribed_scm(elem, fvm, deriv, hvcoord, &
813 : hybrid, dt, tl, nets, nete)
814 : use control_mod, only: tstep_type, qsplit
815 : use derivative_mod, only: derivative_t
816 : use dimensions_mod, only: np, nlev
817 : use element_mod, only: element_t
818 : use hybvcoord_mod, only: hvcoord_t
819 : use hybrid_mod, only: hybrid_t
820 : use se_dyn_time_mod, only: TimeLevel_t, timelevel_qdp
821 : use fvm_control_volume_mod, only: fvm_struct
822 : implicit none
823 :
824 : type (element_t), intent(inout), target :: elem(:)
825 : type(fvm_struct) , intent(inout) :: fvm(:)
826 : type (derivative_t) , intent(in) :: deriv
827 : type (hvcoord_t) :: hvcoord
828 : type (hybrid_t) , intent(in) :: hybrid
829 : real (kind=r8), intent(in) :: dt
830 : type (TimeLevel_t) , intent(in) :: tl
831 : integer , intent(in) :: nets
832 : integer , intent(in) :: nete
833 :
834 : ! Local
835 : integer :: ie,nm1,n0,np1,k,qn0,qnp1,p
836 : real(kind=r8) :: eta_dot_dpdn(np,np,nlev+1)
837 :
838 :
839 0 : nm1 = tl%nm1
840 0 : n0 = tl%n0
841 0 : np1 = tl%np1
842 :
843 0 : call TimeLevel_Qdp(tl, qsplit, qn0, qnp1) ! compute current Qdp() timelevel
844 :
845 0 : do ie=nets,nete
846 0 : do k=1,nlev
847 0 : eta_dot_dpdn(:,:,k)=elem(ie)%derived%omega(:,:,k)
848 : enddo
849 0 : eta_dot_dpdn(:,:,nlev+1) = eta_dot_dpdn(:,:,nlev)
850 :
851 0 : do k=1,nlev
852 0 : elem(ie)%state%dp3d(:,:,k,np1) = elem(ie)%state%dp3d(:,:,k,n0) &
853 0 : + dt*(eta_dot_dpdn(:,:,k+1) - eta_dot_dpdn(:,:,k))
854 : enddo
855 :
856 0 : do k=1,nlev
857 0 : elem(ie)%state%T(:,:,k,np1) = elem(ie)%state%T(:,:,k,n0)
858 : enddo
859 :
860 0 : do p=1,qsize
861 0 : do k=1,nlev
862 0 : elem(ie)%state%Qdp(:,:,k,p,qnp1) = elem(ie)%state%Qdp(:,:,k,p,qn0) &
863 0 : + elem(ie)%state%Qdp(:,:,k,p,qn0)/elem(ie)%state%dp3d(:,:,k,n0) * &
864 0 : dt*(eta_dot_dpdn(:,:,k+1) - eta_dot_dpdn(:,:,k))
865 : enddo
866 : enddo
867 : enddo
868 0 : end subroutine set_prescribed_scm
869 :
870 : end module prim_driver_mod
|
