LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - element_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 27 | 65 | 41.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 15 | 13.3 %

              Line data    Source code
           1             : module element_mod
           2             :
           3             :   use shr_kind_mod,           only: r8=>shr_kind_r8, i8=>shr_kind_i8
           4             :   use coordinate_systems_mod, only: spherical_polar_t, cartesian2D_t, cartesian3D_t, distance
           5             :   use dimensions_mod,         only: np, nc, npsq, nlev, nlevp, qsize_d, max_neigh_edges,ntrac_d
           6             :   use edgetype_mod,           only: edgedescriptor_t
           7             :   use gridgraph_mod,          only: gridvertex_t
           8             :
           9             :   implicit none
          10             :   private
          11             :   integer, public, parameter :: timelevels = 3
          12             :
          13             :
          14             : ! =========== PRIMITIVE-EQUATION DATA-STRUCTURES =====================
          15             :
          16             :   type, public :: elem_state_t
          17             :
          18             :     ! prognostic variables for preqx solver
          19             :
          20             :     ! prognostics must match those in prim_restart_mod.F90
          21             :     ! vertically-lagrangian code advects dp3d instead of ps
          22             :     ! tracers Q, Qdp always use 2 level time scheme
          23             :
          24             :     real (kind=r8) :: v     (np,np,2,nlev,timelevels)            ! velocity
          25             :     real (kind=r8) :: T     (np,np,nlev,timelevels)              ! temperature
          26             :     real (kind=r8) :: dp3d  (np,np,nlev,timelevels)              ! dry delta p on levels
          27             :     real (kind=r8) :: psdry (np,np)                              ! dry surface pressure
          28             :     real (kind=r8) :: phis  (np,np)                              ! surface geopotential (prescribed)
          29             :     real (kind=r8), allocatable :: Qdp(:,:,:,:,:)                ! Tracer mass
          30             :   end type elem_state_t
          31             :
          32             :   !___________________________________________________________________
          33             :   type, public :: derived_state_t
          34             :      !
          35             :      ! storage for subcycling tracers/dynamics
          36             :      !
          37             :     real (kind=r8) :: vn0  (np,np,2,nlev)                      ! velocity for SE tracer advection
          38             :     real (kind=r8) :: dpdiss_biharmonic(np,np,nlev)            ! mean dp dissipation tendency, if nu_p>0
          39             :     real (kind=r8) :: dpdiss_ave(np,np,nlev)                   ! mean dp used to compute psdiss_tens
          40             :
          41             :     ! diagnostics for explicit timestep
          42             :     real (kind=r8) :: phi(np,np,nlev)                          ! geopotential
          43             :     real (kind=r8) :: omega(np,np,nlev)                        ! vertical velocity
          44             :
          45             :     ! tracer advection fields used for consistency and limiters
          46             :     real (kind=r8) :: dp(np,np,nlev)                           ! for dp_tracers at physics timestep
          47             :     real (kind=r8), allocatable :: divdp(:,:,:)                ! divergence of dp
          48             :     real (kind=r8), allocatable :: divdp_proj(:,:,:)           ! DSSed divdp
          49             :     real (kind=r8) :: mass(MAX(qsize_d,ntrac_d)+9)             ! total tracer mass for diagnostics
          50             :
          51             :     ! forcing terms for CAM
          52             :     real (kind=r8), allocatable :: FQ(:,:,:,:)                 ! tracer forcing
          53             :     real (kind=r8) :: FM(np,np,2,nlev)                         ! momentum forcing
          54             :     real (kind=r8), allocatable :: FDP(:,:,:)                  ! save full updated dp right after physics
          55             :     real (kind=r8) :: FT(np,np,nlev)                           ! temperature forcing
          56             :     real (kind=r8) :: etadot_prescribed(np,np,nlevp)           ! prescribed vertical tendency
          57             :     real (kind=r8) :: u_met(np,np,nlev)                        ! zonal component of prescribed meteorology winds
          58             :     real (kind=r8) :: dudt_met(np,np,nlev)                     ! rate of change of zonal component of prescribed meteorology winds
          59             :     real (kind=r8) :: v_met(np,np,nlev)                        ! meridional component of prescribed meteorology winds
          60             :     real (kind=r8) :: dvdt_met(np,np,nlev)                     ! rate of change of meridional component of prescribed meteorology winds
          61             :     real (kind=r8) :: T_met(np,np,nlev)                        ! prescribed meteorology temperature
          62             :     real (kind=r8) :: dTdt_met(np,np,nlev)                     ! rate of change of prescribed meteorology temperature
          63             :     real (kind=r8) :: ps_met(np,np)                            ! surface pressure of prescribed meteorology
          64             :     real (kind=r8) :: dpsdt_met(np,np)                         ! rate of change of surface pressure of prescribed meteorology
          65             :     real (kind=r8) :: nudge_factor(np,np,nlev)                 ! nudging factor (prescribed)
          66             :     real (kind=r8) :: Utnd(npsq,nlev)                          ! accumulated U tendency due to nudging towards prescribed met
          67             :     real (kind=r8) :: Vtnd(npsq,nlev)                          ! accumulated V tendency due to nudging towards prescribed met
          68             :     real (kind=r8) :: Ttnd(npsq,nlev)                          ! accumulated T tendency due to nudging towards prescribed met
          69             :
          70             :     real (kind=r8) :: pecnd(np,np,nlev)                        ! pressure perturbation from condensate
          71             :
          72             :     ! reference profiles
          73             :     real (kind=r8) :: T_ref(np,np,nlev)                        ! reference temperature
          74             :     real (kind=r8) :: dp_ref(np,np,nlev)                       ! reference pressure level thickness
          75             :   end type derived_state_t
          76             :
          77             :   !___________________________________________________________________
          78             :   type, public :: elem_accum_t
          79             :
          80             :
          81             :     ! the "4" timelevels represents data computed at:
          82             :     !  1  t-.5
          83             :     !  2  t+.5   after dynamics
          84             :     !  3  t+.5   after forcing
          85             :     !  4  t+.5   after Robert
          86             :     ! after calling TimeLevelUpdate, all times above decrease by 1.0
          87             :
          88             :
          89             :   end type elem_accum_t
          90             :
          91             :
          92             : ! ============= DATA-STRUCTURES COMMON TO ALL SOLVERS ================
          93             :
          94             :   type, public :: index_t
          95             :      integer :: ia(npsq),ja(npsq)
          96             :      integer :: is,ie
          97             :      integer :: NumUniquePts
          98             :      integer :: UniquePtOffset
          99             :   end type index_t
         100             :
         101             :   !___________________________________________________________________
         102             :   type, public :: element_t
         103             :      integer :: LocalId
         104             :      integer :: GlobalId
         105             :
         106             :      ! Coordinate values of element points
         107             :      type (spherical_polar_t) :: spherep(np,np)                       ! Spherical coords of GLL points
         108             :
         109             :      ! Equ-angular gnomonic projection coordinates
         110             :      type (cartesian2D_t)     :: cartp(np,np)                         ! gnomonic coords of GLL points
         111             :      type (cartesian2D_t)     :: corners(4)                           ! gnomonic coords of element corners
         112             :      real (kind=r8)    :: u2qmap(4,2)                          ! bilinear map from ref element to quad in cubedsphere coordinates
         113             :                                                                       ! SHOULD BE REMOVED
         114             :      ! 3D cartesian coordinates
         115             :      type (cartesian3D_t)     :: corners3D(4)
         116             :
         117             :      ! Element diagnostics
         118             :      real (kind=r8)    :: area                                 ! Area of element
         119             :      real (kind=r8)    :: normDinv                             ! some type of norm of Dinv used for CFL
         120             :      real (kind=r8)    :: dx_short                             ! short length scale in km
         121             :      real (kind=r8)    :: dx_long                              ! long length scale in km
         122             :
         123             :      real (kind=r8)    :: variable_hyperviscosity(np,np)       ! hyperviscosity based on above
         124             :      real (kind=r8)    :: hv_courant                           ! hyperviscosity courant number
         125             :      real (kind=r8)    :: tensorVisc(np,np,2,2)                !og, matrix V for tensor viscosity
         126             :
         127             :      ! Edge connectivity information
         128             : !     integer :: node_numbers(4)
         129             : !     integer :: node_multiplicity(4)                 ! number of elements sharing corner node
         130             :
         131             :      type (GridVertex_t)      :: vertex                               ! element grid vertex information
         132             :      type (EdgeDescriptor_t)  :: desc
         133             :
         134             :      type (elem_state_t)      :: state
         135             :
         136             :      type (derived_state_t)   :: derived
         137             :      ! Metric terms
         138             :      real (kind=r8)    :: met(np,np,2,2)                       ! metric tensor on velocity and pressure grid
         139             :      real (kind=r8)    :: metinv(np,np,2,2)                    ! metric tensor on velocity and pressure grid
         140             :      real (kind=r8)    :: metdet(np,np)                        ! g = SQRT(det(g_ij)) on velocity and pressure grid
         141             :      real (kind=r8)    :: rmetdet(np,np)                       ! 1/metdet on velocity pressure grid
         142             :      real (kind=r8)    :: D(np,np,2,2)                         ! Map covariant field on cube to vector field on the sphere
         143             :      real (kind=r8)    :: Dinv(np,np,2,2)                      ! Map vector field on the sphere to covariant v on cube
         144             :
         145             :
         146             :      ! Mass flux across the sides of each sub-element.
         147             :      ! The storage is redundent since the mass across shared sides
         148             :      ! must be equal in magnitude and opposite in sign.
         149             :      ! The layout is like:
         150             :      !   --------------------------------------------------------------
         151             :      ! ^|    (1,4,3)     |                |              |    (4,4,3) |
         152             :      ! ||                |                |              |            |
         153             :      ! ||(1,4,4)         |                |              |(4,4,4)     |
         154             :      ! ||         (1,4,2)|                |              |     (4,4,2)|
         155             :      ! ||                |                |              |            |
         156             :      ! ||   (1,4,1)      |                |              |  (4,4,1)   |
         157             :      ! |---------------------------------------------------------------
         158             :      ! S|                |                |              |            |
         159             :      ! e|                |                |              |            |
         160             :      ! c|                |                |              |            |
         161             :      ! o|                |                |              |            |
         162             :      ! n|                |                |              |            |
         163             :      ! d|                |                |              |            |
         164             :      !  ---------------------------------------------------------------
         165             :      ! C|                |                |              |            |
         166             :      ! o|                |                |              |            |
         167             :      ! o|                |                |              |            |
         168             :      ! r|                |                |              |            |
         169             :      ! d|                |                |              |            |
         170             :      ! i|                |                |              |            |
         171             :      ! n---------------------------------------------------------------
         172             :      ! a|    (1,1,3)     |                |              |    (4,1,3) |
         173             :      ! t|                |                |              |(4,1,4)     |
         174             :      ! e|(1,1,4)         |                |              |            |
         175             :      !  |         (1,1,2)|                |              |     (4,1,2)|
         176             :      !  |                |                |              |            |
         177             :      !  |    (1,1,1)     |                |              |  (4,1,1)   |
         178             :      !  ---------------------------------------------------------------
         179             :      !          First Coordinate ------->
         180             :      real (kind=r8) :: sub_elem_mass_flux(nc,nc,4,nlev)
         181             :
         182             :      ! Convert vector fields from spherical to rectangular components
         183             :      ! The transpose of this operation is its pseudoinverse.
         184             :      real (kind=r8)    :: vec_sphere2cart(np,np,3,2)
         185             :
         186             :      ! Mass matrix terms for an element on a cube face
         187             :      real (kind=r8)    :: mp(np,np)                            ! mass matrix on v and p grid
         188             :      real (kind=r8)    :: rmp(np,np)                           ! inverse mass matrix on v and p grid
         189             :
         190             :      ! Mass matrix terms for an element on the sphere
         191             :      ! This mass matrix is used when solving the equations in weak form
         192             :      ! with the natural (surface area of the sphere) inner product
         193             :      real (kind=r8)    :: spheremp(np,np)                      ! mass matrix on v and p grid
         194             :      real (kind=r8)    :: rspheremp(np,np)                     ! inverse mass matrix on v and p grid
         195             :
         196             :      integer(i8) :: gdofP(np,np)                     ! global degree of freedom (P-grid)
         197             :
         198             :      real (kind=r8)    :: fcor(np,np)                          ! Coreolis term
         199             :
         200             :      type (index_t) :: idxP
         201             :      type (index_t),pointer :: idxV
         202             :      integer :: FaceNum
         203             :
         204             :      ! force element_t to be a multiple of 8 bytes.
         205             :      ! on BGP, code will crash (signal 7, or signal 15) if 8 byte alignment is off
         206             :      ! check core file for:
         207             :      ! core.63:Generated by interrupt..(Alignment Exception DEAR=0xa1ef671c ESR=0x01800000 CCR0=0x4800a002)
         208             :      integer :: dummy
         209             :   end type element_t
         210             :
         211             :   !___________________________________________________________________
         212             :   public :: element_coordinates
         213             :   public :: element_var_coordinates
         214             :   public :: element_var_coordinates3D
         215             :   public :: GetColumnIdP,GetColumnIdV
         216             :   public :: allocate_element_desc
         217             :   public :: PrintElem
         218             :
         219             : contains
         220             :
         221           0 :   subroutine PrintElem(arr)
         222             :
         223             :     real(kind=r8) :: arr(:,:)
         224             :     integer :: i,j
         225             :
         226           0 :       do j=np,1,-1
         227           0 :          write(6,*) (arr(i,j), i=1,np)
         228             :       enddo
         229             :
         230           0 :   end subroutine PrintElem
         231             : ! ===================== ELEMENT_MOD METHODS ==========================
         232             :
         233           0 :   function GetColumnIdP(elem,i,j) result(col_id)
         234             :
         235             :     ! Get unique identifier for a Physics column on the P-grid
         236             :
         237             :     type(element_t), intent(in) :: elem
         238             :     integer, intent(in) :: i,j
         239             :     integer :: col_id
         240           0 :     col_id = elem%gdofP(i,j)
         241           0 :   end function GetColumnIdP
         242             :
         243             :   !___________________________________________________________________
         244           0 :   function GetColumnIdV(elem,i,j) result(col_id)
         245             :
         246             :     !  Get unique identifier for a Physics column on the V-grid
         247             :
         248             :     type(element_t), intent(in) :: elem
         249             :     integer, intent(in) :: i,j
         250             :     integer :: col_id
         251           0 :     col_id = elem%gdofP(i,j)
         252           0 :   end function GetColumnIdV
         253             :
         254             :   !___________________________________________________________________
         255           0 :   function element_coordinates(start,end,points) result(cart)
         256             :
         257             :     ! Initialize 2D rectilinear element colocation points
         258             :
         259             :     type (cartesian2D_t), intent(in) :: start
         260             :     type (cartesian2D_t), intent(in) :: end
         261             :     real(r8),             intent(in) :: points(:)
         262             :     type (cartesian2D_t)             :: cart(SIZE(points),SIZE(points))
         263             :
         264             :     type (cartesian2D_t) :: length, centroid
         265             :     real(r8)             :: y
         266             :     integer              :: i,j
         267             :
         268           0 :     length%x   = 0.50D0*(end%x-start%x)
         269           0 :     length%y   = 0.50D0*(end%y-start%y)
         270           0 :     centroid%x = 0.50D0*(end%x+start%x)
         271           0 :     centroid%y = 0.50D0*(end%y+start%y)
         272           0 :     do j=1,SIZE(points)
         273           0 :        y = centroid%y + length%y*points(j)
         274           0 :        do i=1,SIZE(points)
         275           0 :           cart(i,j)%x = centroid%x + length%x*points(i)
         276           0 :           cart(i,j)%y = y
         277             :        end do
         278             :     end do
         279           0 :   end function element_coordinates
         280             :
         281             :   !___________________________________________________________________
         282       21600 :   function element_var_coordinates(c,points) result(cart)
         283             :
         284             :     type (cartesian2D_t), intent(in) :: c(4)
         285             :     real(r8),             intent(in) :: points(:)
         286             :     type (cartesian2D_t)             :: cart(SIZE(points),SIZE(points))
         287             :
         288       43200 :     real(r8) :: p(size(points))
         289       43200 :     real(r8) :: q(size(points))
         290             :     integer  :: i,j
         291             :
         292      129600 :     p(:) = (1.0D0-points(:))/2.0D0
         293      108000 :     q(:) = (1.0D0+points(:))/2.0D0
         294             :
         295      108000 :     do j=1,SIZE(points)
         296      453600 :        do i=1,SIZE(points)
         297      691200 :           cart(i,j)%x = p(i)*p(j)*c(1)%x &
         298             :                       + q(i)*p(j)*c(2)%x &
         299             :                       + q(i)*q(j)*c(3)%x &
         300      691200 :                       + p(i)*q(j)*c(4)%x
         301             :           cart(i,j)%y = p(i)*p(j)*c(1)%y &
         302             :                       + q(i)*p(j)*c(2)%y &
         303             :                       + q(i)*q(j)*c(3)%y &
         304      432000 :                       + p(i)*q(j)*c(4)%y
         305             :        end do
         306             :     end do
         307       21600 :   end function element_var_coordinates
         308             :
         309             :   !___________________________________________________________________
         310           0 :   function element_var_coordinates3d(c,points) result(cart)
         311             :
         312             :     type(cartesian3D_t), intent(in) :: c(4)
         313             :     real(r8), intent(in) :: points(:)
         314             :
         315             :     type(cartesian3D_t)  :: cart(SIZE(points),SIZE(points))
         316             :
         317           0 :     real(r8) :: p(size(points))
         318           0 :     real(r8) :: q(size(points)), r
         319             :     integer  :: i,j
         320             :
         321           0 :     p(:) = (1.0D0-points(:))/2.0D0
         322           0 :     q(:) = (1.0D0+points(:))/2.0D0
         323             :
         324           0 :     do j=1,SIZE(points)
         325           0 :        do i=1,SIZE(points)
         326           0 :           cart(i,j)%x = p(i)*p(j)*c(1)%x &
         327             :                       + q(i)*p(j)*c(2)%x &
         328             :                       + q(i)*q(j)*c(3)%x &
         329           0 :                       + p(i)*q(j)*c(4)%x
         330             :           cart(i,j)%y = p(i)*p(j)*c(1)%y &
         331             :                       + q(i)*p(j)*c(2)%y &
         332             :                       + q(i)*q(j)*c(3)%y &
         333           0 :                       + p(i)*q(j)*c(4)%y
         334             :           cart(i,j)%z = p(i)*p(j)*c(1)%z &
         335             :                       + q(i)*p(j)*c(2)%z &
         336             :                       + q(i)*q(j)*c(3)%z &
         337           0 :                       + p(i)*q(j)*c(4)%z
         338             :
         339             :           ! project back to sphere:
         340           0 :           r = distance(cart(i,j))
         341           0 :           cart(i,j)%x = cart(i,j)%x/r
         342           0 :           cart(i,j)%y = cart(i,j)%y/r
         343           0 :           cart(i,j)%z = cart(i,j)%z/r
         344             :        end do
         345             :     end do
         346           0 :   end function element_var_coordinates3d
         347             :
         348             :   !___________________________________________________________________
         349        1536 :   subroutine allocate_element_desc(elem)
         350             :
         351             :     type (element_t), intent(inout)   :: elem(:)
         352             :     integer                           :: num, j,i
         353             :
         354        1536 :     num = SIZE(elem)
         355             :
         356       12336 :     do j=1,num
         357       32400 :        allocate(elem(j)%desc%putmapP(max_neigh_edges))
         358       21600 :        allocate(elem(j)%desc%getmapP(max_neigh_edges))
         359       21600 :        allocate(elem(j)%desc%putmapP_ghost(max_neigh_edges))
         360       21600 :        allocate(elem(j)%desc%getmapP_ghost(max_neigh_edges))
         361       21600 :        allocate(elem(j)%desc%putmapS(max_neigh_edges))
         362       21600 :        allocate(elem(j)%desc%getmapS(max_neigh_edges))
         363       21600 :        allocate(elem(j)%desc%reverse(max_neigh_edges))
         364       21600 :        allocate(elem(j)%desc%globalID(max_neigh_edges))
         365       21600 :        allocate(elem(j)%desc%loc2buf(max_neigh_edges))
         366       98736 :        do i=1,max_neigh_edges
         367       86400 :           elem(j)%desc%loc2buf(i)=i
         368       97200 :           elem(j)%desc%globalID(i)=-1
         369             :        enddo
         370             :
         371             :     end do
         372        1536 :   end subroutine allocate_element_desc
         373             :
         374             :
         375           0 : end module element_mod

Generated by: LCOV version 1.14