scale_atm_phy_tb_dgm_smg.F90

Go to the documentation of this file.

!> module FElib / Atmosphere / Physics turbulence
!!
!! @par Description
!!      Sub-grid scale turbulence process
!!      Smagorinsky-type
!!
!! @author Yuta Kawai, Team SCALE
!!
!! @par Reference
!!  - Brown et al., 1994:
!!    Large-eddy simulation of stable atmospheric boundary layers with a revised stochastic subgrid model.
!!    Roy. Meteor. Soc., 120, 1485-1512
!!  - Scotti et al., 1993:
!!    Generalized Smagorinsky model for anisotropic grids.
!!    Phys. Fluids A, 5, 2306-2308
!!  - Nishizawa et al., 2015:
!!    Influence of grid aspect ratio on planetary boundary layer turbulence in large-eddy simulations
!!    Geosci. Model Dev., 8, 3393–3419
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_tb_dgm_smg
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    eps  => const_eps,     &
    grav => const_grav,    &
    rdry => const_rdry,    &
    cpdry => const_cpdry,  &
    cvdry => const_cvdry,  &
    pres00 => const_pre00, &
    karman  => const_karman, &
    rplanet => const_radius    
 
  use scale_sparsemat  
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
  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
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_phy_tb_dgm_smg_init
  public :: atm_phy_tb_dgm_smg_final
  public :: atm_phy_tb_dgm_smg_cal_grad
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: cal_del_flux_grad
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private, parameter   :: OneOverThree  = 1.0_rp / 3.0_rp
  real(RP), private, parameter   :: twoOverThree  = 2.0_rp / 3.0_rp
  real(RP), private, parameter   :: FourOverThree = 4.0_rp / 3.0_rp
 
  real(RP), private              :: Cs            = 0.13_rp ! Smagorinsky constant (Scotti et al. 1993)
  real(RP), private, parameter   :: PrN           = 0.7_rp  ! Prandtl number in neutral conditions
  real(RP), private, parameter   :: RiC           = 0.25_rp ! critical Richardson number
  real(RP), private, parameter   :: FmC           = 16.0_rp ! fum = sqrt(1 - c*Ri)
  real(RP), private, parameter   :: FhB           = 40.0_rp ! fuh = sqrt(1 - b*Ri)/PrN
  real(RP), private              :: RPrN                    ! 1 / PrN
  real(RP), private              :: RRiC                    ! 1 / RiC
  real(RP), private              :: OnemPrNovRiC            ! PrN / RiC
  
  ! for backscatter
  real(RP), private, parameter   :: CB   = 1.4_rp
  real(RP), private, parameter   :: CBt  = 0.45_rp
  real(RP), private, parameter   :: aN   = 0.47958315233127197_rp ! a_N = sqrt(0.23)
  real(RP), private, parameter   :: atN4 = 0.09_rp                ! a_{\theta N}^4 = 0.3**2
  real(RP), private, parameter   :: C1o  = an**3
  real(RP), private, parameter   :: D1o  = prn * atn4 / an
 
 
  real(RP), private              :: filter_fac    = 2.0_rp
  real(RP), private              :: NU_MAX        = 10000.0_rp
  real(RP), private              :: tke_fac       
 
contains
!OCL SERIAL
  subroutine atm_phy_tb_dgm_smg_init( mesh )
    implicit none    
    class(meshbase3d), intent(in) :: mesh
 
    logical  :: consistent_tke = .true.
 
    namelist / param_atmos_phy_tb_dgm_smg / &
      cs,                                   &
      nu_max,                               &
      filter_fac,                           &
      consistent_tke
    
    integer :: ierr
    !--------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_TB_dgm_smg_setup",*) 'Setup'
    log_info("ATMOS_PHY_TB_dgm_smg_setup",*) 'Smagorinsky-type Eddy Viscocity Model'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_tb_dgm_smg,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_TB_dgm_smg_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_TB_dgm_smg_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_DGM_SMG. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_tb_dgm_smg)
 
    rprn         = 1.0_rp / prn
    rric         = 1.0_rp / ric
    onemprnovric = ( 1.0_rp - prn ) * rric
 
    if ( consistent_tke ) then
      tke_fac = 1.0_rp
    else
      tke_fac = 0.0_rp
    end if
 
    return
  end subroutine atm_phy_tb_dgm_smg_init
 
