scale_atm_dyn_dgm_nonhydro3d_rhot_heve_splitform.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Atmosphere / Dynamics HEVE 
!!
!! @par Description
!!      HEVE DGM scheme for Atmospheric dynamical process. 
!!      To improve the numerical instability due to the aliasing errors, 
!!      the split form based on Gassner et al. (2016, JCP) is used for advection terms. 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_nonhydro3d_rhot_heve_splitform
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    grav => const_grav,  &
    rdry => const_rdry,  &
    cpdry => const_cpdry, &
    cvdry => const_cvdry, &
    pres00 => const_pre00
 
  use scale_sparsemat  
  use scale_element_modalfilter, only: modalfilter
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
  use scale_element_operation_base, only: elementoperationbase3d
  use scale_localmesh_2d, only: localmesh2d  
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmeshfield_base, only: localmeshfield3d
  use scale_meshfield_base, only: meshfield3d
 
  use scale_atm_dyn_dgm_nonhydro3d_common, only: &
    atm_dyn_dgm_nonhydro3d_common_init,                       &
    atm_dyn_dgm_nonhydro3d_common_final,                      &
    dens_vid => prgvar_ddens_id, rhot_vid => prgvar_drhot_id, &
    momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id,   &
    momz_vid => prgvar_momz_id,                               &
    prgvar_num, intrpmat_vpordm1
  
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_init
  public :: atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_final
  public :: atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_cal_tend
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
  real(RP), private, allocatable :: DxT1D_(:,:)
  real(RP), private, allocatable :: DyT1D_(:,:)
  real(RP), private, allocatable :: DzT1D_(:,:)
 
  private :: dx_ab, dy_ab, dz_ab
  private :: dx_abc, dy_abc, dz_abc
 
