LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - fvm_control_volume_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 62 | 68 | 91.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 3 | 5 | 60.0 %

              Line data    Source code
           1             : #ifdef HAVE_CONFIG_H
           2             : #include "config.h"
           3             : #endif
           4             :
           5             : !MODULE FVM_CONTROL_VOLUME_MOD---------------------------------------------CE-for FVM
           6             : ! AUTHOR: Christoph Erath, 11.June 2011                                             !
           7             : ! This module contains everything to initialize the arrival. It also provides the   !
           8             : ! interpolation points for the reconstruction (projection from one face to another  !
           9             : ! when the element is on the cube edge)                                             !
          10             : ! It also intialize the start values, see also fvm_analytic                         !
          11             : !-----------------------------------------------------------------------------------!
          12             : module fvm_control_volume_mod
          13             :   use shr_kind_mod,           only: r8=>shr_kind_r8
          14             :   use coordinate_systems_mod, only: spherical_polar_t
          15             :   use element_mod,            only: element_t
          16             :   use dimensions_mod,         only: nc, nhe, nlev, ntrac_d, qsize_d,ne, np, nhr, ns, nhc
          17             :   use dimensions_mod,         only: fv_nphys, nhe_phys, nhr_phys, ns_phys, nhc_phys,fv_nphys
          18             :   use dimensions_mod,         only: irecons_tracer
          19             :   use cam_abortutils,         only: endrun
          20             :
          21             :   implicit none
          22             :   private
          23             :   integer, parameter, private:: nh = nhr+(nhe-1) ! = 2 (nhr=2; nhe=1)
          24             :                                                  ! = 3 (nhr=2; nhe=2)
          25             :
          26             :   type, public :: fvm_struct
          27             :     ! fvm tracer mixing ratio: (kg/kg)
          28             :     real (kind=r8) :: c(1-nhc:nc+nhc,1-nhc:nc+nhc,nlev,ntrac_d)
          29             :     real (kind=r8) :: se_flux(1-nhe:nc+nhe,1-nhe:nc+nhe,4,nlev)
          30             :
          31             :     real (kind=r8) :: dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,nlev)
          32             :     real (kind=r8) :: dp_ref(nlev)
          33             :     real (kind=r8) :: dp_ref_inverse(nlev)
          34             :     real (kind=r8) :: psc(nc,nc)
          35             :
          36             :     real (kind=r8) :: inv_area_sphere(nc,nc)    ! inverse area_sphere
          37             :     real (kind=r8) :: inv_se_area_sphere(nc,nc) ! inverse area_sphere
          38             :
          39             :     integer                  :: faceno         !face number
          40             :     ! number of south,....,swest and 0 for interior element
          41             :     integer                  :: cubeboundary
          42             :
          43             : #ifdef waccm_debug
          44             :     real (kind=r8) :: CSLAM_gamma(nc,nc,nlev,4)
          45             : #endif
          46             :     real (kind=r8) :: displ_max(1-nhc:nc+nhc,1-nhc:nc+nhc,4)
          47             :     integer        :: flux_vec (2,1-nhc:nc+nhc,1-nhc:nc+nhc,4)
          48             :     !
          49             :     !
          50             :     ! cartesian location of vertices for flux sides
          51             :     !
          52             :     !  x-coordinate of vertex 1: vtx_cart(1,1i,j,1,1)  = fvm%acartx(i)
          53             :     !  y-coordinate of vertex 1: vtx_cart(1,2,i,j,2,1) = fvm%acarty(j)
          54             :     !
          55             :     !  x-coordinate of vertex 2: vtx_cart(2,1,i,j) = fvm%acartx(i+1)
          56             :     !  y-coordinate of vertex 2: vtx_cart(2,2,i,j) = fvm%acarty(j  )
          57             :     !
          58             :     !  x-coordinate of vertex 3: vtx_cart(3,1,i,j) = fvm%acartx(i+1)
          59             :     !  y-coordinate of vertex 3: vtx_cart(3,2,i,j) = fvm%acarty(j+1)
          60             :     !
          61             :     !  x-coordinate of vertex 4: vtx_cart(4,1,i,j) = fvm%acartx(i  )
          62             :     !  y-coordinate of vertex 4: vtx_cart(4,2,i,j) = fvm%acarty(j+1)
          63             :     !
          64             :     real (kind=r8) :: vtx_cart (4,2,1-nhc:nc+nhc,1-nhc:nc+nhc)
          65             :     !
          66             :     ! flux_orient(1,i,j) = panel on which control volume (i,j) is located
          67             :     ! flux_orient(2,i,j) = cshift value for vertex permutation
          68             :     !
          69             :     real (kind=r8) :: flux_orient(2  ,1-nhc:nc+nhc,1-nhc:nc+nhc)
          70             :     !
          71             :     ! i,j: indicator function for non-existent cells (0 for corner halo and 1 elsewhere)
          72             :     !
          73             :     integer                  :: ifct   (1-nhc:nc+nhc,1-nhc:nc+nhc)
          74             :     integer                  :: rot_matrix(2,2,1-nhc:nc+nhc,1-nhc:nc+nhc)
          75             :     !
          76             :     real (kind=r8)           :: dalpha, dbeta             ! central-angle for gnomonic coordinates
          77             :     type (spherical_polar_t) :: center_cart(nc,nc)        ! center of fvm cell in gnomonic coordinates
          78             :     real (kind=r8)           :: area_sphere(nc,nc)        ! spherical area of fvm cell
          79             :     real (kind=r8)           :: spherecentroid(irecons_tracer-1,1-nhc:nc+nhc,1-nhc:nc+nhc) ! centroids
          80             :     !
          81             :     ! pre-computed metric terms (for efficiency)
          82             :     !
          83             :     ! recons_metrics(1,:,:) = spherecentroid(1,:,:)**2 -spherecentroid(3,:,:)
          84             :     ! recons_metrics(2,:,:) = spherecentroid(2,:,:)**2 -spherecentroid(4,:,:)
          85             :     ! recons_metrics(3,:,:) = spherecentroid(1,:,:)*spherecentroid(2,:,:)-spherecentroid(5,:,:)
          86             :     !
          87             :     real (kind=r8)           :: recons_metrics(3,1-nhe:nc+nhe,1-nhe:nc+nhe)
          88             :     !
          89             :     ! recons_metrics_integral(1,:,:) = 2.0_r8*spherecentroid(1,:,:)**2 -spherecentroid(3,:,:)
          90             :     ! recons_metrics_integral(2,:,:) = 2.0_r8*spherecentroid(2,:,:)**2 -spherecentroid(4,:,:)
          91             :     ! recons_metrics_integral(3,:,:) = 2.0_r8*spherecentroid(1,:,:)*spherecentroid(2,:,:)-spherecentroid(5,:,:)
          92             :     !
          93             :     real (kind=r8)           :: recons_metrics_integral(3,1-nhe:nc+nhe,1-nhe:nc+nhe)
          94             :     !
          95             :     integer                  :: jx_min(3), jx_max(3), jy_min(3), jy_max(3) !bounds for computation
          96             :
          97             :     ! provide fixed interpolation points with respect to the arrival grid for
          98             :     ! reconstruction
          99             :     integer                  :: ibase(1-nh:nc+nh,1:nhr,2)
         100             :     real (kind=r8)           :: halo_interp_weight(1:ns,1-nh:nc+nh,1:nhr,2)
         101             :     real (kind=r8)           :: centroid_stretch(7,1-nhe:nc+nhe,1-nhe:nc+nhe) !for finite-difference reconstruction
         102             :     !
         103             :     ! pre-compute weights for reconstruction at cell vertices
         104             :     !
         105             :     !  ! Evaluate constant order terms
         106             :     !  value = fcube(a,b) + &
         107             :     !  ! Evaluate linear order terms
         108             :     !          recons(1,a,b) * (cartx - centroid(1,a,b)) + &
         109             :     !          recons(2,a,b) * (carty - centroid(2,a,b)) + &
         110             :     !  ! Evaluate second order terms
         111             :     !          recons(3,a,b) * (centroid(1,a,b)**2 - centroid(3,a,b)) + &
         112             :     !          recons(4,a,b) * (centroid(2,a,b)**2 - centroid(4,a,b)) + &
         113             :     !          recons(5,a,b) * (centroid(1,a,b) * centroid(2,a,b) - centroid(5,a,b)) + &
         114             :     !
         115             :     !          recons(3,a,b) * (cartx - centroid(1,a,b))**2 + &
         116             :     !          recons(4,a,b) * (carty - centroid(2,a,b))**2 + &
         117             :     !          recons(5,a,b) * (cartx - centroid(1,a,b)) * (carty - centroid(2,a,b))
         118             :     !
         119             :     real (kind=r8)    :: vertex_recons_weights(4,1:irecons_tracer-1,1-nhe:nc+nhe,1-nhe:nc+nhe)
         120             :     !
         121             :     ! for mapping fvm2dyn
         122             :     !
         123             :     real (kind=r8)    :: norm_elem_coord(2,1-nhc:nc+nhc,1-nhc:nc+nhc)
         124             :     !
         125             :     !******************************************
         126             :     !
         127             :     ! separate physics grid variables
         128             :     !
         129             :     !******************************************
         130             :     !
         131             :     real (kind=r8)           , allocatable :: vtx_cart_physgrid(:,:,:,:)
         132             :     real (kind=r8)           , allocatable :: flux_orient_physgrid(:,:,:)
         133             :     integer                  , allocatable :: ifct_physgrid(:,:)
         134             :     integer                  , allocatable :: rot_matrix_physgrid(:,:,:,:)
         135             :     real (kind=r8)           , allocatable :: spherecentroid_physgrid(:,:,:)
         136             :     real (kind=r8)           , allocatable :: recons_metrics_physgrid(:,:,:)
         137             :     real (kind=r8)           , allocatable :: recons_metrics_integral_physgrid(:,:,:)
         138             :     ! centroid_stretch_physgrid for finite-difference reconstruction
         139             :     real (kind=r8)           , allocatable :: centroid_stretch_physgrid       (:,:,:)
         140             :     real (kind=r8)                         :: dalpha_physgrid, dbeta_physgrid             ! central-angle for gnomonic coordinates
         141             :     type (spherical_polar_t) , allocatable :: center_cart_physgrid(:,:)        ! center of fvm cell in gnomonic coordinates
         142             :     real (kind=r8)           , allocatable :: area_sphere_physgrid(:,:)        ! spherical area of fvm cell
         143             :     integer                                :: jx_min_physgrid(3), jx_max_physgrid(3) !bounds for computation
         144             :     integer                                :: jy_min_physgrid(3), jy_max_physgrid(3) !bounds for computation
         145             :     integer                  , allocatable :: ibase_physgrid(:,:,:)
         146             :     real (kind=r8)           , allocatable :: halo_interp_weight_physgrid(:,:,:,:)
         147             :     real (kind=r8)           , allocatable :: vertex_recons_weights_physgrid(:,:,:,:)
         148             :
         149             :     real (kind=r8)           , allocatable :: norm_elem_coord_physgrid(:,:,:)
         150             :     real (kind=r8)           , allocatable :: Dinv_physgrid(:,:,:,:)
         151             :
         152             :     real (kind=r8)           , allocatable :: fc(:,:,:,:)
         153             :     real (kind=r8)           , allocatable :: fc_phys(:,:,:,:)
         154             :     real (kind=r8)           , allocatable :: ft(:,:,:)
         155             :     real (kind=r8)           , allocatable :: fm(:,:,:,:)
         156             :     real (kind=r8)           , allocatable :: dp_phys(:,:,:)
         157             :   end type fvm_struct
         158             :
         159             :   public :: fvm_mesh, fvm_set_cubeboundary, allocate_physgrid_vars
         160             :
         161             :
         162             :   real (kind=r8),parameter, public   :: bignum = 1.0E20_r8
         163             :
         164             : contains
         165       10800 :   subroutine fvm_set_cubeboundary(elem, fvm)
         166             :     implicit none
         167             :     type (element_t) , intent(in)      :: elem
         168             :     type (fvm_struct), intent(inout)   :: fvm
         169             :
         170             :     logical                            :: corner
         171             :     integer                            :: j, mynbr_cnt, mystart
         172             :     integer                            :: nbrsface(8)! store the neighbours in north, south
         173             :
         174       10800 :     fvm%faceno=elem%FaceNum
         175             :     ! write the neighbors in the structure
         176       10800 :     fvm%cubeboundary=0
         177       10800 :     corner=.FALSE.
         178       97200 :     do j=1,8
         179       86400 :        mynbr_cnt = elem%vertex%nbrs_ptr(j+1) - elem%vertex%nbrs_ptr(j) !length of neighbor location
         180       86400 :        mystart = elem%vertex%nbrs_ptr(j)
         181             :        !NOTE: assuming that we do not have multiple corner neighbors (so not a refined mesh)
         182       97200 :        if (mynbr_cnt > 0 ) then
         183       86352 :           nbrsface(j)=elem%vertex%nbrs_face(mystart)
         184             :           ! note that if the element lies on a corner, it will be at j=5,6,7,8
         185       86352 :           if ((nbrsface(j) /= fvm%faceno) .AND. (j<5)) then
         186        1440 :              fvm%cubeboundary=j
         187             :           endif
         188             :        else   ! corner on the cube
         189          48 :           if (.NOT. corner) then
         190          48 :              nbrsface(j)=-1
         191          48 :              fvm%cubeboundary=j
         192          48 :              corner=.TRUE.
         193             :           else
         194           0 :              if ( ne == 0 ) then
         195             :                 ! dont check this condition.  note that we call this code
         196             :                 ! generate phys grid template files, so we need to be able
         197             :                 ! to call create_ari() to create the subcells even though
         198             :                 ! cslam cant run with the unstructed ne=0 case
         199             :              else
         200           0 :                 print *,'Error in fvm_CONTROL_VOLUME_MOD - Subroutine fvm_MESH_ARI: '
         201           0 :                 call endrun('Do not allow one element per face for fvm, please increase ne!')
         202             :              endif
         203             :           endif
         204             :        end if
         205             :     end do
         206       10800 :   end subroutine fvm_set_cubeboundary
         207             :
         208       10800 :   subroutine fvm_mesh(elem, fvm)
         209             :     use fvm_analytic_mod, only : compute_halo_vars
         210             :     use fvm_analytic_mod, only : create_interpolation_points
         211             :     use derivative_mod  , only : subcell_integration
         212             :
         213             :     implicit none
         214             :     type (element_t), intent(in)     :: elem
         215             :     type (fvm_struct), intent(inout) :: fvm
         216             :     integer :: i,j
         217             :     real (kind=r8)            :: tmp(np,np)
         218             :     !
         219             :     ! initialize metric and related terms on panel
         220             :     !
         221             :     call compute_halo_vars(&    !input
         222             :          fvm%faceno,fvm%cubeboundary,nc,nhc,nhe,   &  !input
         223             :          fvm%jx_min,fvm%jx_max,fvm%jy_min,fvm%jy_max,&!output
         224       10800 :          fvm%flux_orient,fvm%ifct,fvm%rot_matrix)     !output
         225       43200 :     do j=1,nc
         226      140400 :       do i=1,nc
         227       97200 :         fvm%norm_elem_coord(1,i,j) = elem%corners(1)%x+(i-0.5_r8)*fvm%dalpha
         228      129600 :         fvm%norm_elem_coord(2,i,j) = elem%corners(1)%y+(j-0.5_r8)*fvm%dalpha
         229             :       end do
         230             :     end do
         231             :
         232             :     !
         233             :     ! overwrite areas for consistency with SE areas (that are O(10E-5) incorrect)
         234             :     !
         235             : !    tmp = 1.0_r8
         236             : !    call subcell_integration(tmp, np, nc, elem%metdet,fvm%area_sphere)
         237             :     !
         238             :     ! do the same for physics grid
         239             :     !
         240             :     call compute_halo_vars(&
         241             :          fvm%faceno,fvm%cubeboundary,fv_nphys,nhc_phys,nhe_phys,&
         242             :          fvm%jx_min_physgrid,fvm%jx_max_physgrid,fvm%jy_min_physgrid,fvm%jy_max_physgrid,&
         243       10800 :          fvm%flux_orient_physgrid,fvm%ifct_physgrid,fvm%rot_matrix_physgrid)
         244       43200 :     do j=1,fv_nphys
         245      140400 :       do i=1,fv_nphys
         246       97200 :         fvm%norm_elem_coord_physgrid(1,i,j) = elem%corners(1)%x+(i-0.5_r8)*fvm%dalpha_physgrid
         247      129600 :         fvm%norm_elem_coord_physgrid(2,i,j) = elem%corners(1)%y+(j-0.5_r8)*fvm%dalpha_physgrid
         248             :       end do
         249             :     end do
         250             :     !
         251             :     ! initialize halo interpolation variables
         252             :     !
         253             :     call create_interpolation_points(elem,&
         254             :          nc,nhc,nhr,ns,nh,fvm%cubeboundary,&
         255       10800 :          fvm%dalpha,fvm%dbeta,fvm%ibase,fvm%halo_interp_weight)
         256             :     call create_interpolation_points(elem,&
         257             :          fv_nphys,nhc_phys,nhr_phys,ns_phys,nhr_phys,fvm%cubeboundary,&
         258       10800 :          fvm%dalpha_physgrid,fvm%dbeta_physgrid,fvm%ibase_physgrid,fvm%halo_interp_weight_physgrid)
         259       10800 :   end subroutine fvm_mesh
         260             :
         261             :
         262        1536 :   subroutine allocate_physgrid_vars(fvm,par)
         263       10800 :     use cam_logfile   , only : iulog
         264             :     use parallel_mod  , only : parallel_t
         265             :     use dimensions_mod, only : nelemd
         266             :     type (fvm_struct), intent(inout) :: fvm(:)
         267             :     type (parallel_t), intent(in)    :: par
         268             :     integer :: ie
         269             :
         270        1536 :     nhc_phys = fv_nphys
         271        1536 :     nhe_phys = 0
         272        1536 :     nhr_phys = 2
         273        1536 :     ns_phys  = MAX(fv_nphys,2)
         274             :
         275        1536 :     if(par%masterproc) then
         276           2 :       write(iulog,*)"allocating physgrid grid vars"
         277           2 :       write(iulog,*)"fv_nphys,nhc_phys,nhe_phys,nhr_phys,ns_phys = ",&
         278           4 :            fv_nphys,nhc_phys,nhe_phys,nhr_phys,ns_phys
         279             :     end if
         280             :
         281       12336 :     do ie=1,nelemd
         282       43200 :       allocate(fvm(ie)%vtx_cart_physgrid      (4,2,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
         283       43200 :       allocate(fvm(ie)%flux_orient_physgrid   (2,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
         284       43200 :       allocate(fvm(ie)%ifct_physgrid         (1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
         285       43200 :       allocate(fvm(ie)%rot_matrix_physgrid   (2,2,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
         286             :
         287           0 :       allocate(fvm(ie)%spherecentroid_physgrid(irecons_tracer-1,&
         288       43200 :            1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
         289       43200 :       allocate(fvm(ie)%recons_metrics_physgrid         (3,1-nhe_phys:fv_nphys+nhe_phys,1-nhe_phys:fv_nphys+nhe_phys))
         290       32400 :       allocate(fvm(ie)%recons_metrics_integral_physgrid(3,1-nhe_phys:fv_nphys+nhe_phys,1-nhe_phys:fv_nphys+nhe_phys))
         291       43200 :       allocate(fvm(ie)%centroid_stretch_physgrid       (7,1-nhe_phys:fv_nphys+nhe_phys,1-nhe_phys:fv_nphys+nhe_phys))
         292       43200 :       allocate(fvm(ie)%center_cart_physgrid(fv_nphys,fv_nphys))
         293       43200 :       allocate(fvm(ie)%area_sphere_physgrid(fv_nphys,fv_nphys))
         294       54000 :       allocate(fvm(ie)%ibase_physgrid(1-nhr_phys:fv_nphys+nhr_phys,1:nhr_phys,2))
         295       64800 :       allocate(fvm(ie)%halo_interp_weight_physgrid(1:ns_phys,1-nhr_phys:fv_nphys+nhr_phys,1:nhr_phys,2))
         296           0 :       allocate(fvm(ie)%vertex_recons_weights_physgrid(4,1:irecons_tracer-1,1-nhe_phys:fv_nphys+nhe_phys,&
         297       43200 :             1-nhe_phys:fv_nphys+nhe_phys))
         298             :
         299       32400 :       allocate(fvm(ie)%norm_elem_coord_physgrid(2,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys    ))
         300       64800 :       allocate(fvm(ie)%Dinv_physgrid           (  1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,2,2))
         301             :
         302       10800 :       allocate(fvm(ie)%fc(nc,nc,nlev,max(ntrac_d,qsize_d)))
         303       64800 :       allocate(fvm(ie)%fc_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev,max(ntrac_d,qsize_d)))
         304       43200 :       allocate(fvm(ie)%ft(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev))
         305       54000 :       allocate(fvm(ie)%fm(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,2,nlev))
         306       33936 :       allocate(fvm(ie)%dp_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev))
         307             :     end do
         308        1536 :   end subroutine allocate_physgrid_vars
         309           0 : end module fvm_control_volume_mod

Generated by: LCOV version 1.14