LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - cube_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 930 | 1049 | 88.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 19 | 25 | 76.0 %

              Line data    Source code
           1             : module cube_mod
           2             :   use shr_kind_mod,           only: r8=>shr_kind_r8
           3             :   use coordinate_systems_mod, only: spherical_polar_t, cartesian3D_t, cartesian2d_t, &
           4             :        projectpoint, cubedsphere2cart, spherical_to_cart, sphere_tri_area,dist_threshold, &
           5             :        change_coordinates
           6             :
           7             :   use physconst,              only: pi, rearth
           8             :   use control_mod,            only: hypervis_scaling, cubed_sphere_map
           9             :   use cam_abortutils,         only: endrun
          10             :
          11             :   implicit none
          12             :   private
          13             :
          14             :   integer,public, parameter :: nfaces = 6          ! number of faces on the cube
          15             :   integer,public, parameter :: nInnerElemEdge = 8  ! number of edges for an interior element
          16             :   integer,public, parameter :: nCornerElemEdge = 4 ! number of corner elements
          17             :
          18             :   real(kind=r8), public, parameter :: cube_xstart = -0.25_R8*PI
          19             :   real(kind=r8), public, parameter :: cube_xend   =  0.25_R8*PI
          20             :   real(kind=r8), public, parameter :: cube_ystart = -0.25_R8*PI
          21             :   real(kind=r8), public, parameter :: cube_yend   =  0.25_R8*PI
          22             :
          23             :
          24             :   type, public :: face_t
          25             :      type (spherical_polar_t) :: sphere0       ! tangent point of face on sphere
          26             :      type (spherical_polar_t) :: sw            ! sw corner of face on sphere
          27             :      type (spherical_polar_t) :: se            ! se corner of face on sphere
          28             :      type (spherical_polar_t) :: ne            ! ne corner of face on sphere
          29             :      type (spherical_polar_t) :: nw            ! nw corner of face on sphere
          30             :      type (cartesian3D_t)     :: P0
          31             :      type (cartesian3D_t)     :: X0
          32             :      type (cartesian3D_t)     :: Y0
          33             :      integer                  :: number
          34             :      integer                  :: padding       ! pad the struct
          35             :   end type face_t
          36             :
          37             :   type, public :: cube_face_coord_t
          38             :      real(r8) :: x             ! x coordinate
          39             :      real(r8) :: y             ! y coordinate
          40             :      type (face_t), pointer :: face     ! face
          41             :   end type cube_face_coord_t
          42             :
          43             :   ! ==========================================
          44             :   ! Public Interfaces
          45             :   ! ==========================================
          46             :
          47             :   public :: CubeTopology
          48             :
          49             :   ! Rotate the North Pole:  used for JW baroclinic test case
          50             :   ! Settings this only changes Coriolis.
          51             :   ! User must also rotate initial condition
          52             :   real (kind=r8), public :: rotate_grid = 0
          53             :
          54             :   ! ===============================
          55             :   ! Public methods for cube
          56             :   ! ===============================
          57             :
          58             :   public  :: cube_init_atomic
          59             :   public  :: convert_gbl_index
          60             :   public  :: vmap,dmap
          61             :   public  :: covariant_rot
          62             :   public  :: contravariant_rot
          63             :   public  :: set_corner_coordinates
          64             :   public  :: assign_node_numbers_to_elem
          65             :
          66             :
          67             :   public  :: CubeEdgeCount
          68             :   public  :: CubeElemCount
          69             :   public  :: CubeSetupEdgeIndex
          70             :   public  :: rotation_init_atomic
          71             :   public  :: ref2sphere
          72             :
          73             :   ! ===============================
          74             :   ! Private methods
          75             :   ! ===============================
          76             :   private :: coordinates_atomic
          77             :   private :: metric_atomic
          78             :   private :: coreolis_init_atomic
          79             :
          80             : contains
          81             :
          82             :   ! =======================================
          83             :   !  cube_init_atomic:
          84             :   !
          85             :   ! Initialize element descriptors for
          86             :   ! cube sphere case for each element ...
          87             :   ! =======================================
          88       21600 :   subroutine cube_init_atomic(elem, gll_points, alpha_in)
          89             :     use element_mod, only : element_t
          90             :     use dimensions_mod, only : np
          91             :
          92             :     type (element_t),   intent(inout) :: elem
          93             :     real(r8),           intent(in)    :: gll_points(np)
          94             :     real(r8), optional, intent(in)    :: alpha_in
          95             :
          96             :     real(r8)                          :: alpha
          97             :
          98       21600 :     if (present(alpha_in)) then
          99       10800 :       alpha = alpha_in
         100             :     else
         101       10800 :       alpha = 1.0_r8
         102             :     end if
         103             :
         104       21600 :     elem%FaceNum=elem%vertex%face_number
         105       21600 :     call coordinates_atomic(elem,gll_points)
         106             :
         107       21600 :     call metric_atomic(elem, gll_points, alpha)
         108             :
         109       21600 :     call coreolis_init_atomic(elem)
         110       21600 :     elem%desc%use_rotation= 0
         111             :
         112       21600 :   end subroutine cube_init_atomic
         113             :
         114             :   ! =======================================
         115             :   ! coordinates_atomic:
         116             :   !
         117             :   ! Initialize element coordinates for
         118             :   ! cube-sphere case ... (atomic)
         119             :   !
         120             :   ! =======================================
         121             :
         122       21600 :   subroutine coordinates_atomic(elem, gll_points)
         123             :     use element_mod,    only: element_t, element_var_coordinates
         124             :     use dimensions_mod, only: np
         125             :
         126             :     type(element_t), intent(inout) :: elem
         127             :     real(r8),        intent(in)    :: gll_points(np)
         128             :
         129             :     real(r8)             :: area1,area2
         130             :     type (cartesian3d_t) :: quad(4)
         131             :     integer              :: face_no,i,j
         132             :
         133       21600 :     face_no = elem%vertex%face_number
         134             :     ! compute the corners in Cartesian coordinates
         135      108000 :     do i=1,4
         136      108000 :        elem%corners3D(i)=cubedsphere2cart(elem%corners(i),face_no)
         137             :     enddo
         138             :
         139             :     ! =========================================
         140             :     ! compute lat/lon coordinates of each GLL point
         141             :     ! =========================================
         142      108000 :     do i=1,np
         143      453600 :     do j=1,np
         144      432000 :        elem%spherep(i,j)=ref2sphere(gll_points(i),gll_points(j),elem%corners3D,cubed_sphere_map,elem%corners,elem%facenum)
         145             :     enddo
         146             :     enddo
         147             :
         148             :     ! also compute the [-pi/2,pi/2] cubed sphere coordinates:
         149      453600 :     elem%cartp=element_var_coordinates(elem%corners,gll_points)
         150             :
         151             :     ! Matrix describing vector conversion to cartesian
         152             :     ! Zonal direction
         153      453600 :     elem%vec_sphere2cart(:,:,1,1) = -SIN(elem%spherep(:,:)%lon)
         154      453600 :     elem%vec_sphere2cart(:,:,2,1) =  COS(elem%spherep(:,:)%lon)
         155      453600 :     elem%vec_sphere2cart(:,:,3,1) =  0.0_r8
         156             :     ! Meridional direction
         157      453600 :     elem%vec_sphere2cart(:,:,1,2) = -SIN(elem%spherep(:,:)%lat)*COS(elem%spherep(:,:)%lon)
         158      453600 :     elem%vec_sphere2cart(:,:,2,2) = -SIN(elem%spherep(:,:)%lat)*SIN(elem%spherep(:,:)%lon)
         159      453600 :     elem%vec_sphere2cart(:,:,3,2) =  COS(elem%spherep(:,:)%lat)
         160             :
         161       21600 :   end subroutine coordinates_atomic
         162             :
         163             :   ! elem_jacobians:
         164             :   !
         165             :   ! Calculate Jacobian associated with mapping
         166             :   ! from arbitrary quadrilateral to [-1,1]^2
         167             :   ! along with its inverse and determinant
         168             :   ! ==========================================
         169             :
         170       21600 :   subroutine elem_jacobians(coords, unif2quadmap)
         171             :
         172             :     use dimensions_mod, only : np
         173             :     type (cartesian2D_t),  dimension(np,np), intent(in) :: coords
         174             :     ! unif2quadmap is the bilinear map from [-1,1]^2 -> arbitrary quadrilateral
         175             :     real (kind=r8), dimension(4,2), intent(out) :: unif2quadmap
         176             :     integer :: ii,jj
         177             :
         178       21600 :     unif2quadmap(1,1)=(coords(1,1)%x+coords(np,1)%x+coords(np,np)%x+coords(1,np)%x)/4.0_r8
         179       21600 :     unif2quadmap(1,2)=(coords(1,1)%y+coords(np,1)%y+coords(np,np)%y+coords(1,np)%y)/4.0_r8
         180       21600 :     unif2quadmap(2,1)=(-coords(1,1)%x+coords(np,1)%x+coords(np,np)%x-coords(1,np)%x)/4.0_r8
         181       21600 :     unif2quadmap(2,2)=(-coords(1,1)%y+coords(np,1)%y+coords(np,np)%y-coords(1,np)%y)/4.0_r8
         182       21600 :     unif2quadmap(3,1)=(-coords(1,1)%x-coords(np,1)%x+coords(np,np)%x+coords(1,np)%x)/4.0_r8
         183       21600 :     unif2quadmap(3,2)=(-coords(1,1)%y-coords(np,1)%y+coords(np,np)%y+coords(1,np)%y)/4.0_r8
         184       21600 :     unif2quadmap(4,1)=(coords(1,1)%x-coords(np,1)%x+coords(np,np)%x-coords(1,np)%x)/4.0_r8
         185       21600 :     unif2quadmap(4,2)=(coords(1,1)%y-coords(np,1)%y+coords(np,np)%y-coords(1,np)%y)/4.0_r8
         186             :
         187       21600 :   end subroutine elem_jacobians
         188             :
         189             :   ! =========================================
         190             :   ! metric_atomic:
         191             :   !
         192             :   ! Initialize cube-sphere metric terms:
         193             :   ! equal angular elements (atomic)
         194             :   ! initialize:
         195             :   !         metdet, rmetdet  (analytic)    = detD, 1/detD
         196             :   !         met                (analytic)    D^t D     (symmetric)
         197             :   !         metdet             (analytic)    = detD
         198             :   !         metinv             (analytic)    Dinv Dinv^t  (symmetic)
         199             :   !         D     (from subroutine vmap)
         200             :   !         Dinv  (computed directly from D)
         201             :   !
         202             :   ! ucontra = Dinv * u  =  metinv * ucov
         203             :   ! ucov    = D^t * u   =  met * ucontra
         204             :   !
         205             :   ! we also compute DE = D*E, where
         206             :   ! E = eigenvectors of metinv as a basis      metinv = E LAMBDA E^t
         207             :   !
         208             :   ! ueig = E^t ucov  = E^t D^t u =  (DE)^t u
         209             :   !
         210             :   !
         211             :   ! so if we want to tweak the mapping by a factor alpha (so he weights add up to 4pi, for example)
         212             :   ! we take:
         213             :   !    NEW       OLD
         214             :   !       D = sqrt(alpha) D  and then rederive all quantities.
         215             :   !    detD = alpha detD
         216             :   !
         217             :   ! where alpha = 4pi/SEMarea, SEMarea = global sum elem(ie)%mv(i,j)*elem(ie)%metdet(i,j)
         218             :   !
         219             :   ! =========================================
         220             :
         221       21600 :   subroutine metric_atomic(elem,gll_points,alpha)
         222             :     use element_mod,    only: element_t
         223             :     use dimensions_mod, only: np
         224             :     use physconst,      only: ra
         225             :
         226             :     type (element_t), intent(inout) :: elem
         227             :     real(r8),         intent(in)    :: alpha
         228             :     real(r8),         intent(in)    :: gll_points(np)
         229             :
         230             :     ! Local variables
         231             :     integer ii
         232             :     integer i,j,nn
         233             :     integer iptr
         234             :
         235             :     real (kind=r8) :: r         ! distance from origin for point on cube tangent to unit sphere
         236             :
         237             :     real (kind=r8) :: const, norm
         238             :     real (kind=r8) :: detD      ! determinant of vector field mapping matrix.
         239             :
         240             :     real (kind=r8) :: x1        ! 1st cube face coordinate
         241             :     real (kind=r8) :: x2        ! 2nd cube face coordinate
         242             :     real (kind=r8) :: tmpD(2,2)
         243             :     real (kind=r8) :: M(2,2),E(2,2),eig(2),DE(2,2),DEL(2,2),V(2,2), nu1, nu2, lamStar1, lamStar2
         244             :     integer :: imaxM(2)
         245             :     real (kind=r8) :: l1, l2, sc,min_svd,max_svd,max_normDinv
         246             :
         247             :
         248             :     ! ==============================================
         249             :     ! Initialize differential mapping operator
         250             :     ! to and from vector fields on the sphere to
         251             :     ! contravariant vector fields on the cube
         252             :     ! i.e. dM/dx^i in Sadourney (1972) and it's
         253             :     ! inverse
         254             :     ! ==============================================
         255             :
         256             :     ! MNL: Calculate Jacobians of bilinear map from cubed-sphere to ref element
         257       21600 :     if (cubed_sphere_map==0) then
         258       21600 :        call elem_jacobians(elem%cartp, elem%u2qmap)
         259             :     endif
         260             :
         261             :     max_svd = 0.0_r8
         262             :     max_normDinv = 0.0_r8
         263             :     min_svd = 1d99
         264      108000 :     do j=1,np
         265      453600 :        do i=1,np
         266      345600 :           x1=gll_points(i)
         267      345600 :           x2=gll_points(j)
         268     2419200 :           call Dmap(elem%D(i,j,:,:),x1,x2,elem%corners3D,cubed_sphere_map,elem%corners,elem%u2qmap,elem%facenum)
         269             :
         270             :
         271             :           ! Numerical metric tensor based on analytic D: met = D^T times D
         272             :           ! (D maps between sphere and reference element)
         273             :           elem%met(i,j,1,1) = elem%D(i,j,1,1)*elem%D(i,j,1,1) + &
         274      345600 :                               elem%D(i,j,2,1)*elem%D(i,j,2,1)
         275             :           elem%met(i,j,1,2) = elem%D(i,j,1,1)*elem%D(i,j,1,2) + &
         276      345600 :                               elem%D(i,j,2,1)*elem%D(i,j,2,2)
         277             :           elem%met(i,j,2,1) = elem%D(i,j,1,1)*elem%D(i,j,1,2) + &
         278      345600 :                               elem%D(i,j,2,1)*elem%D(i,j,2,2)
         279             :           elem%met(i,j,2,2) = elem%D(i,j,1,2)*elem%D(i,j,1,2) + &
         280      345600 :                               elem%D(i,j,2,2)*elem%D(i,j,2,2)
         281             :
         282             :           ! compute D^-1...
         283             :           ! compute determinant of D mapping matrix... if not zero compute inverse
         284             :
         285      345600 :           detD = elem%D(i,j,1,1)*elem%D(i,j,2,2) - elem%D(i,j,1,2)*elem%D(i,j,2,1)
         286             :
         287      345600 :           elem%Dinv(i,j,1,1) =  elem%D(i,j,2,2)/detD
         288      345600 :           elem%Dinv(i,j,1,2) = -elem%D(i,j,1,2)/detD
         289      345600 :           elem%Dinv(i,j,2,1) = -elem%D(i,j,2,1)/detD
         290      345600 :           elem%Dinv(i,j,2,2) =  elem%D(i,j,1,1)/detD
         291             :
         292             :           ! L2 norm = sqrt max eigenvalue of metinv
         293             :           !         = 1/sqrt(min eigenvalue of met)
         294             :           ! l1 and l2 are eigenvalues of met
         295             :           ! (should both be positive, l1 > l2)
         296             :           l1 = (elem%met(i,j,1,1) + elem%met(i,j,2,2) + sqrt(4.0_r8*elem%met(i,j,1,2)*elem%met(i,j,2,1) + &
         297      345600 :               (elem%met(i,j,1,1) - elem%met(i,j,2,2))**2))/2.0_r8
         298             :           l2 = (elem%met(i,j,1,1) + elem%met(i,j,2,2) - sqrt(4.0_r8*elem%met(i,j,1,2)*elem%met(i,j,2,1) + &
         299      345600 :               (elem%met(i,j,1,1) - elem%met(i,j,2,2))**2))/2.0_r8
         300             :           ! Max L2 norm of Dinv is sqrt of max eigenvalue of metinv
         301             :           ! max eigenvalue of metinv is 1/min eigenvalue of met
         302      345600 :           norm = 1.0_r8/sqrt(min(abs(l1),abs(l2)))
         303      345600 :           max_svd = max(norm, max_svd)
         304             :           ! Min L2 norm of Dinv is sqrt of min eigenvalue of metinv
         305             :           ! min eigenvalue of metinv is 1/max eigenvalue of met
         306      345600 :           norm = 1.0_r8/sqrt(max(abs(l1),abs(l2)))
         307      345600 :           min_svd = min(norm, min_svd)
         308             :
         309             :           ! some kind of pseudo-norm of Dinv
         310             :           ! C = 1/sqrt(2) sqrt( |g^x|^2 + |g^y|^2 + 2*|g^x dot g^y|)
         311             :           !   = 1/sqrt(2) sqrt( |g_x|^2 + |g_y|^2 + 2*|g_x dot g_y|) / J
         312             :           ! g^x = Dinv(:,1)    g_x = D(1,:)
         313             :           ! g^y = Dinv(:,2)    g_y = D(2,:)
         314     2419200 :           norm = (2*abs(sum(elem%Dinv(i,j,:,1)*elem%Dinv(i,j,:,2))) + sum(elem%Dinv(i,j,:,1)**2) + sum(elem%Dinv(i,j,:,2)**2))
         315      345600 :           norm = sqrt(norm)
         316             : !          norm = (2*abs(sum(elem%D(1,:,i,j)*elem%D(2,:,i,j))) + sum(elem%D(1,:,i,j)**2) + sum(elem%D(2,:,i,j)**2))
         317             : !          norm = sqrt(norm)/detD
         318             :           max_normDinv = max(norm,max_normDinv)
         319             :
         320             :
         321             :           ! Need inverse of met if not calculated analytically
         322      345600 :           elem%metdet(i,j) = abs(detD)
         323      345600 :           elem%rmetdet(i,j) = 1.0_R8/abs(detD)
         324             :
         325      345600 :           elem%metinv(i,j,1,1) =  elem%met(i,j,2,2)/(detD*detD)
         326      345600 :           elem%metinv(i,j,1,2) = -elem%met(i,j,1,2)/(detD*detD)
         327      345600 :           elem%metinv(i,j,2,1) = -elem%met(i,j,2,1)/(detD*detD)
         328      345600 :           elem%metinv(i,j,2,2) =  elem%met(i,j,1,1)/(detD*detD)
         329             :
         330             :           ! matricies for tensor hyper-viscosity
         331             :           ! compute eigenvectors of metinv (probably same as computed above)
         332     2419200 :           M = elem%metinv(i,j,:,:)
         333             :
         334             :           eig(1) = (M(1,1) + M(2,2) + sqrt(4.0_r8*M(1,2)*M(2,1) + &
         335      345600 :               (M(1,1) - M(2,2))**2))/2.0_r8
         336             :           eig(2) = (M(1,1) + M(2,2) - sqrt(4.0_r8*M(1,2)*M(2,1) + &
         337      345600 :               (M(1,1) - M(2,2))**2))/2.0_r8
         338             :
         339             :           ! use DE to store M - Lambda, to compute eigenvectors
         340      345600 :           DE=M
         341      345600 :           DE(1,1)=DE(1,1)-eig(1)
         342      345600 :           DE(2,2)=DE(2,2)-eig(1)
         343             :
         344     2419200 :           imaxM = maxloc(abs(DE))
         345     2419200 :           if (maxval(abs(DE))==0) then
         346          48 :              E(1,1)=1; E(2,1)=0;
         347      345552 :           elseif ( imaxM(1)==1 .and. imaxM(2)==1 ) then
         348      172784 :              E(2,1)=1; E(1,1) = -DE(2,1)/DE(1,1)
         349      172768 :           else   if ( imaxM(1)==1 .and. imaxM(2)==2 ) then
         350           0 :              E(2,1)=1; E(1,1) = -DE(2,2)/DE(1,2)
         351      172768 :           else   if ( imaxM(1)==2 .and. imaxM(2)==1 ) then
         352         912 :              E(1,1)=1; E(2,1) = -DE(1,1)/DE(2,1)
         353      171856 :           else   if ( imaxM(1)==2 .and. imaxM(2)==2 ) then
         354      171856 :              E(1,1)=1; E(2,1) = -DE(1,2)/DE(2,2)
         355             :           else
         356           0 :              call endrun('Impossible error in cube_mod.F90::metric_atomic()')
         357             :           endif
         358             :
         359             :           ! the other eigenvector is orthgonal:
         360      345600 :           E(1,2)=-E(2,1)
         361      345600 :           E(2,2)= E(1,1)
         362             :
         363             : !normalize columns
         364     1728000 :       E(:,1)=E(:,1)/sqrt(sum(E(:,1)*E(:,1)));
         365     1728000 :       E(:,2)=E(:,2)/sqrt(sum(E(:,2)*E(:,2)));
         366             :
         367             :
         368             : ! OBTAINING TENSOR FOR HV:
         369             :
         370             : ! Instead of the traditional scalar Laplace operator \grad \cdot \grad
         371             : ! we introduce \grad \cdot V \grad
         372             : ! where V = D E LAM LAM^* E^T D^T.
         373             : ! Recall (metric_tensor)^{-1}=(D^T D)^{-1} = E LAM E^T.
         374             : ! Here, LAM = diag( 4/((np-1)dx)^2 , 4/((np-1)dy)^2 ) = diag(  4/(dx_elem)^2, 4/(dy_elem)^2 )
         375             : ! Note that metric tensors and LAM correspondingly are quantities on a unit sphere.
         376             :
         377             : ! This motivates us to use V = D E LAM LAM^* E^T D^T
         378             : ! where LAM^* = diag( nu1, nu2 ) where nu1, nu2 are HV coefficients scaled like (dx)^{hv_scaling/2}, (dy)^{hv_scaling/2}.
         379             : ! (Halves in powers come from the fact that HV consists of two Laplace iterations.)
         380             :
         381             : ! Originally, we took LAM^* = diag(
         382             : !  1/(eig(1)**(hypervis_scaling/4.0_r8))*(rearth**(hypervis_scaling/2.0_r8))
         383             : !  1/(eig(2)**(hypervis_scaling/4.0_r8))*(rearth**(hypervis_scaling/2.0_r8)) ) =
         384             : !  = diag( lamStar1, lamStar2)
         385             : !  \simeq ((np-1)*dx_sphere / 2 )^hv_scaling/2 = SQRT(OPERATOR_HV)
         386             : ! because 1/eig(...) \simeq (dx_on_unit_sphere)^2 .
         387             : ! Introducing the notation OPERATOR = lamStar^2 is useful for conversion formulas.
         388             :
         389             : ! This leads to the following conversion formula: nu_const is nu used for traditional HV on uniform grids
         390             : ! nu_tensor = nu_const * OPERATOR_HV^{-1}, so
         391             : ! nu_tensor = nu_const *((np-1)*dx_sphere / 2 )^{ - hv_scaling} or
         392             : ! nu_tensor = nu_const *(2/( (np-1) * dx_sphere) )^{hv_scaling} .
         393             : ! dx_sphere = 2\pi *rearth/(np-1)/4/NE
         394             : ! [nu_tensor] = [meter]^{4-hp_scaling}/[sec]
         395             :
         396             : ! (1) Later developments:
         397             : ! Apply tensor V only at the second Laplace iteration. Thus, LAM^* should be scaled as (dx)^{hv_scaling}, (dy)^{hv_scaling},
         398             : ! see this code below:
         399             : !          DEL(1:2,1) = (lamStar1**2) *eig(1)*DE(1:2,1)
         400             : !          DEL(1:2,2) = (lamStar2**2) *eig(2)*DE(1:2,2)
         401             :
         402             : ! (2) Later developments:
         403             : ! Bringing [nu_tensor] to 1/[sec]:
         404             : !     lamStar1=1/(eig(1)**(hypervis_scaling/4.0_r8)) *(rearth**2.0_r8)
         405             : !     lamStar2=1/(eig(2)**(hypervis_scaling/4.0_r8)) *(rearth**2.0_r8)
         406             : ! OPERATOR_HV = ( (np-1)*dx_unif_sphere / 2 )^{hv_scaling} * rearth^4
         407             : ! Conversion formula:
         408             : ! nu_tensor = nu_const * OPERATOR_HV^{-1}, so
         409             : ! nu_tensor = nu_const *( 2*rearth /((np-1)*dx))^{hv_scaling} * rearth^{-4.0}.
         410             :
         411             : ! For the baseline coefficient nu=1e15 for NE30,
         412             : ! nu_tensor=7e-8 (BUT RUN TWICE AS SMALL VALUE FOR NOW) for hv_scaling=3.2
         413             : ! and
         414             : ! nu_tensor=1.3e-6 for hv_scaling=4.0.
         415             :
         416             :
         417             : !matrix D*E
         418     1036800 :           DE(1,1)=sum(elem%D(i,j,1,:)*E(:,1))
         419     1036800 :           DE(1,2)=sum(elem%D(i,j,1,:)*E(:,2))
         420     1036800 :           DE(2,1)=sum(elem%D(i,j,2,:)*E(:,1))
         421     1036800 :           DE(2,2)=sum(elem%D(i,j,2,:)*E(:,2))
         422             :
         423      345600 :       lamStar1=1/(eig(1)**(hypervis_scaling/4.0_r8)) *(rearth**2.0_r8)
         424      345600 :       lamStar2=1/(eig(2)**(hypervis_scaling/4.0_r8)) *(rearth**2.0_r8)
         425             :
         426             : !matrix (DE) * Lam^* * Lam , tensor HV when V is applied at each Laplace calculation
         427             : !          DEL(1:2,1) = lamStar1*eig(1)*DE(1:2,1)
         428             : !          DEL(1:2,2) = lamStar2*eig(2)*DE(1:2,2)
         429             :
         430             : !matrix (DE) * (Lam^*)^2 * Lam, tensor HV when V is applied only once, at the last Laplace calculation
         431             : !will only work with hyperviscosity, not viscosity
         432     1036800 :           DEL(1:2,1) = (lamStar1**2) *eig(1)*DE(1:2,1)
         433     1036800 :           DEL(1:2,2) = (lamStar2**2) *eig(2)*DE(1:2,2)
         434             :
         435             : !matrix (DE) * Lam^* * Lam  *E^t *D^t or (DE) * (Lam^*)^2 * Lam  *E^t *D^t
         436     1036800 :           V(1,1)=sum(DEL(1,:)*DE(1,:))
         437     1036800 :           V(1,2)=sum(DEL(1,:)*DE(2,:))
         438     1036800 :           V(2,1)=sum(DEL(2,:)*DE(1,:))
         439     1036800 :           V(2,2)=sum(DEL(2,:)*DE(2,:))
         440             :
         441     2505600 :       elem%tensorVisc(i,j,:,:)=V(:,:)
         442             :
         443             :        end do
         444             :     end do
         445             :
         446             : !    see Paul Ullrich writeup:
         447             : !    max_normDinv might be a tighter bound than max_svd for deformed elements
         448             : !    max_svd >= max_normDinv/sqrt(2), with equality holding if |g^x| = |g^y|
         449             : !    elem%normDinv=max_normDinv/sqrt(2)
         450             :
         451             :     ! this norm is consistent with length scales defined below:
         452       21600 :     elem%normDinv=max_svd
         453             :
         454             :
         455             :     ! compute element length scales, based on SVDs, in km:
         456       21600 :     elem%dx_short = 1.0_r8/(max_svd*0.5_r8*dble(np-1)*ra*1000.0_r8)
         457       21600 :     elem%dx_long  = 1.0_r8/(min_svd*0.5_r8*dble(np-1)*ra*1000.0_r8)
         458             :
         459             :     ! optional noramlization:
         460     1879200 :     elem%D = elem%D * sqrt(alpha)
         461     1879200 :     elem%Dinv = elem%Dinv / sqrt(alpha)
         462      453600 :     elem%metdet = elem%metdet * alpha
         463      453600 :     elem%rmetdet = elem%rmetdet / alpha
         464     1879200 :     elem%met = elem%met * alpha
         465     1879200 :     elem%metinv = elem%metinv / alpha
         466             :
         467       21600 :   end subroutine metric_atomic
         468             :
         469             :
         470             :   ! ========================================
         471             :   ! covariant_rot:
         472             :   !
         473             :   ! 2 x 2 matrix multiply:  Db^T * Da^-T
         474             :   ! for edge rotations: maps face a to face b
         475             :   !
         476             :   ! ========================================
         477             :
         478           0 :   function covariant_rot(Da,Db) result(R)
         479             :
         480             :     real (kind=r8) :: Da(2,2)
         481             :     real (kind=r8) :: Db(2,2)
         482             :     real (kind=r8) :: R(2,2)
         483             :
         484             :     real (kind=r8) :: detDa
         485             :
         486           0 :     detDa = Da(2,2)*Da(1,1) - Da(1,2)*Da(2,1)
         487             :
         488           0 :     R(1,1)=(Da(2,2)*Db(1,1) - Da(1,2)*Db(2,1))/detDa
         489           0 :     R(1,2)=(Da(1,1)*Db(2,1) - Da(2,1)*Db(1,1))/detDa
         490           0 :     R(2,1)=(Da(2,2)*Db(1,2) - Da(1,2)*Db(2,2))/detDa
         491           0 :     R(2,2)=(Da(1,1)*Db(2,2) - Da(2,1)*Db(1,2))/detDa
         492             :
         493           0 :   end function covariant_rot
         494             :
         495             :   ! ========================================
         496             :   ! contravariant_rot:
         497             :   !
         498             :   ! 2 x 2 matrix multiply:  Db^-1 * Da
         499             :   ! that maps a contravariant vector field
         500             :   ! from an edge of cube face a to a contiguous
         501             :   ! edge of cube face b.
         502             :   !
         503             :   ! ========================================
         504             :
         505        8544 :   function contravariant_rot(Da,Db) result(R)
         506             :
         507             :     real(kind=r8), intent(in) :: Da(2,2)
         508             :     real(kind=r8), intent(in) :: Db(2,2)
         509             :     real(kind=r8)             :: R(2,2)
         510             :
         511             :     real(kind=r8)             :: detDb
         512             :
         513        8544 :     detDb = Db(2,2)*Db(1,1) - Db(1,2)*Db(2,1)
         514             :
         515        8544 :     R(1,1)=(Da(1,1)*Db(2,2) - Da(2,1)*Db(1,2))/detDb
         516        8544 :     R(1,2)=(Da(1,2)*Db(2,2) - Da(2,2)*Db(1,2))/detDb
         517        8544 :     R(2,1)=(Da(2,1)*Db(1,1) - Da(1,1)*Db(2,1))/detDb
         518        8544 :     R(2,2)=(Da(2,2)*Db(1,1) - Da(1,2)*Db(2,1))/detDb
         519             :
         520        8544 :   end function contravariant_rot
         521             :
         522             :   ! ========================================================
         523             :   ! Dmap:
         524             :   !
         525             :   ! Initialize mapping that tranforms contravariant
         526             :   ! vector fields on the reference element onto vector fields on
         527             :   ! the sphere.
         528             :   ! ========================================================
         529    48168000 :   subroutine Dmap(D, a,b, corners3D, ref_map, corners, u2qmap, facenum)
         530             :     real (kind=r8), intent(out)  :: D(2,2)
         531             :     real (kind=r8), intent(in)     :: a,b
         532             :     type (cartesian3D_t)   :: corners3D(4)  !x,y,z coords of element corners
         533             :     integer :: ref_map
         534             :     ! only needed for ref_map=0,1
         535             :     type (cartesian2D_t),optional   :: corners(4)    ! gnomonic coords of element corners
         536             :     real (kind=r8),optional  :: u2qmap(4,2)
         537             :     integer,optional  :: facenum
         538             :
         539             :
         540             :
         541    48168000 :     if (ref_map==0) then
         542    48168000 :        if (.not. present ( corners ) ) &
         543           0 :             call endrun('Dmap(): missing arguments for equiangular map')
         544    48168000 :        call dmap_equiangular(D,a,b,corners,u2qmap,facenum)
         545           0 :     else if (ref_map==1) then
         546           0 :        call endrun('equi-distance gnomonic map not yet implemented')
         547           0 :     else if (ref_map==2) then
         548           0 :        call dmap_elementlocal(D,a,b,corners3D)
         549             :     else
         550           0 :        call endrun('bad value of ref_map')
         551             :     endif
         552    48168000 :   end subroutine Dmap
         553             :
         554             :
         555             :
         556             :   ! ========================================================
         557             :   ! Dmap:
         558             :   !
         559             :   ! Equiangular Gnomonic Projection
         560             :   ! Composition of equiangular Gnomonic projection to cubed-sphere face,
         561             :   ! followd by bilinear map to reference element
         562             :   ! ========================================================
         563    48168000 :   subroutine dmap_equiangular(D, a,b, corners,u2qmap,facenum )
         564             :     use dimensions_mod, only : np
         565             :     real (kind=r8), intent(out)  :: D(2,2)
         566             :     real (kind=r8), intent(in)     :: a,b
         567             :     real (kind=r8)    :: u2qmap(4,2)
         568             :     type (cartesian2D_t)     :: corners(4)                          ! gnomonic coords of element corners
         569             :     integer :: facenum
         570             :     ! local
         571             :     real (kind=r8)  :: tmpD(2,2), Jp(2,2),x1,x2,pi,pj,qi,qj
         572             :     real (kind=r8), dimension(4,2) :: unif2quadmap
         573             :
         574             : #if 0
         575             :     ! we shoud get rid of elem%u2qmap() and routine cube_mod.F90::elem_jacobian()
         576             :     ! and replace with this code below:
         577             :     ! but this produces roundoff level changes
         578             :     !unif2quadmap(1,1)=(elem%cartp(1,1)%x+elem%cartp(np,1)%x+elem%cartp(np,np)%x+elem%cartp(1,np)%x)/4.0_r8
         579             :     !unif2quadmap(1,2)=(elem%cartp(1,1)%y+elem%cartp(np,1)%y+elem%cartp(np,np)%y+elem%cartp(1,np)%y)/4.0_r8
         580             :     unif2quadmap(2,1)=(-elem%cartp(1,1)%x+elem%cartp(np,1)%x+elem%cartp(np,np)%x-elem%cartp(1,np)%x)/4.0_r8
         581             :     unif2quadmap(2,2)=(-elem%cartp(1,1)%y+elem%cartp(np,1)%y+elem%cartp(np,np)%y-elem%cartp(1,np)%y)/4.0_r8
         582             :     unif2quadmap(3,1)=(-elem%cartp(1,1)%x-elem%cartp(np,1)%x+elem%cartp(np,np)%x+elem%cartp(1,np)%x)/4.0_r8
         583             :     unif2quadmap(3,2)=(-elem%cartp(1,1)%y-elem%cartp(np,1)%y+elem%cartp(np,np)%y+elem%cartp(1,np)%y)/4.0_r8
         584             :     unif2quadmap(4,1)=(elem%cartp(1,1)%x-elem%cartp(np,1)%x+elem%cartp(np,np)%x-elem%cartp(1,np)%x)/4.0_r8
         585             :     unif2quadmap(4,2)=(elem%cartp(1,1)%y-elem%cartp(np,1)%y+elem%cartp(np,np)%y-elem%cartp(1,np)%y)/4.0_r8
         586             :     Jp(1,1) = unif2quadmap(2,1) + unif2quadmap(4,1)*b
         587             :     Jp(1,2) = unif2quadmap(3,1) + unif2quadmap(4,1)*a
         588             :     Jp(2,1) = unif2quadmap(2,2) + unif2quadmap(4,2)*b
         589             :     Jp(2,2) = unif2quadmap(3,2) + unif2quadmap(4,2)*a
         590             : #else
         591             :     ! input (a,b) shold be a point in the reference element [-1,1]
         592             :     ! compute Jp(a,b)
         593    48168000 :     Jp(1,1) = u2qmap(2,1) + u2qmap(4,1)*b
         594    48168000 :     Jp(1,2) = u2qmap(3,1) + u2qmap(4,1)*a
         595    48168000 :     Jp(2,1) = u2qmap(2,2) + u2qmap(4,2)*b
         596    48168000 :     Jp(2,2) = u2qmap(3,2) + u2qmap(4,2)*a
         597             : #endif
         598             :
         599             :     ! map (a,b) to the [-pi/2,pi/2] equi angular cube face:  x1,x2
         600             :     ! a = gp%points(i)
         601             :     ! b = gp%points(j)
         602    48168000 :     pi = (1-a)/2
         603    48168000 :     pj = (1-b)/2
         604    48168000 :     qi = (1+a)/2
         605    48168000 :     qj = (1+b)/2
         606             :     x1 = pi*pj*corners(1)%x &
         607             :          + qi*pj*corners(2)%x &
         608             :          + qi*qj*corners(3)%x &
         609    48168000 :          + pi*qj*corners(4)%x
         610             :     x2 = pi*pj*corners(1)%y &
         611             :          + qi*pj*corners(2)%y &
         612             :          + qi*qj*corners(3)%y &
         613    48168000 :          + pi*qj*corners(4)%y
         614             :
         615    48168000 :     call vmap(tmpD,x1,x2,facenum)
         616             :
         617             :     ! Include map from element -> ref element in D
         618    48168000 :     D(1,1) = tmpD(1,1)*Jp(1,1) + tmpD(1,2)*Jp(2,1)
         619    48168000 :     D(1,2) = tmpD(1,1)*Jp(1,2) + tmpD(1,2)*Jp(2,2)
         620    48168000 :     D(2,1) = tmpD(2,1)*Jp(1,1) + tmpD(2,2)*Jp(2,1)
         621    48168000 :     D(2,2) = tmpD(2,1)*Jp(1,2) + tmpD(2,2)*Jp(2,2)
         622    48168000 :   end subroutine dmap_equiangular
         623             :
         624             :
         625             :
         626             :   ! ========================================================
         627             :   ! vmap:
         628             :   !
         629             :   ! Initialize mapping that tranforms contravariant
         630             :   ! vector fields on the cube onto vector fields on
         631             :   ! the sphere. This follows Taylor's D matrix
         632             :   !
         633             :   !       | cos(theta)dlambda/dx1  cos(theta)dlambda/dx2 |
         634             :   !   D = |                                              |
         635             :   !       |     dtheta/dx1              dtheta/dx2       |
         636             :   !
         637             :   ! ========================================================
         638             :
         639    48185088 :   subroutine vmap(D, x1, x2, face_no)
         640             :     real(kind=r8), intent(inout) :: D(2,2)
         641             :     real(kind=r8), intent(in)    :: x1
         642             :     real(kind=r8), intent(in)    :: x2
         643             :     integer,       intent(in)    :: face_no
         644             :
         645             :     ! Local variables
         646             :
         647             :     real (kind=r8)    :: poledist ! SQRT(TAN(x1)**2 +TAN(x2)**2)
         648             :     real (kind=r8)    :: r        ! distance from cube point to center of sphere
         649             :
         650             :     real (kind=r8)    :: D11
         651             :     real (kind=r8)    :: D12
         652             :     real (kind=r8)    :: D21
         653             :     real (kind=r8)    :: D22
         654             :     character(len=64) :: errmsg
         655             :
         656    48185088 :     r = SQRT(1.0_r8 + TAN(x1)**2 + TAN(x2)**2)
         657             :
         658    48185088 :     if (face_no >= 1 .and. face_no <= 4) then
         659             :
         660    32123392 :        D11 = 1.0_r8 / (r * COS(x1))
         661    32123392 :        D12 = 0.0_r8
         662    32123392 :        D21 = -TAN(x1)*TAN(x2) / (COS(x1)*r*r)
         663    32123392 :        D22 = 1.0_r8 / (r*r*COS(x1)*COS(x2)*COS(x2))
         664             :
         665    32123392 :        D(1,1) =  D11
         666    32123392 :        D(1,2) =  D12
         667    32123392 :        D(2,1) =  D21
         668    32123392 :        D(2,2) =  D22
         669             :
         670             :
         671    16061696 :     else if (face_no == 6) then
         672     8030848 :        poledist = SQRT( TAN(x1)**2 + TAN(x2)**2)
         673     8030848 :        if (poledist <= DIST_THRESHOLD) then
         674             :
         675             :           ! we set the D transform to the identity matrix
         676             :           ! which works ONLY for swtc1, phi starting at
         677             :           ! 3*PI/2... assumes lon at pole == 0
         678             :
         679          16 :           D(1,1) =  1.0_r8
         680          16 :           D(1,2) =  0.0_r8
         681          16 :           D(2,1) =  0.0_r8
         682          16 :           D(2,2) =  1.0_r8
         683             :
         684             :        else
         685             :
         686     8030832 :           D11 = -TAN(x2)/(poledist*COS(x1)*COS(x1)*r)
         687     8030832 :           D12 =  TAN(x1)/(poledist*COS(x2)*COS(x2)*r)
         688     8030832 :           D21 = -TAN(x1)/(poledist*COS(x1)*COS(x1)*r*r)
         689     8030832 :           D22 = -TAN(x2)/(poledist*COS(x2)*COS(x2)*r*r)
         690             :
         691     8030832 :           D(1,1) =  D11
         692     8030832 :           D(1,2) =  D12
         693     8030832 :           D(2,1) =  D21
         694     8030832 :           D(2,2) =  D22
         695             :
         696             :        end if
         697     8030848 :     else if (face_no == 5) then
         698     8030848 :        poledist = SQRT( TAN(x1)**2 + TAN(x2)**2)
         699     8030848 :        if (poledist <= DIST_THRESHOLD) then
         700             :
         701             :           ! we set the D transform to the identity matrix
         702             :           ! which works ONLY for swtc1, phi starting at
         703             :           ! 3*PI/2... assumes lon at pole == 0, i.e. very specific
         704             :
         705          16 :           D(1,1) =  1.0_r8
         706          16 :           D(1,2) =  0.0_r8
         707          16 :           D(2,1) =  0.0_r8
         708          16 :           D(2,2) =  1.0_r8
         709             :
         710             :        else
         711             :
         712     8030832 :           D11 =  TAN(x2)/(poledist*COS(x1)*COS(x1)*r)
         713     8030832 :           D12 = -TAN(x1)/(poledist*COS(x2)*COS(x2)*r)
         714     8030832 :           D21 =  TAN(x1)/(poledist*COS(x1)*COS(x1)*r*r)
         715     8030832 :           D22 =  TAN(x2)/(poledist*COS(x2)*COS(x2)*r*r)
         716             :
         717     8030832 :           D(1,1) =  D11
         718     8030832 :           D(1,2) =  D12
         719     8030832 :           D(2,1) =  D21
         720     8030832 :           D(2,2) =  D22
         721             :
         722             :        end if
         723             :      else
         724           0 :        write(errmsg, '(a,i0)') 'VMAP: Bad face number, ',face_no
         725           0 :        call endrun(errmsg)
         726             :      end if
         727             :
         728    48185088 :   end subroutine vmap
         729             :
         730             :
         731             :
         732             :
         733             :   ! ========================================================
         734             :   ! Dmap:
         735             :   !
         736             :   ! Initialize mapping that tranforms contravariant
         737             :   ! vector fields on the reference element onto vector fields on
         738             :   ! the sphere.
         739             :   ! For Gnomonic, followed by bilinear, this code uses the old vmap()
         740             :   ! for unstructured grids, this code uses the parametric map that
         741             :   ! maps quads on the sphere directly to the reference element
         742             :   ! ========================================================
         743           0 :   subroutine dmap_elementlocal(D, a,b, corners3D)
         744             :     use element_mod, only : element_t
         745             :
         746             :     type (element_t) :: elem
         747             :     real (kind=r8), intent(out)    :: D(2,2)
         748             :     real (kind=r8), intent(in)     :: a,b
         749             :     type (cartesian3d_t)               ::  corners3D(4)
         750             :
         751             :     type (spherical_polar_t)      :: sphere
         752             :
         753             :     real(kind=r8)               ::  c(3,4), q(4), xx(3), r, lam, th, dd(4,2)
         754             :     real(kind=r8)               ::  sinlam, sinth, coslam, costh
         755             :     real(kind=r8)               ::  D1(2,3), D2(3,3), D3(3,2), D4(3,2)
         756             :     integer :: i,j
         757             :
         758           0 :     sphere = ref2sphere(a,b,corners3D,2) ! use element local map, ref_map=2
         759             :
         760           0 :     c(1,1)=corners3D(1)%x;  c(2,1)=corners3D(1)%y;  c(3,1)=corners3D(1)%z;
         761           0 :     c(1,2)=corners3D(2)%x;  c(2,2)=corners3D(2)%y;  c(3,2)=corners3D(2)%z;
         762           0 :     c(1,3)=corners3D(3)%x;  c(2,3)=corners3D(3)%y;  c(3,3)=corners3D(3)%z;
         763           0 :     c(1,4)=corners3D(4)%x;  c(2,4)=corners3D(4)%y;  c(3,4)=corners3D(4)%z;
         764             :
         765           0 :     q(1)=(1-a)*(1-b); q(2)=(1+a)*(1-b); q(3)=(1+a)*(1+b); q(4)=(1-a)*(1+b);
         766           0 :     q=q/4.0_r8;
         767             :
         768           0 :     do i=1,3
         769           0 :       xx(i)=sum(c(i,:)*q(:))
         770             :     enddo
         771             :
         772           0 :     r=sqrt(xx(1)**2+xx(2)**2+xx(3)**2)
         773             :
         774           0 :     lam=sphere%lon; th=sphere%lat;
         775           0 :     sinlam=sin(lam); sinth=sin(th);
         776           0 :     coslam=cos(lam); costh=cos(th);
         777             :
         778           0 :     D1(1,1)=-sinlam; D1(1,2)=coslam; D1(1,3)=0.0_r8;
         779           0 :     D1(2,1)=0.0_r8;  D1(2,2)=0.0_r8;    D1(2,3)=1.0_r8;
         780             :
         781           0 :     D2(1,1)=(sinlam**2)*(costh**2)+sinth**2; D2(1,2)=-sinlam*coslam*(costh**2); D2(1,3)=-coslam*sinth*costh;
         782           0 :     D2(2,1)=-sinlam*coslam*(costh**2); D2(2,2)=(coslam**2)*(costh**2)+sinth**2; D2(2,3)=-sinlam*sinth*costh;
         783           0 :     D2(3,1)=-coslam*sinth;           D2(3,2)=-sinlam*sinth;               D2(3,3)=costh;
         784             :
         785           0 :     dd(1,1)=-1+b; dd(1,2)=-1+a;
         786           0 :     dd(2,1)=1-b; dd(2,2)=-1-a;
         787           0 :     dd(3,1)=1+b; dd(3,2)=1+a;
         788           0 :     dd(4,1)=-1-b; dd(4,2)=1-a;
         789             :
         790           0 :     dd=dd/4.0_r8
         791             :
         792           0 :     do i=1,3
         793           0 :       do j=1,2
         794           0 :     D3(i,j)=sum(c(i,:)*dd(:,j))
         795             :       enddo
         796             :     enddo
         797             :
         798           0 :     do i=1,3
         799           0 :       do j=1,2
         800           0 :     D4(i,j)=sum(D2(i,:)*D3(:,j))
         801             :       enddo
         802             :     enddo
         803             :
         804           0 :     do i=1,2
         805           0 :       do j=1,2
         806           0 :     D(i,j)=sum(D1(i,:)*D4(:,j))
         807             :       enddo
         808             :     enddo
         809             :
         810           0 :     D=D/r
         811           0 :   end subroutine dmap_elementlocal
         812             :
         813             :
         814             :
         815             :
         816             :
         817             :   ! ========================================
         818             :   ! coreolis_init_atomic:
         819             :   !
         820             :   ! Initialize coreolis term ...
         821             :   !
         822             :   ! ========================================
         823             :
         824       21600 :   subroutine coreolis_init_atomic(elem)
         825             :     use element_mod,    only: element_t
         826             :     use dimensions_mod, only: np
         827             :     use physconst,      only: omega
         828             :
         829             :     type (element_t) :: elem
         830             :
         831             :     ! Local variables
         832             :
         833             :     integer                  :: i,j
         834             :     real (kind=r8) :: lat,lon,rangle
         835             :
         836       21600 :     rangle = rotate_grid * PI / 180._r8
         837      108000 :     do j=1,np
         838      453600 :        do i=1,np
         839      432000 :              if ( rotate_grid /= 0) then
         840           0 :                 lat = elem%spherep(i,j)%lat
         841           0 :                 lon = elem%spherep(i,j)%lon
         842             :                 elem%fcor(i,j)= 2*omega* &
         843           0 :                      (-cos(lon)*cos(lat)*sin(rangle) + sin(lat)*cos(rangle))
         844             :              else
         845      345600 :                 elem%fcor(i,j) = 2.0_r8*omega*SIN(elem%spherep(i,j)%lat)
         846             :              endif
         847             :        end do
         848             :     end do
         849             :
         850       21600 :   end subroutine coreolis_init_atomic
         851             :
         852             :   ! =========================================
         853             :   ! rotation_init_atomic:
         854             :   !
         855             :   ! Initialize cube rotation terms resulting
         856             :   ! from changing cube face coordinate systems
         857             :   !
         858             :   ! =========================================
         859             :
         860             :
         861       10800 :  subroutine rotation_init_atomic(elem, rot_type)
         862             :     use element_mod, only : element_t
         863             :     use dimensions_mod, only : np
         864             :     use control_mod, only : north, south, east, west, neast, seast, swest, nwest
         865             :
         866             :     type (element_t) :: elem
         867             :     character(len=*) rot_type
         868             :
         869             :     ! =======================================
         870             :     ! Local variables
         871             :     ! =======================================
         872             :
         873             :     integer :: myface_no        ! current element face number
         874             :     integer :: nbrface_no       ! neighbor element face number
         875             :     integer :: inbr
         876             :     integer :: nrot,irot
         877             :     integer :: ii,i,j,k
         878             :     integer :: ir,jr
         879             :     integer :: start, cnt
         880             :
         881             :     real (kind=r8) :: Dloc(2,2,np)
         882             :     real (kind=r8) :: Drem(2,2,np)
         883             :     real (kind=r8) :: x1,x2
         884             :
         885             :
         886       10800 :     myface_no = elem%vertex%face_number
         887             :
         888       10800 :     nrot   = 0
         889             :
         890       97200 :     do inbr=1,8
         891       86400 :         cnt = elem%vertex%nbrs_ptr(inbr+1) -  elem%vertex%nbrs_ptr(inbr)
         892             :         start =  elem%vertex%nbrs_ptr(inbr)
         893             :
         894      183552 :         do k = 0, cnt-1
         895       86352 :            nbrface_no = elem%vertex%nbrs_face(start+k)
         896      172752 :            if (myface_no /= nbrface_no) nrot=nrot+1
         897             :         end do
         898             :
         899             :     end do
         900             :
         901       10800 :     if(associated(elem%desc%rot)) then
         902           0 :        if (size(elem%desc%rot) > 0) then
         903             :           !         deallocate(elem%desc%rot)
         904           0 :           NULLIFY(elem%desc%rot)
         905             :        endif
         906             :     end if
         907             :
         908             :     ! =====================================================
         909             :     ! If there are neighbors on other cube faces, allocate
         910             :     ! an array of rotation matrix structs.
         911             :     ! =====================================================
         912             :
         913       10800 :     if (nrot > 0) then
         914        8400 :        allocate(elem%desc%rot(nrot))
         915        1392 :        elem%desc%use_rotation=1
         916        1392 :        irot=0
         917             :
         918       12528 :        do inbr=1,8
         919       11136 :           cnt = elem%vertex%nbrs_ptr(inbr+1) -  elem%vertex%nbrs_ptr(inbr)
         920             :           start =  elem%vertex%nbrs_ptr(inbr)
         921             :
         922       23616 :           do k= 0, cnt-1
         923             :
         924       11088 :              nbrface_no = elem%vertex%nbrs_face(start+k)
         925             :              ! The cube edge (myface_no,nbrface_no) and inbr defines
         926             :              ! a unique rotation given by (D^-1) on myface_no x (D) on nbrface_no
         927             :
         928       22224 :              if (myface_no /= nbrface_no .and. elem%vertex%nbrs(start+k) /= -1 ) then
         929        4224 :                 irot=irot+1
         930             :
         931        4224 :                 if (inbr <= 4) then
         932        1440 :                    allocate(elem%desc%rot(irot)%R(2,2,np))  ! edge
         933             :                 else
         934        2784 :                    allocate(elem%desc%rot(irot)%R(2,2,1 ))   ! corner
         935             :                 end if
         936             :                 ! Initialize Dloc and Drem for no-rotation possibilities
         937       21120 :                 Dloc(1,1,:) = 1.0_r8
         938       21120 :                 Dloc(1,2,:) = 0.0_r8
         939       21120 :                 Dloc(2,1,:) = 0.0_r8
         940       21120 :                 Dloc(2,2,:) = 1.0_r8
         941       21120 :                 Drem(1,1,:) = 1.0_r8
         942       21120 :                 Drem(1,2,:) = 0.0_r8
         943       21120 :                 Drem(2,1,:) = 0.0_r8
         944       21120 :                 Drem(2,2,:) = 1.0_r8
         945             :
         946             :                 ! must compute Dloc on my face, Drem on neighbor face,
         947             :                 ! for each point on edge or corner.
         948             :
         949             :                 ! ====================================
         950             :                 ! Equatorial belt east/west neighbors
         951             :                 ! ====================================
         952             :
         953        4224 :                 if (nbrface_no <= 4 .and. myface_no <= 4) then
         954             :
         955             :                    if (inbr == west) then
         956        1200 :                       do j=1,np
         957         960 :                          x1 = elem%cartp(1,j)%x
         958         960 :                          x2 = elem%cartp(1,j)%y
         959             :                          call Vmap(Dloc(1,1,j), x1,x2,myface_no)
         960        1200 :                          call Vmap(Drem(1,1,j),-x1,x2,nbrface_no)
         961             :                       end do
         962             :                    else if (inbr == east) then
         963        1200 :                       do j=1,np
         964         960 :                          x1 = elem%cartp(np,j)%x
         965         960 :                          x2 = elem%cartp(np,j)%y
         966         960 :                          call Vmap(Dloc(1,1,j), x1,x2,myface_no)
         967        1200 :                          call Vmap(Drem(1,1,j),-x1,x2,nbrface_no)
         968             :                       end do
         969             :                    else if (inbr == swest ) then
         970         232 :                       x1 = elem%cartp(1,1)%x
         971         232 :                       x2 = elem%cartp(1,1)%y
         972         232 :                       call Vmap(Dloc(1,1,1),x1,x2,myface_no)
         973         232 :                       call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
         974             :                    else if (inbr == nwest ) then
         975         232 :                       x1 = elem%cartp(1,np)%x
         976         232 :                       x2 = elem%cartp(1,np)%y
         977         232 :                       call Vmap(Dloc(1,1,1), x1,x2,myface_no)
         978         232 :                       call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
         979             :                    else if (inbr == seast ) then
         980         232 :                       x1 = elem%cartp(np,1)%x
         981         232 :                       x2 = elem%cartp(np,1)%y
         982         232 :                       call Vmap(Dloc(1,1,1), x1,x2,myface_no)
         983         232 :                       call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
         984             :                    else if (inbr == neast ) then
         985         232 :                       x1 = elem%cartp(np,np)%x
         986         232 :                       x2 = elem%cartp(np,np)%y
         987         232 :                       call Vmap(Dloc(1,1,1), x1,x2,myface_no)
         988         232 :                       call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
         989             :                    end if
         990             :
         991             :                 end if
         992             :
         993             :                 ! Northern Neighbors of Equatorial Belt
         994             :
         995        4224 :                 if ( myface_no <= 4 .and. nbrface_no == 6 ) then
         996         704 :                    if (inbr == north) then
         997        1200 :                       do i=1,np
         998         960 :                          ir=np+1-i
         999         960 :                          x1 = elem%cartp(i,np)%x
        1000         960 :                          x2 = elem%cartp(i,np)%y
        1001         960 :                          if ( myface_no == 1) then
        1002         240 :                             call Vmap(Dloc(1,1,i), x1,x2,myface_no)
        1003         240 :                             call Vmap(Drem(1,1,i),x1,-x2,nbrface_no)
        1004             :                          end if
        1005        1200 :                          if ( myface_no == 2) then
        1006         240 :                             call Vmap(Dloc(1,1,i),x1,x2,myface_no)
        1007             :                             call Vmap(Drem(1,1,i),x2,x1,nbrface_no)
        1008             :
        1009             :                          end if
        1010         960 :                          if ( myface_no == 3) then
        1011         240 :                             call Vmap(Dloc(1,1,ir), x1,x2,myface_no)
        1012         240 :                             call Vmap(Drem(1,1,ir),-x1,x2,nbrface_no)
        1013             :                          end if
        1014        1200 :                          if ( myface_no == 4) then
        1015         240 :                             call Vmap(Dloc(1,1,ir), x1,x2,myface_no)
        1016         240 :                             call Vmap(Drem(1,1,ir),-x2,-x1,nbrface_no)
        1017             :                          end if
        1018             :                       end do
        1019         464 :                    else if (inbr == nwest) then
        1020         232 :                       x1 = elem%cartp(1,np)%x
        1021         232 :                       x2 = elem%cartp(1,np)%y
        1022         232 :                       call Vmap(Dloc(1,1,1), x1,x2,myface_no)
        1023         232 :                       if ( myface_no == 1) call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
        1024         232 :                       if ( myface_no == 2) call Vmap(Drem(1,1,1),x2, x1,nbrface_no)
        1025         232 :                       if ( myface_no == 3) call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
        1026         232 :                       if ( myface_no == 4) call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
        1027         232 :                    else if (inbr == neast) then
        1028         232 :                       x1 = elem%cartp(np,np)%x
        1029         232 :                       x2 = elem%cartp(np,np)%y
        1030         232 :                       call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1031         232 :                       if ( myface_no == 1) call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
        1032         232 :                       if ( myface_no == 2) call Vmap(Drem(1,1,1),x2, x1,nbrface_no)
        1033         232 :                       if ( myface_no == 3) call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
        1034         232 :                       if ( myface_no == 4) call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
        1035             :                    end if
        1036             :
        1037             :                 end if
        1038             :
        1039             :                 ! Southern Neighbors of Equatorial Belt
        1040             :
        1041        4224 :                 if ( myface_no <= 4 .and. nbrface_no == 5 ) then
        1042         704 :                    if (inbr == south) then
        1043        1200 :                       do i=1,np
        1044         960 :                          ir=np+1-i
        1045         960 :                          x1 = elem%cartp(i,1)%x
        1046         960 :                          x2 = elem%cartp(i,1)%y
        1047         960 :                          if ( myface_no == 1) then
        1048         240 :                             call Vmap(Dloc(1,1,i), x1, x2,myface_no)
        1049         240 :                             call Vmap(Drem(1,1,i), x1,-x2,nbrface_no)
        1050             :                          end if
        1051         960 :                          if ( myface_no == 2) then
        1052         240 :                             call Vmap(Dloc(1,1,ir),x1,x2,myface_no)
        1053         240 :                             call Vmap(Drem(1,1,ir),-x2,-x1,nbrface_no)
        1054             :                          end if
        1055         960 :                          if ( myface_no == 3) then
        1056         240 :                             call Vmap(Dloc(1,1,ir), x1,x2,myface_no)
        1057         240 :                             call Vmap(Drem(1,1,ir),-x1,x2,nbrface_no)
        1058             :                          end if
        1059        1440 :                          if ( myface_no == 4) then
        1060         240 :                             call Vmap(Dloc(1,1,i), x1,x2,myface_no)
        1061             :                             call Vmap(Drem(1,1,i), x2,x1,nbrface_no)
        1062             :                          end if
        1063             :                       end do
        1064         464 :                    else if (inbr == swest) then
        1065         232 :                       x1 = elem%cartp(1,1)%x
        1066         232 :                       x2 = elem%cartp(1,1)%y
        1067         232 :                       call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1068             :
        1069             :
        1070         232 :                       if ( myface_no == 1) call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
        1071         232 :                       if ( myface_no == 2) call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
        1072         232 :                       if ( myface_no == 3) call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
        1073         232 :                       if ( myface_no == 4) call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
        1074             :
        1075         232 :                    else if (inbr == seast) then
        1076         232 :                       x1 = elem%cartp(np,1)%x
        1077         232 :                       x2 = elem%cartp(np,1)%y
        1078         232 :                       call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1079         232 :                       if ( myface_no == 1) call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
        1080         232 :                       if ( myface_no == 2) call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
        1081         232 :                       if ( myface_no == 3) call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
        1082         232 :                       if ( myface_no == 4) call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
        1083             :                    end if
        1084             :
        1085             :                 end if
        1086             :
        1087             :                 ! Neighbors of Northern Capping Face Number 6
        1088             :
        1089        4224 :                 if ( myface_no == 6 ) then
        1090         704 :                    if (nbrface_no == 1) then
        1091         176 :                       if (inbr == south) then
        1092         300 :                          do i=1,np
        1093         240 :                             x1 = elem%cartp(i,1)%x
        1094         240 :                             x2 = elem%cartp(i,1)%y
        1095         240 :                             call Vmap(Dloc(1,1,i),x1,x2,myface_no)
        1096         300 :                             call Vmap(Drem(1,1,i),x1,-x2,nbrface_no)
        1097             :                          end do
        1098         116 :                       else if (inbr == swest) then
        1099          58 :                          x1 = elem%cartp(1,1)%x
        1100          58 :                          x2 = elem%cartp(1,1)%y
        1101          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1102          58 :                          call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
        1103          58 :                       else if (inbr == seast) then
        1104          58 :                          x1 = elem%cartp(np,1)%x
        1105          58 :                          x2 = elem%cartp(np,1)%y
        1106          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1107          58 :                          call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
        1108             :                       end if
        1109         528 :                    else if (nbrface_no == 2) then
        1110         176 :                       if (inbr == east) then
        1111         300 :                          do j=1,np
        1112         240 :                             x1 = elem%cartp(np,j)%x
        1113         240 :                             x2 = elem%cartp(np,j)%y
        1114         240 :                             call Vmap(Dloc(1,1,j),x1,x2,myface_no)
        1115         300 :                             call Vmap(Drem(1,1,j),x2,x1,nbrface_no)
        1116             :                          end do
        1117         116 :                       else if (inbr == seast) then
        1118          58 :                          x1 = elem%cartp(np,1)%x
        1119          58 :                          x2 = elem%cartp(np,1)%y
        1120          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1121          58 :                          call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
        1122          58 :                       else if (inbr == neast) then
        1123          58 :                          x1 = elem%cartp(np,np)%x
        1124          58 :                          x2 = elem%cartp(np,np)%y
        1125          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1126          58 :                          call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
        1127             :                       end if
        1128         352 :                    else if (nbrface_no == 3) then
        1129         176 :                       if (inbr == north) then
        1130         300 :                          do i=1,np
        1131         240 :                             ir =np+1-i
        1132         240 :                             x1 = elem%cartp(i,np)%x
        1133         240 :                             x2 = elem%cartp(i,np)%y
        1134         240 :                             call Vmap(Dloc(1,1,ir),x1,x2,myface_no)
        1135         300 :                             call Vmap(Drem(1,1,ir),-x1,x2,nbrface_no)
        1136             :                          end do
        1137         116 :                       else if (inbr == nwest) then
        1138          58 :                          x1 = elem%cartp(1,np)%x
        1139          58 :                          x2 = elem%cartp(1,np)%y
        1140          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1141          58 :                          call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
        1142          58 :                       else if (inbr == neast) then
        1143          58 :                          x1 = elem%cartp(np,np)%x
        1144          58 :                          x2 = elem%cartp(np,np)%y
        1145          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1146          58 :                          call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
        1147             :                       end if
        1148         176 :                    else if (nbrface_no == 4) then
        1149         176 :                       if (inbr == west) then
        1150         300 :                          do j=1,np
        1151         240 :                             jr=np+1-j
        1152         240 :                             x1 = elem%cartp(1,j)%x
        1153         240 :                             x2 = elem%cartp(1,j)%y
        1154         240 :                             call Vmap(Dloc(1,1,jr), x1, x2,myface_no )
        1155         300 :                             call Vmap(Drem(1,1,jr),-x2,-x1,nbrface_no)
        1156             :                          end do
        1157         116 :                       else if (inbr == swest) then
        1158          58 :                          x1 = elem%cartp(1,1)%x
        1159          58 :                          x2 = elem%cartp(1,1)%y
        1160          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1161          58 :                          call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
        1162          58 :                       else if (inbr == nwest) then
        1163          58 :                          x1 = elem%cartp(1,np)%x
        1164          58 :                          x2 = elem%cartp(1,np)%y
        1165          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1166          58 :                          call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
        1167             :                       end if
        1168             :                    end if
        1169             :                 end if
        1170             :
        1171             :                 ! Neighbors of South Capping Face Number 5
        1172             :
        1173        4224 :                 if ( myface_no == 5 ) then
        1174         704 :                    if (nbrface_no == 1) then
        1175         176 :                       if (inbr == north) then
        1176         300 :                          do i=1,np
        1177         240 :                             x1 = elem%cartp(i,np)%x
        1178         240 :                             x2 = elem%cartp(i,np)%y
        1179         240 :                             call Vmap(Dloc(1,1,i),x1,x2,myface_no)
        1180         300 :                             call Vmap(Drem(1,1,i),x1,-x2,nbrface_no)
        1181             :                          end do
        1182         116 :                       else if (inbr == nwest) then
        1183          58 :                          x1 = elem%cartp(1,np)%x
        1184          58 :                          x2 = elem%cartp(1,np)%y
        1185          58 :                          call Vmap(Dloc(:,:,1),x1,x2,myface_no)
        1186          58 :                          call Vmap(Drem(:,:,1),x1,-x2,nbrface_no)
        1187          58 :                       else if (inbr == neast) then
        1188          58 :                          x1 = elem%cartp(np,np)%x
        1189          58 :                          x2 = elem%cartp(np,np)%y
        1190          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1191          58 :                          call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
        1192             :                       end if
        1193         528 :                    else if (nbrface_no == 2) then
        1194         176 :                       if (inbr == east) then
        1195         300 :                          do j=1,np
        1196         240 :                             jr=np+1-j
        1197         240 :                             x1 = elem%cartp(np,j)%x
        1198         240 :                             x2 = elem%cartp(np,j)%y
        1199         240 :                             call Vmap(Dloc(1,1,jr),x1,  x2,myface_no)
        1200         300 :                             call Vmap(Drem(1,1,jr),-x2,-x1,nbrface_no)
        1201             :                          end do
        1202         116 :                       else if (inbr == seast) then
        1203          58 :                          x1 = elem%cartp(np,1)%x
        1204          58 :                          x2 = elem%cartp(np,1)%y
        1205          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1206          58 :                          call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
        1207          58 :                       else if (inbr == neast) then
        1208          58 :                          x1 = elem%cartp(np,np)%x
        1209          58 :                          x2 = elem%cartp(np,np)%y
        1210          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1211          58 :                          call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
        1212             :                       end if
        1213         352 :                    else if (nbrface_no == 3) then
        1214         176 :                       if (inbr == south) then
        1215         300 :                          do i=1,np
        1216         240 :                             ir=np+1-i
        1217         240 :                             x1 = elem%cartp(i,1)%x
        1218         240 :                             x2 = elem%cartp(i,1)%y
        1219         240 :                             call Vmap(Dloc(1,1,ir),x1,x2,myface_no)
        1220         300 :                             call Vmap(Drem(1,1,ir),-x1,x2,nbrface_no)
        1221             :                          end do
        1222         116 :                       else if (inbr == swest) then
        1223          58 :                          x1 = elem%cartp(1,1)%x
        1224          58 :                          x2 = elem%cartp(1,1)%y
        1225          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1226          58 :                          call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
        1227          58 :                       else if (inbr == seast) then
        1228          58 :                          x1 = elem%cartp(np,1)%x
        1229          58 :                          x2 = elem%cartp(np,1)%y
        1230          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1231          58 :                          call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
        1232             :                       end if
        1233         176 :                    else if (nbrface_no == 4) then
        1234         176 :                       if (inbr == west) then
        1235         300 :                          do j=1,np
        1236         240 :                             x1 = elem%cartp(1,j)%x
        1237         240 :                             x2 = elem%cartp(1,j)%y
        1238             :                             call Vmap(Dloc(1,1,j),x1,x2,myface_no)
        1239         300 :                             call Vmap(Drem(1,1,j),x2,x1,nbrface_no)
        1240             :                          end do
        1241         116 :                       else if (inbr == swest) then
        1242          58 :                          x1 = elem%cartp(1,1)%x
        1243          58 :                          x2 = elem%cartp(1,1)%y
        1244          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1245          58 :                          call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
        1246          58 :                       else if (inbr == nwest) then
        1247          58 :                          x1 = elem%cartp(1,np)%x
        1248          58 :                          x2 = elem%cartp(1,np)%y
        1249          58 :                          call Vmap(Dloc(1,1,1),x1,x2,myface_no)
        1250          58 :                          call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
        1251             :                       end if
        1252             :                    end if
        1253             :                 end if
        1254             :
        1255        4224 :                 elem%desc%rot(irot)%nbr = inbr
        1256        4224 :                 if (rot_type == "covariant") then
        1257           0 :                    do i=1,SIZE(elem%desc%rot(irot)%R(:,:,:),3)
        1258           0 :                       elem%desc%rot(irot)%R(:,:,i)=covariant_rot(Dloc(:,:,i),Drem(:,:,i))
        1259             :                    end do
        1260        4224 :                 else if (rot_type == "contravariant") then
        1261       12768 :                    do i=1,SIZE(elem%desc%rot(irot)%R(:,:,:),3)
        1262       64032 :                       elem%desc%rot(irot)%R(:,:,i)=contravariant_rot(Dloc(:,:,i),Drem(:,:,i))
        1263             :                    end do
        1264             :                 end if
        1265             :
        1266             :              end if ! end of a unique rotation
        1267             :           end do !k loop over neighbors in that direction
        1268             :        end do !inbr loop
        1269             :     end if !nrot > 0
        1270             :
        1271       10800 :   end subroutine rotation_init_atomic
        1272             :
        1273             :
        1274       21600 :   subroutine set_corner_coordinates(elem)
        1275             :     use element_mod,    only : element_t
        1276             :     use dimensions_mod, only : ne
        1277             :
        1278             :     type (element_t) :: elem
        1279             :
        1280             :     ! Local variables
        1281             :     integer  i,ie,je,face_no,nn
        1282             :     real (kind=r8)  :: dx,dy, startx, starty
        1283             :
        1284       10800 :     if (0==ne) call endrun('Error in set_corner_coordinates: ne is zero')
        1285             :
        1286             :     ! ========================================
        1287             :     ! compute cube face coordinates of element
        1288             :     ! =========================================
        1289             :
        1290       10800 :     call convert_gbl_index(elem%vertex%number,ie,je,face_no)
        1291             :
        1292       10800 :     elem%vertex%face_number = face_no
        1293       10800 :     dx = (cube_xend-cube_xstart)/ne
        1294       10800 :     dy = (cube_yend-cube_ystart)/ne
        1295             :
        1296       10800 :     startx = cube_xstart+ie*dx
        1297       10800 :     starty = cube_ystart+je*dy
        1298             :
        1299       10800 :     elem%corners(1)%x = startx
        1300       10800 :     elem%corners(1)%y = starty
        1301       10800 :     elem%corners(2)%x = startx+dx
        1302       10800 :     elem%corners(2)%y = starty
        1303       10800 :     elem%corners(3)%x = startx+dx
        1304       10800 :     elem%corners(3)%y = starty+dy
        1305       10800 :     elem%corners(4)%x = startx
        1306       10800 :     elem%corners(4)%y = starty+dy
        1307             :
        1308             : #if 0
        1309             :     do i=1,4
        1310             :        elem%node_multiplicity(i) = 4
        1311             :     end do
        1312             :     ie = ie + 1
        1313             :     je = je + 1
        1314             :     if      (ie ==  1 .and. je ==  1) then
        1315             :        elem%node_multiplicity(1) = 3
        1316             :     else if (ie == ne .and. je ==  1) then
        1317             :        elem%node_multiplicity(2) = 3
        1318             :     else if (ie == ne .and. je == ne) then
        1319             :        elem%node_multiplicity(3) = 3
        1320             :     else if (ie ==  1 .and. je == ne) then
        1321             :        elem%node_multiplicity(4) = 3
        1322             :     end if
        1323             : #endif
        1324       10800 :   end subroutine set_corner_coordinates
        1325             :
        1326             :
        1327        1536 :   subroutine assign_node_numbers_to_elem(elements, GridVertex)
        1328             :     use dimensions_mod, only : ne
        1329             :     use element_mod,    only : element_t
        1330             :     use control_mod,    only : north, south, east, west, neast, seast, swest, nwest
        1331             :     use gridgraph_mod,  only : GridVertex_t
        1332             :     implicit none
        1333             :     type (element_t), intent(inout)    :: elements(:)
        1334             :     type (GridVertex_t), intent(in)    :: GridVertex(:)
        1335             :
        1336             :     type (GridVertex_t)                :: vertex
        1337        3072 :     integer                            :: connectivity(6*ne*ne, 4)
        1338             :     integer                            :: nn(4), en(4)
        1339             :     integer el, i, n, direction
        1340             :     integer current_node_num, tot_ne
        1341             :     integer :: start, cnt
        1342             :
        1343        1536 :     current_node_num = 0
        1344        1536 :     tot_ne = 6*ne*ne
        1345             :
        1346           0 :     if (0==ne)      call endrun('Error in assign_node_numbers_to_elem: ne is zero')
        1347        1536 :     if (tot_ne /= SIZE(GridVertex)) call endrun('Error in assign_node_numbers_to_elem: GridVertex not correct length')
        1348             :
        1349    33185280 :     connectivity = 0
        1350             :
        1351     8295936 :     do el = 1,tot_ne
        1352     8294400 :        vertex = GridVertex(el)
        1353     8294400 :        en = 0
        1354    74649600 :        do direction = 1,8
        1355    66355200 :           cnt = vertex%nbrs_ptr(direction+1) -  vertex%nbrs_ptr(direction)
        1356             :           start =  vertex%nbrs_ptr(direction)
        1357             :
        1358   140967936 :           do i=0, cnt-1
        1359    66318336 :              n = vertex%nbrs(start+i)
        1360   132673536 :              if (n /= -1) then
        1361   331591680 :                 nn = connectivity(n,:)
        1362     8294400 :                 select case (direction)
        1363             :                 case (north)
        1364     8294400 :                    if (nn(1)/=0) en(4) = nn(1)
        1365     8294400 :                    if (nn(2)/=0) en(3) = nn(2)
        1366             :                 case (south)
        1367     8294400 :                    if (nn(4)/=0) en(1) = nn(4)
        1368     8294400 :                    if (nn(3)/=0) en(2) = nn(3)
        1369             :                 case (east)
        1370     8294400 :                    if (nn(1)/=0) en(2) = nn(1)
        1371     8294400 :                    if (nn(4)/=0) en(3) = nn(4)
        1372             :                 case (west)
        1373     8294400 :                    if (nn(2)/=0) en(1) = nn(2)
        1374     8294400 :                    if (nn(3)/=0) en(4) = nn(3)
        1375             :                 case (neast)
        1376     8285184 :                    if (nn(1)/=0) en(3) = nn(1)
        1377             :                 case (seast)
        1378     8285184 :                    if (nn(4)/=0) en(2) = nn(4)
        1379             :                 case (swest)
        1380     8285184 :                    if (nn(3)/=0) en(1) = nn(3)
        1381             :                 case (nwest)
        1382    66318336 :                    if (nn(2)/=0) en(4) = nn(2)
        1383             :                 end select
        1384             :              end if
        1385             :           end do
        1386             :        end do !direction
        1387             :
        1388    41472000 :        do i=1,4
        1389    41472000 :           if (en(i) == 0) then
        1390     8297472 :              current_node_num = current_node_num + 1
        1391     8297472 :              en(i) = current_node_num
        1392             :           end if
        1393             :        end do
        1394    41473536 :        connectivity(el,:) = en
        1395             :     end do
        1396             :
        1397        1536 :     if (current_node_num /= (6*ne*ne+2)) then
        1398           0 :        call endrun('Error in assignment of node numbers: Failed Euler test')
        1399             :     end if
        1400             : !    do el = 1,SIZE(elements)
        1401             : !      elements(el)%node_numbers = connectivity(elements(el)%vertex%number, :)
        1402             : !    end do
        1403        1536 :   end subroutine assign_node_numbers_to_elem
        1404             :
        1405             :
        1406             :   ! ================================================
        1407             :   ! convert_gbl_index:
        1408             :   !
        1409             :   ! Convert global element index to cube index
        1410             :   ! ================================================
        1411             :
        1412       10800 :   subroutine convert_gbl_index(number,ie,je,face_no)
        1413             :     use dimensions_mod, only : ne
        1414             :     integer, intent(in)  :: number
        1415             :     integer, intent(out) :: ie,je,face_no
        1416             :
        1417       10800 :     if (0==ne) call endrun('Error in cube_mod:convert_gbl_index: ne is zero')
        1418             :
        1419             :     !  inverse of the function:      number = 1 + ie + ne*je + ne*ne*(face_no-1)
        1420       10800 :     face_no=((number-1)/(ne*ne))+1
        1421       10800 :     ie=MODULO(number-1,ne)
        1422       10800 :     je=(number-1)/ne - (face_no-1)*ne
        1423             :
        1424       10800 :   end subroutine convert_gbl_index
        1425             :
        1426        1536 :   subroutine CubeTopology(GridEdge, GridVertex)
        1427             :     use gridgraph_mod,  only: GridEdge_t, GridVertex_t, initgridedge
        1428             :     use gridgraph_mod,  only: allocate_gridvertex_nbrs, deallocate_gridvertex_nbrs
        1429             :     use dimensions_mod, only: np, ne
        1430             :     use spacecurve_mod, only: IsFactorable, genspacecurve
        1431             :     use control_mod,    only: north, south, east, west, neast, seast, swest, nwest
        1432             :     !-----------------------
        1433             :
        1434             :     ! Since GridVertex fields must be allocated before calling this, it
        1435             :     ! must be intent(inout).
        1436             : !og: is 'target' here necessary?
        1437             : !GridEdge : changed its 'out' attribute to 'inout'
        1438             :     type (GridEdge_t),   intent(inout),target     :: GridEdge(:)
        1439             :     type (GridVertex_t), intent(inout),target     :: GridVertex(:)
        1440             :
        1441             :
        1442        1536 :     integer,allocatable       :: Mesh(:,:)
        1443        1536 :     integer,allocatable       :: Mesh2(:,:),Mesh2_map(:,:,:),sfcij(:,:)
        1444        1536 :     type (GridVertex_t),allocatable        :: GridElem(:,:,:)
        1445             :     integer                   :: i,j,k,ll,number,irev,ne2,i2,j2,sfc_index
        1446             :     integer                   :: EdgeWgtP,CornerWgt
        1447             :     integer                   :: ielem, nedge
        1448             :     integer                   :: offset, ierr, loc
        1449        1536 :     logical, allocatable      :: nbrs_used(:,:,:,:)
        1450             :
        1451             :
        1452        1536 :     if (0==ne) call endrun('Error in CubeTopology: ne is zero')
        1453             :
        1454     8587776 :     allocate(GridElem(ne,ne,nfaces),stat=ierr)
        1455       10752 :     do k = 1, nfaces
        1456      287232 :        do j = 1, ne
        1457     8580096 :           do i = 1, ne
        1458     8570880 :              call allocate_gridvertex_nbrs(GridElem(i,j,k))
        1459             :           end do
        1460             :        end do
        1461             :     end do
        1462             :
        1463        1536 :     if(ierr/=0) then
        1464           0 :        call endrun('error in allocation of GridElem structure')
        1465             :     end if
        1466             :
        1467        9216 :     allocate(nbrs_used(ne,ne,nfaces,8))
        1468    68654592 :     nbrs_used = .false.
        1469             :
        1470             :
        1471             :     number=1
        1472             :     EdgeWgtP   = np
        1473             :     CornerWgt = 1
        1474       10752 :     do k=1,nfaces
        1475      287232 :        do j=1,ne
        1476     8580096 :           do i=1,ne
        1477             :              ! ====================================
        1478             :              ! Number elements
        1479             :              ! ====================================
        1480    74649600 :              GridElem(i,j,k)%nbrs(:)=0
        1481    74649600 :              GridElem(i,j,k)%nbrs_wgt(:)=0
        1482    82944000 :              GridElem(i,j,k)%nbrs_ptr(:)=0
        1483    74649600 :              GridElem(i,j,k)%nbrs_wgt_ghost(:)=1  ! always this value
        1484     8294400 :              GridElem(i,j,k)%SpaceCurve=0
        1485     8294400 :              GridElem(i,j,k)%number=number
        1486     8570880 :              number=number+1
        1487             :
        1488             :           end do
        1489             :        end do
        1490             :     end do
        1491             :
        1492        6144 :     allocate(Mesh(ne,ne))
        1493        1536 :     if(IsFactorable(ne)) then
        1494        1536 :        call GenspaceCurve(Mesh)
        1495             :     else
        1496             :        ! find the smallest ne2 which is a power of 2 and ne2>ne
        1497           0 :        ne2 = 2**ceiling(log(real(ne)) / log(2.0_r8))
        1498           0 :        if (ne2 < ne) then
        1499           0 :          call endrun('Fatal SFC error')
        1500             :        end if
        1501             :
        1502           0 :        allocate(Mesh2(ne2,ne2))
        1503           0 :        allocate(Mesh2_map(ne2,ne2,2))
        1504           0 :        allocate(sfcij(0:ne2*ne2,2))
        1505             :
        1506           0 :        call GenspaceCurve(Mesh2)  ! SFC partition for ne2
        1507             :
        1508             :        ! associate every element on the ne x ne mesh (Mesh)
        1509             :        ! with its closest element on the ne2 x ne2 mesh (Mesh2)
        1510             :        ! Store this as a map from Mesh2 -> Mesh in Mesh2_map.
        1511             :        ! elements in Mesh2 which are not mapped get assigned a value of 0
        1512           0 :        Mesh2_map=0
        1513           0 :        do j=1,ne
        1514           0 :           do i=1,ne
        1515             :              ! map this element to an (i2,j2) element
        1516             :              ! [ (i-.5)/ne , (j-.5)/ne ]  = [ (i2-.5)/ne2 , (j2-.5)/ne2 ]
        1517           0 :              i2=nint( ((i-0.5_r8)/ne)*ne2 + 0.5_r8 )
        1518           0 :              j2=nint( ((j-0.5_r8)/ne)*ne2 + 0.5_r8 )
        1519           0 :              if (i2<1) i2=1
        1520           0 :              if (i2>ne2) i2=ne2
        1521           0 :              if (j2<1) j2=1
        1522           0 :              if (j2>ne2) j2=ne2
        1523           0 :              Mesh2_map(i2,j2,1)=i
        1524           0 :              Mesh2_map(i2,j2,2)=j
        1525             :           enddo
        1526             :        enddo
        1527             :
        1528             :        ! create a reverse index array for Mesh2
        1529             :        ! k = Mesh2(i,j)
        1530             :        ! (i,j) = (sfcij(k,1),sfci(k,2))
        1531           0 :        do j=1,ne2
        1532           0 :           do i=1,ne2
        1533           0 :              k=Mesh2(i,j)
        1534           0 :              sfcij(k,1)=i
        1535           0 :              sfcij(k,2)=j
        1536             :           enddo
        1537             :        enddo
        1538             :
        1539             :        ! generate a SFC for Mesh with the same ordering as the
        1540             :        ! elements in Mesh2 which map to Mesh.
        1541             :        sfc_index=0
        1542           0 :        do k=0,ne2*ne2-1
        1543           0 :           i2=sfcij(k,1)
        1544           0 :           j2=sfcij(k,2)
        1545           0 :           i=Mesh2_map(i2,j2,1)
        1546           0 :           j=Mesh2_map(i2,j2,2)
        1547           0 :           if (i/=0) then
        1548             :              ! (i2,j2) element maps to (i,j) element
        1549           0 :              Mesh(i,j)=sfc_index
        1550           0 :              sfc_index=sfc_index+1
        1551             :           endif
        1552             :        enddo
        1553             :
        1554           0 :        deallocate(Mesh2)
        1555           0 :        deallocate(Mesh2_map)
        1556           0 :        deallocate(sfcij)
        1557             :     endif
        1558             :
        1559             :     ! -------------------------------------------
        1560             :     !  Setup the space-filling curve for face 1
        1561             :     ! -------------------------------------------
        1562        1536 :     offset=0
        1563       47616 :     do j=1,ne
        1564     1430016 :        do i=1,ne
        1565     1428480 :           GridElem(i,j,1)%SpaceCurve = offset + Mesh(i,ne-j+1)
        1566             :        enddo
        1567             :     enddo
        1568             :
        1569             :     ! -------------------------------------------
        1570             :     !  Setup the space-filling curve for face 2
        1571             :     ! -------------------------------------------
        1572        1536 :     offset = offset + ne*ne
        1573       47616 :     do j=1,ne
        1574     1430016 :        do i=1,ne
        1575     1428480 :           GridElem(i,j,2)%SpaceCurve = offset + Mesh(i,ne-j+1)
        1576             :        enddo
        1577             :     enddo
        1578             :
        1579             :     ! -------------------------------------------
        1580             :     !  Setup the space-filling curve for face 6
        1581             :     ! -------------------------------------------
        1582        1536 :     offset = offset + ne*ne
        1583       47616 :     do j=1,ne
        1584     1430016 :        do i=1,ne
        1585     1428480 :           GridElem(i,j,6)%SpaceCurve = offset + Mesh(ne-i+1,ne-j+1)
        1586             :        enddo
        1587             :     enddo
        1588             :
        1589             :     ! -------------------------------------------
        1590             :     !  Setup the space-filling curve for face 4
        1591             :     ! -------------------------------------------
        1592        1536 :     offset = offset + ne*ne
        1593       47616 :     do j=1,ne
        1594     1430016 :        do i=1,ne
        1595     1428480 :           GridElem(i,j,4)%SpaceCurve = offset + Mesh(ne-j+1,i)
        1596             :        enddo
        1597             :     enddo
        1598             :
        1599             :     ! -------------------------------------------
        1600             :     !  Setup the space-filling curve for face 5
        1601             :     ! -------------------------------------------
        1602        1536 :     offset = offset + ne*ne
        1603       47616 :     do j=1,ne
        1604     1430016 :        do i=1,ne
        1605     1428480 :           GridElem(i,j,5)%SpaceCurve = offset + Mesh(i,j)
        1606             :        enddo
        1607             :     enddo
        1608             :
        1609             :
        1610             :     ! -------------------------------------------
        1611             :     !  Setup the space-filling curve for face 3
        1612             :     ! -------------------------------------------
        1613        1536 :     offset = offset + ne*ne
        1614       47616 :     do j=1,ne
        1615     1430016 :        do i=1,ne
        1616     1428480 :           GridElem(i,j,3)%SpaceCurve = offset + Mesh(i,j)
        1617             :        enddo
        1618             :     enddo
        1619             :
        1620             :     ! ==================
        1621             :     ! face interiors
        1622             :     ! ==================
        1623       10752 :     do k=1,6
        1624             :        ! setup  SOUTH, WEST, SW neighbors
        1625      276480 :        do j=2,ne
        1626     8027136 :           do i=2,ne
        1627     7750656 :              nbrs_used(i,j,k,west) = .true.
        1628     7750656 :              nbrs_used(i,j,k,south) = .true.
        1629     7750656 :              nbrs_used(i,j,k,swest) = .true.
        1630             :
        1631             :
        1632     7750656 :              GridElem(i,j,k)%nbrs(west)  = GridElem(i-1,j,k)%number
        1633     7750656 :              GridElem(i,j,k)%nbrs_face(west)  = k
        1634     7750656 :              GridElem(i,j,k)%nbrs_wgt(west)       = EdgeWgtP
        1635     7750656 :              GridElem(i,j,k)%nbrs(south) = GridElem(i,j-1,k)%number
        1636     7750656 :              GridElem(i,j,k)%nbrs_face(south) = k
        1637     7750656 :              GridElem(i,j,k)%nbrs_wgt(south)      = EdgeWgtP
        1638     7750656 :              GridElem(i,j,k)%nbrs(swest) = GridElem(i-1,j-1,k)%number
        1639     7750656 :              GridElem(i,j,k)%nbrs_face(swest) = k
        1640     8017920 :              GridElem(i,j,k)%nbrs_wgt(swest)      = CornerWgt
        1641             :           end do
        1642             :        end do
        1643             :
        1644             :        !  setup EAST, NORTH, NE neighbors
        1645      276480 :        do j=1,ne-1
        1646     8027136 :           do i=1,ne-1
        1647     7750656 :              nbrs_used(i,j,k,east) = .true.
        1648     7750656 :              nbrs_used(i,j,k,north) = .true.
        1649     7750656 :              nbrs_used(i,j,k,neast) = .true.
        1650             :
        1651     7750656 :              GridElem(i,j,k)%nbrs(east)   = GridElem(i+1,j,k)%number
        1652     7750656 :              GridElem(i,j,k)%nbrs_face(east)   = k
        1653     7750656 :              GridElem(i,j,k)%nbrs_wgt(east)        = EdgeWgtP
        1654     7750656 :              GridElem(i,j,k)%nbrs(north)  = GridElem(i,j+1,k)%number
        1655     7750656 :              GridElem(i,j,k)%nbrs_face(north)  = k
        1656     7750656 :              GridElem(i,j,k)%nbrs_wgt(north)       = EdgeWgtP
        1657     7750656 :              GridElem(i,j,k)%nbrs(neast) = GridElem(i+1,j+1,k)%number
        1658     7750656 :              GridElem(i,j,k)%nbrs_face(neast)  = k
        1659     8017920 :              GridElem(i,j,k)%nbrs_wgt(neast)       = CornerWgt
        1660             :           end do
        1661             :        end do
        1662             :
        1663             :        ! Setup the remaining SOUTH, EAST, and SE neighbors
        1664      276480 :        do j=2,ne
        1665     8027136 :           do i=1,ne-1
        1666     7750656 :              nbrs_used(i,j,k,south) = .true.
        1667     7750656 :              nbrs_used(i,j,k,east) = .true.
        1668     7750656 :              nbrs_used(i,j,k,seast) = .true.
        1669             :
        1670             :
        1671             :
        1672     7750656 :              GridElem(i,j,k)%nbrs(south)  = GridElem(i,j-1,k)%number
        1673     7750656 :              GridElem(i,j,k)%nbrs_face(south)  = k
        1674     7750656 :              GridElem(i,j,k)%nbrs_wgt(south)       = EdgeWgtP
        1675     7750656 :              GridElem(i,j,k)%nbrs(east)   = GridElem(i+1,j,k)%number
        1676     7750656 :              GridElem(i,j,k)%nbrs_face(east)   = k
        1677     7750656 :              GridElem(i,j,k)%nbrs_wgt(east)        = EdgeWgtP
        1678     7750656 :              GridElem(i,j,k)%nbrs(seast)  = GridElem(i+1,j-1,k)%number
        1679     7750656 :              GridElem(i,j,k)%nbrs_face(seast)  = k
        1680     8017920 :              GridElem(i,j,k)%nbrs_wgt(seast)       = CornerWgt
        1681             :           enddo
        1682             :        enddo
        1683             :
        1684             :        ! Setup the remaining NORTH, WEST, and NW neighbors
        1685      278016 :        do j=1,ne-1
        1686     8027136 :           do i=2,ne
        1687     7750656 :              nbrs_used(i,j,k,north) = .true.
        1688     7750656 :              nbrs_used(i,j,k,west) = .true.
        1689     7750656 :              nbrs_used(i,j,k,nwest) = .true.
        1690             :
        1691             :
        1692             :
        1693     7750656 :              GridElem(i,j,k)%nbrs(north)  = GridElem(i,j+1,k)%number
        1694     7750656 :              GridElem(i,j,k)%nbrs_face(north)  = k
        1695     7750656 :              GridElem(i,j,k)%nbrs_wgt(north)       = EdgeWgtP
        1696     7750656 :              GridElem(i,j,k)%nbrs(west)   = GridElem(i-1,j,k)%number
        1697     7750656 :              GridElem(i,j,k)%nbrs_face(west)   = k
        1698     7750656 :              GridElem(i,j,k)%nbrs_wgt(west)        = EdgeWgtP
        1699     7750656 :              GridElem(i,j,k)%nbrs(nwest)  = GridElem(i-1,j+1,k)%number
        1700     7750656 :              GridElem(i,j,k)%nbrs_face(nwest)  = k
        1701     8017920 :              GridElem(i,j,k)%nbrs_wgt(nwest)       = CornerWgt
        1702             :           enddo
        1703             :        enddo
        1704             :     end do
        1705             :
        1706             :     ! ======================
        1707             :     ! west/east "belt" edges
        1708             :     ! ======================
        1709             :
        1710        7680 :     do k=1,4
        1711      192000 :        do j=1,ne
        1712      184320 :           nbrs_used(1,j,k,west) = .true.
        1713      184320 :           nbrs_used(ne,j,k,east) = .true.
        1714             :
        1715             :
        1716      184320 :           GridElem(1 ,j,k)%nbrs(west) = GridElem(ne,j,MODULO(2+k,4)+1)%number
        1717      184320 :           GridElem(1 ,j,k)%nbrs_face(west) = MODULO(2+k,4)+1
        1718      184320 :           GridElem(1 ,j,k)%nbrs_wgt(west)  = EdgeWgtP
        1719      184320 :           GridElem(ne,j,k)%nbrs(east) = GridElem(1 ,j,MODULO(k  ,4)+1)%number
        1720      184320 :           GridElem(ne,j,k)%nbrs_face(east) = MODULO(k  ,4)+1
        1721      184320 :           GridElem(ne,j,k)%nbrs_wgt(east)  = EdgeWgtP
        1722             :
        1723             :           !  Special rules for corner 'edges'
        1724      184320 :           if( j /= 1) then
        1725      178176 :              nbrs_used(1,j,k,swest) = .true.
        1726      178176 :              nbrs_used(ne,j,k,seast) = .true.
        1727             :
        1728             :
        1729      178176 :              GridElem(1 ,j,k)%nbrs(swest)   = GridElem(ne,j-1,MODULO(2+k,4)+1)%number
        1730      178176 :              GridElem(1 ,j,k)%nbrs_face(swest)   = MODULO(2+k,4)+1
        1731      178176 :              GridElem(1 ,j,k)%nbrs_wgt(swest)        = CornerWgt
        1732      178176 :              GridElem(ne,j,k)%nbrs(seast)   = GridElem(1 ,j-1,MODULO(k  ,4)+1)%number
        1733      178176 :              GridElem(ne,j,k)%nbrs_face(seast)   = MODULO(k  ,4)+1
        1734      178176 :              GridElem(ne,j,k)%nbrs_wgt(seast)        = CornerWgt
        1735             :           endif
        1736      190464 :           if( j /= ne) then
        1737      178176 :              nbrs_used(1,j,k,nwest) = .true.
        1738      178176 :              nbrs_used(ne,j,k,neast) = .true.
        1739             :
        1740             :
        1741      178176 :              GridElem(1 ,j,k)%nbrs(nwest)   = GridElem(ne,j+1,MODULO(2+k,4)+1)%number
        1742      178176 :              GridElem(1 ,j,k)%nbrs_face(nwest)   = MODULO(2+k,4)+1
        1743      178176 :              GridElem(1 ,j,k)%nbrs_wgt(nwest)        = CornerWgt
        1744      178176 :              GridElem(ne,j,k)%nbrs(neast)   = GridElem(1 ,j+1,MODULO(k  ,4)+1)%number
        1745      178176 :              GridElem(ne,j,k)%nbrs_face(neast)   = MODULO(k  ,4)+1
        1746      178176 :              GridElem(ne,j,k)%nbrs_wgt(neast)        = CornerWgt
        1747             :           endif
        1748             :        end do
        1749             :     end do
        1750             :
        1751             :
        1752             :     ! ==================================
        1753             :     ! south edge of 1 / north edge of 5
        1754             :     ! ==================================
        1755             :
        1756       47616 :     do i=1,ne
        1757       46080 :        nbrs_used(i,1,1,south) = .true.
        1758       46080 :        nbrs_used(i,ne,5,north) = .true.
        1759             :
        1760       46080 :        GridElem(i,1 ,1)%nbrs(south) = GridElem(i,ne,5)%number
        1761       46080 :        GridElem(i,1 ,1)%nbrs_face(south) = 5
        1762       46080 :        GridElem(i,1 ,1)%nbrs_wgt(south)      = EdgeWgtP
        1763       46080 :        GridElem(i,ne,5)%nbrs(north) = GridElem(i,1 ,1)%number
        1764       46080 :        GridElem(i,ne,5)%nbrs_face(north) = 1
        1765       46080 :        GridElem(i,ne,5)%nbrs_wgt(north)      = EdgeWgtP
        1766             :
        1767             :        !  Special rules for corner 'edges'
        1768       46080 :        if( i /= 1) then
        1769       44544 :           nbrs_used(i,1,1,swest) = .true.
        1770       44544 :           nbrs_used(i,ne,5,nwest) = .true.
        1771             :
        1772       44544 :           GridElem(i,1 ,1)%nbrs(swest)    = GridElem(i-1,ne,5)%number
        1773       44544 :           GridElem(i,1 ,1)%nbrs_face(swest)    = 5
        1774       44544 :           GridElem(i,1 ,1)%nbrs_wgt(swest)         = CornerWgt
        1775       44544 :           GridElem(i,ne,5)%nbrs(nwest)    = GridElem(i-1,1 ,1)%number
        1776       44544 :           GridElem(i,ne,5)%nbrs_face(nwest)    = 1
        1777       44544 :           GridElem(i,ne,5)%nbrs_wgt(nwest)         = CornerWgt
        1778             :        endif
        1779       47616 :        if( i /= ne) then
        1780       44544 :           nbrs_used(i,1,1,seast) = .true.
        1781       44544 :           nbrs_used(i,ne,5,neast) = .true.
        1782             :
        1783       44544 :           GridElem(i,1 ,1)%nbrs(seast)    = GridElem(i+1,ne,5)%number
        1784       44544 :           GridElem(i,1 ,1)%nbrs_face(seast)    = 5
        1785       44544 :           GridElem(i,1 ,1)%nbrs_wgt(seast)         = CornerWgt
        1786       44544 :           GridElem(i,ne,5)%nbrs(neast)    = GridElem(i+1,1 ,1)%number
        1787       44544 :           GridElem(i,ne,5)%nbrs_face(neast)    = 1
        1788       44544 :           GridElem(i,ne,5)%nbrs_wgt(neast)         = CornerWgt
        1789             :        endif
        1790             :
        1791             :     end do
        1792             :
        1793             :     ! ==================================
        1794             :     ! south edge of 2 / east edge of 5
        1795             :     ! ==================================
        1796             :
        1797       47616 :     do i=1,ne
        1798       46080 :        irev=ne+1-i
        1799       46080 :        nbrs_used(i,1,2,south) = .true.
        1800       46080 :        nbrs_used(ne,i,5,east) = .true.
        1801             :
        1802             :
        1803       46080 :        GridElem(i,1 ,2)%nbrs(south) = GridElem(ne,irev,5)%number
        1804       46080 :        GridElem(i,1 ,2)%nbrs_face(south) = 5
        1805       46080 :        GridElem(i,1 ,2)%nbrs_wgt(south)      = EdgeWgtP
        1806       46080 :        GridElem(ne,i,5)%nbrs(east)  = GridElem(irev,1 ,2)%number
        1807       46080 :        GridElem(ne,i,5)%nbrs_face(east)  = 2
        1808       46080 :        GridElem(ne,i,5)%nbrs_wgt(east)       = EdgeWgtP
        1809             :
        1810             :        !  Special rules for corner 'edges'
        1811       46080 :        if( i /= 1) then
        1812       44544 :           nbrs_used(i,1,2,swest) = .true.
        1813       44544 :           nbrs_used(ne,i,5,seast) = .true.
        1814             :
        1815             :
        1816       44544 :           GridElem(i,1 ,2)%nbrs(swest) = GridElem(ne,irev+1,5)%number
        1817       44544 :           GridElem(i,1 ,2)%nbrs_face(swest) = 5
        1818       44544 :           GridElem(i,1 ,2)%nbrs_wgt(swest)      = CornerWgt
        1819       44544 :           GridElem(ne,i,5)%nbrs(seast) = GridElem(irev+1,1 ,2)%number
        1820       44544 :           GridElem(ne,i,5)%nbrs_face(seast) = 2
        1821       44544 :           GridElem(ne,i,5)%nbrs_wgt(seast)      = CornerWgt
        1822             :        endif
        1823       47616 :        if(i /= ne) then
        1824       44544 :           nbrs_used(i,1,2,seast) = .true.
        1825       44544 :           nbrs_used(ne,i,5,neast) = .true.
        1826             :
        1827             :
        1828       44544 :           GridElem(i,1 ,2)%nbrs(seast)   = GridElem(ne,irev-1,5)%number
        1829       44544 :           GridElem(i,1 ,2)%nbrs_face(seast)   = 5
        1830       44544 :           GridElem(i,1 ,2)%nbrs_wgt(seast)        = CornerWgt
        1831       44544 :           GridElem(ne,i,5)%nbrs(neast)   = GridElem(irev-1,1 ,2)%number
        1832       44544 :           GridElem(ne,i,5)%nbrs_face(neast)   = 2
        1833       44544 :           GridElem(ne,i,5)%nbrs_wgt(neast)        = CornerWgt
        1834             :        endif
        1835             :     enddo
        1836             :     ! ==================================
        1837             :     ! south edge of 3 / south edge of 5
        1838             :     ! ==================================
        1839             :
        1840       47616 :     do i=1,ne
        1841       46080 :        irev=ne+1-i
        1842       46080 :        nbrs_used(i,1,3,south) = .true.
        1843       46080 :        nbrs_used(i,1,5,south) = .true.
        1844             :
        1845       46080 :        GridElem(i,1,3)%nbrs(south) = GridElem(irev,1,5)%number
        1846       46080 :        GridElem(i,1,3)%nbrs_face(south) = 5
        1847       46080 :        GridElem(i,1,3)%nbrs_wgt(south)      = EdgeWgtP
        1848       46080 :        GridElem(i,1,5)%nbrs(south) = GridElem(irev,1,3)%number
        1849       46080 :        GridElem(i,1,5)%nbrs_face(south) = 3
        1850       46080 :        GridElem(i,1,5)%nbrs_wgt(south)      = EdgeWgtP
        1851             :
        1852             :        !  Special rules for corner 'edges'
        1853       46080 :        if( i /= 1) then
        1854       44544 :           nbrs_used(i,1,3,swest) = .true.
        1855       44544 :           nbrs_used(i,1,5,swest) = .true.
        1856             :
        1857             :
        1858       44544 :           GridElem(i,1,3)%nbrs(swest) = GridElem(irev+1,1,5)%number
        1859       44544 :           GridElem(i,1,3)%nbrs_face(swest) = 5
        1860       44544 :           GridElem(i,1,3)%nbrs_wgt(swest)      = CornerWgt
        1861       44544 :           GridElem(i,1,5)%nbrs(swest) = GridElem(irev+1,1,3)%number
        1862       44544 :           GridElem(i,1,5)%nbrs_face(swest) = 3
        1863       44544 :           GridElem(i,1,5)%nbrs_wgt(swest)      = CornerWgt
        1864             :        endif
        1865       47616 :        if(i /= ne) then
        1866       44544 :           nbrs_used(i,1,3,seast) = .true.
        1867       44544 :           nbrs_used(i,1,5,seast) = .true.
        1868             :
        1869       44544 :           GridElem(i,1,3)%nbrs(seast)    = GridElem(irev-1,1,5)%number
        1870       44544 :           GridElem(i,1,3)%nbrs_face(seast)    = 5
        1871       44544 :           GridElem(i,1,3)%nbrs_wgt(seast)         = CornerWgt
        1872       44544 :           GridElem(i,1,5)%nbrs(seast)    = GridElem(irev-1,1,3)%number
        1873       44544 :           GridElem(i,1,5)%nbrs_face(seast)    = 3
        1874       44544 :           GridElem(i,1,5)%nbrs_wgt(seast)         = CornerWgt
        1875             :        endif
        1876             :     end do
        1877             :
        1878             :     ! ==================================
        1879             :     ! south edge of 4 / west edge of 5
        1880             :     ! ==================================
        1881             :
        1882       47616 :     do i=1,ne
        1883       46080 :        irev=ne+1-i
        1884       46080 :        nbrs_used(i,1,4,south) = .true.
        1885       46080 :        nbrs_used(1,i,5,west) = .true.
        1886             :
        1887       46080 :        GridElem(i,1,4)%nbrs(south) = GridElem(1,i,5)%number
        1888       46080 :        GridElem(i,1,4)%nbrs_face(south) = 5
        1889       46080 :        GridElem(i,1,4)%nbrs_wgt(south)      = EdgeWgtP
        1890       46080 :        GridElem(1,i,5)%nbrs(west)  = GridElem(i,1,4)%number
        1891       46080 :        GridElem(1,i,5)%nbrs_face(west)  = 4
        1892       46080 :        GridElem(1,i,5)%nbrs_wgt(west)       = EdgeWgtP
        1893             :        !  Special rules for corner 'edges'
        1894       46080 :        if( i /= 1) then
        1895       44544 :           nbrs_used(i,1,4,swest) = .true.
        1896       44544 :           nbrs_used(1,i,5,swest) = .true.
        1897             :
        1898       44544 :           GridElem(i,1,4)%nbrs(swest)    = GridElem(1,i-1,5)%number
        1899       44544 :           GridElem(i,1,4)%nbrs_face(swest)    = 5
        1900       44544 :           GridElem(i,1,4)%nbrs_wgt(swest)         = CornerWgt
        1901       44544 :           GridElem(1,i,5)%nbrs(swest)    = GridElem(i-1,1,4)%number
        1902       44544 :           GridElem(1,i,5)%nbrs_face(swest)    = 4
        1903       44544 :           GridElem(1,i,5)%nbrs_wgt(swest)         = CornerWgt
        1904             :        endif
        1905       47616 :        if( i /= ne) then
        1906       44544 :           nbrs_used(i,1,4,seast) = .true.
        1907       44544 :           nbrs_used(1,i,5,nwest) = .true.
        1908             :
        1909       44544 :           GridElem(i,1,4)%nbrs(seast) = GridElem(1,i+1,5)%number
        1910       44544 :           GridElem(i,1,4)%nbrs_face(seast) = 5
        1911       44544 :           GridElem(i,1,4)%nbrs_wgt(seast)      = CornerWgt
        1912       44544 :           GridElem(1,i,5)%nbrs(nwest) = GridElem(i+1,1,4)%number
        1913       44544 :           GridElem(1,i,5)%nbrs_face(nwest) = 4
        1914       44544 :           GridElem(1,i,5)%nbrs_wgt(nwest)      = CornerWgt
        1915             :        endif
        1916             :     end do
        1917             :
        1918             :     ! ==================================
        1919             :     ! north edge of 1 / south edge of 6
        1920             :     ! ==================================
        1921             :
        1922       47616 :     do i=1,ne
        1923       46080 :        nbrs_used(i,ne,1,north) = .true.
        1924       46080 :        nbrs_used(i,1,6,south) = .true.
        1925             :
        1926             :
        1927       46080 :        GridElem(i,ne,1)%nbrs(north) = GridElem(i,1 ,6)%number
        1928       46080 :        GridElem(i,ne,1)%nbrs_face(north) = 6
        1929       46080 :        GridElem(i,ne,1)%nbrs_wgt(north)      = EdgeWgtP
        1930       46080 :        GridElem(i,1 ,6)%nbrs(south) = GridElem(i,ne,1)%number
        1931       46080 :        GridElem(i,1 ,6)%nbrs_face(south) = 1
        1932       46080 :        GridElem(i,1 ,6)%nbrs_wgt(south)      = EdgeWgtP
        1933             :        !  Special rules for corner 'edges'
        1934       46080 :        if( i /= 1) then
        1935       44544 :           nbrs_used(i,ne,1,nwest) = .true.
        1936       44544 :           nbrs_used(i,1,6,swest) = .true.
        1937             :
        1938       44544 :           GridElem(i,ne,1)%nbrs(nwest) = GridElem(i-1,1 ,6)%number
        1939       44544 :           GridElem(i,ne,1)%nbrs_face(nwest) = 6
        1940       44544 :           GridElem(i,ne,1)%nbrs_wgt(nwest)      = CornerWgt
        1941       44544 :           GridElem(i,1 ,6)%nbrs(swest) = GridElem(i-1,ne,1)%number
        1942       44544 :           GridElem(i,1 ,6)%nbrs_face(swest) = 1
        1943       44544 :           GridElem(i,1 ,6)%nbrs_wgt(swest)      = CornerWgt
        1944             :        endif
        1945       47616 :        if( i /= ne) then
        1946       44544 :           nbrs_used(i,ne,1,neast) = .true.
        1947       44544 :           nbrs_used(i,1,6,seast) = .true.
        1948             :
        1949             :
        1950       44544 :           GridElem(i,ne,1)%nbrs(neast) = GridElem(i+1,1 ,6)%number
        1951       44544 :           GridElem(i,ne,1)%nbrs_face(neast) = 6
        1952       44544 :           GridElem(i,ne,1)%nbrs_wgt(neast)      = CornerWgt
        1953       44544 :           GridElem(i,1 ,6)%nbrs(seast) = GridElem(i+1,ne,1)%number
        1954       44544 :           GridElem(i,1 ,6)%nbrs_face(seast) = 1
        1955       44544 :           GridElem(i,1 ,6)%nbrs_wgt(seast)      = CornerWgt
        1956             :        endif
        1957             :     end do
        1958             :
        1959             :     ! ==================================
        1960             :     ! north edge of 2 / east edge of 6
        1961             :     ! ==================================
        1962             :
        1963       47616 :     do i=1,ne
        1964       46080 :        nbrs_used(i,ne,2,north) = .true.
        1965       46080 :        nbrs_used(ne,i,6,east ) = .true.
        1966             :
        1967       46080 :        GridElem(i,ne,2)%nbrs(north) = GridElem(ne,i,6)%number
        1968       46080 :        GridElem(i,ne,2)%nbrs_face(north) = 6
        1969       46080 :        GridElem(i,ne,2)%nbrs_wgt(north)      = EdgeWgtP
        1970       46080 :        GridElem(ne,i,6)%nbrs(east)  = GridElem(i,ne,2)%number
        1971       46080 :        GridElem(ne,i,6)%nbrs_face(east)  = 2
        1972       46080 :        GridElem(ne,i,6)%nbrs_wgt(east)       = EdgeWgtP
        1973             :        !  Special rules for corner 'edges'
        1974       46080 :        if( i /= 1) then
        1975       44544 :           nbrs_used(i,ne,2,nwest) = .true.
        1976       44544 :           nbrs_used(ne,i,6,seast) = .true.
        1977             :
        1978       44544 :           GridElem(i,ne,2)%nbrs(nwest)    = GridElem(ne,i-1,6)%number
        1979       44544 :           GridElem(i,ne,2)%nbrs_face(nwest)    = 6
        1980       44544 :           GridElem(i,ne,2)%nbrs_wgt(nwest)         = CornerWgt
        1981       44544 :           GridElem(ne,i,6)%nbrs(seast)    = GridElem(i-1,ne,2)%number
        1982       44544 :           GridElem(ne,i,6)%nbrs_face(seast)    = 2
        1983       44544 :           GridElem(ne,i,6)%nbrs_wgt(seast)         = CornerWgt
        1984             :        endif
        1985       47616 :        if( i /= ne) then
        1986       44544 :           nbrs_used(i,ne,2,neast) = .true.
        1987       44544 :           nbrs_used(ne,i,6,neast) = .true.
        1988             :
        1989             :
        1990       44544 :           GridElem(i,ne,2)%nbrs(neast) = GridElem(ne,i+1,6)%number
        1991       44544 :           GridElem(i,ne,2)%nbrs_face(neast) = 6
        1992       44544 :           GridElem(i,ne,2)%nbrs_wgt(neast)      = CornerWgt
        1993       44544 :           GridElem(ne,i,6)%nbrs(neast) = GridElem(i+1,ne,2)%number
        1994       44544 :           GridElem(ne,i,6)%nbrs_face(neast) = 2
        1995       44544 :           GridElem(ne,i,6)%nbrs_wgt(neast)      = CornerWgt
        1996             :        endif
        1997             :     end do
        1998             :
        1999             :     ! ===================================
        2000             :     ! north edge of 3 / north edge of 6
        2001             :     ! ===================================
        2002             :
        2003       47616 :     do i=1,ne
        2004       46080 :        irev=ne+1-i
        2005       46080 :        nbrs_used(i,ne,3,north) = .true.
        2006       46080 :        nbrs_used(i,ne,6,north) = .true.
        2007             :
        2008       46080 :        GridElem(i,ne,3)%nbrs(north) = GridElem(irev,ne,6)%number
        2009       46080 :        GridElem(i,ne,3)%nbrs_face(north) = 6
        2010       46080 :        GridElem(i,ne,3)%nbrs_wgt(north)      = EdgeWgtP
        2011       46080 :        GridElem(i,ne,6)%nbrs(north) = GridElem(irev,ne,3)%number
        2012       46080 :        GridElem(i,ne,6)%nbrs_face(north) = 3
        2013       46080 :        GridElem(i,ne,6)%nbrs_wgt(north)      = EdgeWgtP
        2014             :        !  Special rules for corner 'edges'
        2015       46080 :        if( i /= 1) then
        2016       44544 :           nbrs_used(i,ne,3,nwest) = .true.
        2017       44544 :           nbrs_used(i,ne,6,nwest) = .true.
        2018             :
        2019       44544 :           GridElem(i,ne,3)%nbrs(nwest) = GridElem(irev+1,ne,6)%number
        2020       44544 :           GridElem(i,ne,3)%nbrs_face(nwest) = 6
        2021       44544 :           GridElem(i,ne,3)%nbrs_wgt(nwest)      = CornerWgt
        2022       44544 :           GridElem(i,ne,6)%nbrs(nwest) = GridElem(irev+1,ne,3)%number
        2023       44544 :           GridElem(i,ne,6)%nbrs_face(nwest) = 3
        2024       44544 :           GridElem(i,ne,6)%nbrs_wgt(nwest)      = CornerWgt
        2025             :        endif
        2026       47616 :        if( i /= ne) then
        2027       44544 :           nbrs_used(i,ne,3,neast) = .true.
        2028       44544 :           nbrs_used(i,ne,6,neast) = .true.
        2029             :
        2030       44544 :           GridElem(i,ne,3)%nbrs(neast) = GridElem(irev-1,ne,6)%number
        2031       44544 :           GridElem(i,ne,3)%nbrs_face(neast) = 6
        2032       44544 :           GridElem(i,ne,3)%nbrs_wgt(neast)      = CornerWgt
        2033       44544 :           GridElem(i,ne,6)%nbrs(neast) = GridElem(irev-1,ne,3)%number
        2034       44544 :           GridElem(i,ne,6)%nbrs_face(neast) = 3
        2035       44544 :           GridElem(i,ne,6)%nbrs_wgt(neast)      = CornerWgt
        2036             :        endif
        2037             :     end do
        2038             :
        2039             :     ! ===================================
        2040             :     ! north edge of 4 / west edge of 6
        2041             :     ! ===================================
        2042             :
        2043       47616 :     do i=1,ne
        2044       46080 :        irev=ne+1-i
        2045       46080 :        nbrs_used(i,ne,4,north) = .true.
        2046       46080 :        nbrs_used(1,i,6,west) = .true.
        2047             :
        2048       46080 :        GridElem(i,ne,4)%nbrs(north) = GridElem(1,irev,6)%number
        2049       46080 :        GridElem(i,ne,4)%nbrs_face(north) = 6
        2050       46080 :        GridElem(i,ne,4)%nbrs_wgt(north)      = EdgeWgtP
        2051       46080 :        GridElem(1,i,6)%nbrs(west)   = GridElem(irev,ne,4)%number
        2052       46080 :        GridElem(1,i,6)%nbrs_face(west)   = 4
        2053       46080 :        GridElem(1,i,6)%nbrs_wgt(west)        = EdgeWgtP
        2054             :        !  Special rules for corner 'edges'
        2055       46080 :        if( i /= 1) then
        2056       44544 :           nbrs_used(i,ne,4,nwest) = .true.
        2057       44544 :           nbrs_used(1,i,6,swest) = .true.
        2058             :
        2059       44544 :           GridElem(i,ne,4)%nbrs(nwest) = GridElem(1,irev+1,6)%number
        2060       44544 :           GridElem(i,ne,4)%nbrs_face(nwest) = 6
        2061       44544 :           GridElem(i,ne,4)%nbrs_wgt(nwest)      = CornerWgt
        2062       44544 :           GridElem(1,i,6)%nbrs(swest)  = GridElem(irev+1,ne,4)%number
        2063       44544 :           GridElem(1,i,6)%nbrs_face(swest)  = 4
        2064       44544 :           GridElem(1,i,6)%nbrs_wgt(swest)       = CornerWgt
        2065             :        endif
        2066       47616 :        if( i /= ne) then
        2067       44544 :           nbrs_used(i,ne,4,neast) = .true.
        2068       44544 :           nbrs_used(1,i,6,nwest) = .true.
        2069             :
        2070       44544 :           GridElem(i,ne,4)%nbrs(neast) = GridElem(1,irev-1,6)%number
        2071       44544 :           GridElem(i,ne,4)%nbrs_face(neast) = 6
        2072       44544 :           GridElem(i,ne,4)%nbrs_wgt(neast)      = CornerWgt
        2073       44544 :           GridElem(1,i,6)%nbrs(nwest)  = GridElem(irev-1,ne,4)%number
        2074       44544 :           GridElem(1,i,6)%nbrs_face(nwest)  = 4
        2075       44544 :           GridElem(1,i,6)%nbrs_wgt(nwest)       = CornerWgt
        2076             :        endif
        2077             :     end do
        2078             :
        2079             :
        2080             :     ielem = 1                       ! Element counter
        2081       10752 :     do k=1,6
        2082      287232 :        do j=1,ne
        2083     8580096 :           do i=1,ne
        2084     8294400 :              GridVertex(ielem)%nbrs_ptr(1) = 1
        2085    74649600 :              do ll=1,8
        2086    66355200 :                 loc =  GridVertex(ielem)%nbrs_ptr(ll)
        2087    74649600 :                 if (nbrs_used(i,j,k,ll)) then
        2088    66318336 :                    GridVertex(ielem)%nbrs(loc)       = GridElem(i,j,k)%nbrs(ll)
        2089    66318336 :                    GridVertex(ielem)%nbrs_face(loc)  = GridElem(i,j,k)%nbrs_face(ll)
        2090    66318336 :                    GridVertex(ielem)%nbrs_wgt(loc)       = GridElem(i,j,k)%nbrs_wgt(ll)
        2091    66318336 :                    GridVertex(ielem)%nbrs_wgt_ghost(loc) = GridElem(i,j,k)%nbrs_wgt_ghost(ll)
        2092             :
        2093    66318336 :                    GridVertex(ielem)%nbrs_ptr(ll+1) = GridVertex(ielem)%nbrs_ptr(ll)+1
        2094             :                 else
        2095       36864 :                    GridVertex(ielem)%nbrs_ptr(ll+1) = GridVertex(ielem)%nbrs_ptr(ll)
        2096             :                 end if
        2097             :              end do
        2098     8294400 :              GridVertex(ielem)%number     = GridElem(i,j,k)%number
        2099     8294400 :              GridVertex(ielem)%processor_number  = 0
        2100     8294400 :              GridVertex(ielem)%SpaceCurve = GridElem(i,j,k)%SpaceCurve
        2101     8570880 :              ielem=ielem+1
        2102             :           end do
        2103             :        end do
        2104             :     end do
        2105             :
        2106        1536 :     DEALLOCATE(Mesh)
        2107       10752 :      do k = 1, nfaces
        2108      287232 :        do j = 1, ne
        2109     8580096 :           do i = 1, ne
        2110     8570880 :              call deallocate_gridvertex_nbrs(GridElem(i,j,k))
        2111             :           end do
        2112             :        end do
        2113             :     end do
        2114        1536 :     DEALLOCATE(GridElem)
        2115        1536 :     DEALLOCATE(nbrs_used)
        2116             :
        2117             :     ! =======================================
        2118             :     ! Generate cube graph...
        2119             :     ! =======================================
        2120             :
        2121             :     ! ============================================
        2122             :     !  Setup the Grid edges (topology independent)
        2123             :     ! ============================================
        2124        1536 :     call initgridedge(GridEdge,GridVertex)
        2125             :
        2126             :     ! ============================================
        2127             :     !  Setup the Grid edge Indirect addresses
        2128             :     !          (topology dependent)
        2129             :     ! ============================================
        2130        1536 :     nedge = SIZE(GridEdge)
        2131    66319872 :     do i=1,nedge
        2132    66319872 :        call CubeSetupEdgeIndex(GridEdge(i))
        2133             :     enddo
        2134             :
        2135        1536 :   end subroutine CubeTopology
        2136             :
        2137             :   ! ===================================================================
        2138             :   ! CubeEdgeCount:
        2139             :   !
        2140             :   !  Determine the number of Grid Edges
        2141             :   !
        2142             :   ! ===================================================================
        2143             :
        2144        1536 :   function CubeEdgeCount()  result(nedge)
        2145             :     use dimensions_mod, only     : ne
        2146             :     implicit none
        2147             :     integer                     :: nedge
        2148             :
        2149        1536 :     if (0==ne) call endrun('Error in CubeEdgeCount: ne is zero')
        2150        1536 :     nedge = nfaces*(ne*ne*nInnerElemEdge - nCornerElemEdge)
        2151             :
        2152        1536 :   end function CubeEdgeCount
        2153             :
        2154             :   ! ===================================================================
        2155             :   ! CubeElemCount:
        2156             :   !
        2157             :   !  Determine the number of Grid Elem
        2158             :   !
        2159             :   ! ===================================================================
        2160             :
        2161        1536 :   function CubeElemCount()  result(nelem)
        2162             :
        2163             :     use dimensions_mod, only     : ne
        2164             :
        2165             :     implicit none
        2166             :     integer                     :: nelem
        2167        1536 :     if (0==ne) call endrun('Error in CubeElemCount: ne is zero')
        2168             :
        2169        1536 :     nelem = nfaces*ne*ne
        2170        1536 :   end function CubeElemCount
        2171             :
        2172    66318336 :   subroutine CubeSetupEdgeIndex(Edge)
        2173             :     use gridgraph_mod, only : gridedge_t
        2174             :     use dimensions_mod, only : np
        2175             :     use control_mod, only : north, south, east, west, neast, seast, swest, nwest
        2176             :     type (GridEdge_t),target           :: Edge
        2177             :
        2178             :     integer                            :: np0,sFace,dFace
        2179             :     logical                            :: reverse
        2180             :     integer,allocatable                :: forwardV(:), forwardP(:)
        2181             :     integer,allocatable                :: backwardV(:), backwardP(:)
        2182             :
        2183    66318336 :     sFace = Edge%tail_face
        2184    66318336 :     dFace = Edge%head_face
        2185             :     ! Do not reverse the indices
        2186    66318336 :     reverse=.FALSE.
        2187             :
        2188             :     ! Under special conditions use index reversal
        2189             :     if(       (SFace == south .AND. dFace == east)  &
        2190             :          .OR. (sFace == east  .AND. dFace == south) &
        2191             :          .OR. (sFace == north .AND. dFace == west)  &
        2192             :          .OR. (sFace == west  .AND. dFace == north) &
        2193             :          .OR. (sFace == south .AND. dFace == south) &
        2194             :          .OR. (sFace == north .AND. dFace == north) &
        2195             :          .OR. (sFace == east  .AND. dFace == east ) &
        2196    66318336 :          .OR. (sFace == west  .AND. dFace == west ) ) then
        2197      368640 :        reverse=.TRUE.
        2198      368640 :        Edge%reverse=.TRUE.
        2199             :     endif
        2200             :
        2201             :
        2202    66318336 :   end subroutine CubeSetupEdgeIndex
        2203             :
        2204             : !
        2205             : !  HOMME mapping from sphere (or other manifold) to reference element
        2206             : !  one should be able to add any mapping here.  For each new map,
        2207             : !  an associated dmap() routine (which computes the map derivative matrix)
        2208             : !  must also be written
        2209             : !  Note that for conservation, the parameterization of element edges must be
        2210             : !  identical for adjacent elements.  (this is violated with HOMME's default
        2211             : !  equi-angular cubed-sphere mapping for non-cubed sphere grids, hence the
        2212             : !  need for a new map)
        2213             : !
        2214      345600 :   function ref2sphere(a,b, corners3D, ref_map, corners, facenum) result(sphere)
        2215             :     real(kind=r8)    :: a,b
        2216             :     type (spherical_polar_t)      :: sphere
        2217             :     type (cartesian3d_t)            :: corners3D(4)
        2218             :     integer :: ref_map
        2219             :     ! only needed for gnominic maps
        2220             :     type (cartesian2d_t), optional  :: corners(4)
        2221             :     integer, optional               :: facenum
        2222             :
        2223             :
        2224      345600 :     if (ref_map==0) then
        2225      345600 :        if (.not. present(corners) ) &
        2226           0 :             call endrun('ref2sphere(): missing arguments for equiangular map')
        2227      345600 :        sphere = ref2sphere_equiangular(a,b,corners,facenum)
        2228           0 :     elseif (ref_map==1) then
        2229           0 :        call endrun('gnomonic map not yet coded')
        2230           0 :     elseif (ref_map==2) then
        2231           0 :        sphere = ref2sphere_elementlocal(a,b,corners3D)
        2232             :     else
        2233           0 :        call endrun('ref2sphere(): bad value of ref_map')
        2234             :     endif
        2235      345600 :   end function ref2sphere
        2236             :
        2237             : !
        2238             : ! map a point in the referece element to the sphere
        2239             : !
        2240      345600 :   function ref2sphere_equiangular(a,b, corners, face_no) result(sphere)
        2241             :     implicit none
        2242             :     real(kind=r8)    :: a,b
        2243             :     integer,intent(in)            :: face_no
        2244             :     type (spherical_polar_t)      :: sphere
        2245             :     type (cartesian2d_t)          :: corners(4)
        2246             :     ! local
        2247             :     real(kind=r8)               :: pi,pj,qi,qj
        2248             :     type (cartesian2d_t)                 :: cart
        2249             :
        2250             :     ! map (a,b) to the [-pi/2,pi/2] equi angular cube face:  x1,x2
        2251             :     ! a = gp%points(i)
        2252             :     ! b = gp%points(j)
        2253      345600 :     pi = (1-a)/2
        2254      345600 :     pj = (1-b)/2
        2255      345600 :     qi = (1+a)/2
        2256      345600 :     qj = (1+b)/2
        2257             :     cart%x = pi*pj*corners(1)%x &
        2258             :          + qi*pj*corners(2)%x &
        2259             :          + qi*qj*corners(3)%x &
        2260      345600 :          + pi*qj*corners(4)%x
        2261             :     cart%y = pi*pj*corners(1)%y &
        2262             :          + qi*pj*corners(2)%y &
        2263             :          + qi*qj*corners(3)%y &
        2264      345600 :          + pi*qj*corners(4)%y
        2265             :     ! map from [pi/2,pi/2] equ angular cube face to sphere:
        2266      345600 :     sphere=projectpoint(cart,face_no)
        2267             :
        2268      345600 :   end function ref2sphere_equiangular
        2269             :
        2270             : !-----------------------------------------------------------------------------------------
        2271             : ! ELEMENT LOCAL MAP (DOES NOT USE CUBE FACES)
        2272             : ! unlike gnomonic equiangular map, this map will map all straight lines to
        2273             : ! great circle arcs
        2274             : !
        2275             : ! map a point in the referece element to the quad on the sphere by a
        2276             : ! general map, without using faces the map works this way: first, fix
        2277             : ! a coordinate (say, X). Map 4 corners of the ref element (corners are
        2278             : ! (-1,-1),(-1,1),(1,1), and (1,-1)) into 4 X-components of the quad in
        2279             : ! physical space via a bilinear map. Do so for Y and Z components as
        2280             : ! well. It produces a map: Ref element (\xi, \eta) ---> A quad in XYZ
        2281             : ! (ess, a piece of a twisted plane) with vertices of our target quad.  though
        2282             : ! the quad lies in a plane and not on the sphere manifold, its
        2283             : ! vertices belong to the sphere (by initial conditions). The last step
        2284             : ! is to utilize a map (X,Y,X) --> (X,Y,Z)/SQRT(X**2+Y**2+Z**2) to
        2285             : ! project the quad to the unit sphere.
        2286             : ! -----------------------------------------------------------------------------------------
        2287           0 :   function ref2sphere_elementlocal(a,b, corners3D) result(sphere)
        2288             :     use element_mod, only : element_t
        2289             :     implicit none
        2290             :     real(kind=r8)    :: a,b
        2291             :     type (cartesian3d_t)          :: corners3D(4)
        2292             :     type (spherical_polar_t)      :: sphere
        2293             :     real(kind=r8)               ::  q(4) ! local
        2294             :
        2295           0 :     q(1)=(1-a)*(1-b); q(2)=(1+a)*(1-b); q(3)=(1+a)*(1+b); q(4)=(1-a)*(1+b);
        2296           0 :     q=q/4.0_r8;
        2297           0 :     sphere = ref2sphere_elementlocal_q(q,corners3D)
        2298           0 :   end function ref2sphere_elementlocal
        2299             :
        2300           0 :   function ref2sphere_elementlocal_q(q, corners) result(sphere)
        2301             :     implicit none
        2302             :     real(kind=r8)          :: q(4)
        2303             :     type (spherical_polar_t)      :: sphere
        2304             :     type (cartesian3d_t)          :: corners(4)
        2305             :     ! local
        2306             :     type (cartesian3d_t)                 :: cart
        2307             :     real(kind=r8)               ::  c(3,4),  xx(3), r
        2308             :     integer :: i
        2309             :
        2310             : !3D corners fo the quad
        2311           0 :     c(1,1)=corners(1)%x;  c(2,1)=corners(1)%y;  c(3,1)=corners(1)%z;
        2312           0 :     c(1,2)=corners(2)%x;  c(2,2)=corners(2)%y;  c(3,2)=corners(2)%z;
        2313           0 :     c(1,3)=corners(3)%x;  c(2,3)=corners(3)%y;  c(3,3)=corners(3)%z;
        2314           0 :     c(1,4)=corners(4)%x;  c(2,4)=corners(4)%y;  c(3,4)=corners(4)%z;
        2315             :
        2316             : !physical point on a plane (sliced), not yet on the sphere
        2317           0 :     do i=1,3
        2318           0 :       xx(i)=sum(c(i,:)*q(:))
        2319             :     end do
        2320             :
        2321             : !distance from the plane point to the origin
        2322           0 :     r = sqrt(xx(1)**2+xx(2)**2+xx(3)**2)
        2323             :
        2324             : !projecting the plane point to the sphere
        2325           0 :     cart%x=xx(1)/r; cart%y=xx(2)/r; cart%z=xx(3)/r;
        2326             :
        2327             : !XYZ coords of the point to lon/lat
        2328           0 :     sphere=change_coordinates(cart)
        2329             :
        2330           0 :   end function ref2sphere_elementlocal_q
        2331             :
        2332           0 : end module cube_mod

Generated by: LCOV version 1.14