!OCL SERIAL
  subroutine atm_phy_tb_dgm_smg_final()
    implicit none
    !--------------------------------------------------------------------
 
    return
  end subroutine atm_phy_tb_dgm_smg_final
 
!> Calculate parameterized stress tensor and eddy heat flux with turbulent model
!!
!OCL SERIAL  
  subroutine atm_phy_tb_dgm_smg_cal_grad( &
    T11, T12, T13, T21, T22, T23, T31, T32, T33,                & ! (out)
    df1, df2, df3,                                              & ! (out)
    tke, nu, kh,                                                & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd,    & ! (in)
    pres, pt,                                                   & ! (in)
    dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d, & ! (in)
    is_bound                                                    ) ! (in)
 
    use scale_atm_phy_tb_dgm_common, only: &
      atm_phy_tb_dgm_common_calc_lambda
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(out) :: t11(elem%np,lmesh%nea)      !< (1,1) component of stress tensor
    real(rp), intent(out) :: t12(elem%np,lmesh%nea)      !< (1,2) component of stress tensor
    real(rp), intent(out) :: t13(elem%np,lmesh%nea)      !< (1,3) component of stress tensor
    real(rp), intent(out) :: t21(elem%np,lmesh%nea)      !< (2,1) component of stress tensor
    real(rp), intent(out) :: t22(elem%np,lmesh%nea)      !< (2,2) component of stress tensor
    real(rp), intent(out) :: t23(elem%np,lmesh%nea)      !< (2,3) component of stress tensor
    real(rp), intent(out) :: t31(elem%np,lmesh%nea)      !< (3,1) component of stress tensor
    real(rp), intent(out) :: t32(elem%np,lmesh%nea)      !< (3,2) component of stress tensor
    real(rp), intent(out) :: t33(elem%np,lmesh%nea)      !< (3,3) component of stress tensor
    real(rp), intent(out)  :: df1(elem%np,lmesh%nea)     !< Diffusive heat flux in x1 direction / density
    real(rp), intent(out)  :: df2(elem%np,lmesh%nea)     !< Diffusive heat flux in x2 direction / density
    real(rp), intent(out)  :: df3(elem%np,lmesh%nea)     !< Diffusive heat flux in x3 direction / density
    real(rp), intent(out) :: tke(elem%np,lmesh%nea)      !< Parameterized turbulent kinetic energy
    real(rp), intent(out) :: nu(elem%np,lmesh%nea)       !< Eddy viscosity
    real(rp), intent(out) :: kh(elem%np,lmesh%nea)       !< Eddy diffusivity
    real(rp), intent(in)  :: ddens_(elem%np,lmesh%nea)   !< Density perturbation
    real(rp), intent(in)  :: momx_ (elem%np,lmesh%nea)   !< Momentum in x1 direction
    real(rp), intent(in)  :: momy_ (elem%np,lmesh%nea)   !< Momentum in x2 direction
    real(rp), intent(in)  :: momz_ (elem%np,lmesh%nea)   !< Momentum in x3 direction
    real(rp), intent(in)  :: drhot_(elem%np,lmesh%nea)   !< Density x potential temperature perturbation
    real(rp), intent(in)  :: dens_hyd(elem%np,lmesh%nea) !< Reference pressure in hydrostatic balance
    real(rp), intent(in)  :: pres_hyd(elem%np,lmesh%nea) !< Reference density in hydrostatic balance
    real(rp), intent(in)  :: pres(elem%np,lmesh%nea)     !< Pressure
    real(rp), intent(in)  :: pt(elem%np,lmesh%nea)       !< Potential temperature
    type(sparsemat), intent(in) :: dx, dy, dz             !< Differential matrix managed by sparse matrix type
    type(sparsemat), intent(in) :: sx, sy, sz             !< Stiffness matrix managed by sparse matrix type
    type(sparsemat), intent(in) :: lift                   !< Lifting matrix managed by sparse matrix type
    logical, intent(in) :: is_bound(elem%nfptot,lmesh%ne) !< Flag whether nodes are located at domain boundaries
 
    real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
    real(rp) :: dens(elem%np), rdens(elem%np), rhot(elem%np), q(elem%np)
    real(rp) :: ddensdxi(elem%np,3)
    real(rp) :: dveldxi(elem%np,3,3)
    real(rp) :: del_flux_rho (elem%nfptot,lmesh%ne,3)
    real(rp) :: del_flux_mom (elem%nfptot,lmesh%ne,3,3)
    real(rp) :: del_flux_rhot(elem%nfptot,lmesh%ne,3)
 
    real(rp) :: s11(elem%np), s12(elem%np), s22(elem%np), s23(elem%np), s31(elem%np), s33(elem%np)
    real(rp) :: skkovthree
    real(rp) :: tkemultwoovthree
    real(rp) :: coef
 
    real(rp) :: ri ! local gradient Richardson number
    real(rp) :: s2 ! (2SijSij)^1/2
    real(rp) :: fm ! factor in eddy viscosity which represents the stability dependence of 
                   ! the Brown et al (1994)'s subgrid model 
    real(rp) :: pr ! Parandtl number (=Nu/Kh= fm/fh)
 
    real(rp) :: lambda  (elem%np,lmesh%ne) ! basic mixing length
    real(rp) :: lambda_r(elem%np)          ! characteristic subgrid length scale 
    real(rp) :: e(elem%np)  ! subgrid kinetic energy 
    real(rp) :: c1(elem%np) ! factor in the relation with energy disspation rate, lambda_r, and E
 
    integer :: ke
    integer :: p
    !--------------------------------------------------------------------
 
    call cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot,        & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, pt,            & ! (in)
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP,                                               & ! (in)
      lmesh, elem, is_bound )                                                   ! (in)
 
    call atm_phy_tb_dgm_common_calc_lambda( lambda, & ! (out)
      cs, filter_fac, lmesh, elem, lmesh2d, elem2d  ) ! (in)
  
    !$omp parallel do private( &
    !$omp Fx, Fy, Fz, LiftDelFlx,                  &
    !$omp DENS, RHOT, RDENS, Q, DdensDxi, DVelDxi, &
    !$omp S11, S12, S22, S23, S31, S33,            &
    !$omp coef, SkkOvThree, TKEMulTwoOvThree,      &
    !$omp p, Ri, S2, fm, Pr, lambda_r, E, C1       )
    do ke=lmesh%NeS, lmesh%NeE
      !---
      dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
      rdens(:) = 1.0_rp / dens(:)
      rhot(:) = dens(:) * pt(:,ke)
 
      ! gradient of density
      call sparsemat_matmul( dx, dens, fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,1), liftdelflx )
      ddensdxi(:,1) = lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:)
 
      call sparsemat_matmul( dy, dens, fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,2), liftdelflx )
      ddensdxi(:,2) = lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:)
 
      call sparsemat_matmul( dz, dens, fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,3), liftdelflx )
      ddensdxi(:,3) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
      
      ! gradient of u
      q(:) = momx_(:,ke) * rdens(:)
 
      call sparsemat_matmul( dx, momx_(:,ke), fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,1), liftdelflx )
      dveldxi(:,1,1) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
 
      call sparsemat_matmul( dy, momx_(:,ke), fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,1), liftdelflx )
      dveldxi(:,2,1) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
 
      call sparsemat_matmul( dz, momx_(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,1), liftdelflx )
      dveldxi(:,3,1) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
 
      ! gradient of v
      q(:) = momy_(:,ke) * rdens(:)
 
      call sparsemat_matmul( dx, momy_(:,ke), fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,2), liftdelflx )
      dveldxi(:,1,2) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
 
      call sparsemat_matmul( dy, momy_(:,ke), fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,2), liftdelflx )
      dveldxi(:,2,2) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
 
      call sparsemat_matmul( dz, momy_(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,2), liftdelflx )
      dveldxi(:,3,2) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
 
      ! gradient of w
      q(:) = momz_(:,ke) * rdens(:)
 
      call sparsemat_matmul( dx, momz_(:,ke), fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,3), liftdelflx )
      dveldxi(:,1,3) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
 
      call sparsemat_matmul( dy, momz_(:,ke), fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,3), liftdelflx )
      dveldxi(:,2,3) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
 
      call sparsemat_matmul( dz, momz_(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,3), liftdelflx )
      dveldxi(:,3,3) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
 
      ! gradient of pt
      q(:) = rhot(:) * rdens(:)
 
      call sparsemat_matmul( dx, rhot, fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,1), liftdelflx )
      df1(:,ke) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
 
      call sparsemat_matmul( dy, rhot, fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,2), liftdelflx )
      df2(:,ke) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
 
      call sparsemat_matmul( dz, rhot, fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,3), liftdelflx )
      df3(:,ke) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
 
 
      ! Calculate the component of strain velocity tensor
      s11(:) = dveldxi(:,1,1)
      s12(:) = 0.5_rp * ( dveldxi(:,1,2) + dveldxi(:,2,1) )
      s22(:) = dveldxi(:,2,2)
      s23(:) = 0.5_rp * ( dveldxi(:,2,3) + dveldxi(:,3,2) )
      s31(:) = 0.5_rp * ( dveldxi(:,1,3) + dveldxi(:,3,1) )
      s33(:) = dveldxi(:,3,3)
 
      ! Caclulate eddy viscosity & eddy diffusivity
      
      do p=1, elem%Np
        s2 = 2.0_rp * ( s11(p)**2 + s22(p)**2 + s33(p)**2 ) &
           + 4.0_rp * ( s31(p)**2 + s12(p)**2 + s23(p)**2 )
        
        ri = grav / pt(p,ke) * df3(p,ke) / max( s2, eps )
 
        ! The Stability functions fm and fh are given by the appendix A of Brown et al. (1994). 
        if (ri < 0.0_rp ) then ! unstable
          fm = sqrt( 1.0_rp - fmc * ri )
          nu(p,ke) = lambda(p,ke)**2 * sqrt( s2 ) * fm
          pr = fm / sqrt( 1.0_rp - fhb * ri ) * prn
        else if ( ri < ric ) then ! stable
          fm = ( 1.0_rp - ri * rric )**4
          nu(p,ke) = lambda(p,ke)**2 * sqrt( s2 ) * fm
          pr = prn / ( 1.0_rp - onemprnovric * ri )
        else ! strongly stable
          fm = 0.0_rp
          nu(p,ke) = 0.0_rp
          kh(p,ke) = 0.0_rp
          pr = 1.0_rp
        end if
 
        if ( ri < ric ) then
          kh(p,ke) = max( min( nu(p,ke) / pr, nu_max ), eps )
          nu(p,ke) = max( min( nu(p,ke), nu_max ), eps )
          pr = nu(p,ke) / kh(p,ke)
          lambda_r(p) = lambda(p,ke) * sqrt( fm / sqrt( 1.0_rp - ri/pr ) )
        else
          lambda_r(p) = 0.0_rp
        end if
      end do
      