contains
  subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_init( mesh )
    implicit none
    class(meshbase3d), intent(in) :: mesh
 
    integer :: p1, p2, p3, p_
    type(elementbase3d), pointer :: elem
    !--------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_init( mesh )
    elem => mesh%refElem3D
 
    allocate( dxt1d_(elem%Nnode_h1D,elem%Nnode_h1D) )
    allocate( dyt1d_(elem%Nnode_h1D,elem%Nnode_h1D) )
    allocate( dzt1d_(elem%Nnode_v,elem%Nnode_v) )
 
    do p1=1, elem%Nnode_h1D
      dxt1d_(:,p1) = elem%Dx1(p1,1:elem%Nnode_h1D)
    end do
 
    do p2=1, elem%Nnode_h1D
    do p_=1, elem%Nnode_h1D
      dyt1d_(p_,p2) = elem%Dx2(1+(p2-1)*elem%Nnode_h1D,1+(p_-1)*elem%Nnode_h1D)
    end do
    end do
 
    do p3=1, elem%Nnode_v
      dzt1d_(:,p3) = elem%Dx3(elem%Colmask(p3,1),elem%Colmask(:,1))
    end do
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_init
 
 
  subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_final()
    implicit none
    !--------------------------------------------
 
    deallocate( dxt1d_, dyt1d_, dzt1d_ )
    call atm_dyn_dgm_nonhydro3d_common_final()
    
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_final  
 
  !-------------------------------
 
  subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_cal_tend( &
    DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt,         & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dpres_,         & ! (in)
    dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd,         & ! (in)
    coriolis, rtot, cvtot, cptot, dphyddx, dphyddy,      & ! (in)
    element3d_operation, dx, dy, dz, sx, sy, sz, lift,   & ! (in)
    lmesh, elem, lmesh2d, elem2d )                         ! (in)
 
    use scale_atm_dyn_dgm_nonhydro3d_rhot_heve_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalvc_asis
    
    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)  :: therm_hyd(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) :: fx_sp(elem%np), fy_sp(elem%np), fz_sp(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) :: gsqrtdens_(elem%np), rdens_(elem%np), rhot_(elem%np)
    real(rp) :: u_(elem%np), v_(elem%np), w_(elem%np), wt_(elem%np), pot_(elem%np)
    real(rp) :: drho(elem%np), cori(elem%np)
    real(rp) :: gsqrtv(elem%np), rgsqrtv(elem%np)
 
    integer :: ke, ke2d
 
    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%GI3(:,:,1), lmesh%GI3(:,:,2),                        & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP,                                               & ! (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
 
    !$omp parallel do private( ke2d, Cori,                          &
    !$omp RHOT_, GsqrtDens_, rdens_, u_, v_, w_, wt_, pot_,         &
    !$omp DPRES_hyd, GradPhyd_x, GradPhyd_y, drho,                  &
    !$omp GsqrtV, RGsqrtV,                                          &
    !$omp Fx, Fy, Fz, Fx_sp, Fy_sp, Fz_sp, LiftDelFlx               )
    do ke = lmesh%NeS, lmesh%NeE
      !--
      ke2d = lmesh%EMap3Dto2D(ke)
      cori(:) = coriolis(elem%IndexH2Dto3D(:),ke2d)
 
      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)
 
      gsqrtdens_(:) = lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) )
      rdens_(:) = 1.0_rp / gsqrtdens_(:)
      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_(:) 
      pot_(:) = rhot_(:) * rdens_(:)
 
      ke2d = lmesh%EMap3Dto2D(ke)
      cori(:) = coriolis(elem%IndexH2Dto3D(:),ke2d)
 
      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 dx_ab( dxt1d_, gsqrtdens_(:),     u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_ab( dyt1d_, gsqrtdens_(:),     v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
      call dz_ab( dzt1d_, gsqrtdens_(:), wt_(:), elem%Nnode_h1D, elem%Nnode_v, fz_sp )
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
 
      dens_dt(:,ke) = - ( &
            lmesh%Escale(:,ke,1,1) * fx_sp(:) &
          + lmesh%Escale(:,ke,2,2) * fy_sp(:) &
          + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
          + liftdelflx(:) )
      
      !-- MOMX
      call dx_abc( dxt1d_, gsqrtdens_, u_,  u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_abc( dyt1d_, gsqrtdens_, u_,  v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
      call dz_abc( dzt1d_, gsqrtdens_, u_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * dpres_(:,ke)              , fx)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
 
      momx_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * ( fx_sp(:) + fx(:) )  &
            + lmesh%Escale(:,ke,2,2) * fy_sp(:)              &
            + lmesh%Escale(:,ke,3,3) * fz_sp(:)              &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)            &
          - gradphyd_x(:) * rgsqrtv(:)                       &
          + cori(:) * momy_(:,ke)
 
      !-- MOMY
      call dx_abc( dxt1d_, gsqrtdens_, v_,  u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_abc( dyt1d_, gsqrtdens_, v_,  v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )  
      call dz_abc( dzt1d_, gsqrtdens_, v_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * dpres_(:,ke)              , fy)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
 
      momy_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx_sp(:)             &
            + lmesh%Escale(:,ke,2,2) * ( fy_sp(:) + fy(:) ) &
            + lmesh%Escale(:,ke,3,3) * fz_sp(:)             &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)           &
          - gradphyd_y(:) * rgsqrtv(:)                      &
          - cori(:) * momx_(:,ke)
 
      !-- MOMZ
      call dx_abc( dxt1d_, gsqrtdens_, w_,  u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_abc( dyt1d_, gsqrtdens_, w_,  v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
      call dz_abc( dzt1d_, gsqrtdens_, w_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,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_sp(:)             &
          + lmesh%Escale(:,ke,2,2) * fy_sp(:)             &
          + lmesh%Escale(:,ke,3,3) * ( fz_sp(:) + fz(:) ) &
          + liftdelflx(:)                               ) &
          - grav * drho(:)  
 
      !-- RHOT
      call dx_abc( dxt1d_, gsqrtdens_, pot_,      u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_abc( dyt1d_, gsqrtdens_, pot_,      v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )        
      call dz_abc( dzt1d_, gsqrtdens_, pot_,  wt_(:), elem%Nnode_h1D, elem%Nnode_v, fz_sp )        
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
      
      rhot_dt(:,ke) =  - (  &
            lmesh%Escale(:,ke,1,1) * fx_sp(:) &
          + lmesh%Escale(:,ke,2,2) * fy_sp(:) &
          + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
          + liftdelflx(:)                     ) 
    end do
    call prof_rapend( 'cal_dyn_tend_interior', 3)
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_splitform_cal_tend
 
!-----------------------------------------
 
  subroutine dx_ab(DxT1D, a, b, Nnode_h1D, Nnode_v, fx)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dxt1d(nnode_h1d,nnode_h1d)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fx(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fx(i,j,k) = 0.5_rp * sum( dxt1d(:,i) * (a(i,j,k) + a(:,j,k)) * (b(i,j,k) + b(:,j,k)) )
    end do
    end do
    end do
 
    return
  end subroutine dx_ab
 
  subroutine dy_ab(DyT1D, a, b, Nnode_h1D, Nnode_v, fy)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dyt1d(nnode_h1d,nnode_h1d)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fy(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fy(i,j,k) = 0.5_rp * sum( dyt1d(:,j) * (a(i,j,k) + a(i,:,k)) * (b(i,j,k) + b(i,:,k)) )
    end do
    end do
    end do
 
    return
  end subroutine dy_ab
  
  subroutine dz_ab(DzT1D, a, b, Nnode_h1D, Nnode_v, fz)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dzt1d(nnode_v,nnode_v)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fz(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fz(i,j,k) = 0.5_rp * sum( dzt1d(:,k) * (a(i,j,k) + a(i,j,:)) * (b(i,j,k) + b(i,j,:)) )
    end do
    end do
    end do
 
    return
  end subroutine dz_ab
 
  subroutine dx_abc(DxT1D, a, b, c, Nnode_h1D, Nnode_v, fx)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dxt1d(nnode_h1d,nnode_h1d)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: c(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fx(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fx(i,j,k) = 0.25_rp * sum( dxt1d(:,i) * (a(i,j,k) + a(:,j,k)) * (b(i,j,k) + b(:,j,k)) * (c(i,j,k) + c(:,j,k)) )
    end do
    end do
    end do
 
    return
  end subroutine dx_abc
 
  subroutine dy_abc(DyT1D, a, b, c, Nnode_h1D, Nnode_v, fy)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dyt1d(nnode_h1d,nnode_h1d)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: c(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fy(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fy(i,j,k) = 0.25_rp * sum( dyt1d(:,j) * (a(i,j,k) + a(i,:,k)) * (b(i,j,k) + b(i,:,k)) * (c(i,j,k) + c(i,:,k)) )
    end do
    end do
    end do
 
    return
  end subroutine dy_abc
 
  subroutine dz_abc(DzT1D, a, b, c, Nnode_h1D, Nnode_v, fz)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dzt1d(nnode_v,nnode_v)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: c(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fz(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fz(i,j,k) = 0.25_rp * sum( dzt1d(:,k) * (a(i,j,k) + a(i,j,:)) * (b(i,j,k) + b(i,j,:)) * (c(i,j,k) + c(i,j,:)) )
    end do
    end do
    end do
 
    return
  end subroutine dz_abc
 
end module scale_atm_dyn_dgm_nonhydro3d_rhot_heve_splitform