scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve Module Reference

module FElib / Fluid dyn solver / Atmosphere / Global nonhydrostatic model / HEVE More...

Functions/Subroutines
subroutine, public	atm_dyn_dgm_globalnonhydro3d_rhot_heve_init (mesh)
subroutine, public	atm_dyn_dgm_globalnonhydro3d_rhot_heve_final ()
subroutine, public	atm_dyn_dgm_globalnh3d_rhot_heve_cal_tend_shallow_atm_asis (dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, pres_hyd_ref, coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, element3d_operation, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)
subroutine, public	atm_dyn_dgm_globalnonhydro3d_rhot_heve_cal_tend_shallow_atm (dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, pres_hyd_ref, coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, element3d_operation, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)
subroutine, public	atm_dyn_dgm_globalnonhydro3d_rhot_heve_cal_tend_deep_atm (dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, pres_hyd_ref, coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, element3d_operation, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)

Detailed Description

module FElib / Fluid dyn solver / Atmosphere / Global nonhydrostatic model / HEVE

Description

HEVE DGM scheme for Global Atmospheric Dynamical process. The governing equations is a fully compressible nonhydrostatic equations, which consist of mass, momentum, and thermodynamics (density * potential temperature conservation) equations.

Author

Yuta Kawai, Team SCALE

Function/Subroutine Documentation

atm_dyn_dgm_globalnonhydro3d_rhot_heve_init()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve::atm_dyn_dgm_globalnonhydro3d_rhot_heve_init

(

class(meshbase3d), intent(in)

mesh

)

Definition at line 76 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve.F90.

    implicit none
    class(MeshBase3D), intent(in) :: mesh
    !--------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_init( mesh )
 
    return

References scale_atm_dyn_dgm_nonhydro3d_common::atm_dyn_dgm_nonhydro3d_common_init().

atm_dyn_dgm_globalnonhydro3d_rhot_heve_final()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve::atm_dyn_dgm_globalnonhydro3d_rhot_heve_final

Definition at line 87 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve.F90.

    implicit none
    !--------------------------------------------
    
    call atm_dyn_dgm_nonhydro3d_common_final()    
    return

References scale_atm_dyn_dgm_nonhydro3d_common::atm_dyn_dgm_nonhydro3d_common_final().

atm_dyn_dgm_globalnh3d_rhot_heve_cal_tend_shallow_atm_asis()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve::atm_dyn_dgm_globalnh3d_rhot_heve_cal_tend_shallow_atm_asis	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	dens_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momx_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momy_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momz_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	rhot_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	ddens_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momx_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momy_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momz_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	drhot_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dpres_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd_ref,
		real(rp), dimension(elem2d%np,lmesh2d%nea), intent(in)	coriolis,
		real(rp), dimension (elem%np,lmesh%nea), intent(in)	rtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cvtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cptot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dphyddx,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dphyddy,
		class(elementoperationbase3d), intent(in)	element3d_operation,
		type(sparsemat), intent(in)	dx,
		type(sparsemat), intent(in)	dy,
		type(sparsemat), intent(in)	dz,
		type(sparsemat), intent(in)	sx,
		type(sparsemat), intent(in)	sy,
		type(sparsemat), intent(in)	sz,
		type(sparsemat), intent(in)	lift,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(localmesh2d), intent(in)	lmesh2d,
		class(elementbase2d), intent(in)	elem2d )