!      Nu(:,ke) = 0.0_RP
 
      ! if ( backscatter ) then
      ! else 
        e(:) = nu(:,ke)**3 / ( lambda_r(:)**4 + eps )
        c1(:) = c1o
      ! end if
 
      ! TKE
      tke(:,ke) = ( e(:) * lambda_r(:) / c1(:) )**twooverthree
 
      !---
 
      do p=1, elem%Np
        tkemultwoovthree = twooverthree * tke(p,ke) * tke_fac
        skkovthree = ( s11(p) + s22(p) + s33(p) ) * oneoverthree
        coef = 2.0_rp * nu(p,ke)
 
        t11(p,ke) = dens(p) * ( coef * ( s11(p) - skkovthree ) - tkemultwoovthree )
        t12(p,ke) = dens(p) * coef * s12(p)
        t13(p,ke) = dens(p) * coef * s31(p)
 
        t21(p,ke) = dens(p) * coef * s12(p)
        t22(p,ke) = dens(p) * ( coef * ( s22(p) - skkovthree ) - tkemultwoovthree )
        t23(p,ke) = dens(p) * coef * s23(p)
 
        t31(p,ke) = dens(p) * coef * s31(p)
        t32(p,ke) = dens(p) * coef * s23(p)
        t33(p,ke) = dens(p) * ( coef * ( s33(p) - skkovthree ) - tkemultwoovthree )
      end do
 
      df1(:,ke) = kh(:,ke) * df1(:,ke)
      df2(:,ke) = kh(:,ke) * df2(:,ke)
      df3(:,ke) = kh(:,ke) * df3(:,ke)
    end do
  
    return
  end subroutine atm_phy_tb_dgm_smg_cal_grad
 
