Line data Source code
1 : module global_norms_mod
2 :
3 : use shr_kind_mod, only: r8=>shr_kind_r8
4 : use cam_logfile, only: iulog
5 : use edgetype_mod, only: EdgeBuffer_t
6 :
7 : implicit none
8 : private
9 : save
10 :
11 : public :: l1_snorm
12 : public :: l2_snorm
13 : public :: linf_snorm
14 :
15 : public :: l1_vnorm
16 : public :: l2_vnorm
17 : public :: linf_vnorm
18 :
19 : public :: print_cfl
20 : public :: global_integral
21 : public :: global_integrals_general
22 : public :: wrap_repro_sum
23 :
24 : private :: global_maximum
25 : type (EdgeBuffer_t), private :: edgebuf
26 :
27 : interface global_integral
28 : module procedure global_integral_elem
29 : module procedure global_integral_fvm
30 : end interface global_integral
31 :
32 : contains
33 :
34 :
35 : subroutine global_integrals(elem,fld,hybrid,npts,num_flds,nets,nete,I_sphere)
36 : use hybrid_mod, only: hybrid_t
37 : use element_mod, only: element_t
38 : use dimensions_mod, only: np
39 : use physconst, only: pi
40 : use parallel_mod, only: global_shared_buf, global_shared_sum
41 :
42 : type(element_t) , intent(in) :: elem(:)
43 : integer , intent(in) :: npts,nets,nete,num_flds
44 : real (kind=r8), intent(in) :: fld(npts,npts,num_flds,nets:nete)
45 : type (hybrid_t) , intent(in) :: hybrid
46 :
47 : real (kind=r8) :: I_sphere(num_flds)
48 : !
49 : ! Local variables
50 : !
51 : integer :: ie,j,i,q
52 :
53 : real (kind=r8) :: da
54 : real (kind=r8) :: J_tmp(nets:nete,num_flds)
55 : !
56 : ! This algorithm is independent of thread count and task count.
57 : ! This is a requirement of consistancy checking in cam.
58 : !
59 : J_tmp = 0.0_r8
60 :
61 : do ie=nets,nete
62 : do q=1,num_flds
63 : do j=1,np
64 : do i=1,np
65 : da = elem(ie)%mp(i,j)*elem(ie)%metdet(i,j)
66 : J_tmp(ie,q) = J_tmp(ie,q) + da*fld(i,j,q,ie)
67 : end do
68 : end do
69 : end do
70 : end do
71 : do ie=nets,nete
72 : global_shared_buf(ie,1:num_flds) = J_tmp(ie,:)
73 : enddo
74 : call wrap_repro_sum(nvars=num_flds, comm=hybrid%par%comm)
75 : I_sphere(:) =global_shared_sum(1:num_flds) /(4.0_r8*PI)
76 : end subroutine global_integrals
77 :
78 3072 : subroutine global_integrals_general(fld,hybrid,npts,da,num_flds,nets,nete,I_sphere)
79 : use hybrid_mod, only: hybrid_t
80 : use physconst, only: pi
81 : use parallel_mod, only: global_shared_buf, global_shared_sum
82 :
83 : integer, intent(in) :: npts,nets,nete,num_flds
84 : real (kind=r8), intent(in) :: fld(npts,npts,num_flds,nets:nete)
85 : type (hybrid_t), intent(in) :: hybrid
86 : real (kind=r8), intent(in) :: da(npts,npts,nets:nete)
87 :
88 : real (kind=r8) :: I_sphere(num_flds)
89 : !
90 : ! Local variables
91 : !
92 : integer :: ie,j,i,q
93 :
94 6144 : real (kind=r8) :: J_tmp(nets:nete,num_flds)
95 : !
96 : ! This algorithm is independent of thread count and task count.
97 : ! This is a requirement of consistancy checking in cam.
98 : !
99 89424 : J_tmp = 0.0_r8
100 :
101 24672 : do ie=nets,nete
102 100272 : do q=1,num_flds
103 345600 : do j=1,npts
104 1155600 : do i=1,npts
105 1080000 : J_tmp(ie,q) = J_tmp(ie,q) + da(i,j,ie)*fld(i,j,q,ie)
106 : end do
107 : end do
108 : end do
109 : end do
110 24672 : do ie=nets,nete
111 100272 : global_shared_buf(ie,1:num_flds) = J_tmp(ie,:)
112 : enddo
113 3072 : call wrap_repro_sum(nvars=num_flds, comm=hybrid%par%comm)
114 16896 : I_sphere(:) =global_shared_sum(1:num_flds) /(4.0_r8*PI)
115 3072 : end subroutine global_integrals_general
116 :
117 :
118 : ! ================================
119 : ! global_integral:
120 : !
121 : ! eq 81 in Williamson, et. al. p 218
122 : ! for spectral elements
123 : !
124 : ! ================================
125 : ! --------------------------
126 2304 : function global_integral_elem(elem,fld,hybrid,npts,nets,nete) result(I_sphere)
127 : use hybrid_mod, only: hybrid_t
128 : use element_mod, only: element_t
129 : use dimensions_mod, only: np
130 : use physconst, only: pi
131 : use parallel_mod, only: global_shared_buf, global_shared_sum
132 :
133 : type(element_t) , intent(in) :: elem(:)
134 : integer , intent(in) :: npts,nets,nete
135 : real (kind=r8), intent(in) :: fld(npts,npts,nets:nete)
136 : type (hybrid_t) , intent(in) :: hybrid
137 :
138 : real (kind=r8) :: I_sphere
139 :
140 : ! Local variables
141 :
142 : integer :: ie,j,i
143 : real(kind=r8) :: I_tmp(1)
144 :
145 : real (kind=r8) :: da
146 2304 : real (kind=r8) :: J_tmp(nets:nete)
147 : !
148 : ! This algorithm is independent of thread count and task count.
149 : ! This is a requirement of consistancy checking in cam.
150 : !
151 18504 : J_tmp = 0.0_r8
152 :
153 18504 : do ie=nets,nete
154 83304 : do j=1,np
155 340200 : do i=1,np
156 259200 : da = elem(ie)%mp(i,j)*elem(ie)%metdet(i,j)
157 324000 : J_tmp(ie) = J_tmp(ie) + da*fld(i,j,ie)
158 : end do
159 : end do
160 : end do
161 18504 : do ie=nets,nete
162 18504 : global_shared_buf(ie,1) = J_tmp(ie)
163 : enddo
164 2304 : call wrap_repro_sum(nvars=1, comm=hybrid%par%comm)
165 4608 : I_tmp = global_shared_sum(1)
166 2304 : I_sphere = I_tmp(1)/(4.0_r8*PI)
167 :
168 2304 : end function global_integral_elem
169 :
170 0 : function global_integral_fvm(fvm,fld,hybrid,npts,nets,nete) result(I_sphere)
171 : use hybrid_mod, only: hybrid_t
172 : use fvm_control_volume_mod, only: fvm_struct
173 : use physconst, only: pi
174 : use parallel_mod, only: global_shared_buf, global_shared_sum
175 :
176 : type (fvm_struct) , intent(in) :: fvm(:)
177 : integer , intent(in) :: npts,nets,nete
178 : real (kind=r8), intent(in) :: fld(npts,npts,nets:nete)
179 : type (hybrid_t) , intent(in) :: hybrid
180 :
181 : real (kind=r8) :: I_sphere
182 :
183 : ! Local variables
184 :
185 : integer :: ie,j,i
186 : real(kind=r8) :: I_tmp(1)
187 :
188 : real (kind=r8) :: da
189 0 : real (kind=r8) :: J_tmp(nets:nete)
190 : !
191 : ! This algorithm is independent of thread count and task count.
192 : ! This is a requirement of consistancy checking in cam.
193 : !
194 0 : J_tmp = 0.0_r8
195 0 : do ie=nets,nete
196 0 : do j=1,npts
197 0 : do i=1,npts
198 0 : da = fvm(ie)%area_sphere(i,j)
199 0 : J_tmp(ie) = J_tmp(ie) + da*fld(i,j,ie)
200 : end do
201 : end do
202 : end do
203 0 : do ie=nets,nete
204 0 : global_shared_buf(ie,1) = J_tmp(ie)
205 : enddo
206 0 : call wrap_repro_sum(nvars=1, comm=hybrid%par%comm)
207 0 : I_tmp = global_shared_sum(1)
208 0 : I_sphere = I_tmp(1)/(4.0_r8*PI)
209 :
210 0 : end function global_integral_fvm
211 :
212 : !------------------------------------------------------------------------------------
213 :
214 : ! ================================
215 : ! print_cfl:
216 : !
217 : ! Calculate / output CFL info
218 : ! (both advective and based on
219 : ! viscosity or hyperviscosity)
220 : !
221 : ! ================================
222 :
223 1536 : subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,&
224 : dt_remap_actual,dt_tracer_fvm_actual,dt_tracer_se_actual,&
225 : dt_dyn_actual,dt_dyn_visco_actual,dt_dyn_del2_actual,dt_tracer_visco_actual,dt_phys)
226 : !
227 : ! estimate various CFL limits
228 : ! also, for variable resolution viscosity coefficient, make sure
229 : ! worse viscosity CFL (given by dtnu) is not violated by reducing
230 : ! viscosity coefficient in regions where CFL is violated
231 : !
232 : use hybrid_mod, only: hybrid_t
233 : use element_mod, only: element_t
234 : use dimensions_mod, only: np,ne,nelem,nc,nhe,use_cslam,nlev,large_Courant_incr
235 : use dimensions_mod, only: nu_scale_top,nu_div_lev,nu_lev,nu_t_lev
236 :
237 : use quadrature_mod, only: gausslobatto, quadrature_t
238 :
239 : use reduction_mod, only: ParallelMin,ParallelMax
240 : use physconst, only: ra, rearth, pi
241 : use control_mod, only: nu, nu_div, nu_q, nu_p, nu_t, nu_top, fine_ne, max_hypervis_courant
242 : use control_mod, only: tstep_type, hypervis_power, hypervis_scaling
243 : use control_mod, only: sponge_del4_nu_div_fac, sponge_del4_nu_fac, sponge_del4_lev
244 : use cam_abortutils, only: endrun
245 : use parallel_mod, only: global_shared_buf, global_shared_sum
246 : use edge_mod, only: initedgebuffer, FreeEdgeBuffer, edgeVpack, edgeVunpack
247 : use bndry_mod, only: bndry_exchange
248 : use mesh_mod, only: MeshUseMeshFile
249 : use dimensions_mod, only: ksponge_end, kmvis_ref, kmcnd_ref,rho_ref
250 : use physconst, only: cpair
251 : use std_atm_profile,only: std_atm_height
252 : type(element_t) , intent(inout) :: elem(:)
253 : integer , intent(in) :: nets,nete
254 : type (hybrid_t) , intent(in) :: hybrid
255 : real (kind=r8), intent(in) :: dtnu, ptop, pmid(nlev)
256 : !
257 : ! actual time-steps
258 : !
259 : real (kind=r8), intent(in) :: dt_remap_actual,dt_tracer_fvm_actual,dt_tracer_se_actual,&
260 : dt_dyn_actual,dt_dyn_visco_actual,dt_dyn_del2_actual, &
261 : dt_tracer_visco_actual, dt_phys
262 :
263 : ! Element statisics
264 : real (kind=r8) :: max_min_dx,min_min_dx,min_max_dx,max_unif_dx ! used for normalizing scalar HV
265 : real (kind=r8) :: max_normDinv, min_normDinv ! used for CFL
266 : real (kind=r8) :: min_area, max_area,max_ratio !min/max element area
267 : real (kind=r8) :: avg_area, avg_min_dx,tot_area,tot_area_rad
268 : real (kind=r8) :: min_hypervis, max_hypervis, avg_hypervis, stable_hv
269 : real (kind=r8) :: normDinv_hypervis
270 : real (kind=r8) :: x, y, noreast, nw, se, sw
271 3072 : real (kind=r8), dimension(np,np,nets:nete) :: zeta
272 : real (kind=r8) :: lambda_max, lambda_vis, min_gw, lambda,umax, ugw
273 : real (kind=r8) :: scale1, max_laplace,z(nlev)
274 : integer :: ie, i, j, rowind, colind, k
275 : type (quadrature_t) :: gp
276 : character(LEN=256) :: rk_str
277 :
278 : real (kind=r8) :: s_laplacian, s_hypervis, s_rk, s_rk_tracer !Stability region
279 : real (kind=r8) :: dt_max_adv, dt_max_gw, dt_max_tracer_se, dt_max_tracer_fvm
280 : real (kind=r8) :: dt_max_hypervis, dt_max_hypervis_tracer, dt_max_laplacian_top
281 :
282 : real(kind=r8) :: I_sphere, nu_max, nu_div_max
283 3072 : real(kind=r8) :: fld(np,np,nets:nete)
284 :
285 : logical :: top_000_032km, top_032_042km, top_042_090km, top_090_140km, top_140_600km ! model top location ranges
286 : logical :: nu_set,div_set,lev_set
287 :
288 : ! Eigenvalues calculated by folks at UMich (Paul U & Jared W)
289 : select case (np)
290 : case (2)
291 : lambda_max = 0.5_r8
292 : lambda_vis = 0.0_r8 ! need to compute this
293 : case (3)
294 : lambda_max = 1.5_r8
295 : lambda_vis = 12.0_r8
296 : case (4)
297 1536 : lambda_max = 2.74_r8
298 1536 : lambda_vis = 30.0_r8
299 : case (5)
300 : lambda_max = 4.18_r8
301 : lambda_vis = 91.6742_r8
302 : case (6)
303 : lambda_max = 5.86_r8
304 : lambda_vis = 190.1176_r8
305 : case (7)
306 : lambda_max = 7.79_r8
307 : lambda_vis = 374.7788_r8
308 : case (8)
309 : lambda_max = 10.0_r8
310 : lambda_vis = 652.3015_r8
311 : case DEFAULT
312 : lambda_max = 0.0_r8
313 : lambda_vis = 0.0_r8
314 : end select
315 :
316 1536 : if ((lambda_max.eq.0_r8).and.(hybrid%masterthread)) then
317 : print*, "lambda_max not calculated for NP = ",np
318 : print*, "Estimate of gravity wave timestep will be incorrect"
319 : end if
320 : if ((lambda_vis.eq.0_r8).and.(hybrid%masterthread)) then
321 : print*, "lambda_vis not calculated for NP = ",np
322 : print*, "Estimate of viscous CFLs will be incorrect"
323 : end if
324 :
325 12336 : do ie=nets,nete
326 228336 : elem(ie)%variable_hyperviscosity = 1.0_r8
327 : end do
328 :
329 1536 : gp=gausslobatto(np)
330 9216 : min_gw = minval(gp%weights)
331 : !
332 : !******************************************************************************************
333 : !
334 : ! compute some local and global grid metrics
335 : !
336 : !******************************************************************************************
337 : !
338 228336 : fld(:,:,nets:nete)=1.0_r8
339 : ! Calculate surface area by integrating 1.0_r8 over sphere and dividing by 4*PI (Should be 1)
340 1536 : I_sphere = global_integral(elem, fld(:,:,nets:nete),hybrid,np,nets,nete)
341 :
342 1536 : min_normDinv = 1E99_r8
343 1536 : max_normDinv = 0
344 1536 : min_max_dx = 1E99_r8
345 1536 : min_min_dx = 1E99_r8
346 1536 : max_min_dx = 0
347 1536 : min_area = 1E99_r8
348 1536 : max_area = 0
349 1536 : max_ratio = 0
350 12336 : do ie=nets,nete
351 10800 : max_normDinv = max(max_normDinv,elem(ie)%normDinv)
352 10800 : min_normDinv = min(min_normDinv,elem(ie)%normDinv)
353 10800 : min_min_dx = min(min_min_dx,elem(ie)%dx_short)
354 10800 : max_min_dx = max(max_min_dx,elem(ie)%dx_short)
355 10800 : min_max_dx = min(min_max_dx,elem(ie)%dx_long)
356 :
357 226800 : elem(ie)%area = sum(elem(ie)%spheremp(:,:))
358 10800 : min_area = min(min_area,elem(ie)%area)
359 10800 : max_area = max(max_area,elem(ie)%area)
360 10800 : max_ratio = max(max_ratio,elem(ie)%dx_long/elem(ie)%dx_short)
361 :
362 10800 : global_shared_buf(ie,1) = elem(ie)%area
363 12336 : global_shared_buf(ie,2) = elem(ie)%dx_short
364 : enddo
365 1536 : call wrap_repro_sum(nvars=2, comm=hybrid%par%comm)
366 1536 : avg_area = global_shared_sum(1)/dble(nelem)
367 1536 : tot_area_rad = global_shared_sum(1)
368 1536 : avg_min_dx = global_shared_sum(2)/dble(nelem)
369 :
370 1536 : min_area = ParallelMin(min_area,hybrid)
371 1536 : max_area = ParallelMax(max_area,hybrid)
372 1536 : min_normDinv = ParallelMin(min_normDinv,hybrid)
373 1536 : max_normDinv = ParallelMax(max_normDinv,hybrid)
374 1536 : min_min_dx = ParallelMin(min_min_dx,hybrid)
375 1536 : max_min_dx = ParallelMax(max_min_dx,hybrid)
376 1536 : min_max_dx = ParallelMin(min_max_dx,hybrid)
377 1536 : max_ratio = ParallelMax(max_ratio,hybrid)
378 : ! Physical units for area (unit sphere to Earth sphere)
379 1536 : min_area = min_area*rearth*rearth/1000000._r8 !m2 (rearth is in units of km)
380 1536 : max_area = max_area*rearth*rearth/1000000._r8 !m2 (rearth is in units of km)
381 1536 : avg_area = avg_area*rearth*rearth/1000000._r8 !m2 (rearth is in units of km)
382 1536 : tot_area = tot_area_rad*rearth*rearth/1000000._r8!m2 (rearth is in units of km)
383 1536 : if (hybrid%masterthread) then
384 2 : write(iulog,* )""
385 2 : write(iulog,* )"Running Global Integral Diagnostic..."
386 2 : write(iulog,*)"Area of unit sphere is",I_sphere
387 2 : write(iulog,*)"Should be 1.0 to round off..."
388 2 : write(iulog,'(a,f9.3)') 'Element area: max/min',(max_area/min_area)
389 2 : write(iulog,'(a,E23.15)') 'Total Grid area: ',(tot_area)
390 2 : write(iulog,'(a,E23.15)') 'Total Grid area rad^2: ',(tot_area_rad)
391 2 : if (.not.MeshUseMeshFile) then
392 2 : write(iulog,'(a,f6.3,f8.2)') "Average equatorial node spacing (deg, km) = ", &
393 4 : dble(90)/dble(ne*(np-1)), PI*rearth/(2000.0_r8*dble(ne*(np-1)))
394 : end if
395 2 : write(iulog,'(a,2f9.3)') 'norm of Dinv (min, max): ', min_normDinv, max_normDinv
396 2 : write(iulog,'(a,1f8.2)') 'Max Dinv-based element distortion: ', max_ratio
397 2 : write(iulog,'(a,3f8.2)') 'dx based on Dinv svd: ave,min,max = ', avg_min_dx, min_min_dx, max_min_dx
398 2 : write(iulog,'(a,3f8.2)') "dx based on sqrt element area: ave,min,max = ", &
399 4 : sqrt(avg_area)/(np-1),sqrt(min_area)/(np-1),sqrt(max_area)/(np-1)
400 : end if
401 :
402 :
403 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
404 : ! SCALAR, RESOLUTION-AWARE HYPERVISCOSITY
405 : ! this block of code initializes the variable_hyperviscsoity() array
406 : ! based on largest length scale in each element and user specified scaling
407 : ! it then limits the coefficient if the user specifed a max CFL
408 : ! this limiting is based on the smallest length scale of each element
409 : ! since that controls the CFL.
410 : ! Mike Levy
411 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
412 1536 : if (hypervis_power /= 0) then
413 :
414 0 : min_hypervis = 1d99
415 0 : max_hypervis = 0
416 0 : avg_hypervis = 0
417 :
418 :
419 0 : max_unif_dx = min_max_dx ! use this for average resolution, unless:
420 : ! viscosity in namelist specified for smallest element:
421 0 : if (fine_ne>0) then
422 : ! viscosity in namelist specified for regions with a resolution
423 : ! equivilant to a uniform grid with ne=fine_ne
424 : if (np /= 4 ) call endrun('ERROR: setting fine_ne only supported with NP=4')
425 0 : max_unif_dx = (111.28_r8*30)/dble(fine_ne) ! in km
426 : endif
427 :
428 : !
429 : ! note: if L = eigenvalue of metinv, then associated length scale (km) is
430 : ! dx = 1.0_r8/( sqrt(L)*0.5_r8*dble(np-1)*ra*1000.0_r8)
431 : !
432 : ! for viscosity *tensor*, we take at each point:
433 : ! nu1 = nu*(dx1/max_unif_dx)**3.2 dx1 associated with eigenvalue 1
434 : ! nu2 = nu*(dx2/max_unif_dx)**3.2 dx2 associated with eigenvalue 2
435 : ! with this approach:
436 : ! - with this formula, no need to adjust for CFL violations
437 : ! - if nu comes from a 3.2 scaling that is stable for coarse and fine resolutions,
438 : ! this formulat will be stable.
439 : ! - gives the correct answer in long skinny rectangles:
440 : ! large viscosity in the long direction, small viscosity in the short direction
441 : !
442 0 : normDinv_hypervis = 0
443 0 : do ie=nets,nete
444 : ! variable viscosity based on map from ulatlon -> ucontra
445 :
446 : ! dx_long
447 0 : elem(ie)%variable_hyperviscosity = sqrt((elem(ie)%dx_long/max_unif_dx) ** hypervis_power)
448 : elem(ie)%hv_courant = dtnu*(elem(ie)%variable_hyperviscosity(1,1)**2) * &
449 0 : (lambda_vis**2) * ((ra*elem(ie)%normDinv)**4)
450 :
451 : ! Check to see if this is stable
452 0 : if (elem(ie)%hv_courant.gt.max_hypervis_courant) then
453 : stable_hv = sqrt( max_hypervis_courant / &
454 0 : ( dtnu * (lambda_vis)**2 * (ra*elem(ie)%normDinv)**4 ) )
455 :
456 : #if 0
457 : ! Useful print statements for debugging the adjustments to hypervis
458 : print*, "Adjusting hypervis on elem ", elem(ie)%GlobalId
459 : print*, "From ", nu*elem(ie)%variable_hyperviscosity(1,1)**2, " to ", nu*stable_hv
460 : print*, "Difference = ", nu*(/elem(ie)%variable_hyperviscosity(1,1)**2-stable_hv/)
461 : print*, "Factor of ", elem(ie)%variable_hyperviscosity(1,1)**2/stable_hv
462 : print*, " "
463 : #endif
464 : ! make sure that: elem(ie)%hv_courant <= max_hypervis_courant
465 0 : elem(ie)%variable_hyperviscosity = stable_hv
466 0 : elem(ie)%hv_courant = dtnu*(stable_hv**2) * (lambda_vis)**2 * (ra*elem(ie)%normDinv)**4
467 : end if
468 0 : normDinv_hypervis = max(normDinv_hypervis, elem(ie)%hv_courant/dtnu)
469 :
470 0 : min_hypervis = min(min_hypervis, elem(ie)%variable_hyperviscosity(1,1))
471 0 : max_hypervis = max(max_hypervis, elem(ie)%variable_hyperviscosity(1,1))
472 0 : global_shared_buf(ie,1) = elem(ie)%variable_hyperviscosity(1,1)
473 : end do
474 :
475 0 : min_hypervis = ParallelMin(min_hypervis, hybrid)
476 0 : max_hypervis = ParallelMax(max_hypervis, hybrid)
477 0 : call wrap_repro_sum(nvars=1, comm=hybrid%par%comm)
478 0 : avg_hypervis = global_shared_sum(1)/dble(nelem)
479 :
480 0 : normDinv_hypervis = ParallelMax(normDinv_hypervis, hybrid)
481 :
482 : ! apply DSS (aka assembly procedure) to variable_hyperviscosity (makes continuous)
483 0 : call initEdgeBuffer(hybrid%par,edgebuf,elem,1)
484 0 : do ie=nets,nete
485 0 : zeta(:,:,ie) = elem(ie)%variable_hyperviscosity(:,:)*elem(ie)%spheremp(:,:)
486 0 : call edgeVpack(edgebuf,zeta(1,1,ie),1,0,ie)
487 : end do
488 0 : call bndry_exchange(hybrid,edgebuf,location='print_cfl #1')
489 0 : do ie=nets,nete
490 0 : call edgeVunpack(edgebuf,zeta(1,1,ie),1,0,ie)
491 0 : elem(ie)%variable_hyperviscosity(:,:) = zeta(:,:,ie)*elem(ie)%rspheremp(:,:)
492 : end do
493 0 : call FreeEdgeBuffer(edgebuf)
494 :
495 : ! replace hypervis w/ bilinear based on continuous corner values
496 0 : do ie=nets,nete
497 0 : noreast = elem(ie)%variable_hyperviscosity(np,np)
498 0 : nw = elem(ie)%variable_hyperviscosity(1,np)
499 0 : se = elem(ie)%variable_hyperviscosity(np,1)
500 0 : sw = elem(ie)%variable_hyperviscosity(1,1)
501 0 : do i=1,np
502 0 : x = gp%points(i)
503 0 : do j=1,np
504 0 : y = gp%points(j)
505 0 : elem(ie)%variable_hyperviscosity(i,j) = 0.25_r8*( &
506 : (1.0_r8-x)*(1.0_r8-y)*sw + &
507 : (1.0_r8-x)*(y+1.0_r8)*nw + &
508 : (x+1.0_r8)*(1.0_r8-y)*se + &
509 0 : (x+1.0_r8)*(y+1.0_r8)*noreast)
510 : end do
511 : end do
512 : end do
513 1536 : else if (hypervis_scaling/=0) then
514 : ! tensorHV. New eigenvalues are the eigenvalues of the tensor V
515 : ! formulas here must match what is in cube_mod.F90
516 : ! for tensorHV, we scale out the rearth dependency
517 0 : lambda = max_normDinv**2
518 : normDinv_hypervis = (lambda_vis**2) * (max_normDinv**4) * &
519 0 : (lambda**(-hypervis_scaling/2) )
520 : else
521 : ! constant coefficient formula:
522 1536 : normDinv_hypervis = (lambda_vis**2) * (ra*max_normDinv)**4
523 : endif
524 :
525 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
526 : ! TENSOR, RESOLUTION-AWARE HYPERVISCOSITY
527 : ! The tensorVisc() array is computed in cube_mod.F90
528 : ! this block of code will DSS it so the tensor if C0
529 : ! and also make it bilinear in each element.
530 : ! Oksana Guba
531 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
532 1536 : if (hypervis_scaling /= 0) then
533 :
534 0 : call initEdgeBuffer(hybrid%par,edgebuf,elem,1)
535 0 : do rowind=1,2
536 0 : do colind=1,2
537 0 : do ie=nets,nete
538 0 : zeta(:,:,ie) = elem(ie)%tensorVisc(:,:,rowind,colind)*elem(ie)%spheremp(:,:)
539 0 : call edgeVpack(edgebuf,zeta(1,1,ie),1,0,ie)
540 : end do
541 :
542 0 : call bndry_exchange(hybrid,edgebuf)
543 0 : do ie=nets,nete
544 0 : call edgeVunpack(edgebuf,zeta(1,1,ie),1,0,ie)
545 0 : elem(ie)%tensorVisc(:,:,rowind,colind) = zeta(:,:,ie)*elem(ie)%rspheremp(:,:)
546 : end do
547 : enddo !rowind
548 : enddo !colind
549 0 : call FreeEdgeBuffer(edgebuf)
550 :
551 : !IF BILINEAR MAP OF V NEEDED
552 :
553 0 : do rowind=1,2
554 0 : do colind=1,2
555 : ! replace hypervis w/ bilinear based on continuous corner values
556 0 : do ie=nets,nete
557 0 : noreast = elem(ie)%tensorVisc(np,np,rowind,colind)
558 0 : nw = elem(ie)%tensorVisc(1,np,rowind,colind)
559 0 : se = elem(ie)%tensorVisc(np,1,rowind,colind)
560 0 : sw = elem(ie)%tensorVisc(1,1,rowind,colind)
561 0 : do i=1,np
562 0 : x = gp%points(i)
563 0 : do j=1,np
564 0 : y = gp%points(j)
565 0 : elem(ie)%tensorVisc(i,j,rowind,colind) = 0.25_r8*( &
566 : (1.0_r8-x)*(1.0_r8-y)*sw + &
567 : (1.0_r8-x)*(y+1.0_r8)*nw + &
568 : (x+1.0_r8)*(1.0_r8-y)*se + &
569 0 : (x+1.0_r8)*(y+1.0_r8)*noreast)
570 : end do
571 : end do
572 : end do
573 : enddo !rowind
574 : enddo !colind
575 : endif
576 1536 : deallocate(gp%points)
577 1536 : deallocate(gp%weights)
578 :
579 1536 : call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_p ,1.0_r8 ,'_p ')
580 1536 : call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu ,1.0_r8,' ')
581 1536 : call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_div,2.5_r8 ,'_div')
582 :
583 1536 : if (nu_q<0) nu_q = nu_p ! necessary for consistency
584 1536 : if (nu_t<0) nu_t = nu_p ! temperature damping is always equal to nu_p
585 :
586 41472 : nu_div_lev(:) = nu_div
587 41472 : nu_lev(:) = nu
588 41472 : nu_t_lev(:) = nu_p
589 :
590 : !
591 : ! sponge layer strength needed for stability depends on model top location
592 : !
593 1536 : top_000_032km = .false.
594 1536 : top_032_042km = .false.
595 1536 : top_042_090km = .false.
596 1536 : top_090_140km = .false.
597 1536 : top_140_600km = .false.
598 1536 : nu_set = sponge_del4_nu_fac < 0
599 1536 : div_set = sponge_del4_nu_div_fac < 0
600 1536 : lev_set = sponge_del4_lev < 0
601 1536 : if (ptop>1000.0_r8) then
602 : !
603 : ! low top; usually idealized test cases
604 : !
605 0 : top_000_032km = .true.
606 0 : if (hybrid%masterthread) write(iulog,* )"Model top damping configuration: top_000_032km"
607 1536 : else if (ptop>100.0_r8) then
608 : !
609 : ! CAM6 top (~225 Pa) or CAM7 low top
610 : !
611 1536 : top_032_042km = .true.
612 1536 : if (hybrid%masterthread) write(iulog,* )"Model top damping configuration: top_032_042km"
613 0 : else if (ptop>1e-1_r8) then
614 : !
615 : ! CAM7 top (~4.35e-1 Pa)
616 : !
617 0 : top_042_090km = .true.
618 0 : if (hybrid%masterthread) write(iulog,* )"Model top damping configuration: top_042_090km"
619 0 : else if (ptop>1E-4_r8) then
620 : !
621 : ! WACCM top (~4.5e-4 Pa)
622 : !
623 0 : top_090_140km = .true.
624 0 : if (hybrid%masterthread) write(iulog,* )"Model top damping configuration: top_090_140km"
625 : else
626 : !
627 : ! WACCM-x - geospace (~4e-7 Pa)
628 : !
629 0 : top_140_600km = .true.
630 0 : if (hybrid%masterthread) write(iulog,* )"Model top damping configuration: top_140_600km"
631 : end if
632 : !
633 : ! Logging text for sponge layer configuration
634 : !
635 1536 : if (hybrid%masterthread .and. (nu_set .or. div_set .or. lev_set)) then
636 2 : write(iulog,* )""
637 2 : write(iulog,* )"Sponge layer del4 coefficient defaults based on model top location:"
638 : end if
639 : !
640 : ! if user or namelist is not specifying sponge del4 settings here are best guesses (empirically determined)
641 : !
642 1536 : if (top_042_090km) then
643 0 : if (sponge_del4_lev <0) sponge_del4_lev = 4
644 0 : if (sponge_del4_nu_fac <0) sponge_del4_nu_fac = 3.375_r8 !max value without having to increase subcycling of div4
645 0 : if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 3.375_r8 !max value without having to increase subcycling of div4
646 1536 : else if (top_090_140km.or.top_140_600km) then ! defaults for waccm(x)
647 0 : if (sponge_del4_lev <0) sponge_del4_lev = 20
648 0 : if (sponge_del4_nu_fac <0) sponge_del4_nu_fac = 5.0_r8
649 0 : if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 10.0_r8
650 : else
651 1536 : if (sponge_del4_lev <0) sponge_del4_lev = 1
652 1536 : if (sponge_del4_nu_fac <0) sponge_del4_nu_fac = 1.0_r8
653 1536 : if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 1.0_r8
654 : end if
655 :
656 : ! set max wind speed for diagnostics
657 1536 : umax = 120.0_r8
658 1536 : if (top_042_090km) then
659 0 : umax = 240._r8
660 1536 : else if (top_090_140km) then
661 0 : umax = 300._r8
662 1536 : else if (top_140_600km) then
663 0 : umax = 800._r8
664 : end if
665 : !
666 : ! Log sponge layer configuration
667 : !
668 1536 : if (hybrid%masterthread) then
669 2 : if (nu_set) then
670 2 : write(iulog, '(a,e9.2)') ' sponge_del4_nu_fac = ',sponge_del4_nu_fac
671 : end if
672 :
673 2 : if (div_set) then
674 2 : write(iulog, '(a,e9.2)') ' sponge_del4_nu_div_fac = ',sponge_del4_nu_div_fac
675 : end if
676 :
677 2 : if (lev_set) then
678 2 : write(iulog, '(a,i0)') ' sponge_del4_lev = ',sponge_del4_lev
679 : end if
680 2 : write(iulog,* )""
681 : end if
682 :
683 1536 : nu_max = sponge_del4_nu_fac*nu_p
684 1536 : nu_div_max = sponge_del4_nu_div_fac*nu_p
685 41472 : do k=1,nlev
686 : ! Vertical profile from FV dycore (see Lauritzen et al. 2012 DOI:10.1177/1094342011410088)
687 39936 : scale1 = 0.5_r8*(1.0_r8+tanh(2.0_r8*log(pmid(sponge_del4_lev)/pmid(k))))
688 39936 : if (sponge_del4_nu_div_fac /= 1.0_r8) then
689 0 : nu_div_lev(k) = (1.0_r8-scale1)*nu_div+scale1*nu_div_max
690 : end if
691 41472 : if (sponge_del4_nu_fac /= 1.0_r8) then
692 0 : nu_lev(k) = (1.0_r8-scale1)*nu +scale1*nu_max
693 0 : nu_t_lev(k) = (1.0_r8-scale1)*nu_p +scale1*nu_max
694 : end if
695 : end do
696 :
697 1536 : if (hybrid%masterthread)then
698 2 : write(iulog,*) "z computed from barometric formula (using US std atmosphere)"
699 2 : call std_atm_height(pmid(:),z(:))
700 2 : write(iulog,*) "k,pmid_ref,z,nu_lev,nu_t_lev,nu_div_lev"
701 54 : do k=1,nlev
702 54 : write(iulog,'(i3,5e11.4)') k,pmid(k),z(k),nu_lev(k),nu_t_lev(k),nu_div_lev(k)
703 : end do
704 2 : if (nu_top>0) then
705 2 : write(iulog,*) ": ksponge_end = ",ksponge_end
706 2 : write(iulog,*) ": sponge layer Laplacian damping"
707 2 : write(iulog,*) "k, p, z, nu_scale_top, nu (actual Laplacian damping coefficient)"
708 :
709 8 : do k=1,ksponge_end
710 8 : write(iulog,'(i3,4e11.4)') k,pmid(k),z(k),nu_scale_top(k),nu_scale_top(k)*nu_top
711 : end do
712 : end if
713 : end if
714 :
715 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
716 : !
717 : ! time-step information
718 : !
719 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
720 : !
721 : ! S=time-step stability region (i.e. advection w/leapfrog: S=1, viscosity w/forward Euler: S=2)
722 : !
723 1536 : if (tstep_type==1) then
724 0 : S_rk = 2.0_r8
725 0 : rk_str = ' * RK2-SSP 3 stage (same as tracers)'
726 1536 : elseif (tstep_type==2) then
727 0 : S_rk = 2.0_r8
728 0 : rk_str = ' * classic RK3'
729 1536 : elseif (tstep_type==3) then
730 0 : S_rk = 2.0_r8
731 0 : rk_str = ' * Kinnmark&Gray RK4'
732 1536 : elseif (tstep_type==4) then
733 1536 : S_rk = 3.0_r8
734 1536 : rk_str = ' * Kinnmark&Gray RK3 5 stage (3rd order)'
735 : end if
736 1536 : if (hybrid%masterthread) then
737 2 : write(iulog,'(a,f12.8,a)') 'Model top is ',ptop,'Pa'
738 2 : write(iulog,'(a)') ' '
739 2 : write(iulog,'(a)') 'Timestepping methods used in dynamical core:'
740 2 : write(iulog,'(a)')
741 2 : write(iulog,*) rk_str
742 2 : write(iulog,'(a)') ' * Spectral-element advection uses SSP preservation RK3'
743 2 : write(iulog,'(a)') ' * Viscosity operators use forward Euler'
744 : end if
745 1536 : S_laplacian = 2.0_r8 !using forward Euler for sponge diffusion
746 1536 : S_hypervis = 2.0_r8 !using forward Euler for hyperviscosity
747 1536 : S_rk_tracer = 2.0_r8
748 :
749 1536 : ugw = 342.0_r8 !max gravity wave speed
750 :
751 1536 : dt_max_adv = S_rk/(umax*max_normDinv*lambda_max*ra)
752 1536 : dt_max_gw = S_rk/(ugw*max_normDinv*lambda_max*ra)
753 1536 : dt_max_tracer_se = S_rk_tracer*min_gw/(umax*max_normDinv*ra)
754 1536 : if (use_cslam) then
755 1536 : if (large_Courant_incr) then
756 1536 : dt_max_tracer_fvm = dble(nhe)*(4.0_r8*pi*Rearth/dble(4.0_r8*ne*nc))/umax
757 : else
758 0 : dt_max_tracer_fvm = dble(nhe)*(2.0_r8*pi*Rearth/dble(4.0_r8*ne*nc))/umax
759 : end if
760 : else
761 0 : dt_max_tracer_fvm = -1.0_r8
762 : end if
763 125952 : nu_max = MAX(MAXVAL(nu_div_lev(:)),MAXVAL(nu_lev(:)),MAXVAL(nu_t_lev(:)))
764 1536 : dt_max_hypervis = s_hypervis/(nu_max*normDinv_hypervis)
765 1536 : dt_max_hypervis_tracer = s_hypervis/(nu_q*normDinv_hypervis)
766 :
767 84480 : max_laplace = MAX(MAXVAL(nu_scale_top(:))*nu_top,MAXVAL(kmvis_ref(:)/rho_ref(:)))
768 43008 : max_laplace = MAX(max_laplace,MAXVAL(kmcnd_ref(:)/(cpair*rho_ref(:))))
769 1536 : dt_max_laplacian_top = 1.0_r8/(max_laplace*((ra*max_normDinv)**2)*lambda_vis)
770 :
771 1536 : if (hybrid%masterthread) then
772 2 : write(iulog,'(a,f10.2,a)') ' '
773 2 : write(iulog,'(a,f10.2,a)') 'Estimates for maximum stable and actual time-steps for different aspects of algorithm:'
774 2 : write(iulog,'(a,f12.8,a)') '(assume max wind is ',umax,'m/s)'
775 2 : write(iulog,'(a)') '(assume max gravity wave speed is 342m/s)'
776 2 : write(iulog,'(a,f10.2,a)') ' '
777 2 : write(iulog,'(a,f10.2,a,f10.2,a)') '* dt_dyn (time-stepping dycore ; u,v,T,dM) < ',&
778 4 : MIN(dt_max_adv,dt_max_gw),'s ',dt_dyn_actual,'s'
779 2 : if (dt_dyn_actual>MIN(dt_max_adv,dt_max_gw)) write(iulog,*) 'WARNING: dt_dyn theoretically unstable'
780 :
781 2 : write(iulog,'(a,f10.2,a,f10.2,a)') '* dt_dyn_vis (hyperviscosity) ; u,v,T,dM) < ',dt_max_hypervis,&
782 4 : 's ',dt_dyn_visco_actual,'s'
783 2 : if (dt_dyn_visco_actual>dt_max_hypervis) write(iulog,*) 'WARNING: dt_dyn_vis theoretically unstable'
784 2 : if (.not.use_cslam) then
785 0 : write(iulog,'(a,f10.2,a,f10.2,a)') '* dt_tracer_se (time-stepping tracers ; q ) < ',dt_max_tracer_se,'s ',&
786 0 : dt_tracer_se_actual,'s'
787 0 : if (dt_tracer_se_actual>dt_max_tracer_se) write(iulog,*) 'WARNING: dt_tracer_se theoretically unstable'
788 0 : write(iulog,'(a,f10.2,a,f10.2,a)') '* dt_tracer_vis (hyperviscosity tracers; q ) < ',dt_max_hypervis_tracer,'s',&
789 0 : dt_tracer_visco_actual,'s'
790 0 : if (dt_tracer_visco_actual>dt_max_hypervis_tracer) write(iulog,*) 'WARNING: dt_tracer_hypervis theoretically unstable'
791 : end if
792 2 : if (use_cslam) then
793 2 : write(iulog,'(a,f10.2,a,f10.2,a)') '* dt_tracer_fvm (time-stepping tracers ; q ) < ',dt_max_tracer_fvm,&
794 4 : 's ',dt_tracer_fvm_actual
795 2 : if (dt_tracer_fvm_actual>dt_max_tracer_fvm) write(iulog,*) 'WARNING: dt_tracer_fvm theortically unstable'
796 : end if
797 2 : write(iulog,'(a,f10.2)') '* dt_remap (vertical remap dt) ',dt_remap_actual
798 8 : do k=1,ksponge_end
799 6 : max_laplace = MAX(nu_scale_top(k)*nu_top,kmvis_ref(k)/rho_ref(k))
800 6 : max_laplace = MAX(max_laplace,kmcnd_ref(k)/(cpair*rho_ref(k)))
801 6 : dt_max_laplacian_top = 1.0_r8/(max_laplace*((ra*max_normDinv)**2)*lambda_vis)
802 :
803 6 : write(iulog,'(a,f10.2,a,f10.2,a)') '* dt (del2 sponge ; u,v,T,dM) < ',&
804 12 : dt_max_laplacian_top,'s',dt_dyn_del2_actual,'s'
805 8 : if (dt_dyn_del2_actual>dt_max_laplacian_top) then
806 0 : if (k==1) then
807 0 : write(iulog,*) 'WARNING: theoretically unstable in sponge; increase se_hypervis_subcycle_sponge',&
808 0 : ' (this WARNING can sometimes be ignored in level 1)'
809 : else
810 0 : write(iulog,*) 'WARNING: theoretically unstable in sponge; increase se_hypervis_subcycle_sponge'
811 : endif
812 : end if
813 : end do
814 2 : write(iulog,*) ' '
815 2 : if (hypervis_power /= 0) then
816 0 : write(iulog,'(a,3e11.4)')'Scalar hyperviscosity (dynamics): ave,min,max = ', &
817 0 : nu*(/avg_hypervis**2,min_hypervis**2,max_hypervis**2/)
818 : end if
819 2 : write(iulog,*) 'tstep_type = ',tstep_type
820 : end if
821 1536 : end subroutine print_cfl
822 :
823 : !
824 : ! ============================
825 : ! global_maximum:
826 : !
827 : ! Find global maximum on sphere
828 : !
829 : ! ================================
830 :
831 0 : function global_maximum(fld,hybrid,npts,nets,nete) result(Max_sphere)
832 :
833 1536 : use hybrid_mod, only : hybrid_t
834 : use reduction_mod, only : red_max, pmax_mt
835 :
836 : integer , intent(in) :: npts,nets,nete
837 : real (kind=r8), intent(in) :: fld(npts,npts,nets:nete)
838 : type (hybrid_t) , intent(in) :: hybrid
839 :
840 : real (kind=r8) :: Max_sphere
841 :
842 : ! Local variables
843 :
844 : real (kind=r8) :: redp(1)
845 :
846 0 : Max_sphere = MAXVAL(fld(:,:,nets:nete))
847 :
848 0 : redp(1) = Max_sphere
849 0 : call pmax_mt(red_max,redp,1,hybrid)
850 0 : Max_sphere = red_max%buf(1)
851 :
852 0 : end function global_maximum
853 :
854 : ! ==========================================================
855 : ! l1_snorm:
856 : !
857 : ! computes the l1 norm per Williamson et al, p. 218 eq(8)
858 : ! for a scalar quantity
859 : ! ===========================================================
860 :
861 0 : function l1_snorm(elem,fld,fld_exact,hybrid,npts,nets,nete) result(l1)
862 :
863 : use element_mod, only : element_t
864 : use hybrid_mod, only : hybrid_t
865 :
866 : type(element_t) , intent(in) :: elem(:)
867 : integer , intent(in) :: npts,nets,nete
868 : real (kind=r8), intent(in) :: fld(npts,npts,nets:nete) ! computed soln
869 : real (kind=r8), intent(in) :: fld_exact(npts,npts,nets:nete) ! true soln
870 : type (hybrid_t) , intent(in) :: hybrid
871 : real (kind=r8) :: l1
872 :
873 : ! Local variables
874 :
875 0 : real (kind=r8) :: dfld_abs(npts,npts,nets:nete)
876 0 : real (kind=r8) :: fld_exact_abs(npts,npts,nets:nete)
877 : real (kind=r8) :: dfld_abs_int
878 : real (kind=r8) :: fld_exact_abs_int
879 : integer i,j,ie
880 :
881 0 : do ie=nets,nete
882 0 : do j=1,npts
883 0 : do i=1,npts
884 0 : dfld_abs(i,j,ie) = ABS(fld(i,j,ie)-fld_exact(i,j,ie))
885 0 : fld_exact_abs(i,j,ie) = ABS(fld_exact(i,j,ie))
886 : end do
887 : end do
888 : end do
889 :
890 0 : dfld_abs_int = global_integral(elem, dfld_abs(:,:,nets:nete),hybrid,npts,nets,nete)
891 0 : fld_exact_abs_int = global_integral(elem, fld_exact_abs(:,:,nets:nete),hybrid,npts,nets,nete)
892 :
893 0 : l1 = dfld_abs_int/fld_exact_abs_int
894 :
895 0 : end function l1_snorm
896 :
897 : ! ===========================================================
898 : ! l1_vnorm:
899 : !
900 : ! computes the l1 norm per Williamson et al, p. 218 eq(97),
901 : ! for a contravariant vector quantity on the velocity grid.
902 : !
903 : ! ===========================================================
904 :
905 0 : function l1_vnorm(elem, v,vt,hybrid,npts,nets,nete) result(l1)
906 : use element_mod, only : element_t
907 : use hybrid_mod, only : hybrid_t
908 :
909 : type(element_t) , intent(in), target :: elem(:)
910 : integer , intent(in) :: npts,nets,nete
911 : real (kind=r8), intent(in) :: v(npts,npts,2,nets:nete) ! computed soln
912 : real (kind=r8), intent(in) :: vt(npts,npts,2,nets:nete) ! true soln
913 : type (hybrid_t) , intent(in) :: hybrid
914 : real (kind=r8) :: l1
915 :
916 : ! Local variables
917 :
918 0 : real (kind=r8), dimension(:,:,:,:), pointer :: met
919 0 : real (kind=r8) :: dvsq(npts,npts,nets:nete)
920 0 : real (kind=r8) :: vtsq(npts,npts,nets:nete)
921 0 : real (kind=r8) :: dvco(npts,npts,2) ! covariant velocity
922 0 : real (kind=r8) :: vtco(npts,npts,2) ! covariant velocity
923 : real (kind=r8) :: dv1,dv2
924 : real (kind=r8) :: vt1,vt2
925 : real (kind=r8) :: dvsq_int
926 : real (kind=r8) :: vtsq_int
927 :
928 : integer i,j,ie
929 :
930 0 : do ie=nets,nete
931 0 : met => elem(ie)%met
932 0 : do j=1,npts
933 0 : do i=1,npts
934 :
935 0 : dv1 = v(i,j,1,ie)-vt(i,j,1,ie)
936 0 : dv2 = v(i,j,2,ie)-vt(i,j,2,ie)
937 :
938 0 : vt1 = vt(i,j,1,ie)
939 0 : vt2 = vt(i,j,2,ie)
940 :
941 0 : dvco(i,j,1) = met(i,j,1,1)*dv1 + met(i,j,1,2)*dv2
942 0 : dvco(i,j,2) = met(i,j,2,1)*dv1 + met(i,j,2,2)*dv2
943 :
944 0 : vtco(i,j,1) = met(i,j,1,1)*vt1 + met(i,j,1,2)*vt2
945 0 : vtco(i,j,2) = met(i,j,2,1)*vt1 + met(i,j,2,2)*vt2
946 :
947 0 : dvsq(i,j,ie) = SQRT(dvco(i,j,1)*dv1 + dvco(i,j,2)*dv2)
948 0 : vtsq(i,j,ie) = SQRT(vtco(i,j,1)*vt1 + vtco(i,j,2)*vt2)
949 :
950 : end do
951 : end do
952 : end do
953 :
954 0 : dvsq_int = global_integral(elem, dvsq(:,:,nets:nete),hybrid,npts,nets,nete)
955 0 : vtsq_int = global_integral(elem, vtsq(:,:,nets:nete),hybrid,npts,nets,nete)
956 :
957 0 : l1 = dvsq_int/vtsq_int
958 :
959 0 : end function l1_vnorm
960 :
961 : ! ==========================================================
962 : ! l2_snorm:
963 : !
964 : ! computes the l2 norm per Williamson et al, p. 218 eq(83)
965 : ! for a scalar quantity on the pressure grid.
966 : !
967 : ! ===========================================================
968 :
969 0 : function l2_snorm(elem,fld,fld_exact,hybrid,npts,nets,nete) result(l2)
970 : use element_mod, only : element_t
971 : use hybrid_mod, only : hybrid_t
972 :
973 : type(element_t), intent(in) :: elem(:)
974 : integer , intent(in) :: npts,nets,nete
975 : real (kind=r8), intent(in) :: fld(npts,npts,nets:nete) ! computed soln
976 : real (kind=r8), intent(in) :: fld_exact(npts,npts,nets:nete) ! true soln
977 : type (hybrid_t) , intent(in) :: hybrid
978 : real (kind=r8) :: l2
979 :
980 : ! Local variables
981 :
982 0 : real (kind=r8) :: dh2(npts,npts,nets:nete)
983 0 : real (kind=r8) :: fld_exact2(npts,npts,nets:nete)
984 : real (kind=r8) :: dh2_int
985 : real (kind=r8) :: fld_exact2_int
986 : integer i,j,ie
987 :
988 0 : do ie=nets,nete
989 0 : do j=1,npts
990 0 : do i=1,npts
991 0 : dh2(i,j,ie)=(fld(i,j,ie)-fld_exact(i,j,ie))**2
992 0 : fld_exact2(i,j,ie)=fld_exact(i,j,ie)**2
993 : end do
994 : end do
995 : end do
996 :
997 0 : dh2_int = global_integral(elem,dh2(:,:,nets:nete),hybrid,npts,nets,nete)
998 0 : fld_exact2_int = global_integral(elem,fld_exact2(:,:,nets:nete),hybrid,npts,nets,nete)
999 :
1000 0 : l2 = SQRT(dh2_int)/SQRT(fld_exact2_int)
1001 :
1002 0 : end function l2_snorm
1003 :
1004 : ! ==========================================================
1005 : ! l2_vnorm:
1006 : !
1007 : ! computes the l2 norm per Williamson et al, p. 219 eq(98)
1008 : ! for a contravariant vector quantity on the velocity grid.
1009 : !
1010 : ! ===========================================================
1011 :
1012 0 : function l2_vnorm(elem, v,vt,hybrid,npts,nets,nete) result(l2)
1013 : use element_mod, only : element_t
1014 : use hybrid_mod, only : hybrid_t
1015 :
1016 : type(element_t) , intent(in), target :: elem(:)
1017 : integer , intent(in) :: npts,nets,nete
1018 : real (kind=r8), intent(in) :: v(npts,npts,2,nets:nete) ! computed soln
1019 : real (kind=r8), intent(in) :: vt(npts,npts,2,nets:nete) ! true soln
1020 : type (hybrid_t) , intent(in) :: hybrid
1021 : real (kind=r8) :: l2
1022 :
1023 : ! Local variables
1024 :
1025 0 : real (kind=r8), dimension(:,:,:,:), pointer :: met
1026 0 : real (kind=r8) :: dvsq(npts,npts,nets:nete)
1027 0 : real (kind=r8) :: vtsq(npts,npts,nets:nete)
1028 0 : real (kind=r8) :: dvco(npts,npts,2) ! covariant velocity
1029 0 : real (kind=r8) :: vtco(npts,npts,2) ! covariant velocity
1030 : real (kind=r8) :: dv1,dv2
1031 : real (kind=r8) :: vt1,vt2
1032 : real (kind=r8) :: dvsq_int
1033 : real (kind=r8) :: vtsq_int
1034 : integer i,j,ie
1035 :
1036 0 : do ie=nets,nete
1037 0 : met => elem(ie)%met
1038 0 : do j=1,npts
1039 0 : do i=1,npts
1040 :
1041 0 : dv1 = v(i,j,1,ie)-vt(i,j,1,ie)
1042 0 : dv2 = v(i,j,2,ie)-vt(i,j,2,ie)
1043 :
1044 0 : vt1 = vt(i,j,1,ie)
1045 0 : vt2 = vt(i,j,2,ie)
1046 :
1047 0 : dvco(i,j,1) = met(i,j,1,1)*dv1 + met(i,j,1,2)*dv2
1048 0 : dvco(i,j,2) = met(i,j,2,1)*dv1 + met(i,j,2,2)*dv2
1049 :
1050 0 : vtco(i,j,1) = met(i,j,1,1)*vt1 + met(i,j,1,2)*vt2
1051 0 : vtco(i,j,2) = met(i,j,2,1)*vt1 + met(i,j,2,2)*vt2
1052 :
1053 0 : dvsq(i,j,ie) = dvco(i,j,1)*dv1 + dvco(i,j,2)*dv2
1054 0 : vtsq(i,j,ie) = vtco(i,j,1)*vt1 + vtco(i,j,2)*vt2
1055 :
1056 : end do
1057 : end do
1058 : end do
1059 :
1060 0 : dvsq_int = global_integral(elem, dvsq(:,:,nets:nete),hybrid,npts,nets,nete)
1061 0 : vtsq_int = global_integral(elem, vtsq(:,:,nets:nete),hybrid,npts,nets,nete)
1062 :
1063 0 : l2 = SQRT(dvsq_int)/SQRT(vtsq_int)
1064 :
1065 0 : end function l2_vnorm
1066 :
1067 : ! ==========================================================
1068 : ! linf_snorm:
1069 : !
1070 : ! computes the l infinity norm per Williamson et al, p. 218 eq(84)
1071 : ! for a scalar quantity on the pressure grid...
1072 : !
1073 : ! ===========================================================
1074 :
1075 0 : function linf_snorm(fld,fld_exact,hybrid,npts,nets,nete) result(linf)
1076 : use hybrid_mod, only : hybrid_t
1077 : integer , intent(in) :: npts,nets,nete
1078 : real (kind=r8), intent(in) :: fld(npts,npts,nets:nete) ! computed soln
1079 : real (kind=r8), intent(in) :: fld_exact(npts,npts,nets:nete) ! true soln
1080 : type (hybrid_t) , intent(in) :: hybrid
1081 : real (kind=r8) :: linf
1082 :
1083 : ! Local variables
1084 :
1085 0 : real (kind=r8) :: dfld_abs(npts,npts,nets:nete)
1086 0 : real (kind=r8) :: fld_exact_abs(npts,npts,nets:nete)
1087 : real (kind=r8) :: dfld_abs_max
1088 : real (kind=r8) :: fld_exact_abs_max
1089 : integer i,j,ie
1090 :
1091 0 : do ie=nets,nete
1092 0 : do j=1,npts
1093 0 : do i=1,npts
1094 0 : dfld_abs(i,j,ie)=ABS(fld(i,j,ie)-fld_exact(i,j,ie))
1095 0 : fld_exact_abs(i,j,ie)=ABS(fld_exact(i,j,ie))
1096 : end do
1097 : end do
1098 : end do
1099 :
1100 0 : dfld_abs_max = global_maximum(dfld_abs(:,:,nets:nete),hybrid,npts,nets,nete)
1101 0 : fld_exact_abs_max = global_maximum(fld_exact_abs(:,:,nets:nete),hybrid,npts,nets,nete)
1102 :
1103 0 : linf = dfld_abs_max/fld_exact_abs_max
1104 :
1105 0 : end function linf_snorm
1106 :
1107 :
1108 : ! ==========================================================
1109 : ! linf_vnorm:
1110 : !
1111 : ! computes the linf norm per Williamson et al, p. 218 eq(99),
1112 : ! for a contravariant vector quantity on the velocity grid.
1113 : !
1114 : ! ===========================================================
1115 :
1116 0 : function linf_vnorm(elem,v,vt,hybrid,npts,nets,nete) result(linf)
1117 : use hybrid_mod, only : hybrid_t
1118 : use element_mod, only : element_t
1119 :
1120 : type(element_t) , intent(in), target :: elem(:)
1121 : integer , intent(in) :: npts,nets,nete
1122 : real (kind=r8), intent(in) :: v(npts,npts,2,nets:nete) ! computed soln
1123 : real (kind=r8), intent(in) :: vt(npts,npts,2,nets:nete) ! true soln
1124 : type (hybrid_t) , intent(in) :: hybrid
1125 : real (kind=r8) :: linf
1126 :
1127 : ! Local variables
1128 :
1129 0 : real (kind=r8), dimension(:,:,:,:), pointer :: met
1130 0 : real (kind=r8) :: dvsq(npts,npts,nets:nete)
1131 0 : real (kind=r8) :: vtsq(npts,npts,nets:nete)
1132 0 : real (kind=r8) :: dvco(npts,npts,2) ! covariant velocity
1133 0 : real (kind=r8) :: vtco(npts,npts,2) ! covariant velocity
1134 : real (kind=r8) :: dv1,dv2
1135 : real (kind=r8) :: vt1,vt2
1136 : real (kind=r8) :: dvsq_max
1137 : real (kind=r8) :: vtsq_max
1138 : integer i,j,ie
1139 :
1140 0 : do ie=nets,nete
1141 0 : met => elem(ie)%met
1142 :
1143 0 : do j=1,npts
1144 0 : do i=1,npts
1145 :
1146 0 : dv1 = v(i,j,1,ie)-vt(i,j,1,ie)
1147 0 : dv2 = v(i,j,2,ie)-vt(i,j,2,ie)
1148 :
1149 0 : vt1 = vt(i,j,1,ie)
1150 0 : vt2 = vt(i,j,2,ie)
1151 :
1152 0 : dvco(i,j,1) = met(i,j,1,1)*dv1 + met(i,j,1,2)*dv2
1153 0 : dvco(i,j,2) = met(i,j,2,1)*dv1 + met(i,j,2,2)*dv2
1154 :
1155 0 : vtco(i,j,1) = met(i,j,1,1)*vt1 + met(i,j,1,2)*vt2
1156 0 : vtco(i,j,2) = met(i,j,2,1)*vt1 + met(i,j,2,2)*vt2
1157 :
1158 0 : dvsq(i,j,ie) = SQRT(dvco(i,j,1)*dv1 + dvco(i,j,2)*dv2)
1159 0 : vtsq(i,j,ie) = SQRT(vtco(i,j,1)*vt1 + vtco(i,j,2)*vt2)
1160 :
1161 : end do
1162 : end do
1163 : end do
1164 :
1165 0 : dvsq_max = global_maximum(dvsq(:,:,nets:nete),hybrid,npts,nets,nete)
1166 0 : vtsq_max = global_maximum(vtsq(:,:,nets:nete),hybrid,npts,nets,nete)
1167 :
1168 0 : linf = dvsq_max/vtsq_max
1169 :
1170 0 : end function linf_vnorm
1171 :
1172 6912 : subroutine wrap_repro_sum (nvars, comm, nsize)
1173 : use dimensions_mod, only: nelemd
1174 : use shr_reprosum_mod, only: repro_sum => shr_reprosum_calc
1175 : use cam_abortutils, only: endrun
1176 : use parallel_mod, only: global_shared_buf, global_shared_sum, nrepro_vars
1177 :
1178 : integer :: nvars ! number of variables to be summed (cannot exceed nrepro_vars)
1179 : integer :: comm ! mpi communicator
1180 : integer, optional :: nsize ! local buffer size (defaults to nelemd - number of elements in mpi task)
1181 :
1182 : integer nsize_use
1183 :
1184 6912 : if (present(nsize)) then
1185 0 : nsize_use = nsize
1186 : else
1187 6912 : nsize_use = nelemd
1188 : endif
1189 6912 : if (nvars .gt. nrepro_vars) call endrun('ERROR: repro_sum_buffer_size exceeded')
1190 :
1191 : ! Repro_sum contains its own OpenMP, so only one thread should call it (AAM)
1192 :
1193 : !$OMP BARRIER
1194 : !$OMP MASTER
1195 :
1196 6912 : call repro_sum(global_shared_buf, global_shared_sum, nsize_use, nelemd, nvars, commid=comm)
1197 :
1198 :
1199 : !$OMP END MASTER
1200 : !$OMP BARRIER
1201 :
1202 6912 : end subroutine wrap_repro_sum
1203 :
1204 4608 : subroutine automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu,factor,str)
1205 6912 : use physconst, only: rearth
1206 : use control_mod, only: hypervis_scaling,hypervis_power
1207 : use hybrid_mod, only: hybrid_t
1208 : use cam_abortutils, only: endrun
1209 :
1210 : type (hybrid_t), intent(in) :: hybrid
1211 : integer , intent(in) :: ne
1212 : real (kind=r8), intent(in) :: max_min_dx,min_min_dx,factor
1213 : real (kind=r8), intent(inout) :: nu
1214 : character(len=4), intent(in) :: str
1215 :
1216 : real(r8) :: uniform_res_hypervis_scaling,nu_fac
1217 : real(kind=r8) :: nu_min, nu_max
1218 : !
1219 : !************************************************************************************************************
1220 : !
1221 : ! automatically set viscosity coefficients
1222 : !
1223 : !
1224 : ! Use scaling from
1225 : !
1226 : ! - Boville, B. A., 1991: Sensitivity of simulated climate to
1227 : ! model resolution. J. Climate, 4, 469-485.
1228 : !
1229 : ! - TAKAHASHI ET AL., 2006: GLOBAL SIMULATION OF MESOSCALE SPECTRUM
1230 : !
1231 4608 : uniform_res_hypervis_scaling = 1.0_r8/log10(2.0_r8)
1232 : !
1233 : ! compute factor so that at ne30 resolution nu=1E15
1234 : ! scale so that scaling works for other planets
1235 : !
1236 : ! grid spacing in meters = max_min_dx*1000.0_r8
1237 : !
1238 4608 : nu_fac = (rearth/6.37122E6_r8)*1.0E15_r8/(110000.0_r8**uniform_res_hypervis_scaling)
1239 :
1240 4608 : if (nu < 0) then
1241 4608 : if (ne <= 0) then
1242 0 : if (hypervis_power/=0) then
1243 0 : call endrun('ERROR: Automatic scaling of scalar viscosity not implemented')
1244 0 : else if (hypervis_scaling/=0) then
1245 0 : nu_min = factor*nu_fac*(max_min_dx*1000.0_r8)**uniform_res_hypervis_scaling
1246 0 : nu_max = factor*nu_fac*(min_min_dx*1000.0_r8)**uniform_res_hypervis_scaling
1247 0 : nu = factor*nu_min
1248 0 : if (hybrid%masterthread) then
1249 0 : write(iulog,'(a,a)') "Automatically setting nu",TRIM(str)
1250 0 : write(iulog,'(a,2e9.2,a,2f9.2)') "Value at min/max grid spacing: ",nu_min,nu_max,&
1251 0 : " Max/min grid spacing (km) = ",max_min_dx,min_min_dx
1252 : end if
1253 0 : nu = nu_min*(2.0_r8*rearth/(3.0_r8*max_min_dx*1000.0_r8))**hypervis_scaling/(rearth**4)
1254 0 : if (hybrid%masterthread) &
1255 0 : write(iulog,'(a,a,a,e9.3)') "Nu_tensor",TRIM(str)," = ",nu
1256 : end if
1257 : else
1258 4608 : nu = factor*nu_fac*((30.0_r8/ne)*110000.0_r8)**uniform_res_hypervis_scaling
1259 4608 : if (hybrid%masterthread) then
1260 6 : write(iulog,'(a,a,a,e9.2)') "Automatically setting nu",TRIM(str)," =",nu
1261 : end if
1262 : end if
1263 : end if
1264 4608 : end subroutine automatically_set_viscosity_coefficients
1265 : end module global_norms_mod
|