Definition at line 98 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve.F90.

 
    use scale_atm_dyn_dgm_nonhydro3d_rhot_heve_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc_asis
    use scale_const, only: &
      ohm => const_ohm
    implicit none
 
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    class(LocalMesh2D), intent(in) :: lmesh2D
    class(ElementBase2D), intent(in) :: elem2D
    class(ElementOperationBase3D), intent(in) :: element3D_operation
    type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift
    real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: RHOT_dt(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DDENS_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMX_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMY_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMZ_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DRHOT_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DPRES_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DENS_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: PRES_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: PRES_hyd_ref(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: CORIOLIS(elem2D%Np,lmesh2D%NeA)
    real(RP), intent(in)  :: Rtot (elem%Np,lmesh%NeA)    
    real(RP), intent(in)  :: CVtot(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: CPtot(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DPhydDx(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DPhydDy(elem%Np,lmesh%NeA)
 
    real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)
    real(RP) :: DPRES_hyd(elem%Np), GradPhyd_x(elem%Np), GradPhyd_y(elem%Np)
    real(RP) :: del_flux(elem%NfpTot,lmesh%Ne,PRGVAR_NUM)
    real(RP) :: del_flux_hyd(elem%NfpTot,lmesh%Ne,2)
    real(RP) :: RHOT_(elem%Np)
    real(RP) :: rdens_(elem%Np), u_(elem%Np), v_(elem%Np), w_(elem%Np), wt_(elem%np), drho(elem%Np)
 
    real(RP) :: G11(elem%Np), G12(elem%Np), G22(elem%Np)
    real(RP) :: GsqrtV(elem%Np), RGsqrtV(elem%Np)
    real(RP) :: X2D(elem2D%Np,lmesh2D%Ne), Y2D(elem2D%Np,lmesh2D%Ne)
    real(RP) :: X(elem%Np), Y(elem%Np), twoOVdel2(elem%Np)
    real(RP) :: CORI(elem%Np,2)
    logical :: is_panel1to4
    real(RP) :: s
 
    integer :: ke, ke2d
    integer :: p, p12, p3
 
    real(RP) :: gamm, rgamm    
    real(RP) :: rP0
    real(RP) :: RovP0, P0ovR
    !------------------------------------------------------------------------
 
    call prof_rapstart('cal_dyn_tend_bndflux', 3)
    call get_ebnd_flux( &
      del_flux, del_flux_hyd,                                                  & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd,         & ! (in)
      rtot, cvtot, cptot,                                                      & ! (in)
      lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
      lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2),             & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3),  & ! (in)
      lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd,                         & ! (in)
      lmesh, elem, lmesh2d, elem2d                                             ) ! (in)
    call prof_rapend('cal_dyn_tend_bndflux', 3)
 
    !-----
    call prof_rapstart('cal_dyn_tend_interior', 3)
    gamm  = cpdry / cvdry
    rgamm = cvdry / cpdry
    rp0   = 1.0_rp / pres00
    rovp0 = rdry * rp0
    p0ovr = pres00 / rdry
 
    s = 1.0_rp 
    is_panel1to4 = .true.
    if ( lmesh%panelID == 5 ) then
      is_panel1to4 = .false.
    else if ( lmesh%panelID == 6 ) then
      is_panel1to4 = .false.
      s = - 1.0_rp
    end if
 
    !$omp parallel private(                        &
    !$omp RHOT_, rdens_, u_, v_, w_, wt_,          &
    !$omp Fx, Fy, Fz, LiftDelFlx,                  &
    !$omp drho, DPRES_hyd, GradPhyd_x, GradPhyd_y, &
    !$omp G11, G12, G22, GsqrtV, RGsqrtV,          &
    !$omp X, Y, twoOVdel2,                         &
    !$omp CORI, ke, ke2D                           )
 
    !$omp do
    do ke2d = lmesh2d%NeS, lmesh2d%NeE
      x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
      y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
    end do
 
    !$omp do
    do ke = lmesh%NeS, lmesh%NeE
      !--
      ke2d = lmesh%EMap3Dto2D(ke)
      g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1)
      g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2)
      g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2)
      gsqrtv(:)  = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
      rgsqrtv(:) = 1.0_rp / gsqrtv(:)
 
      !--
      rhot_(:) = p0ovr * ( pres_hyd(:,ke) * rp0 )**rgamm + drhot_(:,ke)
      ! DPRES_(:) = PRES00 * ( Rtot(:,ke) * rP0 * RHOT_(:) )**( CPtot(:,ke) / CVtot(:,ke) ) &
      !           - PRES_hyd(:,ke)
 
      rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
      u_(:) = momx_(:,ke) * rdens_(:)
      v_(:) = momy_(:,ke) * rdens_(:)
      w_(:) = momz_(:,ke) * rdens_(:)
      wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:) 
 
      x(:) = x2d(elem%IndexH2Dto3D,ke2d)
      y(:) = y2d(elem%IndexH2Dto3D,ke2d)
      twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
 
      cori(:,1) = s * ohm * twoovdel2(:) * ( - x(:) * y(:)          * momx_(:,ke) + ( 1.0_rp + y(:)**2 ) * momy_(:,ke) )
      cori(:,2) = s * ohm * twoovdel2(:) * ( - ( 1.0_rp + x(:)**2 ) * momx_(:,ke) +  x(:) * y(:)         * momy_(:,ke) )
      if ( is_panel1to4 ) then
        cori(:,1) = s * y(:) * cori(:,1)
        cori(:,2) = s * y(:) * cori(:,2)
      end if
 
      drho(:) = matmul(intrpmat_vpordm1, ddens_(:,ke))
 
      !-- Gradient hydrostatic pressure
      
      dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
 
      call sparsemat_matmul(dx, gsqrtv(:) * dpres_hyd(:), fx)
      call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,1) * dpres_hyd(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
      gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
                    + lmesh%Escale(:,ke,3,3) * fz(:) &
                    + liftdelflx(:)
 
      call sparsemat_matmul(dy, gsqrtv(:) * dpres_hyd(:), fy)
      call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,2) * dpres_hyd(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
      gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
                    + lmesh%Escale(:,ke,3,3) * fz(:) &
                    + liftdelflx(:)
      
      !-- DENS
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * momx_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * momy_(:,ke), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
 
      dens_dt(:,ke) = - ( &
            lmesh%Escale(:,ke,1,1) * fx(:)     &
          + lmesh%Escale(:,ke,2,2) * fy(:)     &
          + lmesh%Escale(:,ke,3,3) * fz(:)     &
          + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
      
      !-- MOMX
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momx_(:,ke)                              &
                                                    + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
 
      momx_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)                                         &
            + lmesh%Escale(:,ke,2,2) * fy(:)                                         &
            + lmesh%Escale(:,ke,3,3) * fz(:)                                         &
            + liftdelflx(:)                   ) / lmesh%Gsqrt(:,ke)                  &
          - twoovdel2(:) * y(:) *                                                    &
            ( x(:) * y(:) * u_(:) - (1.0_rp + y(:)**2) * v_(:) ) * momx_(:,ke)       &
          - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:)         &
          + cori(:,1)              
 
      !-- MOMY
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momy_(:,ke)                              &
                                                    + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
 
      momy_dt(:,ke) = &
            - ( lmesh%Escale(:,ke,1,1) * fx(:)                                         &
              + lmesh%Escale(:,ke,2,2) * fy(:)                                         &
              + lmesh%Escale(:,ke,3,3) * fz(:)                                         &
              + liftdelflx(:)                  ) / lmesh%Gsqrt(:,ke)                   &
            - twoovdel2(:) * x(:) *                                                    &
              ( - (1.0_rp + x(:)**2) * u_(:) + x(:) * y(:) * v_(:) ) * momy_(:,ke)     &
            - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:)         &
            + cori(:,2) 
 
      !-- MOMZ
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) *   u_(:) * momz_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) *   v_(:) * momz_(:,ke), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momz_(:,ke) + rgsqrtv(:) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
      
      momz_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)       &
            + lmesh%Escale(:,ke,2,2) * fy(:)       &
            + lmesh%Escale(:,ke,3,3) * fz(:)       &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)  &
          - grav * drho(:)
 
      !-- RHOT
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * rhot_(:), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * rhot_(:), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * wt_(:) * rhot_(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
      
      rhot_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)      &
            + lmesh%Escale(:,ke,2,2) * fy(:)      &
            + lmesh%Escale(:,ke,3,3) * fz(:)      &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) 
    
    end do
    !$omp end do
    !$omp end parallel
    call prof_rapend('cal_dyn_tend_interior', 3)
 
    return