!-- private --------------------------------------------------------
 
!OCL SERIAL  
  subroutine cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot,  & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, pt_,           & ! (in)
    nx, ny, nz, vmapm, vmapp, lmesh, elem, is_bound                         ) ! (in)
 
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem  
    real(rp), intent(out) ::  del_flux_rho(elem%nfptot*lmesh%ne,3)
    real(rp), intent(out) ::  del_flux_mom(elem%nfptot*lmesh%ne,3,3)
    real(rp), intent(out) ::  del_flux_rhot(elem%nfptot*lmesh%ne,3)
    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) ::  dens_hyd(elem%np*lmesh%nea)
    real(rp), intent(in) ::  pres_hyd(elem%np*lmesh%nea)
    real(rp), intent(in) ::  pt_(elem%np*lmesh%nea)
    real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
    real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
    real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
    logical, intent(in) :: is_bound(elem%nfptot*lmesh%ne)
    
    integer :: i, ip, im
    real(rp) :: densm, densp
    real(rp) :: del
    real(rp) :: facx, facy, facz
    !------------------------------------------------------------------------
    
    !$omp parallel do private ( iM, iP,   &
    !$omp densM, densP,                   &
    !$omp del, facx, facy, facz           )
    do i=1, elem%NfpTot * lmesh%Ne
      im = vmapm(i); ip = vmapp(i)
 
      densm = ddens_(im) + dens_hyd(im)
      densp = ddens_(ip) + dens_hyd(ip)
 
      if ( is_bound(i) ) then
        facx = 1.0_rp
        facy = 1.0_rp 
        facz = 1.0_rp
      else
        ! facx = 1.0_RP - sign(1.0_RP,nx(i))
        ! facy = 1.0_RP - sign(1.0_RP,ny(i))
        ! facz = 1.0_RP - sign(1.0_RP,nz(i))
        facx = 1.0_rp
        facy = 1.0_rp
        facz = 1.0_rp
      end if
 
      del = 0.5_rp * ( densp - densm )
      del_flux_rho(i,1) = facx * del * nx(i)
      del_flux_rho(i,2) = facy * del * ny(i)
      del_flux_rho(i,3) = facz * del * nz(i)
 
      del = 0.5_rp * ( momx_(ip) - momx_(im) )
      del_flux_mom(i,1,1) = facx * del * nx(i)
      del_flux_mom(i,2,1) = facy * del * ny(i)
      del_flux_mom(i,3,1) = facz * del * nz(i)
 
      del = 0.5_rp * ( momy_(ip) - momy_(im) )
      del_flux_mom(i,1,2) = facx * del * nx(i)
      del_flux_mom(i,2,2) = facy * del * ny(i)
      del_flux_mom(i,3,2) = facz * del * nz(i)
 
      del = 0.5_rp * ( momz_(ip) - momz_(im) )
      del_flux_mom(i,1,3) = facx * del * nx(i)
      del_flux_mom(i,2,3) = facy * del * ny(i)
      del_flux_mom(i,3,3) = facz * del * nz(i)
 
      del = 0.5_rp * ( densp * pt_(ip) - densm * pt_(im) )
      del_flux_rhot(i,1) = facx * del * nx(i)
      del_flux_rhot(i,2) = facy * del * ny(i)
      del_flux_rhot(i,3) = facz * del * nz(i)
    end do
 
    return
  end subroutine cal_del_flux_grad
 
end module scale_atm_phy_tb_dgm_smg