Line data Source code
1 : #define OVERLAP 1
2 : module prim_advection_mod
3 : !
4 : ! two formulations. both are conservative
5 : ! u grad Q formulation:
6 : !
7 : ! d/dt[ Q] + U grad Q = 0
8 : !
9 : ! d/dt[ dp/dn ] = div( dp/dn U )
10 : !
11 : ! total divergence formulation:
12 : ! d/dt[dp/dn Q] + div( U dp/dn Q ) = 0
13 : !
14 : ! for convience, rewrite this as dp Q: (since dn does not depend on time or the horizonal):
15 : ! equation is now:
16 : ! d/dt[dp Q] + div( U dp Q ) = 0
17 : !
18 : !
19 : use shr_kind_mod, only: r8=>shr_kind_r8
20 : use dimensions_mod, only: nlev, np, qsize, nc
21 : use physconst, only: cpair
22 : use derivative_mod, only: derivative_t
23 : use element_mod, only: element_t
24 : use fvm_control_volume_mod, only: fvm_struct
25 : use hybvcoord_mod, only: hvcoord_t
26 : use se_dyn_time_mod, only: TimeLevel_t, TimeLevel_Qdp
27 : use control_mod, only: nu_q, nu_p, limiter_option, hypervis_subcycle_q, rsplit
28 : use edge_mod, only: edgevpack, edgevunpack, initedgebuffer, initedgesbuffer
29 :
30 : use edgetype_mod, only: EdgeBuffer_t
31 : use hybrid_mod, only: hybrid_t
32 : use viscosity_mod, only: biharmonic_wk_scalar, neighbor_minmax, &
33 : neighbor_minmax_start, neighbor_minmax_finish
34 : use perf_mod, only: t_startf, t_stopf, t_barrierf
35 : use cam_abortutils, only: endrun
36 : use thread_mod, only: horz_num_threads, vert_num_threads, tracer_num_threads
37 :
38 : implicit none
39 :
40 : private
41 : save
42 :
43 : public :: Prim_Advec_Init1, Prim_Advec_Init2
44 : public :: Prim_Advec_Tracers_remap
45 : public :: prim_advec_tracers_fvm
46 : public :: vertical_remap
47 :
48 : type (EdgeBuffer_t) :: edgeAdv, edgeAdvp1, edgeAdvQminmax, edgeveloc
49 :
50 : integer,parameter :: DSSeta = 1
51 : integer,parameter :: DSSomega = 2
52 : integer,parameter :: DSSdiv_vdp_ave = 3
53 : integer,parameter :: DSSno_var = -1
54 :
55 : real(kind=r8), allocatable :: qmin(:,:,:), qmax(:,:,:)
56 :
57 : !JMD I don't see why this needs to be thread private.
58 : !JMD type (derivative_t), public, allocatable :: deriv(:) ! derivative struct (nthreads)
59 : type (derivative_t), public :: deriv
60 :
61 :
62 : contains
63 :
64 :
65 1536 : subroutine Prim_Advec_Init1(par, elem)
66 : use dimensions_mod, only: nlev, qsize, nelemd,ntrac,use_cslam
67 : use parallel_mod, only: parallel_t, boundaryCommMethod
68 : type(parallel_t) :: par
69 : type (element_t) :: elem(:)
70 : !
71 : ! Shared buffer pointers.
72 : ! Using "=> null()" in a subroutine is usually bad, because it makes
73 : ! the variable have an implicit "save", and therefore shared between
74 : ! threads. But in this case we want shared pointers.
75 : real(kind=r8), pointer :: buf_ptr(:) => null()
76 : real(kind=r8), pointer :: receive_ptr(:) => null()
77 : integer :: advec_remap_num_threads
78 :
79 :
80 : !
81 : ! Set the number of threads used in the subroutine Prim_Advec_tracers_remap()
82 : !
83 1536 : if (use_cslam) then
84 : advec_remap_num_threads = 1
85 : else
86 0 : advec_remap_num_threads = tracer_num_threads
87 : endif
88 : ! this might be called with qsize=0
89 : ! allocate largest one first
90 : ! Currently this is never freed. If it was, only this first one should
91 : ! be freed, as only it knows the true size of the buffer.
92 1536 : if (.not.use_cslam) then
93 : call initEdgeBuffer(par,edgeAdvp1,elem,qsize*nlev + nlev,bndry_type=boundaryCommMethod,&
94 0 : nthreads=horz_num_threads*advec_remap_num_threads)
95 : call initEdgeBuffer(par,edgeAdv,elem,qsize*nlev,bndry_type=boundaryCommMethod, &
96 0 : nthreads=horz_num_threads*advec_remap_num_threads)
97 : ! This is a different type of buffer pointer allocation
98 : ! used for determine the minimum and maximum value from
99 : ! neighboring elements
100 : call initEdgeSBuffer(par,edgeAdvQminmax,elem,qsize*nlev*2,bndry_type=boundaryCommMethod, &
101 0 : nthreads=horz_num_threads*advec_remap_num_threads)
102 : end if
103 1536 : call initEdgeBuffer(par,edgeveloc,elem,2*nlev,bndry_type=boundaryCommMethod)
104 :
105 :
106 : ! Don't actually want these saved, if this is ever called twice.
107 1536 : nullify(buf_ptr)
108 1536 : nullify(receive_ptr)
109 :
110 :
111 : ! this static array is shared by all threads, so dimension for all threads (nelemd), not nets:nete:
112 6144 : allocate (qmin(nlev,qsize,nelemd))
113 4608 : allocate (qmax(nlev,qsize,nelemd))
114 :
115 1536 : end subroutine Prim_Advec_Init1
116 :
117 1536 : subroutine Prim_Advec_Init2(fvm_corners, fvm_points)
118 : use dimensions_mod, only : nc
119 : use derivative_mod, only : derivinit
120 :
121 : real(kind=r8), intent(in) :: fvm_corners(nc+1)
122 : real(kind=r8), intent(in) :: fvm_points(nc)
123 :
124 : ! ==================================
125 : ! Initialize derivative structure
126 : ! ==================================
127 1536 : call derivinit(deriv,fvm_corners, fvm_points)
128 1536 : end subroutine Prim_Advec_Init2
129 :
130 : !
131 : ! fvm driver
132 : !
133 738816 : subroutine Prim_Advec_Tracers_fvm(elem,fvm,hvcoord,hybrid,&
134 : dt,tl,nets,nete,ghostbufQnhc,ghostBufQ1, ghostBufFlux,kmin,kmax)
135 1536 : use fvm_consistent_se_cslam, only: run_consistent_se_cslam
136 : use edgetype_mod, only: edgebuffer_t
137 : implicit none
138 : type (element_t), intent(inout) :: elem(:)
139 : type (fvm_struct), intent(inout) :: fvm(:)
140 : type (hvcoord_t) :: hvcoord
141 : type (hybrid_t), intent(in):: hybrid
142 : type (TimeLevel_t) :: tl
143 :
144 : real(kind=r8) , intent(in) :: dt
145 : integer,intent(in) :: nets,nete,kmin,kmax
146 : type (EdgeBuffer_t), intent(inout):: ghostbufQnhc,ghostBufQ1, ghostBufFlux
147 :
148 738816 : call t_barrierf('sync_prim_advec_tracers_fvm', hybrid%par%comm)
149 738816 : call t_startf('prim_advec_tracers_fvm')
150 :
151 738816 : if (rsplit==0) call endrun('cslam only works for rsplit>0')
152 :
153 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
154 : ! 2D advection step
155 : !
156 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
157 : call run_consistent_se_cslam(elem,fvm,hybrid,dt,tl,nets,nete,hvcoord,&
158 738816 : ghostbufQnhc,ghostBufQ1, ghostBufFlux,kmin,kmax)
159 738816 : call t_stopf('prim_advec_tracers_fvm')
160 738816 : end subroutine Prim_Advec_Tracers_fvm
161 :
162 :
163 :
164 : !=================================================================================================!
165 :
166 0 : subroutine Prim_Advec_Tracers_remap( elem , deriv , hvcoord , hybrid , dt , tl , nets , nete )
167 : implicit none
168 : type (element_t) , intent(inout) :: elem(:)
169 : type (derivative_t) , intent(in ) :: deriv
170 : type (hvcoord_t) , intent(in ) :: hvcoord
171 : type (hybrid_t) , intent(in ) :: hybrid
172 : real(kind=r8) , intent(in ) :: dt
173 : type (TimeLevel_t) , intent(inout) :: tl
174 : integer , intent(in ) :: nets
175 : integer , intent(in ) :: nete
176 :
177 :
178 : !print *,'prim_Advec_Tracers_remap: qsize: ',qsize
179 0 : call Prim_Advec_Tracers_remap_rk2( elem , deriv , hvcoord , hybrid , dt , tl , nets , nete )
180 738816 : end subroutine Prim_Advec_Tracers_remap
181 :
182 :
183 0 : subroutine euler_step_driver(np1_qdp , n0_qdp , dt , elem , hvcoord , hybrid , deriv , nets , nete , DSSopt , rhs_multiplier )
184 :
185 :
186 : integer , intent(in ) :: np1_qdp, n0_qdp
187 : real (kind=r8), intent(in ) :: dt
188 : type (element_t) , intent(inout) :: elem(:)
189 : type (hvcoord_t) , intent(in ) :: hvcoord
190 : type (hybrid_t) , intent(in ) :: hybrid
191 : type (derivative_t) , intent(in ) :: deriv
192 : integer , intent(in ) :: nets
193 : integer , intent(in ) :: nete
194 : integer , intent(in ) :: DSSopt
195 : integer , intent(in ) :: rhs_multiplier
196 :
197 0 : call euler_step( np1_qdp , n0_qdp , dt , elem , hvcoord , hybrid , deriv , nets , nete , DSSopt , rhs_multiplier)
198 :
199 0 : end subroutine euler_step_driver
200 :
201 : !-----------------------------------------------------------------------------
202 : !-----------------------------------------------------------------------------
203 : ! forward-in-time 2 level vertically lagrangian step
204 : ! this code takes a lagrangian step in the horizontal
205 : ! (complete with DSS), and then applies a vertical remap
206 : !
207 : ! This routine may use dynamics fields at timelevel np1
208 : ! In addition, other fields are required, which have to be
209 : ! explicitly saved by the dynamics: (in elem(ie)%derived struct)
210 : !
211 : ! Fields required from dynamics: (in
212 : ! omega it will be DSS'd here, for later use by CAM physics
213 : ! we DSS omega here because it can be done for "free"
214 : ! Consistent mass/tracer-mass advection (used if subcycling turned on)
215 : ! dp() dp at timelevel n0
216 : ! vn0() mean flux < U dp > going from n0 to np1
217 : !
218 : ! 3 stage
219 : ! Euler step from t -> t+.5
220 : ! Euler step from t+.5 -> t+1.0
221 : ! Euler step from t+1.0 -> t+1.5
222 : ! u(t) = u(t)/3 + u(t+2)*2/3
223 : !
224 : !-----------------------------------------------------------------------------
225 : !-----------------------------------------------------------------------------
226 0 : subroutine Prim_Advec_Tracers_remap_rk2( elem , deriv , hvcoord , hybrid , dt , tl , nets , nete )
227 : use derivative_mod, only: divergence_sphere
228 : use control_mod , only: qsplit
229 : use hybrid_mod , only: get_loop_ranges!, PrintHybrid
230 : ! use thread_mod , only : omp_set_num_threads, omp_get_thread_num
231 :
232 : type (element_t) , intent(inout) :: elem(:)
233 : type (derivative_t) , intent(in ) :: deriv
234 : type (hvcoord_t) , intent(in ) :: hvcoord
235 : type (hybrid_t) , intent(in ) :: hybrid
236 : real(kind=r8) , intent(in ) :: dt
237 : type (TimeLevel_t) , intent(inout) :: tl
238 : integer , intent(in ) :: nets
239 : integer , intent(in ) :: nete
240 :
241 : real (kind=r8), dimension(np,np,2 ) :: gradQ
242 : integer :: k,ie
243 : integer :: rkstage,rhs_multiplier
244 : integer :: n0_qdp, np1_qdp
245 : integer :: kbeg,kend,qbeg,qend
246 :
247 : ! call t_barrierf('sync_prim_advec_tracers_remap_k2', hybrid%par%comm)
248 : ! call t_startf('prim_advec_tracers_remap_rk2')
249 : ! call extrae_user_function(1)
250 0 : call TimeLevel_Qdp( tl, qsplit, n0_qdp, np1_qdp) !time levels for qdp are not the same
251 0 : rkstage = 3 ! 3 stage RKSSP scheme, with optimal SSP CFL
252 :
253 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
254 : ! RK2 2D advection step
255 : ! note: stage 3 we take the oppertunity to DSS omega
256 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
257 : ! use these for consistent advection (preserve Q=1)
258 : ! derived%vdp_ave = mean horiz. flux: U*dp
259 : ! derived%omega = advection code will DSS this for the physics, but otherwise
260 : ! it is not needed
261 : ! Also: save a copy of div(U dp) in derived%div(:,:,:,1), which will be DSS'd
262 : ! and a DSS'ed version stored in derived%div(:,:,:,2)
263 :
264 0 : call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
265 :
266 0 : do ie=nets,nete
267 0 : do k=kbeg,kend
268 : ! div( U dp Q),
269 0 : gradQ(:,:,1)=elem(ie)%derived%vn0(:,:,1,k)
270 0 : gradQ(:,:,2)=elem(ie)%derived%vn0(:,:,2,k)
271 : ! elem(ie)%derived%divdp(:,:,k) = divergence_sphere(gradQ,deriv,elem(ie))
272 0 : call divergence_sphere(gradQ,deriv,elem(ie),elem(ie)%derived%divdp(:,:,k))
273 0 : elem(ie)%derived%divdp_proj(:,:,k) = elem(ie)%derived%divdp(:,:,k)
274 : enddo
275 : enddo
276 :
277 :
278 : !rhs_multiplier is for obtaining dp_tracers at each stage:
279 : !dp_tracers(stage) = dp - rhs_multiplier*dt*divdp_proj
280 : ! call t_startf('euler_step')
281 :
282 0 : rhs_multiplier = 0
283 0 : call euler_step_driver( np1_qdp, n0_qdp , dt/2, elem, hvcoord, hybrid, deriv, nets, nete, DSSdiv_vdp_ave, rhs_multiplier )
284 :
285 0 : rhs_multiplier = 1
286 0 : call euler_step_driver( np1_qdp, np1_qdp, dt/2, elem, hvcoord, hybrid, deriv, nets, nete, DSSno_var , rhs_multiplier )
287 :
288 0 : rhs_multiplier = 2
289 0 : call euler_step_driver( np1_qdp, np1_qdp, dt/2, elem, hvcoord, hybrid, deriv, nets, nete, DSSomega , rhs_multiplier )
290 :
291 : ! call t_stopf ('euler_step')
292 :
293 : !to finish the 2D advection step, we need to average the t and t+2 results to get a second order estimate for t+1.
294 0 : call qdp_time_avg( elem , rkstage , n0_qdp , np1_qdp , hybrid, nets , nete )
295 :
296 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
297 : ! Dissipation
298 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
299 0 : if ( limiter_option == 8 ) then
300 : ! dissipation was applied in RHS.
301 : else
302 0 : call advance_hypervis_scalar(edgeadv,elem,hvcoord,hybrid,deriv,tl%np1,np1_qdp,nets,nete,dt)
303 : endif
304 : ! call extrae_user_function(0)
305 :
306 : ! call t_stopf('prim_advec_tracers_remap_rk2')
307 :
308 0 : end subroutine prim_advec_tracers_remap_rk2
309 :
310 : !-----------------------------------------------------------------------------
311 : !-----------------------------------------------------------------------------
312 :
313 0 : subroutine qdp_time_avg( elem , rkstage , n0_qdp , np1_qdp , hybrid , nets , nete )
314 0 : use hybrid_mod, only : hybrid_t, get_loop_ranges
315 : implicit none
316 : type(element_t) , intent(inout) :: elem(:)
317 : integer , intent(in ) :: rkstage , n0_qdp , np1_qdp , nets , nete
318 : type(hybrid_t) :: hybrid
319 : integer :: i,j,ie,q,k
320 : integer :: kbeg,kend,qbeg,qend
321 : real(kind=r8) :: rrkstage
322 :
323 0 : call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
324 :
325 0 : rrkstage=1.0_r8/real(rkstage,kind=r8)
326 0 : do ie=nets,nete
327 0 : do q=qbeg,qend
328 0 : do k=kbeg,kend
329 : !OMP_COLLAPSE_SIMD
330 : !DIR_VECTOR_ALIGNED
331 0 : do j=1,np
332 0 : do i=1,np
333 0 : elem(ie)%state%Qdp(i,j,k,q,np1_qdp) = &
334 0 : rrkstage *( elem(ie)%state%Qdp(i,j,k,q,n0_qdp) + &
335 0 : (rkstage-1)*elem(ie)%state%Qdp(i,j,k,q,np1_qdp) )
336 : enddo
337 : enddo
338 : enddo
339 : enddo
340 : enddo
341 0 : end subroutine qdp_time_avg
342 :
343 : !-----------------------------------------------------------------------------
344 : !-----------------------------------------------------------------------------
345 :
346 0 : subroutine euler_step( np1_qdp , n0_qdp , dt , elem , hvcoord , hybrid , deriv , nets , nete , DSSopt , rhs_multiplier )
347 : ! ===================================
348 : ! This routine is the basic foward
349 : ! euler component used to construct RK SSP methods
350 : !
351 : ! u(np1) = u(n0) + dt2*DSS[ RHS(u(n0)) ]
352 : !
353 : ! n0 can be the same as np1.
354 : !
355 : ! DSSopt = DSSeta or DSSomega: also DSS omega
356 : !
357 : ! ===================================
358 : use dimensions_mod , only: np, nlev
359 : use hybrid_mod , only: hybrid_t!, PrintHybrid
360 : use hybrid_mod , only: get_loop_ranges, threadOwnsTracer
361 : use element_mod , only: element_t
362 : use derivative_mod , only: derivative_t, divergence_sphere, limiter_optim_iter_full
363 : use edge_mod , only: edgevpack, edgevunpack
364 : use bndry_mod , only: bndry_exchange
365 : use hybvcoord_mod , only: hvcoord_t
366 :
367 : integer , intent(in ) :: np1_qdp, n0_qdp
368 : real (kind=r8), intent(in ) :: dt
369 : type (element_t) , intent(inout), target :: elem(:)
370 : type (hvcoord_t) , intent(in ) :: hvcoord
371 : type (hybrid_t) , intent(in ) :: hybrid
372 : type (derivative_t) , intent(in ) :: deriv
373 : integer , intent(in ) :: nets
374 : integer , intent(in ) :: nete
375 : integer , intent(in ) :: DSSopt
376 : integer , intent(in ) :: rhs_multiplier
377 :
378 : ! local
379 : real(kind=r8), dimension(np,np ) :: dpdiss
380 : real(kind=r8), dimension(np,np,nlev) :: dpdissk
381 : real(kind=r8), dimension(np,np,2 ) :: gradQ
382 : real(kind=r8), dimension(np,np,2,nlev ) :: Vstar
383 : real(kind=r8), dimension(np,np ,nlev ) :: Qtens
384 : real(kind=r8), dimension(np,np ,nlev ) :: dp
385 0 : real(kind=r8), dimension(np,np ,nlev,qsize,nets:nete) :: Qtens_biharmonic
386 : real(kind=r8), dimension(np,np) :: div
387 0 : real(kind=r8), pointer, dimension(:,:,:) :: DSSvar
388 : real(kind=r8) :: dp0(nlev)
389 : integer :: ie,q,i,j,k, kptr
390 : integer :: rhs_viss = 0
391 : integer :: kblk,qblk ! The per thead size of the vertical and tracers
392 : integer :: kbeg, kend, qbeg, qend
393 :
394 0 : call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
395 :
396 0 : kblk = kend - kbeg + 1 ! calculate size of the block of vertical levels
397 0 : qblk = qend - qbeg + 1 ! calculate size of the block of tracers
398 :
399 0 : do k = kbeg, kend
400 0 : dp0(k) = ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + &
401 0 : ( hvcoord%hybi(k+1) - hvcoord%hybi(k) )*hvcoord%ps0
402 : enddo
403 :
404 : ! call t_startf('euler_step')
405 :
406 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
407 : ! compute Q min/max values for lim8
408 : ! compute biharmonic mixing term f
409 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
410 0 : rhs_viss = 0
411 0 : if ( limiter_option == 8 ) then
412 : ! when running lim8, we also need to limit the biharmonic, so that term needs
413 : ! to be included in each euler step. three possible algorithms here:
414 : ! 1) most expensive:
415 : ! compute biharmonic (which also computes qmin/qmax) during all 3 stages
416 : ! be sure to set rhs_viss=1
417 : ! cost: 3 biharmonic steps with 3 DSS
418 : !
419 : ! 2) cheapest:
420 : ! compute biharmonic (which also computes qmin/qmax) only on first stage
421 : ! be sure to set rhs_viss=3
422 : ! reuse qmin/qmax for all following stages (but update based on local qmin/qmax)
423 : ! cost: 1 biharmonic steps with 1 DSS
424 : ! main concern: viscosity
425 : !
426 : ! 3) compromise:
427 : ! compute biharmonic (which also computes qmin/qmax) only on last stage
428 : ! be sure to set rhs_viss=3
429 : ! compute qmin/qmax directly on first stage
430 : ! reuse qmin/qmax for 2nd stage stage (but update based on local qmin/qmax)
431 : ! cost: 1 biharmonic steps, 2 DSS
432 : !
433 : ! NOTE when nu_p=0 (no dissipation applied in dynamics to dp equation), we should
434 : ! apply dissipation to Q (not Qdp) to preserve Q=1
435 : ! i.e. laplace(Qdp) ~ dp0 laplace(Q)
436 : ! for nu_p=nu_q>0, we need to apply dissipation to Q * diffusion_dp
437 : !
438 : ! initialize dp, and compute Q from Qdp (and store Q in Qtens_biharmonic)
439 0 : do ie = nets, nete
440 : ! add hyperviscosity to RHS. apply to Q at timelevel n0, Qdp(n0)/dp
441 0 : do k = kbeg, kend
442 : !OMP_COLLAPSE_SIMD
443 : !DIR_VECTOR_ALIGNED
444 0 : do j=1,np
445 0 : do i=1,np
446 0 : dp(i,j,k) = elem(ie)%derived%dp(i,j,k) - rhs_multiplier*dt*elem(ie)%derived%divdp_proj(i,j,k)
447 : enddo
448 : enddo
449 : enddo
450 : !JMD need to update loop based on changes in dungeon21 tag
451 0 : do q = qbeg, qend
452 0 : do k= kbeg, kend
453 0 : Qtens_biharmonic(:,:,k,q,ie) = elem(ie)%state%Qdp(:,:,k,q,n0_qdp)/dp(:,:,k)
454 0 : if ( rhs_multiplier == 1 ) then
455 0 : qmin(k,q,ie)=min(qmin(k,q,ie),minval(Qtens_biharmonic(:,:,k,q,ie)))
456 0 : qmax(k,q,ie)=max(qmax(k,q,ie),maxval(Qtens_biharmonic(:,:,k,q,ie)))
457 : else
458 0 : qmin(k,q,ie)=minval(Qtens_biharmonic(:,:,k,q,ie))
459 0 : qmax(k,q,ie)=maxval(Qtens_biharmonic(:,:,k,q,ie))
460 : endif
461 : enddo
462 : enddo
463 : enddo
464 :
465 : ! compute element qmin/qmax
466 0 : if ( rhs_multiplier == 0 ) then
467 : ! update qmin/qmax based on neighbor data for lim8
468 : ! call t_startf('euler_neighbor_minmax1')
469 0 : call neighbor_minmax(hybrid,edgeAdvQminmax,nets,nete,qmin(:,:,nets:nete),qmax(:,:,nets:nete))
470 : ! call t_stopf('euler_neighbor_minmax1')
471 : endif
472 :
473 : ! get niew min/max values, and also compute biharmonic mixing term
474 0 : if ( rhs_multiplier == 2 ) then
475 0 : rhs_viss = 3
476 : ! two scalings depending on nu_p:
477 : ! nu_p=0: qtens_biharmonic *= dp0 (apply viscsoity only to q)
478 : ! nu_p>0): qtens_biharmonc *= elem()%psdiss_ave (for consistency, if nu_p=nu_q)
479 0 : if ( nu_p > 0 ) then
480 0 : do ie = nets, nete
481 0 : do k = kbeg, kend
482 : !OMP_COLLAPSE_SIMD
483 : !DIR_VECTOR_ALIGNED
484 0 : do j=1,np
485 0 : do i=1,np
486 0 : dpdissk(i,j,k) = elem(ie)%derived%dpdiss_ave(i,j,k)/dp0(k)
487 : enddo
488 : enddo
489 : enddo
490 0 : do q = qbeg,qend
491 0 : do k = kbeg, kend
492 : ! NOTE: divide by dp0 since we multiply by dp0 below
493 : !OMP_COLLAPSE_SIMD
494 : !DIR_VECTOR_ALIGNED
495 0 : do j=1,np
496 0 : do i=1,np
497 0 : Qtens_biharmonic(i,j,k,q,ie)=Qtens_biharmonic(i,j,k,q,ie)*dpdissk(i,j,k)
498 : enddo
499 : enddo
500 : enddo
501 : enddo
502 : enddo
503 : endif
504 :
505 : ! Previous version of biharmonic_wk_scalar_minmax included a min/max
506 : ! calculation into the boundary exchange. This was causing cache issues.
507 : ! Split the single operation into two separate calls
508 : ! call neighbor_minmax()
509 : ! call biharmonic_wk_scalar()
510 : !
511 : #ifdef OVERLAP
512 0 : call neighbor_minmax_start(hybrid,edgeAdvQminmax,nets,nete,qmin(:,:,nets:nete),qmax(:,:,nets:nete))
513 0 : call biharmonic_wk_scalar(elem,qtens_biharmonic,deriv,edgeAdv,hybrid,nets,nete)
514 0 : do ie = nets, nete
515 0 : do q = qbeg, qend
516 0 : do k = kbeg, kend
517 : !OMP_COLLAPSE_SIMD
518 : !DIR_VECTOR_ALIGNED
519 0 : do j=1,np
520 0 : do i=1,np
521 : ! note: biharmonic_wk() output has mass matrix already applied. Un-apply since we apply again below:
522 0 : qtens_biharmonic(i,j,k,q,ie) = &
523 0 : -rhs_viss*dt*nu_q*dp0(k)*Qtens_biharmonic(i,j,k,q,ie) / elem(ie)%spheremp(i,j)
524 : enddo
525 : enddo
526 : enddo
527 : enddo
528 : enddo
529 0 : call neighbor_minmax_finish(hybrid,edgeAdvQminmax,nets,nete,qmin(:,:,nets:nete),qmax(:,:,nets:nete))
530 : #else
531 : call t_startf('euler_neighbor_minmax2')
532 : call neighbor_minmax(hybrid,edgeAdvQminmax,nets,nete,qmin(:,:,nets:nete),qmax(:,:,nets:nete))
533 : call t_stopf('euler_neighbor_minmax2')
534 : call biharmonic_wk_scalar(elem,qtens_biharmonic,deriv,edgeAdv,hybrid,nets,nete)
535 :
536 : do ie = nets, nete
537 : do q = qbeg, qend
538 : do k = kbeg, kend
539 : !OMP_COLLAPSE_SIMD
540 : !DIR_VECTOR_ALIGNED
541 : do j=1,np
542 : do i=1,np
543 : ! note: biharmonic_wk() output has mass matrix already applied. Un-apply since we apply again below:
544 : qtens_biharmonic(i,j,k,q,ie) = &
545 : -rhs_viss*dt*nu_q*dp0(k)*Qtens_biharmonic(i,j,k,q,ie) / elem(ie)%spheremp(i,j)
546 : enddo
547 : enddo
548 : enddo
549 : enddo
550 : enddo
551 : #endif
552 :
553 :
554 : endif
555 : endif ! compute biharmonic mixing term and qmin/qmax
556 :
557 :
558 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
559 : ! 2D Advection step
560 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
561 0 : do ie = nets, nete
562 :
563 :
564 : ! Compute velocity used to advance Qdp
565 0 : do k = kbeg, kend
566 : ! derived variable divdp_proj() (DSS'd version of divdp) will only be correct on 2nd and 3rd stage
567 : ! but that's ok because rhs_multiplier=0 on the first stage:
568 : !OMP_COLLAPSE_SIMD
569 : !DIR_VECTOR_ALIGNED
570 0 : do j=1,np
571 0 : do i=1,np
572 0 : dp(i,j,k) = elem(ie)%derived%dp(i,j,k) - rhs_multiplier * dt * elem(ie)%derived%divdp_proj(i,j,k)
573 0 : Vstar(i,j,1,k) = elem(ie)%derived%vn0(i,j,1,k) / dp(i,j,k)
574 0 : Vstar(i,j,2,k) = elem(ie)%derived%vn0(i,j,2,k) / dp(i,j,k)
575 : enddo
576 : enddo
577 : enddo
578 0 : if ( limiter_option == 8) then
579 : ! Note that the term dpdissk is independent of Q
580 0 : do k = kbeg, kend
581 : ! UN-DSS'ed dp at timelevel n0+1:
582 : !OMP_COLLAPSE_SIMD
583 : !DIR_VECTOR_ALIGNED
584 0 : do j=1,np
585 0 : do i=1,np
586 0 : dpdissk(i,j,k) = dp(i,j,k) - dt * elem(ie)%derived%divdp(i,j,k)
587 : enddo
588 : enddo
589 0 : if ( nu_p > 0 .and. rhs_viss /= 0 ) then
590 : ! add contribution from UN-DSS'ed PS dissipation
591 : ! dpdiss(:,:) = ( hvcoord%hybi(k+1) - hvcoord%hybi(k) ) *
592 : ! elem(ie)%derived%psdiss_biharmonic(:,:)
593 : !OMP_COLLAPSE_SIMD
594 : !DIR_VECTOR_ALIGNED
595 0 : do j=1,np
596 0 : do i=1,np
597 0 : dpdiss(i,j) = elem(ie)%derived%dpdiss_biharmonic(i,j,k)
598 0 : dpdissk(i,j,k) = dpdissk(i,j,k) - rhs_viss * dt * nu_q * dpdiss(i,j) / elem(ie)%spheremp(i,j)
599 : enddo
600 : enddo
601 : endif
602 : ! IMPOSE ZERO THRESHOLD. do this here so it can be turned off for
603 : ! testing
604 0 : do q=qbeg, qend
605 0 : qmin(k,q,ie)=max(qmin(k,q,ie),0.0_r8)
606 : enddo
607 : enddo
608 : endif ! limiter == 8
609 :
610 :
611 : ! advance Qdp
612 0 : do q = qbeg, qend
613 0 : do k = kbeg, kend
614 : ! div( U dp Q),
615 : !OMP_COLLAPSE_SIMD
616 : !DIR_VECTOR_ALIGNED
617 0 : do j=1,np
618 0 : do i=1,np
619 0 : gradQ(i,j,1) = Vstar(i,j,1,k) * elem(ie)%state%Qdp(i,j,k,q,n0_qdp)
620 0 : gradQ(i,j,2) = Vstar(i,j,2,k) * elem(ie)%state%Qdp(i,j,k,q,n0_qdp)
621 : enddo
622 : enddo
623 : ! Qtens(:,:,k) = elem(ie)%state%Qdp(:,:,k,q,n0_qdp) - &
624 : ! dt * divergence_sphere( gradQ , deriv , elem(ie) )
625 0 : call divergence_sphere( gradQ , deriv , elem(ie),div )
626 :
627 : !OMP_COLLAPSE_SIMD
628 : !DIR_VECTOR_ALIGNED
629 0 : do j=1,np
630 0 : do i=1,np
631 0 : Qtens(i,j,k) = elem(ie)%state%Qdp(i,j,k,q,n0_qdp) - dt * div(i,j)
632 : enddo
633 : enddo
634 :
635 : ! optionally add in hyperviscosity computed above:
636 0 : if ( rhs_viss /= 0 ) then
637 : !OMP_COLLAPSE_SIMD
638 : !DIR_VECTOR_ALIGNED
639 0 : do j=1,np
640 0 : do i=1,np
641 0 : Qtens(i,j,k) = Qtens(i,j,k) + Qtens_biharmonic(i,j,k,q,ie)
642 : enddo
643 : enddo
644 : endif
645 : enddo
646 :
647 0 : if ( limiter_option == 8) then
648 : ! apply limiter to Q = Qtens / dp_star
649 0 : call limiter_optim_iter_full( Qtens(:,:,:) , elem(ie)%spheremp(:,:) , qmin(:,q,ie) , &
650 0 : qmax(:,q,ie) , dpdissk, kbeg, kend )
651 : endif
652 :
653 :
654 : ! apply mass matrix, overwrite np1 with solution:
655 : ! dont do this earlier, since we allow np1_qdp == n0_qdp
656 : ! and we dont want to overwrite n0_qdp until we are done using it
657 0 : do k = kbeg, kend
658 : !OMP_COLLAPSE_SIMD
659 : !DIR_VECTOR_ALIGNED
660 0 : do j=1,np
661 0 : do i=1,np
662 0 : elem(ie)%state%Qdp(i,j,k,q,np1_qdp) = elem(ie)%spheremp(i,j) * Qtens(i,j,k)
663 : enddo
664 : enddo
665 : enddo
666 :
667 0 : if ( limiter_option == 4 ) then
668 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
669 : ! sign-preserving limiter, applied after mass matrix
670 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
671 : !JMD !$OMP BARRIER
672 : !JMD !$OMP MASTER
673 0 : call limiter2d_zero(elem(ie)%state%Qdp(:,:,:,q,np1_qdp))
674 : !JMD !$OMP END MASTER
675 : !JMD !$OMP BARRIER
676 : endif
677 :
678 0 : kptr = nlev*(q-1) + kbeg - 1
679 0 : call edgeVpack(edgeAdvp1 , elem(ie)%state%Qdp(:,:,kbeg:kend,q,np1_qdp) , kblk , kptr , ie )
680 : enddo
681 : ! only perform this operation on thread which owns the first tracer
682 0 : if (DSSopt>0) then
683 0 : if (threadOwnsTracer(hybrid,1)) then
684 : ! all zero so we only have to DSS 1:nlev
685 0 : if ( DSSopt == DSSomega ) DSSvar => elem(ie)%derived%omega(:,:,:)
686 0 : if ( DSSopt == DSSdiv_vdp_ave ) DSSvar => elem(ie)%derived%divdp_proj(:,:,:)
687 : ! also DSS extra field
688 0 : do k = kbeg, kend
689 : !OMP_COLLAPSE_SIMD
690 : !DIR_VECTOR_ALIGNED
691 0 : do j=1,np
692 0 : do i=1,np
693 0 : DSSvar(i,j,k) = elem(ie)%spheremp(i,j) * DSSvar(i,j,k)
694 : enddo
695 : enddo
696 : enddo
697 0 : kptr = nlev*qsize + kbeg - 1
698 0 : call edgeVpack( edgeAdvp1 , DSSvar(:,:,kbeg:kend), kblk, kptr, ie)
699 : endif
700 : end if
701 : enddo
702 :
703 0 : call bndry_exchange( hybrid , edgeAdvp1,location='edgeAdvp1')
704 :
705 0 : do ie = nets, nete
706 : ! only perform this operation on thread which owns the first tracer
707 0 : if (DSSopt>0) then
708 0 : if(threadOwnsTracer(hybrid,1)) then
709 0 : if ( DSSopt == DSSomega ) DSSvar => elem(ie)%derived%omega(:,:,:)
710 0 : if ( DSSopt == DSSdiv_vdp_ave ) DSSvar => elem(ie)%derived%divdp_proj(:,:,:)
711 0 : kptr = qsize*nlev + kbeg -1
712 0 : call edgeVunpack( edgeAdvp1 , DSSvar(:,:,kbeg:kend) , kblk , kptr , ie )
713 0 : do k = kbeg, kend
714 : !OMP_COLLAPSE_SIMD
715 : !DIR_VECTOR_ALIGNED
716 0 : do j=1,np
717 0 : do i=1,np
718 0 : DSSvar(i,j,k) = DSSvar(i,j,k) * elem(ie)%rspheremp(i,j)
719 : enddo
720 : enddo
721 : enddo
722 : endif
723 : end if
724 0 : do q = qbeg, qend
725 0 : kptr = nlev*(q-1) + kbeg - 1
726 0 : call edgeVunpack( edgeAdvp1 , elem(ie)%state%Qdp(:,:,kbeg:kend,q,np1_qdp) , kblk , kptr , ie )
727 0 : do k = kbeg, kend
728 : !OMP_COLLAPSE_SIMD
729 : !DIR_VECTOR_ALIGNED
730 0 : do j=1,np
731 0 : do i=1,np
732 0 : elem(ie)%state%Qdp(i,j,k,q,np1_qdp) = elem(ie)%rspheremp(i,j) * elem(ie)%state%Qdp(i,j,k,q,np1_qdp)
733 : enddo
734 : enddo
735 : enddo
736 : enddo
737 : enddo
738 : ! call t_stopf('euler_step')
739 :
740 0 : end subroutine euler_step
741 :
742 :
743 :
744 0 : subroutine limiter2d_zero(Q)
745 : ! mass conserving zero limiter (2D only). to be called just before DSS
746 : !
747 : ! this routine is called inside a DSS loop, and so Q had already
748 : ! been multiplied by the mass matrix. Thus dont include the mass
749 : ! matrix when computing the mass = integral of Q over the element
750 : !
751 : ! ps is only used when advecting Q instead of Qdp
752 : ! so ps should be at one timelevel behind Q
753 : implicit none
754 : real (kind=r8), intent(inout) :: Q(np,np,nlev)
755 :
756 : ! local
757 : ! real (kind=r8) :: dp(np,np)
758 : real (kind=r8) :: mass,mass_new,ml
759 : integer i,j,k
760 :
761 0 : do k = nlev , 1 , -1
762 : mass = 0
763 0 : do j = 1 , np
764 0 : do i = 1 , np
765 : !ml = Q(i,j,k)*dp(i,j)*spheremp(i,j) ! see above
766 0 : ml = Q(i,j,k)
767 0 : mass = mass + ml
768 : enddo
769 : enddo
770 :
771 : ! negative mass. so reduce all postive values to zero
772 : ! then increase negative values as much as possible
773 0 : if ( mass < 0 ) Q(:,:,k) = -Q(:,:,k)
774 : mass_new = 0
775 0 : do j = 1 , np
776 0 : do i = 1 , np
777 0 : if ( Q(i,j,k) < 0 ) then
778 0 : Q(i,j,k) = 0
779 : else
780 0 : ml = Q(i,j,k)
781 0 : mass_new = mass_new + ml
782 : endif
783 : enddo
784 : enddo
785 :
786 : ! now scale the all positive values to restore mass
787 0 : if ( mass_new > 0 ) Q(:,:,k) = Q(:,:,k) * abs(mass) / mass_new
788 0 : if ( mass < 0 ) Q(:,:,k) = -Q(:,:,k)
789 : enddo
790 0 : end subroutine limiter2d_zero
791 :
792 : !-----------------------------------------------------------------------------
793 : !-----------------------------------------------------------------------------
794 :
795 0 : subroutine advance_hypervis_scalar( edgeAdv , elem , hvcoord , hybrid , deriv , nt , nt_qdp , nets , nete , dt2 )
796 : ! hyperviscsoity operator for foward-in-time scheme
797 : ! take one timestep of:
798 : ! Q(:,:,:,np) = Q(:,:,:,np) + dt2*nu*laplacian**order ( Q )
799 : !
800 : ! For correct scaling, dt2 should be the same 'dt2' used in the leapfrog advace
801 : use dimensions_mod , only: np, nlev
802 : use hybrid_mod , only: hybrid_t!, PrintHybrid
803 : use hybrid_mod , only: get_loop_ranges
804 : use element_mod , only: element_t
805 : use derivative_mod , only: derivative_t
806 : use edge_mod , only: edgevpack, edgevunpack
807 : use edgetype_mod , only: EdgeBuffer_t
808 : use bndry_mod , only: bndry_exchange
809 :
810 : implicit none
811 : type (EdgeBuffer_t) , intent(inout) :: edgeAdv
812 : type (element_t) , intent(inout), target :: elem(:)
813 : type (hvcoord_t) , intent(in ) :: hvcoord
814 : type (hybrid_t) , intent(in ) :: hybrid
815 : type (derivative_t) , intent(in ) :: deriv
816 : integer , intent(in ) :: nt
817 : integer , intent(in ) :: nt_qdp
818 : integer , intent(in ) :: nets
819 : integer , intent(in ) :: nete
820 : real (kind=r8), intent(in ) :: dt2
821 :
822 : ! local
823 0 : real (kind=r8), dimension(np,np,nlev,qsize,nets:nete) :: Qtens
824 : real (kind=r8), dimension(np,np,nlev ) :: dp
825 : ! real (kind=r8), dimension( nlev,qsize,nets:nete) :: min_neigh
826 : ! real (kind=r8), dimension( nlev,qsize,nets:nete) :: max_neigh
827 : integer :: k,kptr,ie,ic,q,i,j
828 : integer :: kbeg,kend,qbeg,qend
829 :
830 : ! NOTE: PGI compiler bug: when using spheremp, rspheremp and ps as pointers to elem(ie)% members,
831 : ! data is incorrect (offset by a few numbers actually)
832 : ! removed for now.
833 : ! real (kind=r8), dimension(:,:), pointer :: spheremp,rspheremp
834 : real (kind=r8) :: dt,dp0
835 : integer :: kblk,qblk ! The per thead size of the vertical and tracers
836 :
837 0 : call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
838 :
839 0 : if ( nu_q == 0 ) return
840 : !if ( hypervis_order /= 2 ) return
841 :
842 0 : kblk = kend - kbeg + 1 ! calculate size of the block of vertical levels
843 0 : qblk = qend - qbeg + 1 ! calculate size of the block of tracers
844 :
845 0 : call t_startf('advance_hypervis_scalar')
846 :
847 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
848 : ! hyper viscosity
849 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
850 0 : dt = dt2 / hypervis_subcycle_q
851 :
852 0 : do ic = 1 , hypervis_subcycle_q
853 0 : do ie = nets, nete
854 : ! Qtens = Q/dp (apply hyperviscsoity to dp0 * Q, not Qdp)
855 0 : do k = kbeg, kend
856 : ! various options:
857 : ! 1) biharmonic( Qdp )
858 : ! 2) dp0 * biharmonic( Qdp/dp )
859 : ! 3) dpave * biharmonic(Q/dp)
860 : ! For trace mass / mass consistenciy, we use #2 when nu_p=0
861 : ! and #e when nu_p>0, where dpave is the mean mass flux from the nu_p
862 : ! contribution from dynamics.
863 0 : dp0 = ( hvcoord%hyai(k+1) - hvcoord%hyai(k) ) * hvcoord%ps0 + &
864 0 : ( hvcoord%hybi(k+1) - hvcoord%hybi(k) ) * hvcoord%ps0
865 0 : dp(:,:,k) = elem(ie)%derived%dp(:,:,k) - dt2*elem(ie)%derived%divdp_proj(:,:,k)
866 0 : if (nu_p>0) then
867 0 : do q = qbeg, qend
868 0 : Qtens(:,:,k,q,ie) = elem(ie)%derived%dpdiss_ave(:,:,k)*&
869 0 : elem(ie)%state%Qdp(:,:,k,q,nt_qdp) / dp(:,:,k)
870 : enddo
871 : else
872 0 : do q = qbeg, qend
873 0 : Qtens(:,:,k,q,ie) = dp0*elem(ie)%state%Qdp(:,:,k,q,nt_qdp) / dp(:,:,k)
874 : enddo
875 : endif
876 : enddo
877 : enddo
878 :
879 : ! compute biharmonic operator. Qtens = input and output
880 0 : call biharmonic_wk_scalar( elem , Qtens , deriv , edgeAdv , hybrid , nets , nete )
881 :
882 0 : do ie = nets, nete
883 : !spheremp => elem(ie)%spheremp
884 0 : do q = qbeg, qend
885 0 : do k = kbeg, kend
886 0 : dp0 = ( hvcoord%hyai(k+1) - hvcoord%hyai(k) ) * hvcoord%ps0 + &
887 0 : ( hvcoord%hybi(k+1) - hvcoord%hybi(k) ) * hvcoord%ps0
888 0 : do j = 1 , np
889 0 : do i = 1 , np
890 :
891 : ! advection Qdp. For mass advection consistency:
892 : ! DIFF( Qdp) ~ dp0 DIFF (Q) = dp0 DIFF ( Qdp/dp )
893 0 : elem(ie)%state%Qdp(i,j,k,q,nt_qdp) = elem(ie)%state%Qdp(i,j,k,q,nt_qdp) * elem(ie)%spheremp(i,j) &
894 0 : - dt * nu_q * Qtens(i,j,k,q,ie)
895 : enddo
896 : enddo
897 : enddo
898 :
899 0 : if (limiter_option .ne. 0 ) then
900 : !JMD Only need if threading over the vertical
901 : !JMD!$OMP BARRIER
902 : !JMD!$OMP MASTER
903 : ! smooth some of the negativities introduced by diffusion:
904 0 : call limiter2d_zero( elem(ie)%state%Qdp(:,:,:,q,nt_qdp) )
905 : !JMD!$OMP END MASTER
906 : !JMD!$OMP BARRIER
907 : endif
908 :
909 : enddo
910 0 : do q = qbeg, qend
911 0 : kptr = nlev*(q-1) + kbeg - 1
912 0 : call edgeVpack( edgeAdv , elem(ie)%state%Qdp(:,:,kbeg:kend,q,nt_qdp) , kblk, kptr, ie )
913 : enddo
914 : enddo
915 :
916 0 : call bndry_exchange( hybrid , edgeAdv,location='advance_hypervis_scalar')
917 :
918 0 : do ie = nets, nete
919 0 : do q = qbeg, qend
920 0 : kptr = nlev*(q-1) + kbeg - 1
921 0 : call edgeVunpack( edgeAdv , elem(ie)%state%Qdp(:,:,kbeg:kend,q,nt_qdp) , kblk, kptr, ie )
922 : enddo
923 : !rspheremp => elem(ie)%rspheremp
924 0 : do q = qbeg, qend
925 : ! apply inverse mass matrix
926 0 : do k = kbeg, kend
927 0 : elem(ie)%state%Qdp(:,:,k,q,nt_qdp) = elem(ie)%rspheremp(:,:) * elem(ie)%state%Qdp(:,:,k,q,nt_qdp)
928 : enddo
929 : enddo
930 : enddo
931 :
932 : enddo
933 0 : call t_stopf('advance_hypervis_scalar')
934 0 : end subroutine advance_hypervis_scalar
935 :
936 738816 : subroutine vertical_remap(hybrid,elem,fvm,hvcoord,np1,np1_qdp,nets,nete)
937 : !
938 : ! This routine is called at the end of the vertically Lagrangian
939 : ! dynamics step to compute the vertical flux needed to get back
940 : ! to reference eta levels
941 : !
942 : ! map tracers
943 : ! map velocity components
944 : ! map temperature (either by mapping enthalpy or virtual temperature over log(p)
945 : ! (controlled by vert_remap_uvTq_alg > -20 or <= -20)
946 : !
947 0 : use hybvcoord_mod, only: hvcoord_t
948 : use vertremap_mod, only: remap1
949 : use hybrid_mod, only: hybrid_t, config_thread_region,get_loop_ranges, PrintHybrid
950 : use fvm_control_volume_mod, only: fvm_struct
951 : use dimensions_mod, only: ntrac
952 : use dimensions_mod, only: kord_tr,kord_tr_cslam
953 : use cam_logfile, only: iulog
954 : use physconst, only: pi
955 : use air_composition, only: thermodynamic_active_species_idx_dycore
956 : use cam_thermo, only: get_enthalpy, get_virtual_temp, get_dp, MASS_MIXING_RATIO
957 : use thread_mod, only: omp_set_nested
958 : use control_mod, only: vert_remap_uvTq_alg
959 : type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared)
960 : type(fvm_struct), intent(inout) :: fvm(:)
961 : type (element_t), intent(inout) :: elem(:)
962 : !
963 : real (kind=r8) :: dpc_star(nc,nc,nlev) !Lagrangian levels on CSLAM grid
964 :
965 : type (hvcoord_t) :: hvcoord
966 : integer :: ie,i,j,k,np1,nets,nete,np1_qdp,q, m_cnst
967 : real (kind=r8), dimension(np,np,nlev) :: dp_moist,dp_star_moist, dp_dry,dp_star_dry
968 : real (kind=r8), dimension(np,np,nlev) :: enthalpy_star
969 : real (kind=r8), dimension(np,np,nlev,2):: ttmp
970 : real(r8), parameter :: rad2deg = 180.0_r8/pi
971 : integer :: region_num_threads,qbeg,qend,kord_uvT(1)
972 : type (hybrid_t) :: hybridnew,hybridnew2
973 : real (kind=r8) :: ptop
974 :
975 1477632 : kord_uvT = vert_remap_uvTq_alg
976 :
977 738816 : ptop = hvcoord%hyai(1)*hvcoord%ps0
978 5933616 : do ie=nets,nete
979 : !
980 : ! prepare for mapping of temperature
981 : !
982 5194800 : if (vert_remap_uvTq_alg>-20) then
983 : !
984 : ! compute enthalpy on Lagrangian levels
985 : ! (do it here since qdp is overwritten by remap1)
986 : !
987 10389600 : call get_enthalpy(elem(ie)%state%qdp(:,:,:,1:qsize,np1_qdp), &
988 : elem(ie)%state%t(:,:,:,np1), elem(ie)%state%dp3d(:,:,:,np1), enthalpy_star, &
989 15584400 : active_species_idx_dycore=thermodynamic_active_species_idx_dycore)
990 : else
991 : !
992 : ! map Tv over log(p) following FV and FV3
993 : !
994 0 : call get_virtual_temp(elem(ie)%state%qdp(:,:,:,1:qsize,np1_qdp), enthalpy_star, &
995 0 : dp_dry=elem(ie)%state%dp3d(:,:,:,np1), active_species_idx_dycore=thermodynamic_active_species_idx_dycore)
996 0 : enthalpy_star = enthalpy_star*elem(ie)%state%t(:,:,:,np1)
997 : end if
998 : !
999 : ! update final psdry
1000 : !
1001 : elem(ie)%state%psdry(:,:) = ptop + &
1002 2270127600 : sum(elem(ie)%state%dp3d(:,:,:,np1),3)
1003 : !
1004 : ! compute dry vertical coordinate (Lagrangian and reference levels)
1005 : !
1006 140259600 : do k=1,nlev
1007 2836360800 : dp_star_dry(:,:,k) = elem(ie)%state%dp3d(:,:,k,np1)
1008 135064800 : dp_dry(:,:,k) = ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + &
1009 2971425600 : ( hvcoord%hybi(k+1) - hvcoord%hybi(k) )*elem(ie)%state%psdry(:,:)
1010 2841555600 : elem(ie)%state%dp3d(:,:,k,np1) = dp_dry(:,:,k)
1011 : enddo
1012 : !
1013 5194800 : call get_dp(elem(ie)%state%Qdp(:,:,:,1:qsize,np1_qdp), MASS_MIXING_RATIO,&
1014 10389600 : thermodynamic_active_species_idx_dycore, dp_star_dry, dp_star_moist(:,:,:))
1015 : !
1016 : ! Check if Lagrangian leves have crossed
1017 : !
1018 2841555600 : if (minval(dp_star_moist)<1.0E-12_r8) then
1019 0 : write(iulog,*) "NEGATIVE LAYER THICKNESS DIAGNOSTICS:"
1020 0 : write(iulog,*) " "
1021 0 : do j=1,np
1022 0 : do i=1,np
1023 0 : if (minval(dp_star_moist(i,j,:))<1.0e-12_r8) then
1024 0 : write(iulog,'(A13,2f6.2)') "(lon,lat) = ",&
1025 0 : elem(ie)%spherep(i,j)%lon*rad2deg,elem(ie)%spherep(i,j)%lat*rad2deg
1026 0 : write(iulog,*) " "
1027 0 : do k=1,nlev
1028 0 : write(iulog,'(A21,I5,A1,f16.12,3f10.2)') "k,dp_star_moist,u,v,T: ",k," ",dp_star_moist(i,j,k)/100.0_r8,&
1029 0 : elem(ie)%state%v(i,j,1,k,np1),elem(ie)%state%v(i,j,2,k,np1),elem(ie)%state%T(i,j,k,np1)
1030 : end do
1031 : end if
1032 : end do
1033 : end do
1034 0 : call endrun('negative moist layer thickness. timestep or remap time too large')
1035 : endif
1036 :
1037 5194800 : call remap1(elem(ie)%state%Qdp(:,:,:,1:qsize,np1_qdp),np,1,qsize,qsize,dp_star_dry,dp_dry,ptop,0,.true.,kord_tr)
1038 : !
1039 : ! compute moist reference pressure level thickness
1040 : !
1041 5194800 : call get_dp(elem(ie)%state%Qdp(:,:,:,1:qsize,np1_qdp), MASS_MIXING_RATIO,&
1042 10389600 : thermodynamic_active_species_idx_dycore, dp_dry, dp_moist(:,:,:))
1043 :
1044 : !
1045 : ! Remapping of temperature
1046 : !
1047 5194800 : if (vert_remap_uvTq_alg>-20) then
1048 : !
1049 : ! remap enthalpy energy and back out temperature
1050 : !
1051 5194800 : call remap1(enthalpy_star,np,1,1,1,dp_star_dry,dp_dry,ptop,1,.true.,kord_uvT)
1052 : !
1053 : ! compute sum c^(l)_p*m^(l)*dp on arrival (Eulerian) grid
1054 : !
1055 2841555600 : ttmp(:,:,:,1) = 1.0_r8
1056 5194800 : call get_enthalpy(elem(ie)%state%qdp(:,:,:,1:qsize,np1_qdp), &
1057 : ttmp(:,:,:,1), dp_dry,ttmp(:,:,:,2), &
1058 10389600 : active_species_idx_dycore=thermodynamic_active_species_idx_dycore)
1059 2841555600 : elem(ie)%state%t(:,:,:,np1)=enthalpy_star/ttmp(:,:,:,2)
1060 : else
1061 : !
1062 : ! map Tv over log(p); following FV and FV3
1063 : !
1064 0 : call remap1(enthalpy_star,np,1,1,1,dp_star_moist,dp_moist,ptop,1,.false.,kord_uvT)
1065 0 : call get_virtual_temp(elem(ie)%state%qdp(:,:,:,1:qsize,np1_qdp), ttmp(:,:,:,1), &
1066 0 : dp_dry=dp_dry, active_species_idx_dycore=thermodynamic_active_species_idx_dycore)
1067 : !
1068 : ! convert new Tv to T
1069 : !
1070 0 : elem(ie)%state%t(:,:,:,np1)=enthalpy_star/ttmp(:,:,:,1)
1071 : end if
1072 : !
1073 : ! remap velocity components
1074 : !
1075 5677916400 : call remap1(elem(ie)%state%v(:,:,1,:,np1),np,1,1,1,dp_star_moist,dp_moist,ptop,-1,.false.,kord_uvT)
1076 5678655216 : call remap1(elem(ie)%state%v(:,:,2,:,np1),np,1,1,1,dp_star_moist,dp_moist,ptop,-1,.false.,kord_uvT)
1077 : enddo
1078 :
1079 738816 : if (ntrac>0) then
1080 : !
1081 : ! vertical remapping of CSLAM tracers
1082 : !
1083 5933616 : do ie=nets,nete
1084 1761037200 : dpc_star=fvm(ie)%dp_fvm(1:nc,1:nc,:)
1085 140259600 : do k=1,nlev
1086 545454000 : do j=1,nc
1087 1755842400 : do i=1,nc
1088 : !
1089 : ! new pressure levels on CSLAM grid
1090 : !
1091 4862332800 : fvm(ie)%dp_fvm(i,j,k) = (hvcoord%hyai(k+1) - hvcoord%hyai(k))*hvcoord%ps0 + &
1092 6483110400 : (hvcoord%hybi(k+1) - hvcoord%hybi(k))*fvm(ie)%psc(i,j)
1093 : end do
1094 : end do
1095 : end do
1096 5933616 : if(ntrac>tracer_num_threads) then
1097 5194800 : call omp_set_nested(.true.)
1098 : !$OMP PARALLEL NUM_THREADS(tracer_num_threads), DEFAULT(SHARED), PRIVATE(hybridnew2,qbeg,qend)
1099 5194800 : hybridnew2 = config_thread_region(hybrid,'ctracer')
1100 5194800 : call get_loop_ranges(hybridnew2, qbeg=qbeg, qend=qend)
1101 0 : call remap1(fvm(ie)%c(1:nc,1:nc,:,1:ntrac),nc,qbeg,qend,ntrac,dpc_star, &
1102 12332455200 : fvm(ie)%dp_fvm(1:nc,1:nc,:),ptop,0,.false.,kord_tr_cslam)
1103 : !$OMP END PARALLEL
1104 5194800 : call omp_set_nested(.false.)
1105 : else
1106 0 : call remap1(fvm(ie)%c(1:nc,1:nc,:,1:ntrac),nc,1,ntrac,ntrac,dpc_star, &
1107 0 : fvm(ie)%dp_fvm(1:nc,1:nc,:),ptop,0,.false.,kord_tr_cslam)
1108 : endif
1109 : enddo
1110 : end if
1111 738816 : end subroutine vertical_remap
1112 :
1113 : end module prim_advection_mod
|