References scale_atm_dyn_dgm_nonhydro3d_rhot_heve_numflux::atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc_asis(), and scale_atm_dyn_dgm_nonhydro3d_common::intrpmat_vpordm1.

atm_dyn_dgm_globalnonhydro3d_rhot_heve_cal_tend_shallow_atm()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve::atm_dyn_dgm_globalnonhydro3d_rhot_heve_cal_tend_shallow_atm	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	dens_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momx_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momy_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momz_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	rhot_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	ddens_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momx_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momy_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momz_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	drhot_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dpres_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd_ref,
		real(rp), dimension(elem2d%np,lmesh2d%nea), intent(in)	coriolis,
		real(rp), dimension (elem%np,lmesh%nea), intent(in)	rtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cvtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cptot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dphyddx,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dphyddy,
		class(elementoperationbase3d), intent(in)	element3d_operation,
		type(sparsemat), intent(in)	dx,
		type(sparsemat), intent(in)	dy,
		type(sparsemat), intent(in)	dz,
		type(sparsemat), intent(in)	sx,
		type(sparsemat), intent(in)	sy,
		type(sparsemat), intent(in)	sz,
		type(sparsemat), intent(in)	lift,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(localmesh2d), intent(in)	lmesh2d,
		class(elementbase2d), intent(in)	elem2d )

