LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - derivative_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 366 | 758 | 48.3 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 19 | 41 | 46.3 %

              Line data    Source code
           1             : module derivative_mod
           2             :   use shr_kind_mod,       only: r8=>shr_kind_r8
           3             :   use cam_abortutils,     only: endrun
           4             :   use dimensions_mod, only : np, nc, npdg, nelemd, nlev
           5             :   use quadrature_mod, only : quadrature_t, gauss, gausslobatto,legendre, jacobi
           6             :   ! needed for spherical differential operators:
           7             :   use physconst, only: ra
           8             :   use element_mod, only : element_t
           9             :   use control_mod, only : hypervis_scaling, hypervis_power
          10             :   use perf_mod, only : t_startf, t_stopf
          11             :
          12             : implicit none
          13             : private
          14             :
          15             :   type, public :: derivative_t
          16             :      real (kind=r8) :: Dvv(np,np)
          17             :      real (kind=r8) :: Dvv_diag(np,np)
          18             :      real (kind=r8) :: Dvv_twt(np,np)
          19             :      real (kind=r8) :: Mvv_twt(np,np)  ! diagonal matrix of GLL weights
          20             :      real (kind=r8) :: Mfvm(np,nc+1)
          21             :      real (kind=r8) :: Cfvm(np,nc)
          22             :      real (kind=r8) :: legdg(np,np)
          23             :   end type derivative_t
          24             :
          25             :   type, public :: derivative_stag_t
          26             :      real (kind=r8) :: D(np,np)
          27             :      real (kind=r8) :: M(np,np)
          28             :      real (kind=r8) :: Dpv(np,np)
          29             :      real (kind=r8) :: D_twt(np,np)
          30             :      real (kind=r8) :: M_twt(np,np)
          31             :      real (kind=r8) :: M_t(np,np)
          32             :   end type derivative_stag_t
          33             :
          34             :   real (kind=r8), allocatable :: integration_matrix(:,:)
          35             :   real (kind=r8), allocatable :: integration_matrix_physgrid(:,:)
          36             :   real (kind=r8), allocatable :: boundary_interp_matrix(:,:,:)
          37             :
          38             : ! ======================================
          39             : ! Public Interfaces
          40             : ! ======================================
          41             :
          42             :   public :: subcell_integration
          43             :   public :: subcell_dss_fluxes
          44             :   public :: subcell_div_fluxes
          45             :   public :: subcell_Laplace_fluxes
          46             :   public :: allocate_subcell_integration_matrix_cslam    !for consistent se-cslam algorithm
          47             :   public :: allocate_subcell_integration_matrix_physgrid !for integration se basis functions over physgrid control volumes
          48             :
          49             :   public :: derivinit
          50             :
          51             :   public :: gradient
          52             :   public :: gradient_wk
          53             :   public :: vorticity
          54             :   public :: divergence
          55             :
          56             :   public :: interpolate_gll2fvm_corners
          57             :   public :: interpolate_gll2fvm_points
          58             :   public :: remap_phys2gll
          59             :
          60             :
          61             :   interface divergence
          62             :       module procedure divergence_nonstag
          63             :       module procedure divergence_stag
          64             :   end interface
          65             :
          66             :   interface gradient
          67             :       module procedure gradient_str_nonstag
          68             :       module procedure gradient_str_stag
          69             :   end interface
          70             :
          71             :   interface gradient_wk
          72             :       module procedure gradient_wk_nonstag
          73             :       module procedure gradient_wk_stag
          74             :   end interface
          75             :
          76             :   public :: v2pinit
          77             :
          78             :   private :: dmatinit
          79             :   private :: dvvinit
          80             :   private :: dpvinit
          81             :
          82             : ! these routines compute spherical differential operators as opposed to
          83             : ! the gnomonic coordinate operators above.  Vectors (input or output)
          84             : ! are always expressed in lat-lon coordinates
          85             : !
          86             : ! note that weak derivatives (integrated by parts form) can be defined using
          87             : ! contra or co-variant test functions, so
          88             : !
          89             :   public  :: gradient_sphere
          90             :   public  :: gradient_sphere_wk_testcov
          91             :   public  :: gradient_sphere_wk_testcontra   ! only used for debugging
          92             :   public  :: ugradv_sphere
          93             :   public  :: vorticity_sphere
          94             :   public  :: vorticity_sphere_diag
          95             :   public  :: divergence_sphere
          96             :   public  :: curl_sphere
          97             :   public  :: curl_sphere_wk_testcov
          98             :   public  :: divergence_sphere_wk
          99             :   public  :: laplace_sphere_wk
         100             :   public  :: vlaplace_sphere_wk
         101             :   public  :: vlaplace_sphere_wk_mol
         102             :   public  :: element_boundary_integral
         103             :   public  :: edge_flux_u_cg
         104             :   public  :: gll_to_dgmodal
         105             :   public  :: dgmodal_to_gll
         106             :
         107             :   public  :: limiter_optim_iter_full
         108             :
         109             : contains
         110             :
         111             : ! ==========================================
         112             : ! derivinit:
         113             : !
         114             : ! Initialize the matrices for taking
         115             : ! derivatives and interpolating
         116             : ! ==========================================
         117             :
         118      372480 :   subroutine derivinit(deriv,fvm_corners, fvm_points)
         119             :     type (derivative_t)      :: deriv
         120             :     real (kind=r8),optional :: fvm_corners(nc+1)
         121             :     real (kind=r8),optional :: fvm_points(nc)
         122             :
         123             :     ! Local variables
         124             :     type (quadrature_t) :: gp   ! Quadrature points and weights on pressure grid
         125             :
         126             :     real (kind=r8) :: dmat(np,np)
         127             :     real (kind=r8) :: dpv(np,np)
         128             :     real (kind=r8) :: v2p(np,np)
         129             :     real (kind=r8) :: p2v(np,np)
         130             :     real (kind=r8) :: dvv(np,np)
         131             :     real (kind=r8) :: dvv_diag(np,np)
         132             :     real (kind=r8) :: v2v(np,np)
         133             :     real (kind=r8) :: xnorm
         134             :     integer i,j
         135             :
         136             :     ! ============================================
         137             :     ! initialize matrices in r8 precision
         138             :     ! and transfer results into r8
         139             :     ! floating point precision
         140             :     ! ============================================
         141             :
         142      372480 :     gp=gausslobatto(np)
         143             :
         144             :     ! Legendre polynomials of degree npdg-1, on the np GLL grid:
         145      372480 :     if (npdg>np) call endrun( 'FATAL ERROR: npdg>np')
         146      372480 :     if (npdg>0 .and. npdg<np) then
         147             :        ! in this case, we use a DG basis of Legendre polynomials
         148             :        ! stored at the GLL points.
         149           0 :        do i=1,np
         150           0 :           deriv%legdg(1:npdg,i) = legendre(gp%points(i),npdg-1)
         151             :        end do
         152             :        ! normalize
         153           0 :        do j=1,npdg
         154           0 :           xnorm=sqrt(sum(deriv%legdg(j,:)*deriv%legdg(j,:)*gp%weights(:)))
         155           0 :           deriv%legdg(j,:)=deriv%legdg(j,:)/xnorm
         156             :        enddo
         157             :     endif
         158      372480 :     call dvvinit(dvv,gp)
         159     7822080 :     deriv%Dvv(:,:)   = dvv(:,:)
         160             :
         161     1862400 :     do i=1,np
         162     7822080 :        do j=1,np
         163     7449600 :           if (i.eq.j) then
         164     1489920 :              deriv%dvv_diag(i,j)   = dvv(i,j)
         165             :           else
         166     4469760 :              deriv%dvv_diag(i,j) = 0.0_r8
         167             :           endif
         168             :         end do
         169             :      end do
         170             :
         171      372480 :     v2v = 0.0_r8
         172     1862400 :     do i=1,np
         173     1862400 :        v2v(i,i) = gp%weights(i)
         174             :     end do
         175             :
         176     1862400 :     do i=1,np
         177     7822080 :        do j=1,np
         178     7449600 :           dvv(j,i) = dvv(j,i)*gp%weights(i)
         179             :        end do
         180             :     end do
         181             :
         182     7822080 :     deriv%Dvv_twt = TRANSPOSE(dvv)
         183     7822080 :     deriv%Mvv_twt = v2v
         184      372480 :     if (present(fvm_corners)) &
         185        1536 :          call v2pinit(deriv%Mfvm,gp%points,fvm_corners,np,nc+1)
         186             :
         187      372480 :     if (present(fvm_points)) &
         188        1536 :          call v2pinit(deriv%Cfvm,gp%points,fvm_points,np,nc)
         189             :
         190             :     ! notice we deallocate this memory here even though it was allocated
         191             :     ! by the call to gausslobatto.
         192      372480 :     deallocate(gp%points)
         193      372480 :     deallocate(gp%weights)
         194             :
         195      372480 :   end subroutine derivinit
         196             :
         197             : ! =======================================
         198             : ! dmatinit:
         199             : !
         200             : ! Compute rectangular v->p
         201             : ! derivative matrix (dmat)
         202             : ! =======================================
         203             :
         204             :   subroutine dmatinit(dmat)
         205             :
         206             :     real (kind=r8) :: dmat(np,np)
         207             :
         208             :     ! Local variables
         209             :
         210             :     type (quadrature_t) :: gll
         211             :     type (quadrature_t) :: gs
         212             :
         213             :     integer i,j
         214             :     real(kind=r8)  fact,f1,f2
         215             :     real(kind=r8)  func0,func1
         216             :     real(kind=r8)  dis,c0,c1
         217             :
         218             :     real(kind=r8)  :: leg(np,np)
         219             :     real(kind=r8)  ::  jac(0:np-1)
         220             :     real(kind=r8)  :: djac(0:np-1)
         221             :
         222             :     c0 = 0.0_r8
         223             :     c1 = 1.0_r8
         224             :
         225             :     gll= gausslobatto(np)
         226             :     gs = gauss(np)
         227             :
         228             :     ! =============================================================
         229             :     ! Compute Legendre polynomials on Gauss-Lobatto grid (velocity)
         230             :     ! =============================================================
         231             :
         232             :     do i=1,np
         233             :        leg(:,i) = legendre(gll%points(i),np-1)
         234             :     end do
         235             :
         236             :     ! ================================================================
         237             :     !  Derivatives of velocity cardinal functions on pressure grid
         238             :     !  d(i,j) = D(j,i) = D' (D-transpose) since D(i,j) = dh_j(x_i)/dx
         239             :     ! ================================================================
         240             :
         241             :     fact = np*(np-1)
         242             :
         243             :     do j=1,np
         244             :        call jacobi(np-1,gs%points(j),c0,c0,jac(0:np-1),djac(0:np-1))
         245             :        func0 =  jac(np-1)
         246             :        func1 = djac(np-1)
         247             :        f1 = fact*func0
         248             :        f2 = (c1 - gs%points(j))*(c1 + gs%points(j)) * func1
         249             :        do i = 1, np
         250             :           if ( gs%points(j) /= gll%points(i) ) then
         251             :              dis = gs%points(j) - gll%points(i)
         252             :              dmat(i,j) = func0 / ( leg(np,i)*dis ) + f2 / (fact*leg(np,i)*dis*dis)
         253             :           else
         254             :              dmat(i,j) = c0
         255             :           endif
         256             :        end do
         257             :     end do
         258             :
         259             :     deallocate(gll%points)
         260             :     deallocate(gll%weights)
         261             :
         262             :     deallocate(gs%points)
         263             :     deallocate(gs%weights)
         264             :
         265             : end subroutine dmatinit
         266             :
         267             : ! =======================================
         268             : ! dpvinit:
         269             : !
         270             : ! Compute rectangular p->v
         271             : ! derivative matrix (dmat)
         272             : ! for strong gradients
         273             : ! =======================================
         274             :
         275             : subroutine dpvinit(dmat)
         276             :
         277             : real (kind=r8) :: dmat(np,np)
         278             :
         279             : ! Local variables
         280             :
         281             : type (quadrature_t) :: gll
         282             : type (quadrature_t) :: gs
         283             :
         284             : integer i,j
         285             : real(kind=r8)  dis,c0,c1
         286             :
         287             : real(kind=r8)  :: legv(0:np,np)
         288             : real(kind=r8)  :: dlegv(0:np,np)
         289             :
         290             : real(kind=r8)  :: leg(0:np)
         291             : real(kind=r8)  :: dleg(0:np)
         292             :
         293             : c0 = 0.0_r8
         294             : c1 = 1.0_r8
         295             :
         296             : gll= gausslobatto(np)
         297             : gs = gauss(np)
         298             :
         299             : ! =============================================================
         300             : ! Compute Legendre polynomials on Gauss-Lobatto grid (velocity)
         301             : ! =============================================================
         302             :
         303             : do i=1,np
         304             : call jacobi(np,gll%points(i),c0,c0,legv(0:np,i),dlegv(0:np,i))
         305             : end do
         306             :
         307             : ! ================================================================
         308             : !  Derivatives of velocity cardinal functions on pressure grid
         309             : !  d(i,j) = D(j,i) = D' (D-transpose) since D(i,j) = dh_j(x_i)/dx
         310             : ! ================================================================
         311             :
         312             :     do j=1,np
         313             :        call jacobi(np,gs%points(j),c0,c0,leg(0:np),dleg(0:np))
         314             :        do i = 1, np
         315             :           if ( gs%points(j) /= gll%points(i) ) then
         316             :              dis = gll%points(i) - gs%points(j)
         317             :              dmat(j,i) = dlegv(np,i)/( dleg(np)*dis ) -  legv(np,i)/ (dleg(np)*dis*dis)
         318             :           else
         319             :              dmat(j,i) = c0
         320             :           endif
         321             :        end do
         322             :     end do
         323             :
         324             :     deallocate(gll%points)
         325             :     deallocate(gll%weights)
         326             :
         327             :     deallocate(gs%points)
         328             :     deallocate(gs%weights)
         329             :
         330             :   end subroutine dpvinit
         331             :
         332             : ! =======================================
         333             : ! v2pinit:
         334             : ! Compute interpolation matrix from gll(1:n1) -> gs(1:n2)
         335             : ! =======================================
         336        3072 :   subroutine v2pinit(v2p,gll,gs,n1,n2)
         337             :     integer :: n1,n2
         338             :     real(kind=r8)  ::  v2p(n1,n2)
         339        6144 :     real(kind=r8)  ::  v2p_new(n1,n2)
         340             :     real(kind=r8)  ::  gll(n1),gs(n2)
         341             :     ! Local variables
         342             :
         343             :     integer i,j,k,m,l
         344             :     real(kind=r8)  fact,f1, sum
         345             :     real(kind=r8)  func0,func1
         346             :
         347        6144 :     real(kind=r8)  :: leg(n1,n1)
         348             :     real(kind=r8)  ::  jac(0:n1-1)
         349             :     real(kind=r8)  :: djac(0:n1-1)
         350             :     real(kind=r8)  :: c0,c1
         351             :
         352             :     type(quadrature_t) :: gll_pts
         353        6144 :     real(kind=r8)  :: leg_out(n1,n2)
         354        6144 :     real(kind=r8)  :: gamma(n1)
         355             :
         356        3072 :     c0 = 0.0_r8
         357        3072 :     c1 = 1.0_r8
         358             :
         359             :     ! ==============================================================
         360             :     ! Compute Legendre polynomials on Gauss-Lobatto grid (velocity)
         361             :     ! ==============================================================
         362             :
         363        3072 :     fact = -n1*(n1-1)
         364       15360 :     do i=1,n1
         365       12288 :        leg(:,i) = legendre(gll(i),n1-1)
         366       15360 :        leg(n1,i) = fact * leg(n1,i)
         367             :     end do
         368             :
         369             :     ! ===================================================
         370             :     !  Velocity cardinal functions on pressure grid
         371             :     ! ===================================================
         372             :     ! NEW VERSION, with no division by (gs(j)-gll(i)):
         373             :
         374             :     ! compute legendre polynomials at output points:
         375        3072 :     gll_pts = gausslobatto(n1)
         376             :
         377        3072 :     fact = -n1*(n1-1)
         378       13824 :     do i=1,n2
         379       10752 :        leg_out(:,i) = legendre(gs(i),n1-1)
         380       13824 :        leg_out(n1,i) = fact * leg_out(n1,i)
         381             :     end do
         382             :
         383             :
         384             :     ! compute gamma: (normalization factor for inv(leg)
         385       15360 :     do m=1,n1
         386       12288 :        gamma(m)=0
         387       61440 :        do i=1,n1
         388       61440 :           gamma(m)=gamma(m)+leg(m,i)*leg(m,i)*gll_pts%weights(i)
         389             :        enddo
         390       15360 :        gamma(m)=1/gamma(m)
         391             :     enddo
         392             :
         393             :     ! compute product of leg_out * inv(leg):
         394       13824 :     do j=1,n2   ! this should be fvm points
         395       56832 :        do l=1,n1   ! this should be GLL points
         396       43008 :           sum=0
         397      215040 :           do k=1,n1  ! number of polynomials = number of GLL points
         398      215040 :              sum=sum + leg_out(k,j)*gamma(k)*leg(k,l)
         399             :           enddo
         400       53760 :           v2p_new(l,j) = gll_pts%weights(l)*sum
         401             :        enddo
         402             :     enddo
         403        3072 :     deallocate(gll_pts%points)
         404        3072 :     deallocate(gll_pts%weights)
         405             :
         406       56832 :     v2p=v2p_new
         407        3072 :   end subroutine v2pinit
         408             :
         409             :
         410             :
         411             : ! =======================================
         412             : ! dvvinit:
         413             : !
         414             : ! Compute rectangular v->v
         415             : ! derivative matrix (dvv)
         416             : ! =======================================
         417             :
         418      372480 :   subroutine dvvinit(dvv,gll)
         419             :
         420             :     real(kind=r8)  ::  dvv(np,np)
         421             :     type (quadrature_t)   :: gll
         422             :
         423             :     ! Local variables
         424             :
         425             :     real(kind=r8)  :: leg(np,np)
         426             :     real(kind=r8)  :: c0,c1,c4
         427             :
         428             :     integer i,j
         429             :
         430      372480 :     c0 = 0.0_r8
         431      372480 :     c1 = 1.0_r8
         432      372480 :     c4 = 4.0_r8
         433             :
         434     1862400 :     do i=1,np
         435     1862400 :        leg(:,i) = legendre(gll%points(i),np-1)
         436             :     end do
         437             :
         438     7822080 :     dvv(:,:) = c0
         439     1862400 :     do j=1,np
         440     3724800 :        do i=1,j-1
         441     3724800 :           dvv(j,i) = (c1/(gll%points(i)-gll%points(j)))*leg(np,i)/leg(np,j)
         442             :        end do
         443     1489920 :        dvv(j,j) = c0
         444     4097280 :        do i=j+1,np
         445     3724800 :           dvv(j,i) = (c1/(gll%points(i)-gll%points(j)))*leg(np,i)/leg(np,j)
         446             :        end do
         447             :     end do
         448             :
         449             :
         450      372480 :     dvv(np,np) = + np*(np-1)/c4
         451      372480 :     dvv(1,1)   = - np*(np-1)/c4
         452             :
         453      372480 :   end subroutine dvvinit
         454             :
         455             : !  ================================================
         456             : !  divergence_stag:
         457             : !
         458             : !  Compute divergence (maps v grid -> p grid)
         459             : !  ================================================
         460             :
         461           0 :   subroutine divergence_stag(v,deriv,div)
         462             :
         463             :     real(kind=r8), intent(in) :: v(np,np,2)
         464             :     type (derivative_stag_t), intent(in) :: deriv
         465             :     real(kind=r8), intent(out) :: div(np,np)
         466             :
         467             :     ! Local
         468             :
         469             :     integer i
         470             :     integer j
         471             :     integer l
         472             :
         473             :     real(kind=r8)  sumx00
         474             :     real(kind=r8)  sumy00
         475             :
         476             :     real(kind=r8) :: vtemp(np,np,2)
         477             :
         478             :
         479           0 :      do j=1,np
         480           0 :         do l=1,np
         481             :
         482             :            sumx00=0.0d0
         483             :            sumy00=0.0d0
         484             : !DIR$ UNROLL(NP)
         485           0 :            do i=1,np
         486           0 :               sumx00 = sumx00 + deriv%D(i,l  )*v(i,j  ,1)
         487           0 :               sumy00 = sumy00 + deriv%M(i,l  )*v(i,j  ,2)
         488             :            enddo
         489           0 :           vtemp(j  ,l  ,1) = sumx00
         490           0 :           vtemp(j  ,l  ,2) = sumy00
         491             :        enddo
         492             :     enddo
         493           0 :     do j=1,np
         494           0 :        do i=1,np
         495             :           sumx00=0.0d0
         496             :           sumy00=0.0d0
         497             : !DIR$ UNROLL(NP)
         498           0 :           do l=1,np
         499           0 :              sumx00 = sumx00 +  deriv%M(l,j  )*vtemp(l,i  ,1)
         500           0 :              sumy00 = sumy00 +  deriv%D(l,j  )*vtemp(l,i  ,2)
         501             :           enddo
         502           0 :           div(i  ,j  ) = sumx00 + sumy00
         503             :
         504             :        enddo
         505             :     enddo
         506             :
         507           0 :   end subroutine divergence_stag
         508             :
         509             : !  ================================================
         510             : !  divergence_nonstag:
         511             : !
         512             : !  Compute divergence (maps v->v)
         513             : !  ================================================
         514             :
         515           0 :   subroutine divergence_nonstag(v,deriv,div)
         516             :
         517             :     real(kind=r8), intent(in) :: v(np,np,2)
         518             :     type (derivative_t), intent(in) :: deriv
         519             :
         520             :     real(kind=r8), intent(out) :: div(np,np)
         521             :
         522             :     ! Local
         523             :
         524             :     integer i
         525             :     integer j
         526             :     integer l
         527             :
         528             :     real(kind=r8) ::  dudx00
         529             :     real(kind=r8) ::  dvdy00
         530             :
         531             :     real(kind=r8) ::  vvtemp(np,np)
         532             :
         533           0 :      do j=1,np
         534           0 :         do l=1,np
         535             :            dudx00=0.0d0
         536             :            dvdy00=0.0d0
         537             : !DIR$ UNROLL(NP)
         538           0 :            do i=1,np
         539           0 :               dudx00 = dudx00 + deriv%Dvv(i,l  )*v(i,j  ,1)
         540           0 :               dvdy00 = dvdy00 + deriv%Dvv(i,l  )*v(j  ,i,2)
         541             :            end do
         542             :
         543           0 :            div(l  ,j  ) = dudx00
         544           0 :            vvtemp(j  ,l  ) = dvdy00
         545             :         end do
         546             :     end do
         547           0 :     do j=1,np
         548           0 :        do i=1,np
         549           0 :           div(i,j)=div(i,j)+vvtemp(i,j)
         550             :        end do
         551             :     end do
         552             :
         553           0 :   end subroutine divergence_nonstag
         554             :
         555             : !  ================================================
         556             : !  gradient_wk_stag:
         557             : !
         558             : !  Compute the weak form gradient:
         559             : !  maps scalar field on the pressure grid to the
         560             : !  velocity grid
         561             : !  ================================================
         562             :
         563           0 :   function gradient_wk_stag(p,deriv) result(dp)
         564             :
         565             :     type (derivative_stag_t), intent(in) :: deriv
         566             :     real(kind=r8), intent(in) :: p(np,np)
         567             :
         568             :     real(kind=r8)             :: dp(np,np,2)
         569             :
         570             :     ! Local
         571             :
         572             :     integer i
         573             :     integer j
         574             :     integer l
         575             :
         576             :     real(kind=r8)  sumx00,sumx01
         577             :     real(kind=r8)  sumy00,sumy01
         578             :
         579             :     real(kind=r8)  :: vtempt(np,np,2)
         580             :
         581           0 :     do j=1,np
         582           0 :        do l=1,np
         583             :           sumx00=0.0d0
         584             :           sumy00=0.0d0
         585             : !DIR$ UNROLL(NP)
         586           0 :            do i=1,np
         587           0 :               sumx00 = sumx00 + deriv%D_twt(i,l  )*p(i,j  )
         588           0 :               sumy00 = sumy00 + deriv%M_twt(i,l  )*p(i,j  )
         589             :            enddo
         590           0 :            vtempt(j  ,l  ,1) = sumx00
         591           0 :            vtempt(j  ,l  ,2) = sumy00
         592             :         enddo
         593             :     enddo
         594           0 :     do j=1,np
         595           0 :        do i=1,np
         596             :           sumx00=0.0d0
         597             :           sumy00=0.0d0
         598             : !DIR$ UNROLL(NP)
         599           0 :           do l=1,np
         600           0 :              sumx00 = sumx00 +  deriv%M_twt(l,j  )*vtempt(l,i  ,1)
         601           0 :              sumy00 = sumy00 +  deriv%D_twt(l,j  )*vtempt(l,i  ,2)
         602             :           enddo
         603           0 :           dp(i  ,j  ,1) = sumx00
         604           0 :           dp(i  ,j  ,2) = sumy00
         605             :       enddo
         606             :     enddo
         607             :
         608             :
         609           0 :   end function gradient_wk_stag
         610             :
         611             : !  ================================================
         612             : !  gradient_wk_nonstag:
         613             : !
         614             : !  Compute the weak form gradient:
         615             : !  maps scalar field on the Gauss-Lobatto grid to the
         616             : !  weak gradient on the Gauss-Lobbatto grid
         617             : !  ================================================
         618             :
         619           0 :   function gradient_wk_nonstag(p,deriv) result(dp)
         620             :
         621             :     type (derivative_t), intent(in) :: deriv
         622             :     real(kind=r8), intent(in) :: p(np,np)
         623             :
         624             :     real(kind=r8)             :: dp(np,np,2)
         625             :
         626             :     ! Local
         627             :
         628             :     integer i
         629             :     integer j
         630             :     integer l
         631             :
         632             :     real(kind=r8)  sumx00
         633             :     real(kind=r8)  sumy00
         634             :
         635             :     real(kind=r8)  :: vvtempt(np,np,2)
         636             :
         637           0 :        do j=1,np
         638           0 :           do l=1,np
         639             :              sumx00=0.0d0
         640             :              sumy00=0.0d0
         641             : !DIR$ UNROLL(NP)
         642           0 :              do i=1,np
         643           0 :                 sumx00 = sumx00 + deriv%Dvv_twt(i,l  )*p(i,j  )
         644           0 :                 sumy00 = sumy00 + deriv%Mvv_twt(i,l  )*p(i,j  )
         645             :              end do
         646           0 :              vvtempt(j  ,l  ,1) = sumx00
         647           0 :              vvtempt(j  ,l  ,2) = sumy00
         648             :           end do
         649             :        end do
         650             :
         651           0 :        do j=1,np
         652           0 :           do i=1,np
         653             :              sumx00=0.0d0
         654             :              sumy00=0.0d0
         655             : !DIR$ UNROLL(NP)
         656           0 :              do l=1,np
         657           0 :                 sumx00 = sumx00 +  deriv%Mvv_twt(l,j  )*vvtempt(l,i  ,1)
         658           0 :                 sumy00 = sumy00 +  deriv%Dvv_twt(l,j  )*vvtempt(l,i  ,2)
         659             :              end do
         660           0 :              dp(i  ,j  ,1) = sumx00
         661           0 :              dp(i  ,j  ,2) = sumy00
         662             :           end do
         663             :        end do
         664           0 :   end function gradient_wk_nonstag
         665             :
         666             : !  ================================================
         667             : !  gradient_str_stag:
         668             : !
         669             : !  Compute the *strong* form gradient:
         670             : !  maps scalar field on the pressure grid to the
         671             : !  velocity grid
         672             : !  ================================================
         673             :
         674           0 :   subroutine  gradient_str_stag(p,deriv,dp)
         675             :
         676             :     type (derivative_stag_t), intent(in) :: deriv
         677             :     real(kind=r8), intent(in) :: p(np,np)
         678             :
         679             :     real(kind=r8), intent(out) :: dp(np,np,2)
         680             :
         681             :     ! Local
         682             :
         683             :     integer i
         684             :     integer j
         685             :     integer l
         686             :
         687             :     real(kind=r8)  sumx00
         688             :     real(kind=r8)  sumy00
         689             :
         690             :     real(kind=r8)  :: vtempt(np,np,2)
         691           0 :     do j=1,np
         692           0 :        do l=1,np
         693             :           sumx00=0.0d0
         694             :           sumy00=0.0d0
         695             : !DIR$ UNROLL(NP)
         696           0 :           do i=1,np
         697           0 :              sumx00 = sumx00 + deriv%Dpv(i,l  )*p(i,j  )
         698           0 :              sumy00 = sumy00 + deriv%M_t(i,l  )*p(i,j  )
         699             :           enddo
         700           0 :           vtempt(j  ,l  ,1) = sumx00
         701           0 :           vtempt(j  ,l  ,2) = sumy00
         702             :        enddo
         703             :     enddo
         704           0 :     do j=1,np
         705           0 :        do i=1,np
         706             :           sumx00=0.0d0
         707             :           sumy00=0.0d0
         708             : !DIR$ UNROLL(NP)
         709           0 :           do l=1,np
         710           0 :              sumx00 = sumx00 +  deriv%M_t(l,j  )*vtempt(l,i  ,1)
         711           0 :              sumy00 = sumy00 +  deriv%Dpv(l,j  )*vtempt(l,i  ,2)
         712             :           enddo
         713           0 :           dp(i  ,j  ,1) = sumx00
         714           0 :           dp(i  ,j  ,2) = sumy00
         715             :        enddo
         716             :     enddo
         717             :
         718           0 :   end subroutine gradient_str_stag
         719             :
         720             : !  ================================================
         721             : !  gradient_str_nonstag:
         722             : !
         723             : !  Compute the *strong* gradient on the velocity grid
         724             : !  of a scalar field on the velocity grid
         725             : !  ================================================
         726             :
         727           0 :   subroutine gradient_str_nonstag(s,deriv,ds)
         728             :
         729             :     type (derivative_t), intent(in) :: deriv
         730             :     real(kind=r8), intent(in) :: s(np,np)
         731             :     real(kind=r8), intent(out) :: ds(np,np,2)
         732             :
         733             :     integer i
         734             :     integer j
         735             :     integer l
         736             :     real(kind=r8) ::  dsdx00,dsdx01
         737             :     real(kind=r8) ::  dsdy00,dsdy01
         738           0 :        do j=1,np
         739           0 :           do l=1,np
         740             :              dsdx00=0.0d0
         741             :              dsdy00=0.0d0
         742             : !DIR$ UNROLL(NP)
         743           0 :              do i=1,np
         744           0 :                 dsdx00 = dsdx00 + deriv%Dvv(i,l  )*s(i,j  )
         745           0 :                 dsdy00 = dsdy00 + deriv%Dvv(i,l  )*s(j  ,i)
         746             :              end do
         747           0 :              ds(l  ,j  ,1) = dsdx00
         748           0 :              ds(j  ,l  ,2) = dsdy00
         749             :           end do
         750             :        end do
         751           0 :   end subroutine  gradient_str_nonstag
         752             :
         753             : !  ================================================
         754             : !  vorticity:
         755             : !
         756             : !  Compute the vorticity of the velocity field on the
         757             : !  velocity grid
         758             : !  ================================================
         759             :
         760           0 :   subroutine vorticity(v,deriv,vort)
         761             :
         762             :     type (derivative_t), intent(in) :: deriv
         763             :     real(kind=r8), intent(in) :: v(np,np,2)
         764             :
         765             :     real(kind=r8), intent(out) :: vort(np,np)
         766             :
         767             :     integer i
         768             :     integer j
         769             :     integer l
         770             :
         771             :     real(kind=r8) ::  dvdx00,dvdx01
         772             :     real(kind=r8) ::  dudy00,dudy01
         773             :
         774             :     real(kind=r8)  :: vvtemp(np,np)
         775           0 :     do j=1,np
         776           0 :        do l=1,np
         777             :           dudy00=0.0d0
         778             :           dvdx00=0.0d0
         779             : !DIR$ UNROLL(NP)
         780           0 :           do i=1,np
         781           0 :              dvdx00 = dvdx00 + deriv%Dvv(i,l  )*v(i,j  ,2)
         782           0 :              dudy00 = dudy00 + deriv%Dvv(i,l  )*v(j  ,i,1)
         783             :           enddo
         784           0 :           vort(l  ,j  ) = dvdx00
         785           0 :           vvtemp(j  ,l  ) = dudy00
         786             :        enddo
         787             :     enddo
         788           0 :     do j=1,np
         789           0 :        do i=1,np
         790           0 :           vort(i,j)=vort(i,j)-vvtemp(i,j)
         791             :        end do
         792             :     end do
         793             :
         794           0 :   end subroutine vorticity
         795             :
         796             : !  ================================================
         797             : !  interpolate_gll2fvm_points:
         798             : !
         799             : !  shape funtion interpolation from data on GLL grid to cellcenters on physics grid
         800             : !  Author: Christoph Erath
         801             : !  ================================================
         802           0 :   function interpolate_gll2fvm_points(v,deriv) result(p)
         803             :
         804             :     real(kind=r8), intent(in) :: v(np,np)
         805             :     type (derivative_t)         :: deriv
         806             :     real(kind=r8) :: p(nc,nc)
         807             :
         808             :     ! Local
         809             :     integer i
         810             :     integer j
         811             :     integer l
         812             :
         813             :     real(kind=r8)  sumx00,sumx01
         814             :     real(kind=r8)  sumx10,sumx11
         815             :     real(kind=r8)  vtemp(np,nc)
         816             :
         817           0 :     do j=1,np
         818           0 :        do l=1,nc
         819             :           sumx00=0.0d0
         820             : !DIR$ UNROLL(NP)
         821           0 :           do i=1,np
         822           0 :              sumx00 = sumx00 + deriv%Cfvm(i,l  )*v(i,j  )
         823             :           enddo
         824           0 :           vtemp(j  ,l) = sumx00
         825             :         enddo
         826             :     enddo
         827           0 :     do j=1,nc
         828           0 :        do i=1,nc
         829             :           sumx00=0.0d0
         830             : !DIR$ UNROLL(NP)
         831           0 :           do l=1,np
         832           0 :              sumx00 = sumx00 + deriv%Cfvm(l,j  )*vtemp(l,i)
         833             :           enddo
         834           0 :           p(i  ,j  ) = sumx00
         835             :        enddo
         836             :     enddo
         837           0 :   end function interpolate_gll2fvm_points
         838             : !  ================================================
         839             : !  interpolate_gll2fvm_corners:
         840             : !
         841             : !  shape funtion interpolation from data on GLL grid to physics grid
         842             : !
         843             : !  ================================================
         844           0 :   function interpolate_gll2fvm_corners(v,deriv) result(p)
         845             :
         846             :     real(kind=r8), intent(in) :: v(np,np)
         847             :     type (derivative_t), intent(in) :: deriv
         848             :     real(kind=r8) :: p(nc+1,nc+1)
         849             :
         850             :     ! Local
         851             :     integer i
         852             :     integer j
         853             :     integer l
         854             :
         855             :     real(kind=r8)  sumx00,sumx01
         856             :     real(kind=r8)  sumx10,sumx11
         857             :     real(kind=r8)  vtemp(np,nc+1)
         858             :
         859           0 :     do j=1,np
         860           0 :        do l=1,nc+1
         861             :           sumx00=0.0d0
         862             : !DIR$ UNROLL(NP)
         863           0 :           do i=1,np
         864           0 :              sumx00 = sumx00 + deriv%Mfvm(i,l  )*v(i,j  )
         865             :           enddo
         866           0 :           vtemp(j  ,l) = sumx00
         867             :         enddo
         868             :     enddo
         869           0 :     do j=1,nc+1
         870           0 :        do i=1,nc+1
         871             :           sumx00=0.0d0
         872             : !DIR$ UNROLL(NP)
         873           0 :           do l=1,np
         874           0 :              sumx00 = sumx00 + deriv%Mfvm(l,j  )*vtemp(l,i)
         875             :           enddo
         876           0 :           p(i  ,j  ) = sumx00
         877             :        enddo
         878             :     enddo
         879           0 :   end function interpolate_gll2fvm_corners
         880             :
         881             :
         882             : !  ================================================
         883             : !  remap_phys2gll:
         884             : !
         885             : !  interpolate to an equally spaced (in reference element coordinate system)
         886             : !  "physics" grid to the GLL grid
         887             : !
         888             : !  1st order, monotone, conservative
         889             : !  MT initial version 2013
         890             : !  ================================================
         891           0 :   function remap_phys2gll(pin,nphys) result(pout)
         892             :     integer :: nphys
         893             :     real(kind=r8), intent(in) :: pin(nphys*nphys)
         894             :     real(kind=r8) :: pout(np,np)
         895             :
         896             :     ! Local
         897             :     integer, save  :: nphys_init=0
         898             :     integer, save  :: nintersect
         899             :     real(kind=r8),save,pointer :: acell(:)  ! arrivial cell index of i'th intersection
         900             :     real(kind=r8),save,pointer :: dcell(:)  ! departure cell index of i'th intersection
         901             :     real(kind=r8),save,pointer :: delta(:)  ! length of i'th intersection
         902             :     real(kind=r8),save,pointer :: delta_a(:)  ! length of arrival cells
         903             :     integer in_i,in_j,ia,ja,id,jd,count,i,j
         904             :     logical :: found
         905             :
         906             :     real(kind=r8) :: tol = 1.0e-13_r8
         907             :     real(kind=r8) :: weight,x1,x2,dx
         908           0 :     real(kind=r8) :: gll_edges(np+1),phys_edges(nphys+1)
         909             :     type(quadrature_t) :: gll_pts
         910           0 :     if (nphys_init/=nphys) then
         911             :        ! setup (most be done on masterthread only) since all data is static
         912             :        ! MT: move barrier inside if loop - we dont want a barrier every regular call
         913             : !OMP BARRIER
         914             : !OMP MASTER
         915           0 :        nphys_init=nphys
         916             :        ! find number of intersections
         917           0 :        nintersect = np+nphys-1  ! max number of possible intersections
         918           0 :        allocate(acell(nintersect))
         919           0 :        allocate(dcell(nintersect))
         920           0 :        allocate(delta(nintersect))
         921           0 :        allocate(delta_a(np))
         922             :
         923             :        ! compute phys grid cell edges on [-1,1]
         924           0 :        do i=1,nphys+1
         925           0 :           dx = 2d0/nphys
         926           0 :           phys_edges(i)=-1 + (i-1)*dx
         927             :        enddo
         928             :
         929             :        ! compute GLL cell edges on [-1,1]
         930           0 :        gll_pts = gausslobatto(np)
         931           0 :        gll_edges(1)=-1
         932           0 :        do i=2,np
         933           0 :           gll_edges(i) = gll_edges(i-1) + gll_pts%weights(i-1)
         934             :        enddo
         935           0 :        gll_edges(np+1)=1
         936           0 :        delta_a=gll_pts%weights
         937           0 :        deallocate(gll_pts%points)
         938           0 :        deallocate(gll_pts%weights)
         939             :
         940           0 :        count=0
         941           0 :        x1=-1
         942           0 :        do while ( abs(x1-1) > tol )
         943             :           ! find point x2 closet to x1 and x2>x1:
         944           0 :           x2 = 1.1_r8
         945           0 :           do ia=2,np+1
         946           0 :              if (gll_edges(ia)>x1) then
         947           0 :                 if ( ( gll_edges(ia)-x1) < (x2-x1) ) then
         948           0 :                    x2=gll_edges(ia)
         949             :                 endif
         950             :              endif
         951             :           enddo
         952           0 :           do id=2,nphys+1
         953           0 :              if (phys_edges(id)>x1) then
         954           0 :                 if ( ( phys_edges(id)-x1) < (x2-x1) ) then
         955           0 :                    x2=phys_edges(id)
         956             :                 endif
         957             :              endif
         958             :           enddo
         959           0 :           print *,'x2=',x2
         960           0 :           if (x2>1+tol) call endrun('ERROR: did not find next intersection point')
         961           0 :           if (x2<=x1) call endrun('ERROR: next intersection point did not advance')
         962           0 :           count=count+1
         963           0 :           if (count>nintersect) call endrun('ERROR: search failuer: nintersect was too small')
         964           0 :           delta(count)=x2-x1
         965             :
         966           0 :           found=.false.
         967           0 :           do ia=1,np
         968           0 :              if (gll_edges(ia) <= x1+tol  .and.  x2-tol <= gll_edges(ia+1)) then
         969           0 :                 found=.true.
         970           0 :                 acell(count)=ia
         971             :              endif
         972             :           enddo
         973           0 :           if (.not. found) call endrun('ERROR: interval search problem')
         974             :
         975           0 :           found=.false.
         976           0 :           do id=1,nphys
         977           0 :              if (phys_edges(id) <= x1+tol .and.  x2-tol <= phys_edges(id+1)) then
         978           0 :                 found=.true.
         979           0 :                 dcell(count)=id
         980             :              endif
         981             :           enddo
         982           0 :           if (.not. found) call endrun('ERROR: interval search problem')
         983           0 :           x1=x2
         984             :        enddo
         985             :        ! reset to actual number of intersections
         986           0 :        nintersect=count
         987             : !OMP END MASTER
         988             : !OMP BARRIER
         989             :     endif
         990             :
         991           0 :     pout=0
         992           0 :     do in_i = 1,nintersect
         993           0 :        do in_j = 1,nintersect
         994           0 :           ia = acell(in_i)
         995           0 :           ja = acell(in_j)
         996           0 :           id = dcell(in_i)
         997           0 :           jd = dcell(in_j)
         998             :           ! mass in intersection region:  value*area_intersect
         999             :           ! value_arrival = value*area_intersect/area_arrival
        1000           0 :           weight = (  delta(in_i)*delta(in_j) ) / ( delta_a(ia)*delta_a(ja))
        1001             :           ! accumulate contribution from each intersection region:
        1002           0 :           pout(ia,ja) = pout(ia,ja) + weight*pin(id+(jd-1)*nphys)
        1003             :        enddo
        1004             :     enddo
        1005             :
        1006           0 :     end function remap_phys2gll
        1007             :
        1008             : !----------------------------------------------------------------
        1009             :
        1010             : !DIR$ ATTRIBUTES FORCEINLINE :: gradient_sphere
        1011 57852549000 :   subroutine gradient_sphere(s,deriv,Dinv,ds)
        1012             : !
        1013             : !   input s:  scalar
        1014             : !   output  ds: spherical gradient of s, lat-lon coordinates
        1015             : !
        1016             :
        1017             :     type (derivative_t), intent(in) :: deriv
        1018             :     real(kind=r8), intent(in), dimension(np,np,2,2) :: Dinv
        1019             :     real(kind=r8), intent(in) :: s(np,np)
        1020             :     real(kind=r8), intent(out) :: ds(np,np,2)
        1021             :
        1022             :     integer i
        1023             :     integer j
        1024             :     integer l
        1025             :
        1026             :     real(kind=r8) ::  dsdx00, dsdy00
        1027             :     real(kind=r8) ::  v1(np,np),v2(np,np)
        1028             :
        1029 >28926*10^7 :     do j=1,np
        1030 >12149*10^8 :        do l=1,np
        1031             :           dsdx00=0.0d0
        1032             :           dsdy00=0.0d0
        1033             : !DIR$ UNROLL(NP)
        1034 >46282*10^8 :           do i=1,np
        1035 >37025*10^8 :              dsdx00 = dsdx00 + deriv%Dvv(i,l  )*s(i,j  )
        1036 >46282*10^8 :              dsdy00 = dsdy00 + deriv%Dvv(i,l  )*s(j  ,i)
        1037             :           end do
        1038 >92564*10^7 :           v1(l  ,j  ) = dsdx00*ra
        1039 >11570*10^8 :           v2(j  ,l  ) = dsdy00*ra
        1040             :        end do
        1041             :     end do
        1042             :     ! convert covarient to latlon
        1043 >28926*10^7 :     do j=1,np
        1044 >12149*10^8 :        do i=1,np
        1045 >92564*10^7 :           ds(i,j,1)=Dinv(i,j,1,1)*v1(i,j) + Dinv(i,j,2,1)*v2(i,j)
        1046 >11570*10^8 :           ds(i,j,2)=Dinv(i,j,1,2)*v1(i,j) + Dinv(i,j,2,2)*v2(i,j)
        1047             :        enddo
        1048             :     enddo
        1049             :
        1050 57852549000 :     end subroutine gradient_sphere
        1051             :
        1052             :
        1053  8883108000 :   function curl_sphere_wk_testcov(s,deriv,elem) result(ds)
        1054             : !
        1055             : !   integrated-by-parts gradient, w.r.t. COVARIANT test functions
        1056             : !   input s:  scalar  (assumed to be s*khat)
        1057             : !   output  ds: weak curl, lat/lon coordinates
        1058             : !
        1059             : ! starting with:
        1060             : !   PHIcov1 = (PHI,0)  covariant vector
        1061             : !   PHIcov2 = (0,PHI)  covariant vector
        1062             : !
        1063             : !   ds1 = integral[ PHIcov1 dot curl(s*khat) ]
        1064             : !   ds2 = integral[ PHIcov2 dot curl(s*khat) ]
        1065             : ! integrate by parts:
        1066             : !   ds1 = integral[ vor(PHIcov1) * s ]
        1067             : !   ds2 = integral[ vor(PHIcov1) * s ]
        1068             : !
        1069             : !     PHIcov1 = (PHI^mn,0)
        1070             : !     PHIcov2 = (0,PHI^mn)
        1071             : !  vorticity() acts on covariant vectors:
        1072             : !   ds1 = sum wij g  s_ij 1/g (  (PHIcov1_2)_x  - (PHIcov1_1)_y )
        1073             : !       = -sum wij s_ij  d/dy (PHI^mn )
        1074             : ! for d/dy component, only sum over i=m
        1075             : !       = -sum  w_mj s_mj   d( PHI^n)(j)
        1076             : !           j
        1077             : !
        1078             : !   ds2 = sum wij g  s_ij 1/g (  (PHIcov2_2)_x  - (PHIcov2_1)_y )
        1079             : !       = +sum wij s_ij  d/dx (PHI^mn )
        1080             : ! for d/dx component, only sum over j=n
        1081             : !       = +sum  w_in s_in  d( PHI^m)(i)
        1082             : !           i
        1083             : !
        1084             :     type (derivative_t), intent(in) :: deriv
        1085             :     type (element_t), intent(in) :: elem
        1086             :     real(kind=r8), intent(in) :: s(np,np)
        1087             :
        1088             :     real(kind=r8) :: ds(np,np,2)
        1089             :
        1090             :     integer i,j,l,m,n
        1091             :     real(kind=r8) ::  dscontra(np,np,2)
        1092             :
        1093  8883108000 :     dscontra=0
        1094 44415540000 :     do n=1,np
        1095 >18654*10^7 :        do m=1,np
        1096             : !DIR$ UNROLL(NP)
        1097 >74618*10^7 :           do j=1,np
        1098             :              ! phi(n)_y  sum over second index, 1st index fixed at m
        1099 >56851*10^7 :              dscontra(m,n,1)=dscontra(m,n,1)-(elem%mp(m,j)*s(m,j)*deriv%Dvv(n,j) )*ra
        1100             :              ! phi(m)_x  sum over first index, second index fixed at n
        1101 >71064*10^7 :              dscontra(m,n,2)=dscontra(m,n,2)+(elem%mp(j,n)*s(j,n)*deriv%Dvv(m,j) )*ra
        1102             :           enddo
        1103             :        enddo
        1104             :     enddo
        1105             :
        1106             :     ! convert contra -> latlon
        1107 44415540000 :     do j=1,np
        1108 >18654*10^7 :        do i=1,np
        1109 >14212*10^7 :           ds(i,j,1)=(elem%D(i,j,1,1)*dscontra(i,j,1) + elem%D(i,j,1,2)*dscontra(i,j,2))
        1110 >17766*10^7 :           ds(i,j,2)=(elem%D(i,j,2,1)*dscontra(i,j,1) + elem%D(i,j,2,2)*dscontra(i,j,2))
        1111             :        enddo
        1112             :     enddo
        1113  8883108000 :     end function curl_sphere_wk_testcov
        1114             :
        1115             :
        1116  8883108000 :   function gradient_sphere_wk_testcov(s,deriv,elem) result(ds)
        1117             : !
        1118             : !   integrated-by-parts gradient, w.r.t. COVARIANT test functions
        1119             : !   input s:  scalar
        1120             : !   output  ds: weak gradient, lat/lon coordinates
        1121             : !   ds = - integral[ div(PHIcov) s ]
        1122             : !
        1123             : !     PHIcov1 = (PHI^mn,0)
        1124             : !     PHIcov2 = (0,PHI^mn)
        1125             : !   div() acts on contra components, so convert test function to contra:
        1126             : !     PHIcontra1 =  metinv PHIcov1  = (a^mn,b^mn)*PHI^mn
        1127             : !                                     a = metinv(1,1)  b=metinv(2,1)
        1128             : !
        1129             : !   ds1 = sum wij g  s_ij 1/g ( g a PHI^mn)_x  + ( g b PHI^mn)_y )
        1130             : !       = sum  wij s_ij  ag(m,n)  d/dx( PHI^mn ) + bg(m,n) d/dy( PHI^mn)
        1131             : !          i,j
        1132             : ! for d/dx component, only sum over j=n
        1133             : !       = sum  w_in s_in  ag(m,n)  d( PHI^m)(i)
        1134             : !          i
        1135             : ! for d/dy component, only sum over i=m
        1136             : !       = sum  w_mj s_mj  bg(m,n)  d( PHI^n)(j)
        1137             : !          j
        1138             : !
        1139             : !
        1140             : ! This formula is identical to gradient_sphere_wk_testcontra, except that
        1141             : !    g(m,n) is replaced by a(m,n)*g(m,n)
        1142             : !  and we have two terms for each componet of ds
        1143             : !
        1144             : !
        1145             :     type (derivative_t), intent(in) :: deriv
        1146             :     type (element_t), intent(in) :: elem
        1147             :     real(kind=r8), intent(in) :: s(np,np)
        1148             :
        1149             :     real(kind=r8) :: ds(np,np,2)
        1150             :
        1151             :     integer i,j,l,m,n
        1152             :     real(kind=r8) ::  dscontra(np,np,2)
        1153             :
        1154             :
        1155  8883108000 :     dscontra=0
        1156 44415540000 :     do n=1,np
        1157 >18654*10^7 :        do m=1,np
        1158             : !DIR$ UNROLL(NP)
        1159 >74618*10^7 :           do j=1,np
        1160 >11370*10^8 :              dscontra(m,n,1)=dscontra(m,n,1)-(&
        1161 >56851*10^7 :                   (elem%mp(j,n)*elem%metinv(m,n,1,1)*elem%metdet(m,n)*s(j,n)*deriv%Dvv(m,j) ) +&
        1162             :                   (elem%mp(m,j)*elem%metinv(m,n,2,1)*elem%metdet(m,n)*s(m,j)*deriv%Dvv(n,j) ) &
        1163 >17055*10^8 :                   ) *ra
        1164             :
        1165             :              dscontra(m,n,2)=dscontra(m,n,2)-(&
        1166             :                   (elem%mp(j,n)*elem%metinv(m,n,1,2)*elem%metdet(m,n)*s(j,n)*deriv%Dvv(m,j) ) +&
        1167             :                   (elem%mp(m,j)*elem%metinv(m,n,2,2)*elem%metdet(m,n)*s(m,j)*deriv%Dvv(n,j) ) &
        1168 >71064*10^7 :                   ) *ra
        1169             :           enddo
        1170             :        enddo
        1171             :     enddo
        1172             :     ! convert contra -> latlon
        1173 44415540000 :     do j=1,np
        1174 >18654*10^7 :        do i=1,np
        1175 >14212*10^7 :           ds(i,j,1)=(elem%D(i,j,1,1)*dscontra(i,j,1) + elem%D(i,j,1,2)*dscontra(i,j,2))
        1176 >17766*10^7 :           ds(i,j,2)=(elem%D(i,j,2,1)*dscontra(i,j,1) + elem%D(i,j,2,2)*dscontra(i,j,2))
        1177             :        enddo
        1178             :     enddo
        1179             :
        1180  8883108000 :     end function gradient_sphere_wk_testcov
        1181             :
        1182             :
        1183           0 :   function gradient_sphere_wk_testcontra(s,deriv,elem) result(ds)
        1184             : !
        1185             : !   integrated-by-parts gradient, w.r.t. CONTRA test functions
        1186             : !   input s:  scalar
        1187             : !   output  ds: weak gradient, lat/lon coordinates
        1188             : !
        1189             : !   integral[ div(phivec) s ] = sum  spheremp()* divergence_sphere(phivec) *s
        1190             : !   ds1 = above formual with phivec=(PHI,0) in CONTRA coordinates
        1191             : !   ds2 = above formual with phivec=(0,PHI) in CONTRA coordinates
        1192             : !
        1193             : ! PHI = (phi,0)
        1194             : !   s1 =  sum w_ij s_ij g_ij 1/g_ij ( g_ij PHI^mn )x
        1195             : !      =  sum w_ij s_ij g_mn dx(PHI^mn)_ij
        1196             : !         ij
        1197             : ! because x derivative is zero for j<>n, only have to sum over j=n
        1198             : !   s1(m,n)  =  sum w_i,n g_mn dx(PHI^m)_i,n s_i,n
        1199             : !                i
        1200             : !
        1201             :     type (derivative_t), intent(in) :: deriv
        1202             :     type (element_t), intent(in) :: elem
        1203             :     real(kind=r8), intent(in) :: s(np,np)
        1204             :
        1205             :     real(kind=r8) :: ds(np,np,2)
        1206             :
        1207             :     integer i,j,l,m,n
        1208             :     real(kind=r8) ::  dscov(np,np,2)
        1209             :
        1210           0 :     dscov=0
        1211           0 :     do n=1,np
        1212           0 :        do m=1,np
        1213             : !DIR$ UNROLL(NP)
        1214           0 :           do j=1,np
        1215             :              ! phi(m)_x  sum over first index, second index fixed at n
        1216           0 :              dscov(m,n,1)=dscov(m,n,1)-(elem%mp(j,n)*elem%metdet(m,n)*s(j,n)*deriv%Dvv(m,j) )*ra
        1217             :              ! phi(n)_y  sum over second index, 1st index fixed at m
        1218           0 :              dscov(m,n,2)=dscov(m,n,2)-(elem%mp(m,j)*elem%metdet(m,n)*s(m,j)*deriv%Dvv(n,j) )*ra
        1219             :           enddo
        1220             :        enddo
        1221             :     enddo
        1222             :
        1223             :     ! convert covariant -> latlon
        1224           0 :     ds(:,:,1)=elem%Dinv(:,:,1,1)*dscov(:,:,1) + elem%Dinv(:,:,2,1)*dscov(:,:,2)
        1225           0 :     ds(:,:,2)=elem%Dinv(:,:,1,2)*dscov(:,:,1) + elem%Dinv(:,:,2,2)*dscov(:,:,2)
        1226             :
        1227           0 :     end function gradient_sphere_wk_testcontra
        1228             :
        1229   242562600 :   function ugradv_sphere(u,v,deriv,elem) result(ugradv)
        1230             : !
        1231             : !   input:  vectors u and v  (latlon coordinates)
        1232             : !   output: vector  [ u dot grad ] v  (latlon coordinates)
        1233             : !
        1234             :     type (derivative_t), intent(in) :: deriv
        1235             :     type (element_t), intent(in) :: elem
        1236             :     real(kind=r8), intent(in) :: u(np,np,2)
        1237             :     real(kind=r8), intent(in) :: v(np,np,2)
        1238             :
        1239             :     real(kind=r8) :: ugradv(np,np,2)
        1240             :     real(kind=r8) :: dum_cart(np,np,3)
        1241             :     real(kind=r8) :: temp(np,np,2)
        1242             :
        1243             :     integer :: component
        1244             :
        1245             :     ! latlon -> cartesian
        1246   970250400 :     do component=1,3
        1247             :        ! Summing along the third dimension is a sum over components for each point.
        1248             :        ! (This is just a faster way of doing a dot product for each grid point,
        1249             :        ! since reindexing the inputs to use the intrinsic effectively would be
        1250             :        ! just asking for trouble.)
        1251 38810016000 :        dum_cart(:,:,component)=sum( elem%vec_sphere2cart(:,:,component,:)*v(:,:,:) ,3)
        1252             :     end do
        1253             :
        1254             :     ! Do ugradv on the cartesian components.
        1255   970250400 :     do component=1,3
        1256             :        ! Dot u with the gradient of each component
        1257   727687800 :        call gradient_sphere(dum_cart(:,:,component),deriv,elem%Dinv,temp)
        1258 38810016000 :        dum_cart(:,:,component) = sum( u(:,:,:) * temp,3)
        1259             :     enddo
        1260             :
        1261             :     ! cartesian -> latlon
        1262   727687800 :     do component=1,2
        1263             :        ! vec_sphere2cart is its own pseudoinverse.
        1264 33716201400 :        ugradv(:,:,component) = sum(dum_cart(:,:,:)*elem%vec_sphere2cart(:,:,:,component), 3)
        1265             :     end do
        1266             :
        1267   242562600 :   end function ugradv_sphere
        1268             :
        1269             :
        1270             :
        1271           0 :   function curl_sphere(s,deriv,elem) result(ds)
        1272             : !
        1273             : !   input s:  scalar  (assumed to be  s khat)
        1274             : !   output  curl(s khat) vector in lat-lon coordinates
        1275             : !
        1276             : !   This subroutine can be used to compute divergence free velocity fields,
        1277             : !   since div(ds)=0
        1278             : !
        1279             : !    first compute:
        1280             : !    curl(s khat) = (1/jacobian) ( ds/dy, -ds/dx ) in contra-variant coordinates
        1281             : !    then map to lat-lon
        1282             : !
        1283             :     type (derivative_t), intent(in) :: deriv
        1284             :     type (element_t), intent(in) :: elem
        1285             :     real(kind=r8), intent(in) :: s(np,np)
        1286             :
        1287             :     real(kind=r8) :: ds(np,np,2)
        1288             :
        1289             :     integer i
        1290             :     integer j
        1291             :     integer l
        1292             :
        1293             :     real(kind=r8) ::  dsdx00
        1294             :     real(kind=r8) ::  dsdy00
        1295             :     real(kind=r8) ::  v1(np,np),v2(np,np)
        1296             :
        1297           0 :     do j=1,np
        1298           0 :        do l=1,np
        1299             :           dsdx00=0.0d0
        1300             :           dsdy00=0.0d0
        1301             : !DIR$ UNROLL(NP)
        1302           0 :           do i=1,np
        1303           0 :              dsdx00 = dsdx00 + deriv%Dvv(i,l  )*s(i,j  )
        1304           0 :              dsdy00 = dsdy00 + deriv%Dvv(i,l  )*s(j  ,i)
        1305             :           end do
        1306           0 :           v2(l  ,j  ) = -dsdx00*ra
        1307           0 :           v1(j  ,l  ) =  dsdy00*ra
        1308             :        end do
        1309             :     end do
        1310             :     ! convert contra -> latlon *and* divide by jacobian
        1311           0 :     do j=1,np
        1312           0 :        do i=1,np
        1313           0 :           ds(i,j,1)=(elem%D(i,j,1,1)*v1(i,j) + elem%D(i,j,1,2)*v2(i,j))/elem%metdet(i,j)
        1314           0 :           ds(i,j,2)= (elem%D(i,j,2,1)*v1(i,j) + elem%D(i,j,2,2)*v2(i,j))/elem%metdet(i,j)
        1315             :        enddo
        1316             :     enddo
        1317             :
        1318           0 :     end function curl_sphere
        1319             :
        1320             :
        1321             : !--------------------------------------------------------------------------
        1322             :
        1323             :
        1324             :
        1325 19221591600 :   subroutine divergence_sphere_wk(v,deriv,elem,div)
        1326             : !
        1327             : !   input:  v = velocity in lat-lon coordinates
        1328             : !   ouput:  div(v)  spherical divergence of v, integrated by parts
        1329             : !
        1330             : !   Computes  -< grad(psi) dot v >
        1331             : !   (the integrated by parts version of < psi div(v) > )
        1332             : !
        1333             : !   note: after DSS, divergence_sphere() and divergence_sphere_wk()
        1334             : !   are identical to roundoff, as theory predicts.
        1335             : !
        1336             :     real(kind=r8), intent(in) :: v(np,np,2)  ! in lat-lon coordinates
        1337             :     type (derivative_t), intent(in) :: deriv
        1338             :     type (element_t), intent(in) :: elem
        1339             :     real(kind=r8),intent(out)  :: div(np,np)
        1340             :
        1341             :     ! Local
        1342             :
        1343             :     integer i,j,m,n
        1344             :
        1345             :     real(kind=r8) ::  vtemp(np,np,2)
        1346             :     real(kind=r8) ::  ggtemp(np,np,2)
        1347             :     real(kind=r8) ::  gtemp(np,np,2)
        1348             :     real(kind=r8) ::  psi(np,np)
        1349             :     real(kind=r8) :: xtmp
        1350             :
        1351             :     ! latlon- > contra
        1352 96107958000 :     do j=1,np
        1353 >40365*10^7 :        do i=1,np
        1354 >30754*10^7 :           vtemp(i,j,1)=(elem%Dinv(i,j,1,1)*v(i,j,1) + elem%Dinv(i,j,1,2)*v(i,j,2))
        1355 >38443*10^7 :           vtemp(i,j,2)=(elem%Dinv(i,j,2,1)*v(i,j,1) + elem%Dinv(i,j,2,2)*v(i,j,2))
        1356             :        enddo
        1357             :     enddo
        1358             :
        1359 96107958000 :     do n=1,np
        1360 >40365*10^7 :        do m=1,np
        1361             :
        1362 >30754*10^7 :           div(m,n)=0
        1363             : !DIR$ UNROLL(NP)
        1364 >16146*10^8 :           do j=1,np
        1365 >12301*10^8 :              div(m,n)=div(m,n)-(elem%spheremp(j,n)*vtemp(j,n,1)*deriv%Dvv(m,j) &
        1366             :                               +elem%spheremp(m,j)*vtemp(m,j,2)*deriv%Dvv(n,j)) &
        1367 >27679*10^8 :                               * ra
        1368             :           enddo
        1369             :
        1370             :        end do
        1371             :     end do
        1372             :
        1373 19221591600 :   end subroutine divergence_sphere_wk
        1374             :
        1375             :
        1376             :
        1377           0 :   function element_boundary_integral(v,deriv,elem) result(result)
        1378             : !
        1379             : !   input:  v = velocity in lat-lon coordinates
        1380             : !   ouput:  result(i,j) = contour integral of PHI_ij * v dot normal
        1381             : !           where PHI_ij = cardinal function at i,j GLL point
        1382             : !
        1383             : !   this routine is used just to check spectral element integration by parts identities
        1384             : !
        1385             :     real(kind=r8), intent(in) :: v(np,np,2)  ! in lat-lon coordinates
        1386             :     type (derivative_t), intent(in) :: deriv
        1387             :     type (element_t), intent(in) :: elem
        1388             :     real(kind=r8) :: result(np,np)
        1389             :
        1390             :     ! Local
        1391             :     real(kind=r8) :: ucontra(np,np,2)  ! in lat-lon coordinates
        1392             :     integer i,j
        1393             :
        1394             :     ! latlon->contra
        1395           0 :     do j=1,np
        1396           0 :        do i=1,np
        1397           0 :           ucontra(i,j,1)=(elem%Dinv(i,j,1,1)*v(i,j,1) + elem%Dinv(i,j,1,2)*v(i,j,2))
        1398           0 :           ucontra(i,j,2)=(elem%Dinv(i,j,2,1)*v(i,j,1) + elem%Dinv(i,j,2,2)*v(i,j,2))
        1399             :        enddo
        1400             :     enddo
        1401             :
        1402             :     ! note: GLL weights  weight(i) = Mvv_twt(i,i)
        1403           0 :     result=0
        1404           0 :     j=1
        1405           0 :     do i=1,np
        1406           0 :        result(i,j)=result(i,j)-deriv%Mvv_twt(i,i)*elem%metdet(i,j)*ucontra(i,j,2)*ra
        1407             :     enddo
        1408             :
        1409           0 :     j=np
        1410           0 :     do i=1,np
        1411           0 :        result(i,j)=result(i,j)+deriv%Mvv_twt(i,i)*elem%metdet(i,j)*ucontra(i,j,2)*ra
        1412             :     enddo
        1413             :
        1414           0 :     i=1
        1415           0 :     do j=1,np
        1416           0 :        result(i,j)=result(i,j)-deriv%Mvv_twt(j,j)*elem%metdet(i,j)*ucontra(i,j,1)*ra
        1417             :     enddo
        1418             :
        1419           0 :     i=np
        1420           0 :     do j=1,np
        1421           0 :        result(i,j)=result(i,j)+deriv%Mvv_twt(j,j)*elem%metdet(i,j)*ucontra(i,j,1)*ra
        1422             :     enddo
        1423           0 :   end function element_boundary_integral
        1424             :
        1425             :
        1426             :
        1427           0 :   function edge_flux_u_cg( v,p,pedges, deriv, elem, u_is_contra) result(result)
        1428             : !
        1429             : !
        1430             : !   input:  v = velocity in contra or lat-lon coordinates (CONTINUIOUS)
        1431             : !           p      = scalar on this element
        1432             : !           pedges = scalar edge data from neighbor elements
        1433             : !
        1434             : !   ouput:  result(i,j) = contour integral of PHI_ij * pstar * v dot normal
        1435             : !           where PHI_ij = cardinal function at i,j GLL point
        1436             : !           pstar = centered or other flux
        1437             : !
        1438             :     real(kind=r8), intent(in) :: v(np,np,2)
        1439             :     real(kind=r8), intent(in) :: p(np,np)
        1440             :     real(kind=r8), intent(in) :: pedges(0:np+1,0:np+1)
        1441             :     type (derivative_t), intent(in) :: deriv
        1442             :     type (element_t), intent(in) :: elem
        1443             :     real(kind=r8) :: result(np,np)
        1444             :     logical :: u_is_contra
        1445             :
        1446             :     ! Local
        1447             :     real(kind=r8) :: ucontra(np,np,2)  ! in lat-lon coordinates
        1448             :     real(kind=r8) :: flux,pstar
        1449             :     integer i,j
        1450             :
        1451             :
        1452           0 :     result=0
        1453             :
        1454             :
        1455           0 :     if (u_is_contra) then
        1456           0 :        ucontra=v
        1457             :     else
        1458             :        ! latlon->contra
        1459           0 :        do j=1,np
        1460           0 :           do i=1,np
        1461           0 :              ucontra(i,j,1)=(elem%Dinv(i,j,1,1)*v(i,j,1) + elem%Dinv(i,j,1,2)*v(i,j,2))
        1462           0 :              ucontra(i,j,2)=(elem%Dinv(i,j,2,1)*v(i,j,1) + elem%Dinv(i,j,2,2)*v(i,j,2))
        1463             :           enddo
        1464             :        enddo
        1465             :     endif
        1466             :     ! upwind
        1467           0 :     do i=1,np
        1468           0 :        j=1
        1469           0 :        pstar=p(i,j)
        1470           0 :        if (ucontra(i,j,2)>0) pstar=pedges(i,0)
        1471           0 :        flux = -pstar*ucontra(i,j,2)*( deriv%Mvv_twt(i,i)*elem%metdet(i,j)*ra)
        1472           0 :        result(i,j)=result(i,j)+flux
        1473             :
        1474           0 :        j=np
        1475           0 :        pstar=p(i,j)
        1476           0 :        if (ucontra(i,j,2)<0) pstar=pedges(i,np+1)
        1477           0 :        flux = pstar*ucontra(i,j,2)* ( deriv%Mvv_twt(i,i)*elem%metdet(i,j)*ra)
        1478           0 :        result(i,j)=result(i,j)+flux
        1479             :     enddo
        1480             :
        1481           0 :     do j=1,np
        1482           0 :        i=1
        1483           0 :        pstar=p(i,j)
        1484           0 :        if (ucontra(i,j,1)>0) pstar=pedges(0,j)
        1485           0 :        flux = -pstar*ucontra(i,j,1)* ( deriv%Mvv_twt(j,j)*elem%metdet(i,j)*ra)
        1486           0 :        result(i,j)=result(i,j)+flux
        1487             :
        1488           0 :        i=np
        1489           0 :        pstar=p(i,j)
        1490           0 :        if (ucontra(i,j,1)<0) pstar=pedges(np+1,j)
        1491           0 :        flux = pstar*ucontra(i,j,1)* ( deriv%Mvv_twt(j,j)*elem%metdet(i,j)*ra)
        1492           0 :        result(i,j)=result(i,j)+flux
        1493             :     end do
        1494             :
        1495           0 :   end function edge_flux_u_cg
        1496             :
        1497             :
        1498             : !DIR$ ATTRIBUTES FORCEINLINE :: vorticity_sphere
        1499 16129854000 :   subroutine vorticity_sphere(v,deriv,elem,vort)
        1500             : !
        1501             : !   input:  v = velocity in lat-lon coordinates
        1502             : !   ouput:  spherical vorticity of v
        1503             : !
        1504             :
        1505             :     type (derivative_t), intent(in) :: deriv
        1506             :     type (element_t), intent(in) :: elem
        1507             :     real(kind=r8), intent(in) :: v(np,np,2)
        1508             :
        1509             :     real(kind=r8), intent(out) :: vort(np,np)
        1510             :
        1511             :     integer i
        1512             :     integer j
        1513             :     integer l
        1514             :
        1515             :     real(kind=r8) ::  dvdx00,dudy00
        1516             :     real(kind=r8) ::  vco(np,np,2)
        1517             :     real(kind=r8) ::  vtemp(np,np)
        1518             :
        1519             :     ! convert to covariant form
        1520 80649270000 :     do j=1,np
        1521 >33872*10^7 :        do i=1,np
        1522 >25807*10^7 :           vco(i,j,1)=(elem%D(i,j,1,1)*v(i,j,1) + elem%D(i,j,2,1)*v(i,j,2))
        1523 >32259*10^7 :           vco(i,j,2)=(elem%D(i,j,1,2)*v(i,j,1) + elem%D(i,j,2,2)*v(i,j,2))
        1524             :        enddo
        1525             :     enddo
        1526             :
        1527 80649270000 :     do j=1,np
        1528 >33872*10^7 :        do l=1,np
        1529             :
        1530             :           dudy00=0.0d0
        1531             :           dvdx00=0.0d0
        1532             :
        1533             : !DIR$ UNROLL(NP)
        1534 >12903*10^8 :           do i=1,np
        1535 >10323*10^8 :              dvdx00 = dvdx00 + deriv%Dvv(i,l  )*vco(i,j  ,2)
        1536 >12903*10^8 :              dudy00 = dudy00 + deriv%Dvv(i,l  )*vco(j  ,i,1)
        1537             :           enddo
        1538             :
        1539 >25807*10^7 :           vort(l  ,j  ) = dvdx00
        1540 >32259*10^7 :           vtemp(j  ,l  ) = dudy00
        1541             :        enddo
        1542             :     enddo
        1543             :
        1544 80649270000 :     do j=1,np
        1545 >33872*10^7 :        do i=1,np
        1546 >32259*10^7 :           vort(i,j)=(vort(i,j)-vtemp(i,j))*(elem%rmetdet(i,j)*ra)
        1547             :        end do
        1548             :     end do
        1549             :
        1550 16129854000 :   end subroutine vorticity_sphere
        1551             :
        1552           0 :   function vorticity_sphere_diag(v,deriv,elem) result(vort)
        1553             :   !
        1554             :   !   input:  v = velocity in lat-lon coordinates
        1555             :   !   ouput:  diagonal component of spherical vorticity of v
        1556             :   !
        1557             :
        1558             :       type (derivative_t), intent(in) :: deriv
        1559             :       type (element_t), intent(in) :: elem
        1560             :       real(kind=r8), intent(in) :: v(np,np,2)
        1561             :
        1562             :       real(kind=r8) :: vort(np,np)
        1563             :
        1564             :       integer i
        1565             :       integer j
        1566             :       integer l
        1567             :
        1568             :       real(kind=r8) :: dvdx00,dudy00
        1569             :       real(kind=r8) :: vco(np,np,2)
        1570             :       real(kind=r8) :: vtemp(np,np)
        1571             :       real(kind=r8) :: rdx
        1572             :       real(kind=r8) :: rdy
        1573             :
        1574             :       ! convert to covariant form
        1575             :
        1576           0 :       do j=1,np
        1577           0 :          do i=1,np
        1578           0 :             vco(i,j,1)=(elem%D(i,j,1,1)*v(i,j,1) + elem%D(i,j,2,1)*v(i,j,2))
        1579           0 :             vco(i,j,2)=(elem%D(i,j,1,2)*v(i,j,1) + elem%D(i,j,2,2)*v(i,j,2))
        1580             :          enddo
        1581             :       enddo
        1582             :
        1583             :
        1584           0 :       do j=1,np
        1585           0 :          do l=1,np
        1586             :             dudy00=0.0d0
        1587             :             dvdx00=0.0d0
        1588             : !DIR$ UNROLL(NP)
        1589           0 :             do i=1,np
        1590           0 :                dvdx00 = dvdx00 + deriv%Dvv_diag(i,l)*vco(i,j ,2)
        1591           0 :                dudy00 = dudy00 + deriv%Dvv_diag(i,l)*vco(j ,i,1)
        1592             :             enddo
        1593           0 :             vort(l ,j) = dvdx00
        1594           0 :             vtemp(j ,l) = dudy00
        1595             :          enddo
        1596             :       enddo
        1597             :
        1598           0 :       do j=1,np
        1599           0 :          do i=1,np
        1600           0 :           vort(i,j)=(vort(i,j)-vtemp(i,j))*(elem%rmetdet(i,j)*ra)
        1601             :          end do
        1602             :       end do
        1603             :
        1604           0 :   end function vorticity_sphere_diag
        1605             :
        1606             : !DIR$ ATTRIBUTES FORCEINLINE :: divergence_sphere
        1607 23619162600 :   subroutine divergence_sphere(v,deriv,elem,div)
        1608             : !
        1609             : !   input:  v = velocity in lat-lon coordinates
        1610             : !   ouput:  div(v)  spherical divergence of v
        1611             : !
        1612             :
        1613             :
        1614             :     real(kind=r8), intent(in) :: v(np,np,2)  ! in lat-lon coordinates
        1615             :     type (derivative_t), intent(in) :: deriv
        1616             :     type (element_t), intent(in) :: elem
        1617             :     real(kind=r8), intent(out) :: div(np,np)
        1618             :
        1619             :     ! Local
        1620             :
        1621             :     integer i
        1622             :     integer j
        1623             :     integer l
        1624             :
        1625             :     real(kind=r8) ::  dudx00
        1626             :     real(kind=r8) ::  dvdy00
        1627             :     real(kind=r8) ::  gv(np,np,2),vvtemp(np,np)
        1628             :
        1629             :     ! convert to contra variant form and multiply by g
        1630 >11809*10^7 :     do j=1,np
        1631 >49600*10^7 :        do i=1,np
        1632 >37790*10^7 :           gv(i,j,1)=elem%metdet(i,j)*(elem%Dinv(i,j,1,1)*v(i,j,1) + elem%Dinv(i,j,1,2)*v(i,j,2))
        1633 >47238*10^7 :           gv(i,j,2)=elem%metdet(i,j)*(elem%Dinv(i,j,2,1)*v(i,j,1) + elem%Dinv(i,j,2,2)*v(i,j,2))
        1634             :        enddo
        1635             :     enddo
        1636             :
        1637             :     ! compute d/dx and d/dy
        1638 >11809*10^7 :     do j=1,np
        1639 >49600*10^7 :        do l=1,np
        1640             :           dudx00=0.0d0
        1641             :           dvdy00=0.0d0
        1642             : !DIR$ UNROLL(NP)
        1643 >18895*10^8 :           do i=1,np
        1644 >15116*10^8 :              dudx00 = dudx00 + deriv%Dvv(i,l  )*gv(i,j  ,1)
        1645 >18895*10^8 :              dvdy00 = dvdy00 + deriv%Dvv(i,l  )*gv(j  ,i,2)
        1646             :           end do
        1647 >37790*10^7 :           div(l  ,j  ) = dudx00
        1648 >47238*10^7 :           vvtemp(j  ,l  ) = dvdy00
        1649             :        end do
        1650             :     end do
        1651             :
        1652 >11809*10^7 :     do j=1,np
        1653 >49600*10^7 :        do i=1,np
        1654 >47238*10^7 :           div(i,j)=(div(i,j)+vvtemp(i,j))*(elem%rmetdet(i,j)*ra)
        1655             :        end do
        1656             :     end do
        1657             :
        1658 23619162600 :   end subroutine divergence_sphere
        1659             :
        1660             :
        1661             : !DIR$ ATTRIBUTES FORCEINLINE :: laplace_sphere_wk
        1662 19221591600 :   subroutine laplace_sphere_wk(s,deriv,elem,laplace,var_coef,mol_nu)
        1663             : !
        1664             : !   input:  s = scalar
        1665             : !   ouput:  -< grad(PHI), grad(s) >   = weak divergence of grad(s)
        1666             : !     note: for this form of the operator, grad(s) does not need to be made C0
        1667             : !
        1668             :     real(kind=r8), intent(in) :: s(np,np)
        1669             :     type (derivative_t), intent(in)     :: deriv
        1670             :     type (element_t), intent(in)        :: elem
        1671             :     real(kind=r8)                       :: laplace(np,np)
        1672             :     logical,       intent(in)           :: var_coef
        1673             :     real(kind=r8), intent(in), optional :: mol_nu(np,np)  !variable nu (e.g. molecular diffusion)
        1674             :
        1675             :     real(kind=r8)             :: laplace2(np,np)
        1676             :     integer i,j
        1677             :
        1678             :     ! Local
        1679             :     real(kind=r8) :: grads(np,np,2), oldgrads(np,np,2)
        1680             :
        1681 19221591600 :     call gradient_sphere(s,deriv,elem%Dinv,grads)
        1682             :     !
        1683             :     ! molecular diffusion coefficient
        1684             :     !
        1685 19221591600 :     if (present(mol_nu)) then
        1686           0 :       if (var_coef) &
        1687           0 :            call endrun('ERROR: this option is only for non-grid dependent variable viscosity')
        1688           0 :       grads(:,:,1) = grads(:,:,1)*mol_nu(:,:)
        1689           0 :       grads(:,:,2) = grads(:,:,2)*mol_nu(:,:)
        1690             :     end if
        1691             :
        1692 19221591600 :     if (var_coef) then
        1693 18847566000 :        if (hypervis_power/=0 ) then
        1694             :           ! scalar viscosity with variable coefficient
        1695           0 :           grads(:,:,1) = grads(:,:,1)*elem%variable_hyperviscosity(:,:)
        1696           0 :           grads(:,:,2) = grads(:,:,2)*elem%variable_hyperviscosity(:,:)
        1697 18847566000 :        else if (hypervis_scaling /=0 ) then
        1698             :           ! tensor hv, (3)
        1699           0 :           oldgrads=grads
        1700           0 :           do j=1,np
        1701           0 :              do i=1,np
        1702           0 :                 grads(i,j,1) = oldgrads(i,j,1)*elem%tensorVisc(i,j,1,1) + &
        1703           0 :                                oldgrads(i,j,2)*elem%tensorVisc(i,j,1,2)
        1704             :                 grads(i,j,2) = oldgrads(i,j,1)*elem%tensorVisc(i,j,2,1) + &
        1705           0 :                                oldgrads(i,j,2)*elem%tensorVisc(i,j,2,2)
        1706             :              end do
        1707             :           end do
        1708             :        else
        1709             :           ! do nothing: constant coefficient viscosity
        1710             :        endif
        1711             :     endif
        1712             :
        1713             :     ! note: divergnece_sphere and divergence_sphere_wk are identical *after* bndry_exchange
        1714             :     ! if input is C_0.  Here input is not C_0, so we should use divergence_sphere_wk().
        1715             :     ! laplace=divergence_sphere_wk(grads,deriv,elem)
        1716 19221591600 :     call divergence_sphere_wk(grads,deriv,elem,laplace)
        1717             :
        1718 19221591600 :   end subroutine laplace_sphere_wk
        1719             :
        1720             : !DIR$ ATTRIBUTES FORCEINLINE :: vlaplace_sphere_wk
        1721  8883108000 :   subroutine vlaplace_sphere_wk(v,deriv,elem,undamprrcart,laplace,var_coef,nu_ratio)
        1722             : !
        1723             : !   input:  v = vector in lat-lon coordinates
        1724             : !   ouput:  weak laplacian of v, in lat-lon coordinates
        1725             : !
        1726             : !   logic:
        1727             : !      tensorHV:     requires cartesian
        1728             : !      nu_div/=nu:   requires contra formulatino
        1729             : !
        1730             : !   One combination NOT supported:  tensorHV and nu_div/=nu then abort
        1731             : !
        1732             :     real(kind=r8),           intent(in)  :: v(np,np,2)
        1733             :     type (derivative_t),     intent(in)  :: deriv
        1734             :     type (element_t),        intent(in)  :: elem
        1735             :     logical,                 intent(in)  :: undamprrcart
        1736             :     real(kind=r8),           intent(out) :: laplace(np,np,2)
        1737             :     logical,       optional, intent(in)  :: var_coef
        1738             :     real(kind=r8), optional, intent(in)  :: nu_ratio
        1739             :
        1740             :
        1741  8883108000 :     if (hypervis_scaling/=0 .and. var_coef) then
        1742             :        ! tensorHV is turned on - requires cartesian formulation
        1743           0 :        if (present(nu_ratio)) then
        1744           0 :           if (nu_ratio /= 1._r8) then
        1745           0 :              call endrun('ERROR: tensorHV can not be used with nu_div/=nu')
        1746             :           endif
        1747             :        endif
        1748           0 :        laplace=vlaplace_sphere_wk_cartesian(v,deriv,elem,var_coef,undamprrcart)
        1749             :     else
        1750             :        ! all other cases, use contra formulation:
        1751  8883108000 :        laplace=vlaplace_sphere_wk_contra(v,deriv,elem,var_coef,undamprrcart,nu_ratio)
        1752             :     endif
        1753             :
        1754  8883108000 :   end subroutine vlaplace_sphere_wk
        1755             :   !
        1756             :   ! version of vlaplace_sphere_wk for molecular diffusion
        1757             :   !
        1758           0 :   subroutine vlaplace_sphere_wk_mol(v,deriv,elem,undamprrcart,mol_nu,laplace)
        1759             :     !
        1760             :     !   input:  v = vector in lat-lon coordinates
        1761             :     !   ouput:  weak laplacian of v, in lat-lon coordinates
        1762             :     !
        1763             :     real(kind=r8), intent(in)      :: v(np,np,2)
        1764             :     type (derivative_t), intent(in):: deriv
        1765             :     type (element_t), intent(in)   :: elem
        1766             :     logical,             intent(in) :: undamprrcart
        1767             :     real(kind=r8), intent(in)      :: mol_nu(np,np)
        1768             :     real(kind=r8), intent(out)     :: laplace(np,np,2)
        1769             :
        1770             :     real(kind=r8) :: vor(np,np),div(np,np)
        1771             :
        1772             :     integer :: n,m
        1773             :
        1774           0 :     call divergence_sphere(v,deriv,elem,div)
        1775           0 :     call vorticity_sphere(v,deriv,elem,vor)
        1776             :
        1777           0 :     div = div*mol_nu(:,:)
        1778           0 :     vor = vor*mol_nu(:,:)
        1779             :
        1780             :     laplace = gradient_sphere_wk_testcov(div,deriv,elem) - &
        1781           0 :          curl_sphere_wk_testcov(vor,deriv,elem)
        1782             :
        1783           0 :     if (undamprrcart) then
        1784           0 :       do n=1,np
        1785           0 :         do m=1,np
        1786             :           ! add in correction so we dont damp rigid rotation
        1787           0 :           laplace(m,n,1)=laplace(m,n,1) + 2*elem%spheremp(m,n)*v(m,n,1)*(ra**2)
        1788           0 :           laplace(m,n,2)=laplace(m,n,2) + 2*elem%spheremp(m,n)*v(m,n,2)*(ra**2)
        1789             :         enddo
        1790             :       enddo
        1791             :     end if
        1792             :
        1793             :
        1794           0 :   end subroutine vlaplace_sphere_wk_mol
        1795             :
        1796             :
        1797             :
        1798           0 :   function vlaplace_sphere_wk_cartesian(v,deriv,elem,var_coef,undamprrcart) result(laplace)
        1799             : !
        1800             : !   input:  v = vector in lat-lon coordinates
        1801             : !   ouput:  weak laplacian of v, in lat-lon coordinates
        1802             :
        1803             :     real(kind=r8),       intent(in) :: v(np,np,2)
        1804             :     type (derivative_t), intent(in) :: deriv
        1805             :     type (element_t),    intent(in) :: elem
        1806             :     logical,             intent(in) :: var_coef
        1807             :     logical,             intent(in) :: undamprrcart
        1808             :     real(kind=r8)                   :: laplace(np,np,2)
        1809             :     ! Local
        1810             :
        1811             :     integer component
        1812             :     real(kind=r8) :: dum_cart(np,np,3)
        1813             :     real(kind=r8) :: dum_cart2(np,np)
        1814             :
        1815             :
        1816             :     ! latlon -> cartesian
        1817           0 :     do component=1,3
        1818           0 :        dum_cart2(:,:) = elem%vec_sphere2cart(:,:,component,1)*v(:,:,1) + &
        1819           0 :                                 elem%vec_sphere2cart(:,:,component,2)*v(:,:,2)
        1820             :        ! Do laplace on cartesian comps
        1821           0 :        call laplace_sphere_wk(dum_cart2,deriv,elem,dum_cart(:,:,component),var_coef)
        1822             :     enddo
        1823             :
        1824             :     ! cartesian -> latlon
        1825           0 :     do component=1,2
        1826             :        ! vec_sphere2cart is its own pseudoinverse.
        1827           0 :        laplace(:,:,component) = dum_cart(:,:,1)*elem%vec_sphere2cart(:,:,1,component) + &
        1828             :                                 dum_cart(:,:,2)*elem%vec_sphere2cart(:,:,2,component) + &
        1829           0 :                                 dum_cart(:,:,3)*elem%vec_sphere2cart(:,:,3,component)
        1830             :     end do
        1831             :
        1832           0 :     if (undamprrcart) then
        1833             :       ! add in correction so we dont damp rigid rotation
        1834           0 :       laplace(:,:,1)=laplace(:,:,1) + 2*elem%spheremp(:,:)*v(:,:,1)*(ra**2)
        1835           0 :       laplace(:,:,2)=laplace(:,:,2) + 2*elem%spheremp(:,:)*v(:,:,2)*(ra**2)
        1836             :     end if
        1837             :
        1838           0 :   end function vlaplace_sphere_wk_cartesian
        1839             :
        1840             :
        1841             :
        1842  8883108000 :   function vlaplace_sphere_wk_contra(v,deriv,elem,var_coef,undamprrcart,nu_ratio) result(laplace)
        1843             : !
        1844             : !   input:  v = vector in lat-lon coordinates
        1845             : !   ouput:  weak laplacian of v, in lat-lon coordinates
        1846             : !
        1847             : !   du/dt = laplace(u) = grad(div) - curl(vor)
        1848             : !   < PHI du/dt > = < PHI laplace(u) >        PHI = covariant, u = contravariant
        1849             : !                 = < PHI grad(div) >  - < PHI curl(vor) >
        1850             : !                 = grad_wk(div) - curl_wk(vor)
        1851             : !
        1852             :     real(kind=r8),           intent(in) :: v(np,np,2)
        1853             :     logical,                 intent(in) :: var_coef
        1854             :     type (derivative_t),     intent(in) :: deriv
        1855             :     type (element_t),        intent(in) :: elem
        1856             :     logical,                 intent(in) :: undamprrcart
        1857             :     real(kind=r8), optional, intent(in) :: nu_ratio
        1858             :
        1859             :     real(kind=r8)                       :: laplace(np,np,2)
        1860             :
        1861             :     ! Local
        1862             :
        1863             :     integer i,j,l,m,n
        1864             :     real(kind=r8) :: vor(np,np),div(np,np)
        1865             :     real(kind=r8) :: v1,v2,div1,div2,vor1,vor2,phi_x,phi_y
        1866             :
        1867  8883108000 :     call divergence_sphere(v,deriv,elem,div)
        1868  8883108000 :     call vorticity_sphere(v,deriv,elem,vor)
        1869             :
        1870  8883108000 :     if (var_coef .and. hypervis_power/=0 ) then
        1871             :           ! scalar viscosity with variable coefficient
        1872           0 :           div = div*elem%variable_hyperviscosity(:,:)
        1873           0 :           vor = vor*elem%variable_hyperviscosity(:,:)
        1874             :     endif
        1875             :
        1876 >18654*10^7 :     if (present(nu_ratio)) div = nu_ratio*div
        1877             :
        1878             :     laplace = gradient_sphere_wk_testcov(div,deriv,elem) - &
        1879 >38197*10^7 :          curl_sphere_wk_testcov(vor,deriv,elem)
        1880             :
        1881  8883108000 :     if (undamprrcart) then
        1882 44415540000 :       do n=1,np
        1883 >18654*10^7 :         do m=1,np
        1884             :           ! add in correction so we dont damp rigid rotation
        1885 >14212*10^7 :           laplace(m,n,1)=laplace(m,n,1) + 2*elem%spheremp(m,n)*v(m,n,1)*(ra**2)
        1886 >17766*10^7 :           laplace(m,n,2)=laplace(m,n,2) + 2*elem%spheremp(m,n)*v(m,n,2)*(ra**2)
        1887             :         enddo
        1888             :       enddo
        1889             :     end if
        1890  8883108000 :   end function vlaplace_sphere_wk_contra
        1891             :
        1892           0 :   function gll_to_dgmodal(p,deriv) result(phat)
        1893             : !
        1894             : !   input:  v = velocity in lat-lon coordinates
        1895             : !   ouput:  phat = Legendre coefficients
        1896             : !
        1897             : !   Computes  < g dot p  > = SUM  g(i,j) p(i,j) w(i) w(j)
        1898             : !   (the quadrature approximation on the *reference element* of the integral of p against
        1899             : !    all Legendre polynomials up to degree npdg
        1900             : !
        1901             : !   for npdg < np, this routine gives the (exact) modal expansion of p/spheremp()
        1902             : !
        1903             :     real(kind=r8), intent(in) :: p(np,np)
        1904             :     type (derivative_t), intent(in) :: deriv
        1905             :     real(kind=r8) :: phat(npdg,npdg)
        1906             :
        1907             :     ! Local
        1908             :     integer i,j,m,n
        1909           0 :     real(kind=r8) :: A(np,npdg)
        1910           0 :     A=0
        1911           0 :     phat=0
        1912             :
        1913             :     ! N^3 tensor product formulation:
        1914           0 :     do m=1,npdg
        1915           0 :     do j=1,np
        1916             : !DIR$ UNROLL(NP)
        1917           0 :     do i=1,np
        1918           0 :        A(j,m)=A(j,m)+( p(i,j)*deriv%Mvv_twt(i,i)*deriv%Mvv_twt(j,j)  )*deriv%legdg(m,i)
        1919             :     enddo
        1920             :     enddo
        1921             :     enddo
        1922             :
        1923           0 :     do n=1,npdg
        1924           0 :     do m=1,npdg
        1925             : !DIR$ UNROLL(NP)
        1926           0 :     do j=1,np
        1927           0 :        phat(m,n)=phat(m,n)+A(j,m)*deriv%legdg(n,j)
        1928             :     enddo
        1929             :     enddo
        1930             :     enddo
        1931             :
        1932           0 :   end function
        1933             :
        1934           0 :   function dgmodal_to_gll(phat,deriv) result(p)
        1935             :     !
        1936             :     !   input:  phat = coefficients of Legendre expansion
        1937             :     !   ouput:  p    = sum expansion to evaluate phat at GLL points
        1938             :     !
        1939             :     real(kind=r8) :: p(np,np)
        1940             :     type (derivative_t), intent(in) :: deriv
        1941             :     real(kind=r8) :: phat(npdg,npdg)
        1942             :     ! Local
        1943             :     integer i,j,m,n
        1944           0 :     real(kind=r8) :: A(npdg,np)
        1945             :
        1946           0 :     p(:,:)=0
        1947             :     ! tensor product version
        1948           0 :     A=0
        1949           0 :     do i=1,np
        1950           0 :       do n=1,npdg
        1951           0 :         do m=1,npdg
        1952           0 :           A(n,i)=A(n,i)+phat(m,n)*deriv%legdg(m,i)
        1953             :         enddo
        1954             :       enddo
        1955             :     enddo
        1956           0 :     do j=1,np
        1957           0 :       do i=1,np
        1958           0 :         do n=1,npdg
        1959           0 :           p(i,j) = p(i,j)+A(n,i)*deriv%legdg(n,j)
        1960             :         enddo
        1961             :       enddo
        1962             :     enddo
        1963             :
        1964           0 :   end function dgmodal_to_gll
        1965             :
        1966  7433758800 :   subroutine subcell_dss_fluxes(dss, p, n, metdet, C, fluxes)
        1967             :
        1968             :     integer              , intent(in)  :: p
        1969             :     integer              , intent(in)  :: n
        1970             :     real (kind=r8), intent(in)  :: dss     (p,p)
        1971             :     real (kind=r8), intent(in)  :: metdet  (p,p)
        1972             :     real (kind=r8), intent(in)  :: C       (2,2,2)
        1973             :
        1974             :     real (kind=r8)              :: fluxes  (n,n,4)
        1975             :
        1976 14867517600 :     real (kind=r8)              :: Bp(p,p)
        1977 14867517600 :     real (kind=r8)              :: Tp(p,p)
        1978 14867517600 :     real (kind=r8)              :: Lp(p,p)
        1979 14867517600 :     real (kind=r8)              :: Rp(p,p)
        1980             :
        1981 14867517600 :     real (kind=r8)              :: B(n,n)
        1982 14867517600 :     real (kind=r8)              :: T(n,n)
        1983 14867517600 :     real (kind=r8)              :: L(n,n)
        1984  7433758800 :     real (kind=r8)              :: R(n,n)
        1985             :
        1986             :     integer            :: i,j
        1987             :
        1988 >39398*10^7 :     fluxes  = 0
        1989             :
        1990 >15610*10^7 :     Bp = 0
        1991 >15610*10^7 :     Tp = 0
        1992 >15610*10^7 :     Rp = 0
        1993 >15610*10^7 :     Lp = 0
        1994             :
        1995 44602552800 :     Bp(:,1)  = dss(:,1)  ! bottom
        1996 37168794000 :     Tp(:,p)  = dss(:,p)  ! top
        1997 37168794000 :     Rp(p,:)  = dss(p,:)  ! right
        1998 37168794000 :     Lp(1,:)  = dss(1,:)  ! left
        1999             :
        2000  7433758800 :     Bp(1,1)  = C(1,1,2)
        2001  7433758800 :     Lp(1,1)  = C(1,1,1)
        2002  7433758800 :     Bp(p,1)  = C(2,1,2)
        2003  7433758800 :     Rp(p,1)  = C(2,1,1)
        2004             :
        2005  7433758800 :     Tp(1,p)  = C(1,2,2)
        2006  7433758800 :     Lp(1,p)  = C(1,2,1)
        2007  7433758800 :     Tp(p,p)  = C(2,2,2)
        2008  7433758800 :     Rp(p,p)  = C(2,2,1)
        2009             :
        2010             :
        2011  7433758800 :     call subcell_integration(Bp, p, n, metdet,B)
        2012  7433758800 :     call subcell_integration(Tp, p, n, metdet,T)
        2013  7433758800 :     call subcell_integration(Lp, p, n, metdet,L)
        2014  7433758800 :     call subcell_integration(Rp, p, n, metdet,R)
        2015             :
        2016 29735035200 :     do i = 1,n
        2017 96638864400 :     do j = 1,n
        2018 66903829200 :       if (1<j) T(i,j) = T(i,j) + T(i,j-1)
        2019 89205105600 :       if (1<i) R(i,j) = R(i,j) + R(i-1,j)
        2020             :     end do
        2021             :     end do
        2022             :
        2023 29735035200 :     do i = n,1,-1
        2024 96638864400 :     do j = n,1,-1
        2025 66903829200 :       if (j<n) B(i,j) = B(i,j) + B(i,j+1)
        2026 89205105600 :       if (i<n) L(i,j) = L(i,j) + L(i+1,j)
        2027             :     end do
        2028             :     end do
        2029             :
        2030 29735035200 :     do i = 1,n
        2031 96638864400 :       do j = 1,n
        2032 66903829200 :         if (1==j) fluxes(i,j,1) =  B(i,j)
        2033 66903829200 :         if (n==i) fluxes(i,j,2) =  R(i,j)
        2034 66903829200 :         if (j==n) fluxes(i,j,3) =  T(i,j)
        2035 66903829200 :         if (1==i) fluxes(i,j,4) =  L(i,j)
        2036             :
        2037 66903829200 :         if (1< j) fluxes(i,j,1) =   B(i,j) - T(i,j-1)
        2038 66903829200 :         if (i< n) fluxes(i,j,2) =   R(i,j) - L(i+1,j)
        2039 66903829200 :         if (j< n) fluxes(i,j,3) =   T(i,j) - B(i,j+1)
        2040 89205105600 :         if (1< i) fluxes(i,j,4) =   L(i,j) - R(i-1,j)
        2041             :       end do
        2042             :     end do
        2043             :
        2044  7433758800 :   end subroutine subcell_dss_fluxes
        2045             :
        2046  2898698400 :   subroutine subcell_div_fluxes(u, p, n, metdet,fluxes)
        2047             :
        2048             :     implicit none
        2049             :
        2050             :     integer              , intent(in)  :: p
        2051             :     integer              , intent(in)  :: n
        2052             :     real (kind=r8), intent(in)  :: u(p,p,2)
        2053             :     real (kind=r8), intent(in)  :: metdet(p,p)
        2054             :
        2055  5797396800 :     real (kind=r8)              :: v(p,p,2)
        2056             :     real (kind=r8)              :: fluxes(n,n,4)
        2057  5797396800 :     real (kind=r8)              :: tb(n,p)
        2058  5797396800 :     real (kind=r8)              :: lr(p,n)
        2059  5797396800 :     real (kind=r8)              :: flux_l(n,n)
        2060  5797396800 :     real (kind=r8)              :: flux_r(n,n)
        2061  5797396800 :     real (kind=r8)              :: flux_b(n,n)
        2062  5797396800 :     real (kind=r8)              :: flux_t(n,n)
        2063             :
        2064             :     integer i,j
        2065             :
        2066             :     if (.not.ALLOCATED(integration_matrix)      .or. &
        2067  2898698400 :         SIZE(integration_matrix,1).ne.n .or. &
        2068             :         SIZE(integration_matrix,2).ne.p) then
        2069           0 :       call endrun( 'FATAL ERROR: allocate_subcell_integration_matrix not called')
        2070             :     end if
        2071             :
        2072 60872666400 :     v(:,:,1) = u(:,:,1)*metdet(:,:)
        2073 60872666400 :     v(:,:,2) = u(:,:,2)*metdet(:,:)
        2074             :
        2075  2898698400 :     tb = MATMUL(integration_matrix, v(:,:,2))
        2076  2898698400 :     flux_b(:,:) = MATMUL(tb,TRANSPOSE(boundary_interp_matrix(:,1,:)))
        2077  2898698400 :     flux_t(:,:) = MATMUL(tb,TRANSPOSE(boundary_interp_matrix(:,2,:)))
        2078             :
        2079  2898698400 :     lr = MATMUL(v(:,:,1),TRANSPOSE(integration_matrix))
        2080  2898698400 :     flux_l(:,:) = MATMUL(boundary_interp_matrix(:,1,:),lr)
        2081  2898698400 :     flux_r(:,:) = MATMUL(boundary_interp_matrix(:,2,:),lr)
        2082             :
        2083 37683079200 :     fluxes(:,:,1) = -flux_b(:,:)*ra
        2084 37683079200 :     fluxes(:,:,2) =  flux_r(:,:)*ra
        2085 37683079200 :     fluxes(:,:,3) =  flux_t(:,:)*ra
        2086 37683079200 :     fluxes(:,:,4) = -flux_l(:,:)*ra
        2087             :
        2088  2898698400 :   end subroutine subcell_div_fluxes
        2089             :
        2090  4535060400 :   subroutine subcell_Laplace_fluxes(u, deriv, elem, p, n, fluxes)
        2091             :
        2092             :     implicit none
        2093             :
        2094             :     integer              , intent(in)  :: p
        2095             :     integer              , intent(in)  :: n
        2096             :     type (derivative_t)  , intent(in)  :: deriv
        2097             :     type (element_t)     , intent(in)  :: elem
        2098             :     real (kind=r8), intent(in)  :: u(p,p)
        2099             :
        2100  9070120800 :     real (kind=r8)              :: g(p,p,2)
        2101  9070120800 :     real (kind=r8)              :: v(p,p,2)
        2102  9070120800 :     real (kind=r8)              :: div(p,p,2)
        2103  9070120800 :     real (kind=r8)              :: sub_int(n,n,2)
        2104             :     real (kind=r8)              :: fluxes(n,n,4)
        2105             :
        2106             :     integer i,j
        2107             :
        2108  4535060400 :     call gradient_sphere(u,deriv,elem%Dinv,g)
        2109             :
        2110 95236268400 :     v(:,:,1) = elem%Dinv(:,:,1,1)*g(:,:,1) + elem%Dinv(:,:,1,2)*g(:,:,2)
        2111 95236268400 :     v(:,:,2) = elem%Dinv(:,:,2,1)*g(:,:,1) + elem%Dinv(:,:,2,2)*g(:,:,2)
        2112             :
        2113 22675302000 :     do j=1,p
        2114 95236268400 :     do i=1,p
        2115 >36280*10^7 :        div(i,j,1) = -SUM(elem%spheremp(:,j)*v(:,j,1)*deriv%Dvv(i,:))
        2116 >38094*10^7 :        div(i,j,2) = -SUM(elem%spheremp(i,:)*v(i,:,2)*deriv%Dvv(j,:))
        2117             :     end do
        2118             :     end do
        2119 >19500*10^7 :     div = div * ra
        2120             :
        2121 95236268400 :     div(:,:,1) = div(:,:,1) / elem%spheremp(:,:)
        2122 95236268400 :     div(:,:,2) = div(:,:,2) / elem%spheremp(:,:)
        2123             :
        2124  4535060400 :     call subcell_integration(div(:,:,1), p, n, elem%metdet, sub_int(:,:,1))
        2125  4535060400 :     call subcell_integration(div(:,:,2), p, n, elem%metdet, sub_int(:,:,2))
        2126             :
        2127 18140241600 :     do i=1,n
        2128 45350604000 :     do j=2,n
        2129 27210362400 :       sub_int(j,i,1) = sub_int(j,i,1) + sub_int(j-1,i,1)
        2130 40815543600 :       sub_int(i,j,2) = sub_int(i,j,2) + sub_int(i,j-1,2)
        2131             :     end do
        2132             :     end do
        2133             :
        2134 >24035*10^7 :     fluxes = 0
        2135 18140241600 :     do i=1,n
        2136 58955785200 :     do j=1,n
        2137 40815543600 :       if (1.lt.j) fluxes(i,j,1) = -sub_int(i,  j-1,2)
        2138 40815543600 :       if (i.lt.n) fluxes(i,j,2) =  sub_int(i,  j,  1)
        2139 40815543600 :       if (j.lt.n) fluxes(i,j,3) =  sub_int(i,  j,  2)
        2140 54420724800 :       if (1.lt.i) fluxes(i,j,4) = -sub_int(i-1,j,  1)
        2141             :     end do
        2142             :     end do
        2143             :
        2144  4535060400 :   end subroutine subcell_Laplace_fluxes
        2145             :
        2146             :
        2147             :   ! Given a field defined on the unit element, [-1,1]x[-1,1]
        2148             :   ! sample values, sampled_val, premultiplied by integration weights,
        2149             :   ! and a number, np, of Gauss-Lobatto-Legendre points. Divide
        2150             :   ! the square up into intervals by intervals sub-squares so that
        2151             :   ! there are now intervals**2 sub-cells.  Integrate the
        2152             :   ! function defined by sampled_val over each of these
        2153             :   ! sub-cells and return the integrated values as an
        2154             :   ! intervals by intervals matrix.
        2155             :   !
        2156             :   ! Efficiency is obtained by computing and caching the appropriate
        2157             :   ! integration matrix the first time the function is called.
        2158 39789279000 :   subroutine subcell_integration(sampled_val, np, intervals, metdet,values,separate_physics_grid)
        2159             :
        2160             :     implicit none
        2161             :
        2162             :     integer              , intent(in)  :: np
        2163             :     integer              , intent(in)  :: intervals
        2164             :     real (kind=r8), intent(in)  :: sampled_val(np,np)
        2165             :     real (kind=r8), intent(in)  :: metdet(np,np)
        2166 79578558000 :     real (kind=r8)              :: val(np,np)
        2167             :     real (kind=r8)              :: values(intervals,intervals)
        2168             :
        2169             :     logical, optional :: separate_physics_grid
        2170             :
        2171             :     logical :: ltmp
        2172             :     integer i,j
        2173             :
        2174 39789279000 :     if (present(separate_physics_grid)) then
        2175   978075000 :       ltmp = separate_physics_grid
        2176             :     else
        2177             :       ltmp = .false.
        2178             :     end if
        2179             :
        2180             :
        2181   978075000 :     if (ltmp) then
        2182             :         if (.not.ALLOCATED(integration_matrix_physgrid)      .or. &
        2183           0 :              SIZE(integration_matrix_physgrid,1).ne.intervals .or. &
        2184             :              SIZE(integration_matrix_physgrid,2).ne.np) then
        2185           0 :           call endrun( 'FATAL ERROR: allocate_subcell_integration_matrix_physgrid not called')
        2186             :         end if
        2187           0 :         val = sampled_val * metdet
        2188             :         values = MATMUL(integration_matrix_physgrid, &
        2189           0 :              MATMUL(val,TRANSPOSE(integration_matrix_physgrid)))
        2190             :     else
        2191             :       if (.not.ALLOCATED(integration_matrix)      .or. &
        2192 39789279000 :            SIZE(integration_matrix,1).ne.intervals .or. &
        2193             :            SIZE(integration_matrix,2).ne.np) then
        2194           0 :         call endrun( 'FATAL ERROR: allocate_subcell_integration_matrix not called')
        2195             :       end if
        2196             :       ! Multiply sampled values by spectral element weights
        2197 >83557*10^7 :       val = sampled_val * metdet
        2198             :
        2199             :       ! Multiply the sampled values by the weighted jacobians.
        2200             :       ! Symmetry allows us to write this as J^t V J
        2201             :       ! where J is a vector.
        2202             :
        2203             :       values = MATMUL(integration_matrix, &
        2204 39789279000 :            MATMUL(val,TRANSPOSE(integration_matrix)))
        2205             :
        2206             :     end if
        2207 39789279000 :   end subroutine subcell_integration
        2208             :
        2209             :
        2210             :   ! Helper subroutine that will fill in a matrix needed to
        2211             :   ! integrate a function defined on the GLL points of a unit
        2212             :   ! square on sub-cells.  So np is the number of integration
        2213             :   ! GLL points defined on the unit square (actually [-1,1]x[-1,1])
        2214             :   ! and intervals is the number to cut it up into, say a 3 by 3
        2215             :   ! set of uniform sub-cells.  This function will fill the
        2216             :   ! subcell_integration matrix with the correct coefficients
        2217             :   ! to integrate over each subcell.
        2218        1536 :   subroutine allocate_subcell_integration_matrix_cslam(np, intervals)
        2219             :     !-----------------
        2220             :     !-----------------
        2221             :     use quadrature_mod, only : gausslobatto, quadrature_t
        2222             :
        2223             :     implicit none
        2224             :
        2225             :     integer              , intent(in)  :: np
        2226             :     integer              , intent(in)  :: intervals
        2227             :     real (kind=r8)              :: values(intervals,intervals)
        2228             :
        2229             :
        2230             :     real(kind=r8), parameter :: zero = 0.0_r8, one=1.0_r8, two=2.0_r8
        2231             :
        2232             :
        2233        3072 :     real (kind=r8) :: sub_gll        (intervals,np)
        2234             :
        2235        3072 :     real (kind=r8) :: Lagrange_interp(intervals,np,np)
        2236             :     type (quadrature_t)   :: gll
        2237             :
        2238        3072 :     real (kind=r8) :: legrange_div(np)
        2239             :     real (kind=r8) :: a,b,x,y, x_j, x_i
        2240             :     real (kind=r8) :: r(1)
        2241             :     integer i,j,n,m
        2242             :
        2243        1536 :     if (ALLOCATED(integration_matrix)) deallocate(integration_matrix)
        2244        6144 :     allocate(integration_matrix(intervals,np))
        2245        1536 :     if (ALLOCATED(boundary_interp_matrix)) deallocate(boundary_interp_matrix)
        2246        6144 :     allocate(boundary_interp_matrix(intervals,2,np))
        2247             :
        2248        1536 :     gll = gausslobatto(np)
        2249             :
        2250             :     ! The GLL (Gauss-Lobatto-Legendre) points are from [-1,1],
        2251             :     ! we have a bunch of sub-intervals defined by intervals that
        2252             :     ! go from [a,b] so we need to linearly map [-1,1] -> [a,b]
        2253             :     ! all the  GLL points by  y = (a/2)(1-x) + (b/2)(1+x)
        2254        6144 :     do i=1,intervals
        2255        4608 :       a = -one + (i-one)*two/intervals
        2256        4608 :       b = -one +  i     *two/intervals
        2257       24576 :       sub_gll(i,:) = (a+b)/two + gll%points(:)/intervals
        2258             :     end do
        2259             :
        2260             :     ! Now to interpolate from the values at the input GLL
        2261             :     ! points to the sub-GLL points.  Do this by Lagrange
        2262             :     ! interpolation.  The jth Lagrange interpolating polynomial
        2263             :     ! for points x_i is
        2264             :     !              \prod_{i\ne j} (x-x_i)/(x_j-x_i)
        2265             :     ! These are then multiplied by the sampled values y_i
        2266             :     ! and summed.
        2267             :
        2268             :     ! Save some time by pre-computing the denominitor. I think
        2269             :     ! this is OK since all the points are of order 1 so should
        2270             :     ! be well behaved.
        2271        7680 :     do n = 1,np
        2272        6144 :       x_j = gll%points(n)
        2273        6144 :       x   = one
        2274       30720 :       do m = 1,np
        2275       30720 :         if (m.ne.n) then
        2276       18432 :           x_i = gll%points(m)
        2277       18432 :           x = x * (x_j-x_i)
        2278             :         endif
        2279             :       end do
        2280        7680 :       legrange_div(n)= x
        2281             :     end do
        2282        6144 :     do i=1,intervals
        2283       24576 :       do n=1,np
        2284       18432 :         x = sub_gll(i,n)
        2285       96768 :         do j = 1,np
        2286             :           y = one
        2287      368640 :           do m = 1,np
        2288      368640 :             if (m.ne.j) then
        2289      221184 :               x_i = gll%points(m)
        2290      221184 :               y = y * (x-x_i)
        2291             :             end if
        2292             :           end do
        2293       92160 :           Lagrange_interp(i,n,j) = y/legrange_div(j)
        2294             :         end do
        2295             :       end do
        2296             :     end do
        2297             :
        2298             :     ! Integration is the GLL weights times Jacobians times
        2299             :     ! the interpolated values:
        2300             :     !                   w^t I Y I^t w
        2301             :     ! where
        2302             :     ! w is GLL weights and Jacobians,
        2303             :     ! I is the Lagrange_interp matrix, and
        2304             :     ! Y is the coefficient matrix, sampled_val.
        2305             :     ! This can be written  J Y J^t where
        2306             :     !                       J = w^t I
        2307             :     ! J is integration_matrix
        2308        6144 :     do i=1,intervals
        2309       24576 :       integration_matrix(i,:) = MATMUL(gll%weights(:),Lagrange_interp(i,:,:))
        2310             :     end do
        2311             :
        2312             :     ! There is still the Jacobian to consider.  We are
        2313             :     ! integrating over [a,b] x [c,d] where
        2314             :     !        |b-a| = |d-c| = 2/Intervals
        2315             :     ! Multiply the weights appropriately given that
        2316             :     ! they are defined for a 2x2 square
        2317       26112 :     integration_matrix = integration_matrix/intervals
        2318             :
        2319       26112 :     boundary_interp_matrix(:,1,:) = Lagrange_interp(:,1,:)
        2320       26112 :     boundary_interp_matrix(:,2,:) = Lagrange_interp(:,np,:)
        2321        1536 :   end subroutine allocate_subcell_integration_matrix_cslam
        2322             :
        2323        1536 :   subroutine allocate_subcell_integration_matrix_physgrid(np, intervals)
        2324             :     !-----------------
        2325             :     !-----------------
        2326             :     use quadrature_mod, only : gausslobatto, quadrature_t
        2327             :
        2328             :     implicit none
        2329             :
        2330             :     integer              , intent(in)  :: np
        2331             :     integer              , intent(in)  :: intervals
        2332             :     real (kind=r8)              :: values(intervals,intervals)
        2333             :
        2334             :
        2335             :     real(kind=r8), parameter :: zero = 0.0_r8, one=1.0_r8, two=2.0_r8
        2336             :
        2337             :
        2338        3072 :     real (kind=r8) :: sub_gll        (intervals,np)
        2339             :
        2340        3072 :     real (kind=r8) :: Lagrange_interp(intervals,np,np)
        2341             :     type (quadrature_t)   :: gll
        2342             :
        2343        3072 :     real (kind=r8) :: legrange_div(np)
        2344             :     real (kind=r8) :: a,b,x,y, x_j, x_i
        2345             :     real (kind=r8) :: r(1)
        2346             :     integer i,j,n,m
        2347             :
        2348        1536 :     if (ALLOCATED(integration_matrix_physgrid)) deallocate(integration_matrix_physgrid)
        2349        6144 :     allocate(integration_matrix_physgrid(intervals,np))
        2350             :
        2351        1536 :     gll = gausslobatto(np)
        2352             :
        2353             :     ! The GLL (Gauss-Lobatto-Legendre) points are from [-1,1],
        2354             :     ! we have a bunch of sub-intervals defined by intervals that
        2355             :     ! go from [a,b] so we need to linearly map [-1,1] -> [a,b]
        2356             :     ! all the  GLL points by  y = (a/2)(1-x) + (b/2)(1+x)
        2357        6144 :     do i=1,intervals
        2358        4608 :       a = -one + (i-one)*two/intervals
        2359        4608 :       b = -one +  i     *two/intervals
        2360       24576 :       sub_gll(i,:) = (a+b)/two + gll%points(:)/intervals
        2361             :     end do
        2362             :
        2363             :     ! Now to interpolate from the values at the input GLL
        2364             :     ! points to the sub-GLL points.  Do this by Lagrange
        2365             :     ! interpolation.  The jth Lagrange interpolating polynomial
        2366             :     ! for points x_i is
        2367             :     !              \prod_{i\ne j} (x-x_i)/(x_j-x_i)
        2368             :     ! These are then multiplied by the sampled values y_i
        2369             :     ! and summed.
        2370             :
        2371             :     ! Save some time by pre-computing the denominitor. I think
        2372             :     ! this is OK since all the points are of order 1 so should
        2373             :     ! be well behaved.
        2374        7680 :     do n = 1,np
        2375        6144 :       x_j = gll%points(n)
        2376        6144 :       x   = one
        2377       30720 :       do m = 1,np
        2378       30720 :         if (m.ne.n) then
        2379       18432 :           x_i = gll%points(m)
        2380       18432 :           x = x * (x_j-x_i)
        2381             :         endif
        2382             :       end do
        2383        7680 :       legrange_div(n)= x
        2384             :     end do
        2385        6144 :     do i=1,intervals
        2386       24576 :       do n=1,np
        2387       18432 :         x = sub_gll(i,n)
        2388       96768 :         do j = 1,np
        2389             :           y = one
        2390      368640 :           do m = 1,np
        2391      368640 :             if (m.ne.j) then
        2392      221184 :               x_i = gll%points(m)
        2393      221184 :               y = y * (x-x_i)
        2394             :             end if
        2395             :           end do
        2396       92160 :           Lagrange_interp(i,n,j) = y/legrange_div(j)
        2397             :         end do
        2398             :       end do
        2399             :     end do
        2400        6144 :     do i=1,intervals
        2401       24576 :       integration_matrix_physgrid(i,:) = MATMUL(gll%weights(:),Lagrange_interp(i,:,:))
        2402             :     end do
        2403       26112 :     integration_matrix_physgrid = integration_matrix_physgrid/intervals
        2404        1536 :   end subroutine allocate_subcell_integration_matrix_physgrid
        2405             :
        2406             :
        2407             :
        2408           0 :   subroutine limiter_optim_iter_full(ptens,sphweights,minp,maxp,dpmass,kbeg,kend)
        2409             :     !
        2410             :     !The idea here is the following: We need to find a grid field which is closest
        2411             :     !to the initial field (in terms of weighted sum), but satisfies the min/max constraints.
        2412             :     !So, first we find values which do not satisfy constraints and bring these values
        2413             :     !to a closest constraint. This way we introduce some mass change (addmass),
        2414             :     !so, we redistribute addmass in the way that l2 error is smallest.
        2415             :     !This redistribution might violate constraints thus, we do a few iterations.
        2416             :     !
        2417             :     ! O. Guba ~2012                    Documented in Guba, Taylor & St-Cyr, JCP 2014
        2418             :     ! I. Demeshko & M. Taylor 7/2015:  Removed indirect addressing.
        2419             :     ! N. Lopez & M. Taylor 8/2015:     Mass redistributon tweak which is better at
        2420             :     !                                  linear coorelation preservation
        2421             :     !
        2422             :     use dimensions_mod, only : np, np, nlev
        2423             :
        2424             :     real (kind=r8), dimension(nlev), intent(inout)   :: minp, maxp
        2425             :     real (kind=r8), dimension(np*np,nlev), intent(inout)   :: ptens
        2426             :     real (kind=r8), dimension(np*np,nlev), intent(in), optional  :: dpmass
        2427             :     real (kind=r8), dimension(np*np), intent(in)   :: sphweights
        2428             :     integer, intent(in) :: kbeg, kend
        2429             :
        2430             :     real (kind=r8), dimension(np,np) :: ptens_mass
        2431             :     integer  k1, k, i, j, iter, weightsnum
        2432             :     real (kind=r8) :: addmass, weightssum, mass, sumc
        2433             :     real (kind=r8) :: x(np*np),c(np*np)
        2434             :     integer :: maxiter = np*np-1
        2435             :     real (kind=r8) :: tol_limiter = 5.0e-14_r8
        2436             :
        2437           0 :     do k = kbeg, kend
        2438             :
        2439           0 :      do k1=1,np*np
        2440           0 :        c(k1)=sphweights(k1)*dpmass(k1,k)
        2441           0 :        x(k1)=ptens(k1,k)/dpmass(k1,k)
        2442             :      enddo
        2443             :
        2444           0 :      sumc=sum(c)
        2445           0 :      if (sumc <= 0 ) CYCLE   ! this should never happen, but if it does, dont limit
        2446           0 :      mass=sum(c*x)
        2447             :
        2448             :
        2449             :
        2450             :       ! relax constraints to ensure limiter has a solution:
        2451             :       ! This is only needed if runnign with the SSP CFL>1 or
        2452             :       ! due to roundoff errors
        2453           0 :       if( mass < minp(k)*sumc ) then
        2454           0 :         minp(k) = mass / sumc
        2455             :       endif
        2456           0 :       if( mass > maxp(k)*sumc ) then
        2457           0 :         maxp(k) = mass / sumc
        2458             :       endif
        2459             :
        2460             :
        2461             :
        2462           0 :      do iter=1,maxiter
        2463             :
        2464             :       addmass=0.0d0
        2465             :
        2466           0 :        do k1=1,np*np
        2467           0 :          if((x(k1)>maxp(k))) then
        2468           0 :            addmass=addmass+(x(k1)-maxp(k))*c(k1)
        2469           0 :            x(k1)=maxp(k)
        2470             :          endif
        2471           0 :          if((x(k1)<minp(k))) then
        2472           0 :            addmass=addmass-(minp(k)-x(k1))*c(k1)
        2473           0 :            x(k1)=minp(k)
        2474             :          endif
        2475             :        enddo !k1
        2476             :
        2477           0 :        if(abs(addmass)<=tol_limiter*abs(mass)) exit
        2478             :
        2479           0 :        weightssum=0.0d0
        2480           0 :        if(addmass>0)then
        2481           0 :         do k1=1,np*np
        2482           0 :           if(x(k1)<maxp(k))then
        2483           0 :             weightssum=weightssum+c(k1)
        2484             :           endif
        2485             :         enddo !k1
        2486           0 :         do k1=1,np*np
        2487           0 :           if(x(k1)<maxp(k))then
        2488           0 :               x(k1)=x(k1)+addmass/weightssum
        2489             :           endif
        2490             :         enddo
        2491             :       else
        2492           0 :         do k1=1,np*np
        2493           0 :           if(x(k1)>minp(k))then
        2494           0 :             weightssum=weightssum+c(k1)
        2495             :           endif
        2496             :         enddo
        2497           0 :         do k1=1,np*np
        2498           0 :           if(x(k1)>minp(k))then
        2499           0 :             x(k1)=x(k1)+addmass/weightssum
        2500             :           endif
        2501             :         enddo
        2502             :       endif
        2503             :
        2504             :
        2505             :    enddo!end of iteration
        2506             :
        2507           0 :    do k1=1,np*np
        2508           0 :       ptens(k1,k)=x(k1)
        2509             :    enddo
        2510             :
        2511             :   enddo
        2512             :
        2513           0 :   do k = kbeg, kend
        2514           0 :     do k1=1,np*np
        2515           0 :       ptens(k1,k)=ptens(k1,k)*dpmass(k1,k)
        2516             :     enddo
        2517             :   enddo
        2518             :
        2519           0 :   end subroutine limiter_optim_iter_full
        2520             :
        2521             :
        2522             :
        2523             :
        2524             :
        2525           0 : end module derivative_mod

Generated by: LCOV version 1.14