LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - viscosity_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 72 | 210 | 34.3 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 2 | 14 | 14.3 %

              Line data    Source code
           1             : module viscosity_mod
           2             : !
           3             : !  This module should be renamed "global_deriv_mod.F90"
           4             : !
           5             : !  It is a collection of derivative operators that must be applied to the field
           6             : !  over the sphere (as opposed to derivative operators that can be applied element
           7             : !  by element)
           8             : !
           9             : !
          10             :   use shr_kind_mod,   only: r8=>shr_kind_r8
          11             :   use thread_mod,     only: max_num_threads, omp_get_num_threads
          12             :   use dimensions_mod, only: np, nc, nlev,nlevp, qsize,nelemd
          13             :   use hybrid_mod,     only: hybrid_t, get_loop_ranges, config_thread_region
          14             :   use parallel_mod,   only: parallel_t
          15             :   use element_mod,    only: element_t
          16             :   use derivative_mod, only: derivative_t, laplace_sphere_wk, vlaplace_sphere_wk, vorticity_sphere, derivinit, divergence_sphere
          17             :   use edgetype_mod,   only: EdgeBuffer_t, EdgeDescriptor_t
          18             :   use edge_mod,       only: edgevpack, edgevunpack, edgeVunpackmin, edgeSunpackmin, &
          19             :        edgeVunpackmax, initEdgeBuffer, FreeEdgeBuffer, edgeSunpackmax, edgeSpack
          20             :   use bndry_mod,      only: bndry_exchange, bndry_exchange_start,bndry_exchange_finish
          21             :   use control_mod,    only: hypervis_scaling, nu, nu_div
          22             :   use thread_mod,     only: vert_num_threads
          23             :
          24             :   implicit none
          25             :   save
          26             :
          27             :   public :: biharmonic_wk_scalar
          28             :   public :: biharmonic_wk_omega
          29             :   public :: neighbor_minmax, neighbor_minmax_start,neighbor_minmax_finish
          30             :
          31             :   !
          32             :   ! compute vorticity/divergence and then project to make continious
          33             :   ! high-level routines uses only for I/O
          34             :   public :: compute_zeta_C0
          35             :   public :: compute_div_C0
          36             :
          37             :   interface compute_zeta_C0
          38             :     module procedure compute_zeta_C0_hybrid       ! hybrid version
          39             :     module procedure compute_zeta_C0_par          ! single threaded
          40             :   end interface compute_zeta_C0
          41             :   interface compute_div_C0
          42             :     module procedure compute_div_C0_hybrid
          43             :     module procedure compute_div_C0_par
          44             :   end interface compute_div_C0
          45             :
          46             :   public :: compute_zeta_C0_contra    ! for older versions of sweq which carry
          47             :   public :: compute_div_C0_contra     ! velocity around in contra-coordinates
          48             :
          49             :   type (EdgeBuffer_t)          :: edge1
          50             :
          51             : CONTAINS
          52             :
          53     6649344 : subroutine biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid,nt,nets,nete,kbeg,kend)
          54             :   use derivative_mod, only : subcell_Laplace_fluxes
          55             :   use dimensions_mod, only : use_cslam, nu_div_lev,nu_lev
          56             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          57             :   ! compute weak biharmonic operator
          58             :   !    input:  h,v (stored in elem()%, in lat-lon coordinates
          59             :   !    output: ttens,vtens  overwritten with weak biharmonic of h,v (output in lat-lon coordinates)
          60             :   !
          61             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          62             :   type (hybrid_t)      , intent(in) :: hybrid
          63             :   type (element_t)     , intent(inout), target :: elem(:)
          64             :   integer              , intent(in)  :: nt,nets,nete
          65             :   integer              , intent(in)  :: kbeg, kend
          66             :   real (kind=r8), intent(out), dimension(nc,nc,4,nlev,nets:nete) :: dpflux
          67             :   real (kind=r8), dimension(np,np,2,nlev,nets:nete)  :: vtens
          68             :   real (kind=r8), dimension(np,np,nlev,nets:nete) :: ttens,dptens
          69             :   type (EdgeBuffer_t)  , intent(inout) :: edge3
          70             :   type (derivative_t)  , intent(in) :: deriv
          71             :   ! local
          72             :   integer :: i,j,k,kptr,ie,kblk
          73             : !  real (kind=r8), dimension(:,:), pointer :: rspheremv
          74             :   real (kind=r8), dimension(np,np) :: tmp
          75             :   real (kind=r8), dimension(np,np) :: tmp2
          76             :   real (kind=r8), dimension(np,np,2) :: v
          77             :
          78             :   real (kind=r8), dimension(np,np,nlev) :: lap_p_wk
          79             :   real (kind=r8), dimension(np,np,nlevp) :: T_i
          80             :
          81             :
          82             :   real (kind=r8) :: nu_ratio1, nu_ratio2, dp_thresh
          83             :   logical var_coef1
          84             :
          85     6649344 :   kblk = kend - kbeg + 1
          86             :
          87 >23049*10^7 :   if (use_cslam) dpflux = 0
          88             :   !if tensor hyperviscosity with tensor V is used, then biharmonic operator is (\grad\cdot V\grad) (\grad \cdot \grad)
          89             :   !so tensor is only used on second call to laplace_sphere_wk
          90     6649344 :   var_coef1 = .true.
          91     6649344 :   if(hypervis_scaling > 0)    var_coef1 = .false.
          92     6649344 :   dp_thresh=.025_r8  ! tunable coefficient
          93    53402544 :   do ie=nets,nete
          94             : !$omp parallel do num_threads(vert_num_threads) private(k,tmp)
          95  4394800800 :     do k=kbeg,kend
          96  4348047600 :       nu_ratio1=1
          97  4348047600 :       nu_ratio2=1
          98  4348047600 :       if (nu_div_lev(k)/=nu_lev(k)) then
          99  4348047600 :         if(hypervis_scaling /= 0) then
         100             :           ! we have a problem with the tensor in that we cant seperate
         101             :           ! div and curl components.  So we do, with tensor V:
         102             :           ! nu * (del V del ) * ( nu_ratio * grad(div) - curl(curl))
         103           0 :           nu_ratio1=nu_div_lev(k)/nu_lev(k)
         104             :           nu_ratio2=1
         105             :         else
         106  4348047600 :           nu_ratio1=sqrt(nu_div_lev(k)/nu_lev(k))
         107  4348047600 :           nu_ratio2=sqrt(nu_div_lev(k)/nu_lev(k))
         108             :         endif
         109             :       endif
         110             :
         111 91308999600 :       tmp=elem(ie)%state%T(:,:,k,nt)-elem(ie)%derived%T_ref(:,:,k)
         112  4348047600 :       call laplace_sphere_wk(tmp,deriv,elem(ie),ttens(:,:,k,ie),var_coef=var_coef1)
         113             :
         114 91308999600 :       tmp=elem(ie)%state%dp3d(:,:,k,nt)-elem(ie)%derived%dp_ref(:,:,k)
         115  4348047600 :       call laplace_sphere_wk(tmp,deriv,elem(ie),dptens(:,:,k,ie),var_coef=var_coef1)
         116             :
         117  4348047600 :       call vlaplace_sphere_wk(elem(ie)%state%v(:,:,:,k,nt),deriv,elem(ie),.true.,vtens(:,:,:,k,ie), &
         118  8742848400 :            var_coef=var_coef1,nu_ratio=nu_ratio1)
         119             :     enddo
         120             :
         121    46753200 :     kptr = kbeg - 1
         122    46753200 :     call edgeVpack(edge3,ttens(:,:,kbeg:kend,ie),kblk,kptr,ie)
         123             :
         124    46753200 :     kptr = kbeg - 1 + nlev
         125 91355752800 :     call edgeVpack(edge3,vtens(:,:,1,kbeg:kend,ie),kblk,kptr,ie)
         126             :
         127    46753200 :     kptr = kbeg - 1 + 2*nlev
         128 91355752800 :     call edgeVpack(edge3,vtens(:,:,2,kbeg:kend,ie),kblk,kptr,ie)
         129             :
         130    46753200 :     kptr = kbeg - 1 + 3*nlev
         131    53402544 :     call edgeVpack(edge3,dptens(:,:,kbeg:kend,ie),kblk,kptr,ie)
         132             :   enddo
         133             :
         134     6649344 :   call bndry_exchange(hybrid,edge3,location='biharmonic_wk_dp3d')
         135             :
         136    53402544 :   do ie=nets,nete
         137             : !CLEAN    rspheremv     => elem(ie)%rspheremp(:,:)
         138             :
         139    46753200 :     kptr = kbeg - 1
         140    46753200 :     call edgeVunpack(edge3,ttens(:,:,kbeg:kend,ie),kblk,kptr,ie)
         141             :
         142    46753200 :     kptr = kbeg - 1 + nlev
         143 >18266*10^7 :     call edgeVunpack(edge3,vtens(:,:,1,kbeg:kend,ie),kblk,kptr,ie)
         144             :
         145    46753200 :     kptr = kbeg - 1 + 2*nlev
         146 >18266*10^7 :     call edgeVunpack(edge3,vtens(:,:,2,kbeg:kend,ie),kblk,kptr,ie)
         147             :
         148    46753200 :     kptr = kbeg - 1 + 3*nlev
         149    46753200 :     call edgeVunpack(edge3,dptens(:,:,kbeg:kend,ie),kblk,kptr,ie)
         150             :
         151    46753200 :     if (use_cslam) then
         152  4394800800 :       do k=1,nlev
         153             : !CLEAN        tmp(:,:)= rspheremv(:,:)*dptens(:,:,k,ie)
         154 91308999600 :         tmp(:,:)= elem(ie)%rspheremp(:,:)*dptens(:,:,k,ie)
         155  4394800800 :         call subcell_Laplace_fluxes(tmp, deriv, elem(ie), np, nc,dpflux(:,:,:,k,ie))
         156             :       enddo
         157             :     endif
         158             :
         159             :     ! apply inverse mass matrix, then apply laplace again
         160             :     !$omp parallel do num_threads(vert_num_threads) private(k,v,tmp,tmp2)
         161  4401450144 :     do k=kbeg,kend
         162             : !CLEAN      tmp(:,:)=rspheremv(:,:)*ttens(:,:,k,ie)
         163 91308999600 :       tmp(:,:)=elem(ie)%rspheremp(:,:)*ttens(:,:,k,ie)
         164  4348047600 :       call laplace_sphere_wk(tmp,deriv,elem(ie),ttens(:,:,k,ie),var_coef=.true.)
         165             : !CLEAN      tmp2(:,:)=rspheremv(:,:)*dptens(:,:,k,ie)
         166 91308999600 :       tmp2(:,:)=elem(ie)%rspheremp(:,:)*dptens(:,:,k,ie)
         167  4348047600 :       call laplace_sphere_wk(tmp2,deriv,elem(ie),dptens(:,:,k,ie),var_coef=.true.)
         168             : !CLEAN      v(:,:,1)=rspheremv(:,:)*vtens(:,:,1,k,ie)
         169             : !CLEAN      v(:,:,2)=rspheremv(:,:)*vtens(:,:,2,k,ie)
         170             :
         171 91308999600 :       v(:,:,1)=elem(ie)%rspheremp(:,:)*vtens(:,:,1,k,ie)
         172 91308999600 :       v(:,:,2)=elem(ie)%rspheremp(:,:)*vtens(:,:,2,k,ie)
         173             :       call vlaplace_sphere_wk(v(:,:,:),deriv,elem(ie),.true.,vtens(:,:,:,k,ie), &
         174  4394800800 :            var_coef=.true.,nu_ratio=nu_ratio2)
         175             :
         176             :     enddo
         177             :   enddo
         178             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         179     6649344 : end subroutine biharmonic_wk_dp3d
         180             :
         181             :
         182     1112832 : subroutine biharmonic_wk_omega(elem,ptens,deriv,edge3,hybrid,nets,nete,kbeg,kend)
         183             :   type (hybrid_t)      , intent(in) :: hybrid
         184             :   type (element_t)     , intent(inout), target :: elem(:)
         185             :   integer              , intent(in)  :: nets,nete
         186             :   integer              , intent(in)  :: kbeg, kend
         187             :   real (kind=r8), dimension(np,np,nlev,nets:nete) :: ptens
         188             :   type (EdgeBuffer_t)  , intent(inout) :: edge3
         189             :   type (derivative_t)  , intent(in) :: deriv
         190             :
         191             :   ! local
         192             :   integer :: i,j,k,kptr,ie,kblk
         193     1112832 :   real (kind=r8), dimension(:,:), pointer :: rspheremv
         194             :   real (kind=r8), dimension(np,np) :: tmp
         195             :   real (kind=r8), dimension(np,np) :: tmp2
         196             :   real (kind=r8), dimension(np,np,2) :: v
         197             :   real (kind=r8) :: nu_ratio1, nu_ratio2
         198             :   logical var_coef1
         199             :
         200     1112832 :   kblk = kend - kbeg + 1
         201             :
         202             :   !if tensor hyperviscosity with tensor V is used, then biharmonic operator is (\grad\cdot V\grad) (\grad \cdot \grad)
         203             :   !so tensor is only used on second call to laplace_sphere_wk
         204     1112832 :   var_coef1 = .true.
         205           0 :   if(hypervis_scaling > 0)    var_coef1 = .false.
         206             :
         207     1112832 :   nu_ratio1=1
         208     1112832 :   nu_ratio2=1
         209             :
         210     8937432 :   do ie=nets,nete
         211             :
         212             :     !$omp parallel do num_threads(vert_num_threads) private(k,tmp)
         213   735512400 :     do k=kbeg,kend
         214 15281443800 :       tmp=elem(ie)%derived%omega(:,:,k)
         215   735512400 :       call laplace_sphere_wk(tmp,deriv,elem(ie),ptens(:,:,k,ie),var_coef=var_coef1)
         216             :     enddo
         217             :
         218     7824600 :     kptr = kbeg - 1
         219     8937432 :     call edgeVpack(edge3,ptens(:,:,kbeg:kend,ie),kblk,kptr,ie)
         220             :   enddo
         221             :
         222     1112832 :   call bndry_exchange(hybrid,edge3,location='biharmonic_wk_omega')
         223             :
         224     8937432 :   do ie=nets,nete
         225     7824600 :     rspheremv     => elem(ie)%rspheremp(:,:)
         226             :
         227     7824600 :     kptr = kbeg - 1
         228     7824600 :     call edgeVunpack(edge3,ptens(:,:,kbeg:kend,ie),kblk,kptr,ie)
         229             :
         230             :     ! apply inverse mass matrix, then apply laplace again
         231             :     !$omp parallel do num_threads(vert_num_threads) private(k,tmp)
         232   736625232 :     do k=kbeg,kend
         233 15281443800 :       tmp(:,:)=rspheremv(:,:)*ptens(:,:,k,ie)
         234   735512400 :       call laplace_sphere_wk(tmp,deriv,elem(ie),ptens(:,:,k,ie),var_coef=.true.)
         235             :     enddo
         236             :   enddo
         237             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         238     7762176 : end subroutine biharmonic_wk_omega
         239             :
         240             :
         241           0 : subroutine biharmonic_wk_scalar(elem,qtens,deriv,edgeq,hybrid,nets,nete)
         242             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         243             : ! compute weak biharmonic operator
         244             : !    input:  qtens = Q
         245             : !    output: qtens = weak biharmonic of Q
         246             : !
         247             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         248             : type (hybrid_t)      , intent(in) :: hybrid
         249             : type (element_t)     , intent(inout), target :: elem(:)
         250             : integer :: nets,nete
         251             : real (kind=r8), dimension(np,np,nlev,qsize,nets:nete) :: qtens
         252             : type (EdgeBuffer_t)  , intent(inout) :: edgeq
         253             : type (derivative_t)  , intent(in) :: deriv
         254             :
         255             : ! local
         256             : integer :: k,kptr,i,j,ie,ic,q
         257             : integer :: kbeg,kend,qbeg,qend
         258             : real (kind=r8), dimension(np,np) :: lap_p
         259             : logical var_coef1
         260             : integer :: kblk,qblk   ! The per thead size of the vertical and tracers
         261             :
         262           0 :   call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
         263             :
         264             :    !if tensor hyperviscosity with tensor V is used, then biharmonic operator is (\grad\cdot V\grad) (\grad \cdot \grad)
         265             :    !so tensor is only used on second call to laplace_sphere_wk
         266           0 :    var_coef1 = .true.
         267           0 :    if(hypervis_scaling > 0)    var_coef1 = .false.
         268             :
         269             :
         270           0 :    kblk = kend - kbeg + 1   ! calculate size of the block of vertical levels
         271           0 :    qblk = qend - qbeg + 1   ! calculate size of the block of tracers
         272             :
         273           0 :    do ie=nets,nete
         274           0 :       do q=qbeg,qend
         275           0 :          do k=kbeg,kend
         276           0 :            lap_p(:,:)=qtens(:,:,k,q,ie)
         277           0 :            call laplace_sphere_wk(lap_p,deriv,elem(ie),qtens(:,:,k,q,ie),var_coef=var_coef1)
         278             :          enddo
         279           0 :          kptr = nlev*(q-1) + kbeg - 1
         280           0 :          call edgeVpack(edgeq, qtens(:,:,kbeg:kend,q,ie),kblk,kptr,ie)
         281             :       enddo
         282             :    enddo
         283             :
         284             :
         285           0 :    call bndry_exchange(hybrid,edgeq,location='biharmonic_wk_scalar')
         286             :
         287           0 :    do ie=nets,nete
         288             :
         289             :       ! apply inverse mass matrix, then apply laplace again
         290           0 :       do q=qbeg,qend
         291           0 :         kptr = nlev*(q-1) + kbeg - 1
         292           0 :         call edgeVunpack(edgeq, qtens(:,:,kbeg:kend,q,ie),kblk,kptr,ie)
         293           0 :         do k=kbeg,kend
         294           0 :            lap_p(:,:)=elem(ie)%rspheremp(:,:)*qtens(:,:,k,q,ie)
         295           0 :            call laplace_sphere_wk(lap_p,deriv,elem(ie),qtens(:,:,k,q,ie),var_coef=.true.)
         296             :         enddo
         297             :       enddo
         298             :    enddo
         299             :
         300             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         301           0 : end subroutine biharmonic_wk_scalar
         302             :
         303             :
         304           0 : subroutine make_C0(zeta,elem,hybrid,nets,nete)
         305             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         306             : ! apply DSS (aka assembly procedure) to zeta.
         307             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         308             :
         309             : type (hybrid_t)      , intent(in) :: hybrid
         310             : type (element_t)     , intent(in), target :: elem(:)
         311             : integer :: nets,nete
         312             : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
         313             :
         314             : ! local
         315             : integer :: k,i,j,ie,ic,kptr,nthread_save
         316             :
         317             :
         318           0 :     call initEdgeBuffer(hybrid%par,edge1,elem,nlev)
         319             :
         320           0 : do ie=nets,nete
         321             : #if (defined COLUMN_OPENMP)
         322             : !$omp parallel do num_threads(vert_num_threads) private(k)
         323             : #endif
         324           0 :    do k=1,nlev
         325           0 :       zeta(:,:,k,ie)=zeta(:,:,k,ie)*elem(ie)%spheremp(:,:)
         326             :    enddo
         327           0 :    kptr=0
         328           0 :    call edgeVpack(edge1, zeta(1,1,1,ie),nlev,kptr,ie)
         329             : enddo
         330           0 : call bndry_exchange(hybrid,edge1,location='make_C0')
         331           0 : do ie=nets,nete
         332           0 :    kptr=0
         333           0 :    call edgeVunpack(edge1, zeta(1,1,1,ie),nlev,kptr, ie)
         334             : #if (defined COLUMN_OPENMP)
         335             : !$omp parallel do num_threads(vert_num_threads) private(k)
         336             : #endif
         337           0 :    do k=1,nlev
         338           0 :       zeta(:,:,k,ie)=zeta(:,:,k,ie)*elem(ie)%rspheremp(:,:)
         339             :    enddo
         340             : enddo
         341             :
         342           0 : call FreeEdgeBuffer(edge1)
         343             :
         344           0 : end subroutine
         345             :
         346             :
         347           0 : subroutine make_C0_vector(v,elem,hybrid,nets,nete)
         348             : #if 1
         349             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         350             : ! apply DSS to a velocity vector
         351             : ! this is a low-performance routine used for I/O and analysis.
         352             : ! no need to optimize
         353             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         354             : type (hybrid_t)      , intent(in) :: hybrid
         355             : type (element_t)     , intent(in), target :: elem(:)
         356             : integer :: nets,nete
         357             : real (kind=r8), dimension(np,np,2,nlev,nets:nete) :: v
         358             :
         359             : ! local
         360             : integer :: k,i,j,ie,ic,kptr
         361             : type (EdgeBuffer_t)          :: edge2
         362           0 : real (kind=r8), dimension(np,np,nlev,nets:nete) :: v1
         363             :
         364           0 : v1(:,:,:,:) = v(:,:,1,:,:)
         365           0 : call make_C0(v1,elem,hybrid,nets,nete)
         366           0 : v(:,:,1,:,:) = v1(:,:,:,:)
         367             :
         368           0 : v1(:,:,:,:) = v(:,:,2,:,:)
         369           0 : call make_C0(v1,elem,hybrid,nets,nete)
         370           0 : v(:,:,2,:,:) = v1(:,:,:,:)
         371             : #else
         372             : type (hybrid_t)      , intent(in) :: hybrid
         373             : type (element_t)     , intent(in), target :: elem(:)
         374             : integer :: nets,nete
         375             : real (kind=r8), dimension(np,np,2,nlev,nets:nete) :: v
         376             :
         377             : ! local
         378             : integer :: k,i,j,ie,ic,kptr
         379             : type (EdgeBuffer_t)          :: edge2
         380             : real (kind=r8), dimension(np,np,nlev,nets:nete) :: v1
         381             :
         382             :
         383             :
         384             :     call initEdgeBuffer(hybrid%par,edge2,elem,2*nlev)
         385             :
         386             : do ie=nets,nete
         387             : #if (defined COLUMN_OPENMP)
         388             : !$omp parallel do num_threads(vert_num_threads) private(k)
         389             : #endif
         390             :    do k=1,nlev
         391             :       v(:,:,1,k,ie)=v(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
         392             :       v(:,:,2,k,ie)=v(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
         393             :    enddo
         394             :    kptr=0
         395             :    call edgeVpack(edge2, v(1,1,1,1,ie),2*nlev,kptr,ie)
         396             : enddo
         397             : call bndry_exchange(hybrid,edge2,location='make_C0_vector')
         398             : do ie=nets,nete
         399             :    kptr=0
         400             :    call edgeVunpack(edge2, v(1,1,1,1,ie),2*nlev,kptr,ie)
         401             : #if (defined COLUMN_OPENMP)
         402             : !$omp parallel do num_threads(vert_num_threads) private(k)
         403             : #endif
         404             :    do k=1,nlev
         405             :       v(:,:,1,k,ie)=v(:,:,1,k,ie)*elem(ie)%rspheremp(:,:)
         406             :       v(:,:,2,k,ie)=v(:,:,2,k,ie)*elem(ie)%rspheremp(:,:)
         407             :    enddo
         408             : enddo
         409             :
         410             : call FreeEdgeBuffer(edge2)
         411             : #endif
         412           0 : end subroutine
         413             :
         414             :
         415             :
         416             :
         417             :
         418             :
         419           0 : subroutine compute_zeta_C0_contra(zeta,elem,hybrid,nets,nete,nt)
         420             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         421             : ! compute C0 vorticity.  That is, solve:
         422             : !     < PHI, zeta > = <PHI, curl(elem%state%v >
         423             : !
         424             : !    input:  v (stored in elem()%, in contra-variant coordinates)
         425             : !    output: zeta(:,:,:,:)
         426             : !
         427             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         428             :
         429             : type (hybrid_t)      , intent(in) :: hybrid
         430             : type (element_t)     , intent(in), target :: elem(:)
         431             : integer :: nt,nets,nete
         432             : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
         433             : real (kind=r8), dimension(np,np,2) :: ulatlon
         434             : real (kind=r8), dimension(np,np) :: v1,v2
         435             :
         436             : ! local
         437             : integer :: k,ie
         438             : type (derivative_t)          :: deriv
         439             :
         440           0 : call derivinit(deriv)
         441             :
         442           0 : do k=1,nlev
         443           0 : do ie=nets,nete
         444           0 :     v1 = elem(ie)%state%v(:,:,1,k,nt)
         445           0 :     v2 = elem(ie)%state%v(:,:,2,k,nt)
         446           0 :     ulatlon(:,:,1) = elem(ie)%D(:,:,1,1)*v1 + elem(ie)%D(:,:,1,2)*v2
         447           0 :     ulatlon(:,:,2) = elem(ie)%D(:,:,2,1)*v1 + elem(ie)%D(:,:,2,2)*v2
         448           0 :    call vorticity_sphere(ulatlon,deriv,elem(ie),zeta(:,:,k,ie))
         449             : enddo
         450             : enddo
         451             :
         452           0 : call make_C0(zeta,elem,hybrid,nets,nete)
         453             :
         454           0 : end subroutine
         455             :
         456             :
         457             :
         458           0 : subroutine compute_div_C0_contra(zeta,elem,hybrid,nets,nete,nt)
         459             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         460             : ! compute C0 divergence. That is, solve:
         461             : !     < PHI, zeta > = <PHI, div(elem%state%v >
         462             : !
         463             : !    input:  v (stored in elem()%, in contra-variant coordinates)
         464             : !    output: zeta(:,:,:,:)
         465             : !
         466             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         467             :
         468             : type (hybrid_t)      , intent(in) :: hybrid
         469             : type (element_t)     , intent(in), target :: elem(:)
         470             : integer :: nt,nets,nete
         471             : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
         472             : real (kind=r8), dimension(np,np,2) :: ulatlon
         473             : real (kind=r8), dimension(np,np) :: v1,v2
         474             :
         475             : ! local
         476             : integer :: k,ie
         477             : type (derivative_t)          :: deriv
         478             :
         479           0 : call derivinit(deriv)
         480             :
         481           0 : do k=1,nlev
         482           0 : do ie=nets,nete
         483           0 :     v1 = elem(ie)%state%v(:,:,1,k,nt)
         484           0 :     v2 = elem(ie)%state%v(:,:,2,k,nt)
         485           0 :     ulatlon(:,:,1) = elem(ie)%D(:,:,1,1)*v1 + elem(ie)%D(:,:,1,2)*v2
         486           0 :     ulatlon(:,:,2) = elem(ie)%D(:,:,2,1)*v1 + elem(ie)%D(:,:,2,2)*v2
         487           0 :    call divergence_sphere(ulatlon,deriv,elem(ie),zeta(:,:,k,ie))
         488             : enddo
         489             : enddo
         490             :
         491           0 : call make_C0(zeta,elem,hybrid,nets,nete)
         492             :
         493           0 : end subroutine
         494             :
         495           0 : subroutine compute_zeta_C0_par(zeta,elem,par,nt)
         496             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         497             : ! compute C0 vorticity.  That is, solve:
         498             : !     < PHI, zeta > = <PHI, curl(elem%state%v >
         499             : !
         500             : !    input:  v (stored in elem()%, in lat-lon coordinates)
         501             : !    output: zeta(:,:,:,:)
         502             : !
         503             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         504             : type (parallel_t) :: par
         505             : type (element_t)     , intent(in), target :: elem(:)
         506             : real (kind=r8), dimension(np,np,nlev,nelemd) :: zeta
         507             : integer :: nt
         508             :
         509             : ! local
         510             : type (hybrid_t)              :: hybrid
         511             : integer :: k,i,j,ie,ic
         512             : type (derivative_t)          :: deriv
         513             :
         514             : ! single thread
         515           0 : hybrid = config_thread_region(par,'serial')
         516             :
         517           0 : call compute_zeta_C0_hybrid(zeta,elem,hybrid,1,nelemd,nt)
         518             :
         519           0 : end subroutine
         520             :
         521             :
         522           0 : subroutine compute_div_C0_par(zeta,elem,par,nt)
         523             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         524             : ! compute C0 divergence. That is, solve:
         525             : !     < PHI, zeta > = <PHI, div(elem%state%v >
         526             : !
         527             : !    input:  v (stored in elem()%, in lat-lon coordinates)
         528             : !    output: zeta(:,:,:,:)
         529             : !
         530             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         531             :
         532             : type (parallel_t) :: par
         533             : type (element_t)     , intent(in), target :: elem(:)
         534             : real (kind=r8), dimension(np,np,nlev,nelemd) :: zeta
         535             : integer :: nt
         536             :
         537             : ! local
         538             : type (hybrid_t)              :: hybrid
         539             : integer :: k,i,j,ie,ic
         540             : type (derivative_t)          :: deriv
         541             :
         542             : ! single thread
         543           0 : hybrid = config_thread_region(par,'serial')
         544             :
         545           0 : call compute_div_C0_hybrid(zeta,elem,hybrid,1,nelemd,nt)
         546             :
         547           0 : end subroutine
         548             :
         549             :
         550             :
         551           0 : subroutine compute_zeta_C0_hybrid(zeta,elem,hybrid,nets,nete,nt)
         552             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         553             : ! compute C0 vorticity.  That is, solve:
         554             : !     < PHI, zeta > = <PHI, curl(elem%state%v >
         555             : !
         556             : !    input:  v (stored in elem()%, in lat-lon coordinates)
         557             : !    output: zeta(:,:,:,:)
         558             : !
         559             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         560             :
         561             : type (hybrid_t)      , intent(in) :: hybrid
         562             : type (element_t)     , intent(in), target :: elem(:)
         563             : integer :: nt,nets,nete
         564             : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
         565             :
         566             : ! local
         567             : integer :: k,i,j,ie,ic
         568             : type (derivative_t)          :: deriv
         569             :
         570           0 : call derivinit(deriv)
         571             :
         572           0 : do ie=nets,nete
         573             : #if (defined COLUMN_OPENMP)
         574             : !$omp parallel do num_threads(vert_num_threads) private(k)
         575             : #endif
         576           0 : do k=1,nlev
         577           0 :    call vorticity_sphere(elem(ie)%state%v(:,:,:,k,nt),deriv,elem(ie),zeta(:,:,k,ie))
         578             : enddo
         579             : enddo
         580             :
         581           0 : call make_C0(zeta,elem,hybrid,nets,nete)
         582             :
         583           0 : end subroutine
         584             :
         585             :
         586           0 : subroutine compute_div_C0_hybrid(zeta,elem,hybrid,nets,nete,nt)
         587             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         588             : ! compute C0 divergence. That is, solve:
         589             : !     < PHI, zeta > = <PHI, div(elem%state%v >
         590             : !
         591             : !    input:  v (stored in elem()%, in lat-lon coordinates)
         592             : !    output: zeta(:,:,:,:)
         593             : !
         594             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         595             :
         596             : type (hybrid_t)      , intent(in) :: hybrid
         597             : type (element_t)     , intent(in), target :: elem(:)
         598             : integer :: nt,nets,nete
         599             : real (kind=r8), dimension(np,np,nlev,nets:nete) :: zeta
         600             :
         601             : ! local
         602             : integer :: k,i,j,ie,ic
         603             : type (derivative_t)          :: deriv
         604             :
         605           0 : call derivinit(deriv)
         606             :
         607           0 : do ie=nets,nete
         608             : #if (defined COLUMN_OPENMP)
         609             : !$omp parallel do num_threads(vert_num_threads) private(k)
         610             : #endif
         611           0 : do k=1,nlev
         612           0 :    call divergence_sphere(elem(ie)%state%v(:,:,:,k,nt),deriv,elem(ie),zeta(:,:,k,ie))
         613             : enddo
         614             : enddo
         615             :
         616           0 : call make_C0(zeta,elem,hybrid,nets,nete)
         617             :
         618           0 : end subroutine
         619             :
         620             :
         621             :
         622             :
         623             :
         624             :
         625             :
         626             :
         627           0 : subroutine neighbor_minmax(hybrid,edgeMinMax,nets,nete,min_neigh,max_neigh)
         628             :
         629             :    type (hybrid_t)      , intent(in) :: hybrid
         630             :    type (EdgeBuffer_t)  , intent(inout) :: edgeMinMax
         631             :    integer :: nets,nete
         632             :    real (kind=r8) :: min_neigh(nlev,qsize,nets:nete)
         633             :    real (kind=r8) :: max_neigh(nlev,qsize,nets:nete)
         634             :    integer :: kblk, qblk
         635             :    ! local
         636             :    integer:: ie, q, k, kptr
         637             :    integer:: kbeg, kend, qbeg, qend
         638             :
         639           0 :    call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
         640             :
         641           0 :    kblk = kend - kbeg + 1   ! calculate size of the block of vertical levels
         642           0 :    qblk = qend - qbeg + 1   ! calculate size of the block of tracers
         643             :
         644           0 :    do ie=nets,nete
         645           0 :       do q = qbeg, qend
         646           0 :          kptr = nlev*(q - 1) + kbeg - 1
         647           0 :          call  edgeSpack(edgeMinMax,min_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
         648           0 :          kptr = qsize*nlev + nlev*(q - 1) + kbeg - 1
         649           0 :          call  edgeSpack(edgeMinMax,max_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
         650             :       enddo
         651             :    enddo
         652             :
         653           0 :    call bndry_exchange(hybrid,edgeMinMax,location='neighbor_minmax')
         654             :
         655           0 :    do ie=nets,nete
         656           0 :       do q=qbeg,qend
         657           0 :          kptr = nlev*(q - 1) + kbeg - 1
         658           0 :          call  edgeSunpackMIN(edgeMinMax,min_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
         659           0 :          kptr = qsize*nlev + nlev*(q - 1) + kbeg - 1
         660           0 :          call  edgeSunpackMAX(edgeMinMax,max_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
         661           0 :          do k=kbeg,kend
         662           0 :             min_neigh(k,q,ie) = max(min_neigh(k,q,ie),0.0_r8)
         663             :          enddo
         664             :       enddo
         665             :    enddo
         666             :
         667           0 : end subroutine neighbor_minmax
         668             :
         669             :
         670           0 : subroutine neighbor_minmax_start(hybrid,edgeMinMax,nets,nete,min_neigh,max_neigh)
         671             :
         672             :    type (hybrid_t)      , intent(in) :: hybrid
         673             :    type (EdgeBuffer_t)  , intent(inout) :: edgeMinMax
         674             :    integer :: nets,nete
         675             :    real (kind=r8) :: min_neigh(nlev,qsize,nets:nete)
         676             :    real (kind=r8) :: max_neigh(nlev,qsize,nets:nete)
         677             :    integer :: kblk, qblk
         678             :    integer :: kbeg, kend, qbeg, qend
         679             :
         680             :    ! local
         681             :    integer :: ie,q, k,kptr
         682             :
         683           0 :    call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
         684             :
         685           0 :    kblk = kend - kbeg + 1   ! calculate size of the block of vertical levels
         686           0 :    qblk = qend - qbeg + 1   ! calculate size of the block of tracers
         687             :
         688           0 :    do ie=nets,nete
         689           0 :       do q=qbeg, qend
         690           0 :          kptr = nlev*(q - 1) + kbeg - 1
         691           0 :          call  edgeSpack(edgeMinMax,min_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
         692           0 :          kptr = qsize*nlev + nlev*(q - 1) + kbeg - 1
         693           0 :          call  edgeSpack(edgeMinMax,max_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
         694             :       enddo
         695             :    enddo
         696             :
         697           0 :    call bndry_exchange_start(hybrid,edgeMinMax,location='neighbor_minmax_start')
         698             :
         699           0 : end subroutine neighbor_minmax_start
         700             :
         701           0 : subroutine neighbor_minmax_finish(hybrid,edgeMinMax,nets,nete,min_neigh,max_neigh)
         702             :
         703             :    type (hybrid_t)      , intent(in) :: hybrid
         704             :    type (EdgeBuffer_t)  , intent(inout) :: edgeMinMax
         705             :    integer :: nets,nete
         706             :    real (kind=r8) :: min_neigh(nlev,qsize,nets:nete)
         707             :    real (kind=r8) :: max_neigh(nlev,qsize,nets:nete)
         708             :    integer :: kblk, qblk
         709             :    integer :: ie,q, k,kptr
         710             :    integer :: kbeg, kend, qbeg, qend
         711             :
         712           0 :    call get_loop_ranges(hybrid,kbeg=kbeg,kend=kend,qbeg=qbeg,qend=qend)
         713             :
         714           0 :    kblk = kend - kbeg + 1   ! calculate size of the block of vertical levels
         715           0 :    qblk = qend - qbeg + 1   ! calculate size of the block of tracers
         716             :
         717           0 :    call bndry_exchange_finish(hybrid,edgeMinMax,location='neighbor_minmax_finish')
         718             :
         719           0 :    do ie=nets,nete
         720           0 :       do q=qbeg, qend
         721           0 :          kptr = nlev*(q - 1) + kbeg - 1
         722           0 :          call  edgeSunpackMIN(edgeMinMax,min_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
         723           0 :          kptr = qsize*nlev + nlev*(q - 1) + kbeg - 1
         724           0 :          call  edgeSunpackMAX(edgeMinMax,max_neigh(kbeg:kend,q,ie),kblk,kptr,ie)
         725           0 :          do k=kbeg,kend
         726           0 :             min_neigh(k,q,ie) = max(min_neigh(k,q,ie),0.0_r8)
         727             :          enddo
         728             :       enddo
         729             :    enddo
         730             :
         731           0 : end subroutine neighbor_minmax_finish
         732             :
         733             : end module viscosity_mod

Generated by: LCOV version 1.14