Definition at line 336 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve.F90.

 
    use scale_atm_dyn_dgm_nonhydro3d_rhot_heve_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc
    use scale_const, only: &
      ohm => const_ohm
    implicit none
 
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    class(LocalMesh2D), intent(in) :: lmesh2D
    class(ElementBase2D), intent(in) :: elem2D
    class(ElementOperationBase3D), intent(in) :: element3D_operation
    type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift
    real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: RHOT_dt(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DDENS_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMX_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMY_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMZ_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DRHOT_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DPRES_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DENS_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: PRES_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: PRES_hyd_ref(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: CORIOLIS(elem2D%Np,lmesh2D%NeA)
    real(RP), intent(in)  :: Rtot (elem%Np,lmesh%NeA)    
    real(RP), intent(in)  :: CVtot(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: CPtot(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DPhydDx(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DPhydDy(elem%Np,lmesh%NeA)
 
    real(RP) :: Flux(elem%Np,3,5), DFlux(elem%Np,4,5)
    real(RP) :: del_flux(elem%NfpTot,PRGVAR_NUM,lmesh%Ne)
    real(RP) :: rdens_, u_, v_, w_, pt_
    real(RP) :: drho(elem%Np)
    real(RP) :: GsqrtV(elem%Np), RGsqrtV(elem%Np), RGsqrt(elem%Np)
    real(RP) :: Gsqrt_, GsqrtDPRES_, E11, E22, E33
 
    real(RP) :: G11(elem%Np), G12(elem%Np), G22(elem%Np)
    real(RP) :: X2D(elem2D%Np,lmesh2D%Ne), Y2D(elem2D%Np,lmesh2D%Ne)
    real(RP) :: X(elem%Np), Y(elem%Np), twoOVdel2(elem%Np)
    real(RP) :: CORI(elem%Np,2)
    logical :: is_panel1to4
    real(RP) :: s
 
    integer :: ke, ke2d
    integer :: p, p12, p3
 
    real(RP) :: gamm, rgamm    
    real(RP) :: rP0
    real(RP) :: RovP0, P0ovR
    !------------------------------------------------------------------------
 
    call prof_rapstart('cal_dyn_tend_bndflux', 3)
    call get_ebnd_flux( &
      del_flux,                                                                & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd,         & ! (in)
      rtot, cvtot, cptot,                                                      & ! (in)
      lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
      lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2),             & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3),  & ! (in)
      lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd,                         & ! (in)
      lmesh, elem, lmesh2d, elem2d                                             ) ! (in)
    call prof_rapend('cal_dyn_tend_bndflux', 3)
 
    !-----
    call prof_rapstart('cal_dyn_tend_interior', 3)
    gamm  = cpdry / cvdry
    rgamm = cvdry / cpdry
    rp0   = 1.0_rp / pres00
    rovp0 = rdry * rp0
    p0ovr = pres00 / rdry
 
    s = 1.0_rp 
    is_panel1to4 = .true.
    if ( lmesh%panelID == 5 ) then
      is_panel1to4 = .false.
    else if ( lmesh%panelID == 6 ) then
      is_panel1to4 = .false.
      s = - 1.0_rp
    end if
 
    !$omp parallel private( ke, ke2d, p,                  &
    !$omp rdens_, u_, v_, w_, pt_,                        &
    !$omp drho, CORI,                                     &
    !$omp RGsqrt, GsqrtV, RGsqrtV, Gsqrt_, GsqrtDPRES_,   &
    !$omp G11, G12, G22, X, Y, twoOVdel2,                 &
    !$omp E11, E22, E33, DFlux, Flux    )
 
    !$omp do
    do ke2d = lmesh2d%NeS, lmesh2d%NeE
      x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
      y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
    end do
 
    !$omp do
    do ke = lmesh%NeS, lmesh%NeE
      !--
      ke2d = lmesh%EMap3Dto2D(ke)
 
      do p=1, elem%Np
        g11(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,1,1)
        g12(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,1,2)
        g22(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,2,2)
      end do
      do p=1, elem%Np
        gsqrtv(p)  = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d)
        rgsqrtv(p) = 1.0_rp / gsqrtv(p)
        rgsqrt(p) = 1.0_rp / lmesh%Gsqrt(p,ke)
      end do
 
      !--
 
      do p=1, elem%Np
        gsqrt_ = lmesh%Gsqrt(p,ke)
        flux(p,1,dens_vid) = gsqrt_ * momx_(p,ke)
        flux(p,2,dens_vid) = gsqrt_ * momy_(p,ke)
        flux(p,3,dens_vid) = gsqrt_ * ( &
            momz_(p,ke) * rgsqrtv(p)        &
          + lmesh%GI3(p,ke,1) * momx_(p,ke) &
          + lmesh%GI3(p,ke,2) * momy_(p,ke) )
      end do
      do p=1, elem%Np
        pt_ = ( p0ovr * (pres_hyd(p,ke) * rp0)**rgamm + drhot_(p,ke) ) &
            / ( ddens_(p,ke) + dens_hyd(p,ke) )
 
        flux(p,1,rhot_vid) = flux(p,1,dens_vid) * pt_
        flux(p,2,rhot_vid) = flux(p,2,dens_vid) * pt_
        flux(p,3,rhot_vid) = flux(p,3,dens_vid) * pt_
      end do
 
      do p=1, elem%Np
        w_ = momz_(p,ke) / (ddens_(p,ke) + dens_hyd(p,ke))
 
        flux(p,1,momz_vid) = flux(p,1,dens_vid) * w_
        flux(p,2,momz_vid) = flux(p,2,dens_vid) * w_
        flux(p,3,momz_vid) = flux(p,3,dens_vid) * w_ + lmesh%Gsqrt(p,ke) * rgsqrtv(p) * dpres_(p,ke)
      end do
 
      do p=1, elem%Np
        gsqrtdpres_ = lmesh%Gsqrt(p,ke) * dpres_(p,ke)
 
        rdens_ = 1.0_rp / (ddens_(p,ke) + dens_hyd(p,ke))
        u_ = momx_(p,ke) * rdens_
        v_ = momy_(p,ke) * rdens_
 
        flux(p,1,momx_vid) = flux(p,1,dens_vid) * u_ + g11(p) * gsqrtdpres_
        flux(p,2,momx_vid) = flux(p,2,dens_vid) * u_ + g12(p) * gsqrtdpres_
        flux(p,3,momx_vid) = flux(p,3,dens_vid) * u_ + gsqrtdpres_ * ( g11(p) * lmesh%GI3(p,ke,1) + g12(p) * lmesh%GI3(p,ke,2) ) 
 
        flux(p,1,momy_vid) = flux(p,1,dens_vid) * v_ + g12(p) * gsqrtdpres_
        flux(p,2,momy_vid) = flux(p,2,dens_vid) * v_ + g22(p) * gsqrtdpres_ 
        flux(p,3,momy_vid) = flux(p,3,dens_vid) * v_ + gsqrtdpres_ * ( g12(p) * lmesh%GI3(p,ke,1) + g22(p) * lmesh%GI3(p,ke,2) ) 
      end do
 
      call element3d_operation%Div_var5( &
        flux, del_flux(:,:,ke), & ! (in)
        dflux                   ) ! (out)
 
 
      do p=1, elem%Np
        e11 = lmesh%Escale(p,ke,1,1)
        e22 = lmesh%Escale(p,ke,2,2)
        e33 = lmesh%Escale(p,ke,3,3)
 
        dens_dt(p,ke) = - ( &
              e11 * dflux(p,1,dens_vid)    &
            + e22 * dflux(p,2,dens_vid)    &
            + e33 * dflux(p,3,dens_vid)    &
            + dflux(p,4,dens_vid) ) * rgsqrt(p)
          
        rhot_dt(p,ke) = - ( &
              e11 * dflux(p,1,rhot_vid)    &
            + e22 * dflux(p,2,rhot_vid)    &
            + e33 * dflux(p,3,rhot_vid)    &
            + dflux(p,4,rhot_vid) ) * rgsqrt(p)
      end do
 
      call element3d_operation%VFilterPM1( ddens_(:,ke), & ! (in)
        drho ) ! (out)  
 
      do p=1, elem%Np
        e11 = lmesh%Escale(p,ke,1,1)
        e22 = lmesh%Escale(p,ke,2,2)
        e33 = lmesh%Escale(p,ke,3,3)
            
        momz_dt(p,ke) = - ( &
              e11 * dflux(p,1,momz_vid)     &
            + e22 * dflux(p,2,momz_vid)     &
            + e33 * dflux(p,3,momz_vid)     &
            + dflux(p,4,momz_vid) ) * rgsqrt(p) &
            - grav * drho(p)            
      end do
 
      !--
      do p=1, elem%Np
        x(p) = x2d(elem%IndexH2Dto3D(p),ke2d)
        y(p) = y2d(elem%IndexH2Dto3D(p),ke2d)
        twoovdel2(p) = 2.0_rp / ( 1.0_rp + x(p)**2 + y(p)**2 )
      end do
 
      do p=1, elem%Np
        cori(p,1) = s * ohm * twoovdel2(p) * ( - x(p) * y(p)          * momx_(p,ke) + ( 1.0_rp + y(p)**2 ) * momy_(p,ke) )
        cori(p,2) = s * ohm * twoovdel2(p) * ( - ( 1.0_rp + x(p)**2 ) * momx_(p,ke) +  x(p) * y(p)         * momy_(p,ke) )
      end do
      if ( is_panel1to4 ) then
        do p=1, elem%Np
          cori(p,1) = s * y(p) * cori(p,1)
          cori(p,2) = s * y(p) * cori(p,2)
        end do
      end if
 
      do p=1, elem%Np
        rdens_ = 1.0_rp / (ddens_(p,ke) + dens_hyd(p,ke))
        u_ = momx_(p,ke) * rdens_
        v_ = momy_(p,ke) * rdens_
 
        momx_dt(p,ke) = - ( g11(p) * dphyddx(p,ke) + g12(p) * dphyddy(p,ke) ) &
                        - twoovdel2(p) * y(p) *                                            &
                          ( x(p) * y(p) * u_ - (1.0_rp + y(p)**2) * v_ ) * momx_(p,ke)     &    
                        + cori(p,1)
        
        momy_dt(p,ke) = - ( g12(p) * dphyddx(p,ke) + g22(p) * dphyddy(p,ke) ) &
                        - twoovdel2(p) * x(p) *                                            &
                          ( - (1.0_rp + x(p)**2) * u_ + x(p) * y(p) * v_ ) * momy_(p,ke)   &
                        + cori(p,2)
      end do
 
      do p=1, elem%Np
        e11 = lmesh%Escale(p,ke,1,1)
        e22 = lmesh%Escale(p,ke,2,2)
        e33 = lmesh%Escale(p,ke,3,3)
  
        momx_dt(p,ke) = momx_dt(p,ke) - ( &
              e11 * dflux(p,1,momx_vid)    &
            + e22 * dflux(p,2,momx_vid)    &
            + e33 * dflux(p,3,momx_vid)    &
            + dflux(p,4,momx_vid) ) * rgsqrt(p)
        
        momy_dt(p,ke) = momy_dt(p,ke) - ( &
              e11 * dflux(p,1,momy_vid)    &
            + e22 * dflux(p,2,momy_vid)    &
            + e33 * dflux(p,3,momy_vid)    &
            + dflux(p,4,momy_vid) ) * rgsqrt(p)
      end do
    end do
    !$omp end do
    !$omp end parallel
    call prof_rapend('cal_dyn_tend_interior', 3)
 
    return

References scale_atm_dyn_dgm_nonhydro3d_rhot_heve_numflux::atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc().

atm_dyn_dgm_globalnonhydro3d_rhot_heve_cal_tend_deep_atm()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve::atm_dyn_dgm_globalnonhydro3d_rhot_heve_cal_tend_deep_atm	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	dens_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momx_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momy_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momz_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	rhot_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	ddens_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momx_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momy_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momz_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	drhot_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dpres_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd_ref,
		real(rp), dimension(elem2d%np,lmesh2d%nea), intent(in)	coriolis,
		real(rp), dimension (elem%np,lmesh%nea), intent(in)	rtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cvtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cptot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dphyddx,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dphyddy,
		class(elementoperationbase3d), intent(in)	element3d_operation,
		type(sparsemat), intent(in)	dx,
		type(sparsemat), intent(in)	dy,
		type(sparsemat), intent(in)	dz,
		type(sparsemat), intent(in)	sx,
		type(sparsemat), intent(in)	sy,
		type(sparsemat), intent(in)	sz,
		type(sparsemat), intent(in)	lift,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(localmesh2d), intent(in)	lmesh2d,
		class(elementbase2d), intent(in)	elem2d )

Definition at line 599 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve.F90.

 
    use scale_atm_dyn_dgm_nonhydro3d_rhot_heve_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc_asis
    use scale_const, only: &
      ohm => const_ohm
    implicit none
 
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    class(LocalMesh2D), intent(in) :: lmesh2D
    class(ElementBase2D), intent(in) :: elem2D
    class(ElementOperationBase3D), intent(in) :: element3D_operation
    type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift
    real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: RHOT_dt(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DDENS_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMX_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMY_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMZ_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DRHOT_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DPRES_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DENS_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: PRES_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: PRES_hyd_ref(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: CORIOLIS(elem2D%Np,lmesh2D%NeA)
    real(RP), intent(in)  :: Rtot (elem%Np,lmesh%NeA)    
    real(RP), intent(in)  :: CVtot(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: CPtot(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DPhydDx(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DPhydDy(elem%Np,lmesh%NeA)
 
    real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)
    real(RP) :: DPRES_hyd(elem%Np), GradPhyd_x(elem%Np), GradPhyd_y(elem%Np)
    real(RP) :: del_flux(elem%NfpTot,lmesh%Ne,PRGVAR_NUM)
    real(RP) :: del_flux_hyd(elem%NfpTot,lmesh%Ne,2)
    real(RP) :: RHOT_(elem%Np)
    real(RP) :: rdens_(elem%Np), u_(elem%Np), v_(elem%Np), w_(elem%Np), wt_(elem%np), drho(elem%Np)
 
    real(RP) :: G11(elem%Np), G12(elem%Np), G22(elem%Np)
    real(RP) :: GsqrtV(elem%Np), RGsqrtV(elem%Np), Rgam2(elem%Np)
    real(RP) :: X2D(elem2D%Np,lmesh2D%Ne), Y2D(elem2D%Np,lmesh2D%Ne)
    real(RP) :: X(elem%Np), Y(elem%Np), twoOVdel2(elem%Np)
    real(RP) :: OM1(elem%Np), OM2(elem%Np), OM3(elem%Np), DEL(elem%Np), R(elem%Np)
    logical :: is_panel1to4
    real(RP) :: s
 
    integer :: ke, ke2d
    integer :: p
    
    real(RP) :: rgamm    
    real(RP) :: rP0
    real(RP) :: P0ovR
    !------------------------------------------------------------------------
 
    call prof_rapstart('cal_dyn_tend_bndflux', 3)
    call get_ebnd_flux( &
      del_flux, del_flux_hyd,                                                  & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd,         & ! (in)
      rtot, cvtot, cptot,                                                      & ! (in)
      lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
      lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2),             & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3),  & ! (in)
      lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd,                         & ! (in)
      lmesh, elem, lmesh2d, elem2d                                             ) ! (in)
    call prof_rapend('cal_dyn_tend_bndflux', 3)
 
    !-----
    call prof_rapstart('cal_dyn_tend_interior', 3)
    rgamm = cvdry / cpdry
    rp0   = 1.0_rp / pres00
    p0ovr = pres00 / rdry
 
    s = 2.0_rp * ohm
    is_panel1to4 = .true.
    if ( lmesh%panelID == 5 ) then
      is_panel1to4 = .false.
    else if ( lmesh%panelID == 6 ) then
      is_panel1to4 = .false.
      s = - s
    end if
 
    !$omp parallel private(                           &
    !$omp RHOT_, rdens_, u_, v_, w_, wt_,             &
    !$omp Fx, Fy, Fz, LiftDelFlx,                     &
    !$omp drho, DPRES_hyd, GradPhyd_x, GradPhyd_y,    &
    !$omp G11, G12, G22, Rgam2, GsqrtV, RGsqrtV,      &
    !$omp X, Y, twoOVdel2,                            &
    !$omp OM1, OM2, OM3, DEL, R, ke, ke2D             )
 
    !$omp do
    do ke2d = lmesh2d%NeS, lmesh2d%NeE
      x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
      y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
    end do
 
    !$omp do
    do ke = lmesh%NeS, lmesh%NeE
      !--
      ke2d = lmesh%EMap3Dto2D(ke)
      rgam2(:) = 1.0_rp / lmesh%gam(:,ke)**2
      g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1) * rgam2(:)
      g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2) * rgam2(:)
      g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2) * rgam2(:)
      gsqrtv(:)  = lmesh%Gsqrt(:,ke) * rgam2(:) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
      rgsqrtv(:) = 1.0_rp / gsqrtv(:)
 
      !--
      rhot_(:) = p0ovr * ( pres_hyd(:,ke) * rp0 )**rgamm + drhot_(:,ke)
      ! DPRES_(:) = PRES00 * ( Rtot(:,ke) * rP0 * RHOT_(:) )**( CPtot(:,ke) / CVtot(:,ke) ) &
      !           - PRES_hyd(:,ke)
 
      rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
      u_(:) = momx_(:,ke) * rdens_(:)
      v_(:) = momy_(:,ke) * rdens_(:)
      w_(:) = momz_(:,ke) * rdens_(:)
      wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:) 
 
      x(:) = x2d(elem%IndexH2Dto3D,ke2d)
      y(:) = y2d(elem%IndexH2Dto3D,ke2d)
      del(:) = sqrt( 1.0_rp + x(:)**2 + y(:)**2 )
      twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
 
      r(:) = rplanet * lmesh%gam(:,ke)
 
      !-  pnl=1~4: OM1:                     0, OM2: s del / (r (1+Y^2)) ,   s OM3 : Y /del
      !-  pnl=5,6: OM1: - s X del/(r (1+X^2)), OM2: -  s Y del/(r(1+Y^2)), OM3 : s/del
      if ( is_panel1to4 ) then
        om1(:) = 0.0_rp
        om2(:) = s * del(:) / ( r(:) * ( 1.0_rp + y(:)**2 ) )
        om3(:) = s * y(:) / del(:)
      else
        om1(:) = - s * x(:) * del(:) / ( r(:) * ( 1.0_rp + x(:)**2 ) )
        om2(:) = - s * y(:) * del(:) / ( r(:) * ( 1.0_rp + y(:)**2 ) )
        om3(:) =   s / del(:)
      end if
 
      drho(:) = matmul(intrpmat_vpordm1, ddens_(:,ke))
 
      !-- Gradient hydrostatic pressure
      
      dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
 
      call sparsemat_matmul(dx, gsqrtv(:) * dpres_hyd(:), fx)
      call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,1) * dpres_hyd(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
      gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
                    + lmesh%Escale(:,ke,3,3) * fz(:) &
                    + liftdelflx(:)
 
      call sparsemat_matmul(dy, gsqrtv(:) * dpres_hyd(:), fy)
      call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,2) * dpres_hyd(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
      gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
                    + lmesh%Escale(:,ke,3,3) * fz(:) &
                    + liftdelflx(:)
      
      !-- DENS
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * momx_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * momy_(:,ke), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
 
      dens_dt(:,ke) = - ( &
            lmesh%Escale(:,ke,1,1) * fx(:)     &
          + lmesh%Escale(:,ke,2,2) * fy(:)     &
          + lmesh%Escale(:,ke,3,3) * fz(:)     &
          + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
      
      !-- MOMX
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momx_(:,ke)                              &
                                                    + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
 
      momx_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)                                                &
            + lmesh%Escale(:,ke,2,2) * fy(:)                                                &
            + lmesh%Escale(:,ke,3,3) * fz(:)                                                &
            + liftdelflx(:)                   ) / lmesh%Gsqrt(:,ke)                         &
          - twoovdel2(:) * y(:) *                                                           & !-> metric terms
            ( x(:) * y(:) * u_(:) - ( 1.0_rp + y(:)**2 ) * v_(:) ) * momx_(:,ke)            & !
          - 2.0_rp * u_(:) * momz_(:,ke) / r(:)                                             & 
          - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:)                & !-> gradient hydrostatic pressure
          - lmesh%Gsqrt(:,ke) * (  g11(:) * ( om2(:) * momz_(:,ke) - om3(:) * momy_(:,ke) ) & !-> Coriolis term
                                 - g12(:) * ( om1(:) * momz_(:,ke) - om3(:) * momx_(:,ke) ) )
 
      !-- MOMY
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momy_(:,ke)                              &
                                                    + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
 
      momy_dt(:,ke) = &
            - ( lmesh%Escale(:,ke,1,1) * fx(:)                                                &
              + lmesh%Escale(:,ke,2,2) * fy(:)                                                &
              + lmesh%Escale(:,ke,3,3) * fz(:)                                                &
              + liftdelflx(:)                  ) / lmesh%Gsqrt(:,ke)                          &
            - twoovdel2(:) * x(:) *                                                           & !-> metric terms
              ( - (1.0_rp + x(:)**2) * u_(:) + x(:) * y(:) * v_(:) ) * momy_(:,ke)            & !
            - 2.0_rp * v_(:) * momz_(:,ke) / r(:)                                             & 
            - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:)                & !-> gradient hydrostatic pressure
            - lmesh%Gsqrt(:,ke) * (  g12(:) * ( om2(:) * momz_(:,ke) - om3(:) * momy_(:,ke) ) & !-> Coriolis term
                                   - g22(:) * ( om1(:) * momz_(:,ke) - om3(:) * momx_(:,ke) ) )
 
      !-- MOMZ
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) *   u_(:) * momz_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) *   v_(:) * momz_(:,ke), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momz_(:,ke) + rgsqrtv(:) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
      
      momz_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)                                                                    &
            + lmesh%Escale(:,ke,2,2) * fy(:)                                                                    &
            + lmesh%Escale(:,ke,3,3) * fz(:)                                                                    &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)                                                               &
          - 0.25_rp * r(:) * twoovdel2(:)**2 * ( 1.0_rp * x(:)**2 ) * ( 1.0_rp * y(:)**2 )                      & !-> metric terms
            * ( - ( 1.0_rp + x(:)**2 ) * momx_(:,ke) * u_(:)                                                    & !
                + 2.0_rp * x(:) * y(:) * momx_(:,ke) * v_(:)                                                    & !
                - ( 1.0_rp + y(:)**2 ) * momy_(:,ke) * v_(:)  )                                                 & 
          + 2.0_rp * dpres_(:,ke) / r(:)                                                                        & !-> metric term with gradient of pressure deviaition
          - lmesh%Gsqrt(:,ke) * ( om1(:) * momy_(:,ke) - om2(:) * momx_(:,ke) )                                 & !-> Coriolis term
          - grav * rgam2(:) * drho(:)                                                                             !-> buoyancy term
 
      !-- RHOT
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * rhot_(:), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * rhot_(:), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * wt_(:) * rhot_(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
      
      rhot_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)      &
            + lmesh%Escale(:,ke,2,2) * fy(:)      &
            + lmesh%Escale(:,ke,3,3) * fz(:)      &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) 
    end do
    !$omp end do
    !$omp end parallel
    call prof_rapend('cal_dyn_tend_interior', 3)
 
    return

References scale_atm_dyn_dgm_nonhydro3d_rhot_heve_numflux::atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc_asis(), and scale_atm_dyn_dgm_nonhydro3d_common::intrpmat_vpordm1.