LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - global_norms_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 233 | 484 | 48.1 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 5 | 13 | 38.5 %

              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      144384 :     nu_div_lev(:) = nu_div
         587      144384 :     nu_lev(:)     = nu
         588      144384 :     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           0 :       top_032_042km = .true.
         612           0 :       if (hybrid%masterthread) write(iulog,* )"Model top damping configuration: top_032_042km"
         613        1536 :     else if (ptop>1e-1_r8) then
         614             :       !
         615             :       ! CAM7 top (~4.35e-1 Pa)
         616             :       !
         617        1536 :       top_042_090km = .true.
         618        1536 :       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        1536 :        if (sponge_del4_lev       <0) sponge_del4_lev        = 4
         644        1536 :        if (sponge_del4_nu_fac    <0) sponge_del4_nu_fac     = 3.375_r8 !max value without having to increase subcycling of div4
         645        1536 :        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           0 :     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           0 :       if (sponge_del4_lev       <0) sponge_del4_lev        = 1
         652           0 :       if (sponge_del4_nu_fac    <0) sponge_del4_nu_fac     = 1.0_r8
         653           0 :       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        1536 :        umax = 240._r8
         660           0 :     else if (top_090_140km) then
         661           0 :        umax = 300._r8
         662           0 :     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      144384 :     do k=1,nlev
         686             :       ! Vertical profile from FV dycore (see Lauritzen et al. 2012 DOI:10.1177/1094342011410088)
         687      142848 :       scale1        = 0.5_r8*(1.0_r8+tanh(2.0_r8*log(pmid(sponge_del4_lev)/pmid(k))))
         688      142848 :       if (sponge_del4_nu_div_fac /= 1.0_r8) then
         689      142848 :         nu_div_lev(k) = (1.0_r8-scale1)*nu_div+scale1*nu_div_max
         690             :       end if
         691      144384 :       if (sponge_del4_nu_fac /= 1.0_r8) then
         692      142848 :         nu_lev(k)     = (1.0_r8-scale1)*nu    +scale1*nu_max
         693      142848 :         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         188 :       do k=1,nlev
         702         188 :         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          10 :         do k=1,ksponge_end
         710          10 :            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      434688 :     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      290304 :     max_laplace = MAX(MAXVAL(nu_scale_top(:))*nu_top,MAXVAL(kmvis_ref(:)/rho_ref(:)))
         768      145920 :     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          10 :       do k=1,ksponge_end
         799           8 :         max_laplace = MAX(nu_scale_top(k)*nu_top,kmvis_ref(k)/rho_ref(k))
         800           8 :         max_laplace = MAX(max_laplace,kmcnd_ref(k)/(cpair*rho_ref(k)))
         801           8 :         dt_max_laplacian_top   = 1.0_r8/(max_laplace*((ra*max_normDinv)**2)*lambda_vis)
         802             :
         803           8 :         write(iulog,'(a,f10.2,a,f10.2,a)') '* dt    (del2 sponge           ; u,v,T,dM) < ',&
         804          16 :              dt_max_laplacian_top,'s',dt_dyn_del2_actual,'s'
         805          10 :         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

Generated by: LCOV version 1.14