scale_atm_phy_tb_dgm_common.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / Physics turbulence / Common
!!
!! @par Description
!!      A common modules for atmospheric turbulent parameterization
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_tb_dgm_common
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    eps  => const_eps
  
  use scale_sparsemat
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  use scale_localmesh_2d, only: localmesh2d  
  use scale_localmesh_3d, only: localmesh3d
 
  use scale_model_var_manager, only: &
    modelvarmanager, variableinfo
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_phy_tb_dgm_common_setup_variables
  public :: atm_phy_tb_dgm_common_calc_lambda
 
  public :: atm_phy_tb_dgm_common_cal_grad_qtrc
  public :: atm_phy_tb_dgm_common_cal_tend
  public :: atm_phy_tb_dgm_common_cal_tend_qtrc
  
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public, parameter :: atmos_phy_tb_aux_t13_id      = 1
  integer, public, parameter :: atmos_phy_tb_aux_t23_id      = 2
  integer, public, parameter :: atmos_phy_tb_aux_t33_id      = 3
  integer, public, parameter :: atmos_phy_tb_aux_diffflx3_id = 4
  integer, public, parameter :: atmos_phy_tb_aux_t31_id      = 5
  integer, public, parameter :: atmos_phy_tb_aux_t32_id      = 6
  integer, public, parameter :: atmos_phy_tb_aux_diffflx1_id = 7
  integer, public, parameter :: atmos_phy_tb_aux_diffflx2_id = 8
  integer, public, parameter :: atmos_phy_tb_aux_t11_id      = 9
  integer, public, parameter :: atmos_phy_tb_aux_t12_id      = 10
  integer, public, parameter :: atmos_phy_tb_aux_t21_id      = 11
  integer, public, parameter :: atmos_phy_tb_aux_t22_id      = 12
  integer, public, parameter :: atmos_phy_tb_aux_scalar_num  = 4
  integer, public, parameter :: atmos_phy_tb_aux_hvec_num    = 2
  integer, public, parameter :: atmos_phy_tb_aux_htensor_num = 1
  integer, public, parameter :: atmos_phy_tb_aux_num         = 12
  
  integer, public, parameter :: atmos_phy_tb_diag_tke_id     = 1
  integer, public, parameter :: atmos_phy_tb_diag_nu_id      = 2
  integer, public, parameter :: atmos_phy_tb_diag_kh_id      = 3
  integer, public, parameter :: atmos_phy_tb_diag_num        = 3
 
  integer, public, parameter :: atmos_phy_tb_momx_t_id  = 1
  integer, public, parameter :: atmos_phy_tb_momy_t_id  = 2
  integer, public, parameter :: atmos_phy_tb_momz_t_id  = 3
  integer, public, parameter :: atmos_phy_tb_rhot_t_id  = 4
  integer, public, parameter :: atmos_phy_tb_tends_num1 = 4
  !-
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
  
  private :: cal_bnd_flux
  private :: cal_bnd_flux_grad_qtrc
  private :: cal_bnd_flux_qtrc
 
  private :: fact, p1, p2, p3
  
contains
!OCL SERIAL
  subroutine atm_phy_tb_dgm_common_setup_variables( &
    tends, auxvars, diagvars,                       & ! (inout)
    tends_manager, auxvar_manager, diagvar_manager, & ! (inout)
    tends_num_tot, mesh3d                           ) ! (in)
    use scale_tracer, only: &
      qa, tracer_name, tracer_desc, tracer_unit    
    use scale_mesh_base3d, only: meshbase3d
    use scale_meshfield_base, only: meshfield3d
    implicit none
    integer, intent(in) :: tends_num_tot
    type(meshfield3d), intent(inout) :: tends(tends_num_tot)
    type(meshfield3d), intent(inout) :: auxvars(atmos_phy_tb_aux_num)
    type(meshfield3d), intent(inout) :: diagvars(atmos_phy_tb_diag_num)
    type(modelvarmanager), intent(inout) :: tends_manager
    type(modelvarmanager), intent(inout) :: auxvar_manager
    type(modelvarmanager), intent(inout) :: diagvar_manager
    class(meshbase3d), intent(in) :: mesh3d
 
    ! type(VariableInfo), intent(out) :: auxvar_info(ATMOS_PHY_TB_AUX_NUM)
    ! type(VariableInfo), intent(out) :: diagvar_info(ATMOS_PHY_TB_DIAG_NUM)
    ! type(VariableInfo), intent(out) :: tend_info(ATMOS_PHY_TB_TENDS_NUM1)
 
    type(variableinfo) :: atmos_phy_tb_aux_vinfo(atmos_phy_tb_aux_num)
    DATA atmos_phy_tb_aux_vinfo / &
      variableinfo( atmos_phy_tb_aux_t13_id, 'TB_T13', 'stress tensor (T13)',                         & ! rho x nu x S : [kg/m3] x [m2/s] x [s-1]
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                ), &
      variableinfo( atmos_phy_tb_aux_t23_id, 'TB_T23', 'stress tensor (T23)',                         &
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                ), &
      variableinfo( atmos_phy_tb_aux_t33_id, 'TB_T33', 'stress tensor (T33)',                         &
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                ), &
      variableinfo( atmos_phy_tb_aux_diffflx3_id, 'DIFF_FLX3', 'diffusive heat flux (z) / density',   &
                    'm2/s.K/m',  3, 'XYZ',  ''                                                     ), &                                      
      variableinfo( atmos_phy_tb_aux_t31_id, 'TB_T31', 'stress tensor (T31)',                         &
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                ), &
      variableinfo( atmos_phy_tb_aux_t32_id, 'TB_T32', 'stress tensor (T32)',                         &
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                ), &
      variableinfo( atmos_phy_tb_aux_diffflx1_id, 'DIFF_FLX1', 'diffusive heat flux (x) / density',   &
                    'm2/s.K/m',  3, 'XYZ',  ''                                                     ), &
      variableinfo( atmos_phy_tb_aux_diffflx2_id, 'DIFF_FLX2', 'diffusive heat flux (y) / density',   &
                    'm2/s.K/m',  3, 'XYZ',  ''                                                     ), &
      variableinfo( atmos_phy_tb_aux_t11_id, 'TB_T11', 'stress tensor (T11)',                         &
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                ), &
      variableinfo( atmos_phy_tb_aux_t12_id, 'TB_T12', 'stress tensor (T12)',                         &
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                ), &
      variableinfo( atmos_phy_tb_aux_t21_id, 'TB_T21', 'stress tensor (T21)',                         &
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                ), &
      variableinfo( atmos_phy_tb_aux_t22_id, 'TB_T22', 'stress tensor (T22)',                         &
                    'kg.m-3.m2.s-2',  3, 'XYZ',  ''                                                )  /
 
    type(variableinfo) :: atmos_phy_tb_diag_vinfo(atmos_phy_tb_diag_num)
    DATA atmos_phy_tb_diag_vinfo / &
      variableinfo( atmos_phy_tb_diag_tke_id, 'TKE', 'SGS turbluence kinetic energy',   &
                    'm2/s2',  3, 'XYZ',  ''                                          ), &
      variableinfo( atmos_phy_tb_diag_nu_id, 'NU', 'eddy viscosity',                    &
                    'm2/s',  3, 'XYZ',  ''                                           ), &
      variableinfo( atmos_phy_tb_diag_kh_id, 'KH', 'eddy diffusion',                    &
                    'm2/s',  3, 'XYZ',  ''                                           )  /
                    
    type(variableinfo) :: atmos_phy_tb_tend_vinfo(atmos_phy_tb_tends_num1)
    DATA atmos_phy_tb_tend_vinfo / &
      variableinfo( atmos_phy_tb_momx_t_id, 'TB_MOMX_t', 'tendency of x-momentum in TB process',    &
                    'kg/m2/s2',  3, 'XYZ',  ''                                                   ), &
      variableinfo( atmos_phy_tb_momy_t_id, 'TB_MOMY_t', 'tendency of y-momentum in TB process',    &
                    'kg/m2/s2',  3, 'XYZ',  ''                                                   ), &
      variableinfo( atmos_phy_tb_momz_t_id, 'TB_MOMZ_t', 'tendency of z-momentum in TB process',    &
                    'kg/m2/s2',  3, 'XYZ',  ''                                                   ), &
      variableinfo( atmos_phy_tb_rhot_t_id, 'TB_RHOT_t', 'tendency of rho*PT in TB process',        &
                    'kg/m3.K/s', 3, 'XYZ',  ''                                                   )  / 
 
    type(variableinfo) :: qtrc_vinfo_tmp
 
    integer :: iv
    integer :: iq
    logical :: reg_file_hist
    !----------------------------------------------------------
 
    ! auxvar_info(:) = ATMOS_PHY_TB_AUX_VINFO
    ! diagvar_info(:) = ATMOS_PHY_TB_DIAG_VINFO
    ! tend_info(:) = ATMOS_PHY_TB_TEND_VINFO
 
    !- Initialize an object to manage tendencies with turbulent model
 
    reg_file_hist = .true.    
    do iv = 1, atmos_phy_tb_tends_num1
      call tends_manager%Regist( &
        atmos_phy_tb_tend_vinfo(iv), mesh3d,     & ! (in) 
        tends(iv), reg_file_hist,                & ! (out)
        fill_zero=.true.                         )
    end do
 
    qtrc_vinfo_tmp%ndims    = 3
    qtrc_vinfo_tmp%dim_type = 'XYZ'
    qtrc_vinfo_tmp%STDNAME  = ''
 
    do iq = 1, qa
      iv = atmos_phy_tb_tends_num1 + iq 
      qtrc_vinfo_tmp%keyID = iv
      qtrc_vinfo_tmp%NAME  = 'TB_'//trim(tracer_name(iq))//'_t'
      qtrc_vinfo_tmp%DESC  = 'tendency of '//trim(tracer_desc(iq))//' in TB process'
      qtrc_vinfo_tmp%UNIT  = trim(tracer_unit(iq))//'/s'
 
      call tends_manager%Regist( &
        qtrc_vinfo_tmp, mesh3d,                  & ! (in) 
        tends(iv), reg_file_hist,                & ! (out)
        fill_zero=.true.                         )
    end do
 
    !- Initialize an object to manage auxiliary variables with turbulent model
 
    reg_file_hist = .true.    
    do iv = 1, atmos_phy_tb_aux_num
      call auxvar_manager%Regist( &
        atmos_phy_tb_aux_vinfo(iv), mesh3d, & ! (in) 
        auxvars(iv), reg_file_hist,         & ! (out)
        fill_zero=.true.                    )
    end do
 
    !- Initialize an object to manage diagnostic variables with turbulent model
    
    reg_file_hist = .true.    
    do iv = 1, atmos_phy_tb_diag_num
      call diagvar_manager%Regist( &
        atmos_phy_tb_diag_vinfo(iv), mesh3d, & ! (in) 
        diagvars(iv), reg_file_hist,         & ! (out)
        fill_zero=.true.                     )
    end do
 
    return
  end subroutine atm_phy_tb_dgm_common_setup_variables
 
!OCL SERIAL  
  subroutine atm_phy_tb_dgm_common_calc_lambda( lambda, & ! (out)
    cs, filter_fac, lmesh, elem, lmesh2d, elem2d        ) ! (in)
 
    use scale_const, only: &
      karman  => const_karman
    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) :: lambda(elem%np,lmesh%ne)
    real(rp), intent(in) :: cs
    real(rp), intent(in) :: filter_fac
 
    integer :: ke
 
    real(rp) :: vol
    real(rp) :: he
    real(rp) :: lambda0
    real(rp) :: zs(elem2d%np)
    real(rp) :: z1(elem2d%np)
    real(rp) :: dz(elem  %np)
    real(rp) :: fz
    real(rp) :: elem_aspect_eff
 
    real(rp), parameter :: oneoverthree  = 1.0_rp / 3.0_rp
    !--------------------------------------------------------------------
 
    !$omp parallel do private( &
    !$omp vol, lambda0, Zs, Z1, dz, he, FZ, elem_aspect_eff )
    do ke=lmesh%NeS, lmesh%NeE
      vol = sum( elem%IntWeight_lgl(:) * lmesh%J(:,ke) * lmesh%Gsqrt(:,ke) )
      he = ( vol / ( dble(elem%PolyOrder_h+1)**2 * dble(elem%PolyOrder_v+1) ) )**oneoverthree
 
      fz = lmesh%zlev(elem%Colmask(elem%Nnode_v,1),ke) - lmesh%zlev(elem%Colmask(1,1),ke)      
      elem_aspect_eff = fact( fz/dble(elem%PolyOrder_v+1), sqrt(vol/fz)/dble(elem%PolyOrder_h+1), sqrt(vol/fz)/dble(elem%PolyOrder_h+1) )
      lambda0 = elem_aspect_eff * cs * filter_fac * he
 
      zs(:) = lmesh%zlev(elem%Hslice(:,1),lmesh%EMap3Dto2D(ke))
      z1(:) = lmesh%zlev(elem%Hslice(:,2),lmesh%EMap3Dto2D(ke))
      dz(:) = max( lmesh%zlev(:,ke) - zs(elem%IndexH2Dto3D(:)), he )
 
      !lambda(:,ke) = sqrt( 1.0_RP / (1.0_RP / lambda0**2 + 1.0_RP / ( KARMAN * max( dz(:), EPS ) )**2 ) )
      lambda(:,ke) = sqrt( 1.0_rp / (1.0_rp / lambda0**2 + 1.0_rp / ( karman * ( dz(:) + 1.0e-4_rp ) )**2 ) )    
    end do
 
    return
  end subroutine atm_phy_tb_dgm_common_calc_lambda
 
 
!> Calculate parameterized diffusive mass flux of tracer with turbulent model
!!
!OCL SERIAL  
  subroutine atm_phy_tb_dgm_common_cal_grad_qtrc( &
    DFQ1, DFQ2, DFQ3,                                           & ! (out)
    drdx, drdy, drdz,                                           & ! (inout)
    kh, qtrc, ddens, dens_hyd,                                  & ! (in)
    dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d, & ! (in)
    is_bound, cal_grad_dens                                     ) ! (in)
 
    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) :: dfq1(elem%np,lmesh%nea)      !< Diffusive mass flux in x1 direction / density (Kh dq/dx1)
    real(rp), intent(out) :: dfq2(elem%np,lmesh%nea)      !< Diffusive mass flux in x2 direction / density (Kh dq/dx2)
    real(rp), intent(out) :: dfq3(elem%np,lmesh%nea)      !< Diffusive mass flux in x3 direction / density (Kh dq/dx3)
    real(rp), intent(inout) :: drdx(elem%np,lmesh%nea)    !< Spatial gradient of density in x1 direction
    real(rp), intent(inout) :: drdy(elem%np,lmesh%nea)    !< Spatial gradient of density in x2 direction
    real(rp), intent(inout) :: drdz(elem%np,lmesh%nea)    !< Spatial gradient of density in x3 direction
    real(rp), intent(in)  :: kh(elem%np,lmesh%nea)        !< Eddy diffusivity
    real(rp), intent(in)  :: qtrc(elem%np,lmesh%nea)      !< Mass faction of tracer
    real(rp), intent(in)  :: ddens(elem%np,lmesh%nea)     !< Density perturbation
    real(rp), intent(in)  :: dens_hyd(elem%np,lmesh%nea)  !< Reference desity in hydrostatic state
    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
    logical, intent(in) :: cal_grad_dens                  !< Flag whether spatial gradients of density are calcuated
 
    real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
    real(rp) :: bnd_flux(elem%nfptot,lmesh%ne,3)
    real(rp) :: bnd_flux_rho (elem%nfptot,lmesh%ne,3)    
    
    real(rp) :: dens(elem%np), rdens(elem%np)
    real(rp) :: rhoxqtrc(elem%np)
 
    integer :: ke
    integer :: p
    !--------------------------------------------------------------------
 
    call cal_bnd_flux_grad_qtrc( bnd_flux, bnd_flux_rho,                      & ! (out)
      qtrc, ddens, dens_hyd,                                                  & ! (in)
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP,                                               & ! (in)
      lmesh, elem, is_bound, cal_grad_dens )                                    ! (in)
 
    !$omp parallel private( ke, Fx, Fy, Fz, LiftDelFlx, RHOxQTRC, DENS, RDENS )
 
    ! Calculate gradient of density
    if ( cal_grad_dens ) then
      !$omp do
      do ke=lmesh%NeS, lmesh%NeE
        dens(:) = dens_hyd(:,ke) + ddens(:,ke)
 
        call sparsemat_matmul( dx, dens, fx )
        call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux_rho(:,ke,1), liftdelflx )
        drdx(:,ke) = lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:)
  
        call sparsemat_matmul( dy, dens, fy )
        call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux_rho(:,ke,2), liftdelflx )
        drdy(:,ke) = lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:)
  
        call sparsemat_matmul( dz, dens, fz )
        call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux_rho(:,ke,3), liftdelflx )
        drdz(:,ke) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)        
      end do
    end if
 
    ! Calculate gradient of tracer
    !$omp do
    do ke=lmesh%NeS, lmesh%NeE
      dens(:) = dens_hyd(:,ke) + ddens(:,ke)
      rdens(:) = 1.0_rp / dens(:)
      rhoxqtrc(:) = dens(:) * qtrc(:,ke)
 
      !---
      call sparsemat_matmul( dx, rhoxqtrc(:), fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux(:,ke,1), liftdelflx )
      dfq1(:,ke) = kh(:,ke) * ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - qtrc(:,ke) * drdx(:,ke) ) * rdens(:)
 
      call sparsemat_matmul( dy, rhoxqtrc(:), fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux(:,ke,2), liftdelflx )
      dfq2(:,ke) = kh(:,ke) * ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - qtrc(:,ke) * drdy(:,ke) ) * rdens(:)
 
      call sparsemat_matmul( dz, rhoxqtrc(:), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux(:,ke,3), liftdelflx )
      dfq3(:,ke) = kh(:,ke) * ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - qtrc(:,ke) * drdz(:,ke) ) * rdens(:)
    end do
 
    !$omp end parallel
 
    return
  end subroutine atm_phy_tb_dgm_common_cal_grad_qtrc
 
!> Calculate tendecies with turbulent model
!!
!OCL SERIAL  
  subroutine atm_phy_tb_dgm_common_cal_tend( &
    MOMX_t, MOMY_t, MOMZ_t, RHOT_t,                             & ! (out)
    t11, t12, t13, t21, t22, t23, t31, t32, t33,                & ! (in)
    df1, df2, df3,                                              & ! (in)
    nu, kh,                                                     & ! (in)
    ddens_, momx_, momy_, momz_, drhot_,                        & ! (in)
    dens_hyd, pres_hyd,  pres_, pt_,                            & ! (in)
    dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d, & ! (in)
    is_bound                                                    ) ! (in)
 
    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) :: momx_t(elem%np,lmesh%nea)   !< Tendency of momentum in x1 direction with turbulent model
    real(rp), intent(out) :: momy_t(elem%np,lmesh%nea)   !< Tendency of momentum in x2 direction with turbulent model
    real(rp), intent(out) :: momz_t(elem%np,lmesh%nea)   !< Tendency of momentum in x3 direction with turbulent model
    real(rp), intent(out) :: rhot_t(elem%np,lmesh%nea)   !< Tendency of density x potential temperature with turbulent model
    real(rp), intent(in) :: t11(elem%np,lmesh%nea)       !< (1,1) component of stress tensor
    real(rp), intent(in) :: t12(elem%np,lmesh%nea)       !< (1,2) component of stress tensor
    real(rp), intent(in) :: t13(elem%np,lmesh%nea)       !< (1,3) component of stress tensor
    real(rp), intent(in) :: t21(elem%np,lmesh%nea)       !< (2,1) component of stress tensor
    real(rp), intent(in) :: t22(elem%np,lmesh%nea)       !< (2,2) component of stress tensor
    real(rp), intent(in) :: t23(elem%np,lmesh%nea)       !< (2,3) component of stress tensor
    real(rp), intent(in) :: t31(elem%np,lmesh%nea)       !< (3,1) component of stress tensor
    real(rp), intent(in) :: t32(elem%np,lmesh%nea)       !< (3,2) component of stress tensor
    real(rp), intent(in) :: t33(elem%np,lmesh%nea)       !< (3,3) component of stress tensor
    real(rp), intent(in)  :: df1(elem%np,lmesh%nea)      !< Diffusive heat flux in x1 direction / density
    real(rp), intent(in)  :: df2(elem%np,lmesh%nea)      !< Diffusive heat flux in x2 direction / density
    real(rp), intent(in)  :: df3(elem%np,lmesh%nea)      !< Diffusive heat flux in x3 direction / density
    real(rp), intent(in)  :: nu   (elem%np,lmesh%nea)    !< Eddy viscosity
    real(rp), intent(in)  :: 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
 
    integer :: ke
 
    real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
    real(rp) :: dens(elem%np)
    real(rp) :: bnd_flux_mom(elem%nfptot,lmesh%ne,3)
    real(rp) :: bnd_flux_rhot(elem%nfptot,lmesh%ne)    
    !--------------------------------------------------------------------
 
    call cal_bnd_flux( bnd_flux_mom, bnd_flux_rhot,                           & ! (out)
      t11, t12, t13, t21, t22, t23, t31, t32, t33,                            & ! (in)
      df1, df2, df3,                                                          & ! (in)
      nu, kh,                                                                 & ! (in)
      ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd,                & ! (in)
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP, lmesh, elem, is_bound   )                       ! (in)
 
    !$omp parallel do private( &
    !$omp Fx, Fy, Fz, LiftDelFlx,                    &
    !$omp DENS                                       )
    do ke=lmesh%NeS, lmesh%NeE
      dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
 
      ! MOMX
      call sparsemat_matmul( dx, t11(:,ke), fx )
      call sparsemat_matmul( dy, t12(:,ke), fy )
      call sparsemat_matmul( dz, t13(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux_mom(:,ke,1), liftdelflx )
 
      momx_t(:,ke) = lmesh%Escale(:,ke,1,1) * fx(:) + lmesh%Escale(:,ke,2,2) * fy(:) &
                   + lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
 
      ! MOMY
      call sparsemat_matmul( dx, t21(:,ke), fx )
      call sparsemat_matmul( dy, t22(:,ke), fy )
      call sparsemat_matmul( dz, t23(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux_mom(:,ke,2), liftdelflx )
 
      momy_t(:,ke) = lmesh%Escale(:,ke,1,1) * fx(:) + lmesh%Escale(:,ke,2,2) * fy(:) &
                   + lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
 
      ! MOMZ
      call sparsemat_matmul( dx, t31(:,ke), fx )
      call sparsemat_matmul( dy, t32(:,ke), fy )
      call sparsemat_matmul( dz, t33(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux_mom(:,ke,3), liftdelflx )
 
      momz_t(:,ke) = lmesh%Escale(:,ke,1,1) * fx(:) + lmesh%Escale(:,ke,2,2) * fy(:) &
                   + lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
 
      ! RHOT
      call sparsemat_matmul( dx, dens(:) * df1(:,ke), fx )
      call sparsemat_matmul( dy, dens(:) * df2(:,ke), fy )
      call sparsemat_matmul( dz, dens(:) * df3(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * bnd_flux_rhot(:,ke), liftdelflx )
 
      rhot_t(:,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
 
    return
  end subroutine atm_phy_tb_dgm_common_cal_tend
 
!> Calculate tendecies of tracer density with turbulent model
!!
!OCL SERIAL  
  subroutine atm_phy_tb_dgm_common_cal_tend_qtrc( &
    RHOQ_t,                                                     & ! (out)
    dfq1, dfq2, dfq3,                                           & ! (in)
    kh, ddens_,dens_hyd,                                        & ! (in)
    dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d, & ! (in)
    is_bound                                                    ) ! (in)
 
    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) :: rhoq_t(elem%np,lmesh%nea)    !< Tendency of tracer mass fraction
    real(rp), intent(in)  :: dfq1(elem%np,lmesh%nea)      !< Diffusive mass flux in x1 direction / density (Kh dq/dx1)
    real(rp), intent(in)  :: dfq2(elem%np,lmesh%nea)      !< Diffusive mass flux in x2 direction / density (Kh dq/dx2)
    real(rp), intent(in)  :: dfq3(elem%np,lmesh%nea)      !< Diffusive mass flux in x3 direction / density (Kh dq/dx3)
    real(rp), intent(in)  :: kh   (elem%np,lmesh%nea)     !< Eddy diffusivity
    real(rp), intent(in)  :: ddens_(elem%np,lmesh%nea)    !< Density perturbation
    real(rp), intent(in)  :: dens_hyd(elem%np,lmesh%nea)  !< Reference pressure in hydrostatic state
    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
 
    integer :: ke
 
    real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
    real(rp) :: dens(elem%np), rhot(elem%np)
    real(rp) :: del_flux(elem%nfptot,lmesh%ne)
    !--------------------------------------------------------------------
 
    call cal_bnd_flux_qtrc( del_flux,                                         & ! (out)
      dfq1, dfq2, dfq3,                                                       & ! (in)
      kh, ddens_, dens_hyd,                                                   & ! (in)
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP,                                               & ! (in)
      lmesh, elem, is_bound )                                                   ! (in)
 
    !$omp parallel do private( &
    !$omp Fx, Fy, Fz, LiftDelFlx,            &
    !$omp DENS                               )
    do ke=lmesh%NeS, lmesh%NeE
      dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
 
      ! RHOQ
      call sparsemat_matmul( dx, dens(:) * dfq1(:,ke), fx )
      call sparsemat_matmul( dy, dens(:) * dfq2(:,ke), fy )
      call sparsemat_matmul( dz, dens(:) * dfq3(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux(:,ke), liftdelflx )
 
      rhoq_t(:,ke) = lmesh%Escale(:,ke,1,1) * fx(:) + lmesh%Escale(:,ke,2,2) * fy(:) &
                   + lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
    end do
 
    return
  end subroutine atm_phy_tb_dgm_common_cal_tend_qtrc
 
!-- private ---------------------------------------------
 
 
!OCL SERIAL  
  subroutine cal_bnd_flux_grad_qtrc( bnd_flux, bnd_flux_rho, & ! (out)
    qtrc_, ddens_, dens_hyd_,                                & ! (in)
    nx, ny, nz, vmapm, vmapp, lmesh, elem, is_bound,         & ! (in)
    cal_grad_dens                                            ) ! (in)
 
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem  
    real(rp), intent(out) ::  bnd_flux(elem%nfptot*lmesh%ne,3)
    real(rp), intent(out) ::  bnd_flux_rho(elem%nfptot*lmesh%ne,3)
    real(rp), intent(in) ::  qtrc_(elem%np*lmesh%nea)
    real(rp), intent(in) ::  ddens_(elem%np*lmesh%nea)
    real(rp), intent(in) ::  dens_hyd_(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)
    logical, intent(in) :: cal_grad_dens
    
    integer :: i, ip, im
    real(rp) :: del
    real(rp) :: facx, facy, facz
    real(rp) :: densm, densp
 
    !------------------------------------------------------------------------
    
    !$omp parallel do private ( iM, iP,       &
    !$omp del, facx, facy, facz, densM, densP )
    do i=1, elem%NfpTot * lmesh%Ne
      im = vmapm(i); ip = vmapp(i)
 
 
      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
 
      densm = ddens_(im) + dens_hyd_(im)
      densp = ddens_(ip) + dens_hyd_(ip)
 
      if ( cal_grad_dens ) then
        del = 0.5_rp * ( densp - densm )
        ! del = 0.5_RP * ( sqrt_DENS_Kh_P - sqrt_DENS_Kh_M )
        bnd_flux_rho(i,1) = facx * del * nx(i)
        bnd_flux_rho(i,2) = facy * del * ny(i)
        bnd_flux_rho(i,3) = facz * del * nz(i)
      end if
      
      del = 0.5_rp * ( densp * qtrc_(ip) - densm * qtrc_(im) )
      bnd_flux(i,1) = facx * del * nx(i)
      bnd_flux(i,2) = facy * del * ny(i)
      bnd_flux(i,3) = facz * del * nz(i)
    end do
 
    return
  end subroutine cal_bnd_flux_grad_qtrc
 
!OCL SERIAL  
  subroutine cal_bnd_flux( bnd_flux_mom, bnd_flux_rhot,                & ! (out)
    t11, t12, t13, t21, t22, t23, t31, t32, t33,                       & ! (in)
    df1, df2, df3,                                                     & ! (in)
    nu, kh,                                                            & ! (in)
    ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd,           & ! (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) ::  bnd_flux_mom (elem%nfptot*lmesh%ne,3)
    real(rp), intent(out) ::  bnd_flux_rhot(elem%nfptot*lmesh%ne) 
    real(rp), intent(in)  :: t11(elem%np*lmesh%nea), t12(elem%np*lmesh%nea), t13(elem%np*lmesh%nea)
    real(rp), intent(in)  :: t21(elem%np*lmesh%nea), t22(elem%np*lmesh%nea), t23(elem%np*lmesh%nea)
    real(rp), intent(in)  :: t31(elem%np*lmesh%nea), t32(elem%np*lmesh%nea), t33(elem%np*lmesh%nea)
    real(rp), intent(in)  :: df1(elem%np*lmesh%nea)
    real(rp), intent(in)  :: df2(elem%np*lmesh%nea)
    real(rp), intent(in)  :: df3(elem%np*lmesh%nea)
    real(rp), intent(in)  :: nu   (elem%np*lmesh%nea) ! Eddy viscosity
    real(rp), intent(in)  :: kh   (elem%np*lmesh%nea) ! Eddy diffusivity
    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) :: 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) :: taum_x, taup_x
    real(rp) :: taum_y, taup_y
    real(rp) :: taum_z, taup_z
    real(rp) :: nx_, ny_, nz_
    !------------------------------------------------------------------------
 
    !$omp parallel do private( iM, iP, &
    !$omp densM, densP, TauM_x, TauP_x, TauM_y, TauP_y, TauM_z, TauP_z, nx_, ny_, nz_  )
    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 ( ip > elem%Np * lmesh%Ne .and. abs(nz(i)) > eps ) then ! Tentative implementation for the treatmnet of lower/upper boundary. 
        nx_ = nx(i)
        ny_ = ny(i)
        nz_ = nz(i)
      else
        ! nx_ = ( 1.0_RP + sign(1.0_RP,nx(i)) ) * nx(i)
        ! ny_ = ( 1.0_RP + sign(1.0_RP,ny(i)) ) * ny(i)
        ! nz_ = ( 1.0_RP + sign(1.0_RP,nz(i)) ) * nz(i)
        nx_ = nx(i)
        ny_ = ny(i)
        nz_ = nz(i)
      end if
 
      taum_x = t11(im) * nx_ + t12(im) * ny_ + t13(im) * nz_
      taup_x = t11(ip) * nx_ + t12(ip) * ny_ + t13(ip) * nz_
 
      taum_y = t21(im) * nx_ + t22(im) * ny_ + t23(im) * nz_
      taup_y = t21(ip) * nx_ + t22(ip) * ny_ + t23(ip) * nz_
 
      taum_z = t31(im) * nx_ + t32(im) * ny_ + t33(im) * nz_
      taup_z = t31(ip) * nx_ + t32(ip) * ny_ + t33(ip) * nz_
 
      bnd_flux_mom(i,1) = 0.5_rp * ( taup_x - taum_x )
      bnd_flux_mom(i,2) = 0.5_rp * ( taup_y - taum_y )
      bnd_flux_mom(i,3) = 0.5_rp * ( taup_z - taum_z )
      bnd_flux_rhot(i)  = 0.5_rp * ( densp * ( df1(ip) * nx_ + df2(ip) * ny_ + df3(ip) * nz_ ) &
                                   - densm * ( df1(im) * nx_ + df2(im) * ny_ + df3(im) * nz_ ) )
    end do
 
    return
  end subroutine cal_bnd_flux
 
!OCL SERIAL  
  subroutine cal_bnd_flux_qtrc( bnd_flux,             & ! (out)
    dfq1, dfq2, dfq3,                                 & ! (in)
    kh, ddens_, dens_hyd,                             & ! (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) ::  bnd_flux(elem%nfptot*lmesh%ne) 
    real(rp), intent(in)  :: dfq1(elem%np*lmesh%nea)
    real(rp), intent(in)  :: dfq2(elem%np*lmesh%nea)
    real(rp), intent(in)  :: dfq3(elem%np*lmesh%nea)
    real(rp), intent(in)  :: kh   (elem%np*lmesh%nea) ! Eddy diffusivity
    real(rp), intent(in) ::  ddens_(elem%np*lmesh%nea)
    real(rp), intent(in) ::  dens_hyd(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) :: nx_, ny_, nz_
    !------------------------------------------------------------------------
 
    !$omp parallel do private( iM, iP, &
    !$omp densM, densP, nx_, ny_, nz_  )
    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 ( ip > elem%Np * lmesh%Ne .and. abs(nz(i)) > eps ) then ! Tentative implementation for the treatmnet of lower/upper boundary. 
        nx_ = nx(i)
        ny_ = ny(i)
        nz_ = nz(i)
      else
        ! nx_ = ( 1.0_RP + sign(1.0_RP,nx(i)) ) * nx(i)
        ! ny_ = ( 1.0_RP + sign(1.0_RP,ny(i)) ) * ny(i)
        ! nz_ = ( 1.0_RP + sign(1.0_RP,nz(i)) ) * nz(i)
        nx_ = nx(i)
        ny_ = ny(i)
        nz_ = nz(i)
      end if
 
      if ( is_bound(i) ) then
        bnd_flux(i)  = - densm * ( dfq1(im) * nx_ + dfq2(im) * ny_ + dfq3(im) * nz_ )
      else
        bnd_flux(i)  = 0.5_rp * ( densp * ( dfq1(ip) * nx_ + dfq2(ip) * ny_ + dfq3(ip) * nz_ ) &
                                - densm * ( dfq1(im) * nx_ + dfq2(im) * ny_ + dfq3(im) * nz_ ) )
      end if
    end do
 
    return
  end subroutine cal_bnd_flux_qtrc
 
!OCL SERIAL
  elemental function fact(dz, dx, dy)
    implicit none
    real(rp), intent(in) :: dz
    real(rp), intent(in) :: dx
    real(rp), intent(in) :: dy
    real(rp) :: fact ! (out)
 
    real(rp), parameter :: oot = -1.0_rp/3.0_rp
    real(rp), parameter :: fot =  5.0_rp/3.0_rp
    real(rp), parameter :: eot = 11.0_rp/3.0_rp
    real(rp), parameter :: tof = -3.0_rp/4.0_rp
    real(rp) :: a1, a2, b1, b2, dmax
    !--------------------------------------------------------------------
 
    dmax = max(dz, dx, dy)
    if ( dz == dmax ) then
       a1 = dx / dmax
       a2 = dy / dmax
    else if ( dx == dmax ) then
       a1 = dz / dmax
       a2 = dy / dmax
    else ! dy == dmax
       a1 = dz / dmax
       a2 = dx / dmax
    end if
    b1 = atan( a1/a2 )
    b2 = atan( a2/a1 )
 
   fact = 1.736_rp * (a1*a2)**oot &
         * ( 4.0_rp*p1(b1)*a1**oot + 0.222_rp*p2(b1)*a1**fot + 0.077*p3(b1)*a1**eot - 3.0_rp*b1 &
           + 4.0_rp*p1(b2)*a2**oot + 0.222_rp*p2(b2)*a2**fot + 0.077*p3(b2)*a2**eot - 3.0_rp*b2 &
           )**tof
   return
  end function fact
!OCL SERIAL
  elemental function p1(z)
    implicit none
    real(rp), intent(in) :: z
    real(rp) :: p1 ! (out)
 
    real(rp), parameter :: twooverthree  = 2.0_rp / 3.0_rp
    !--------------------------------------------------------------------
 
    p1 = 2.5_rp * p2(z) - 1.5_rp * sin(z) * cos(z)**twooverthree
    return
  end function p1
!OCL SERIAL
  elemental function p2(z)
    implicit none
    real(rp), intent(in) :: z
    real(rp) :: p2 ! (out)
    !--------------------------------------------------------------------
 
    p2 = 0.986_rp * z + 0.073_rp * z**2 - 0.418_rp * z**3 + 0.120_rp * z**4
    return
  end function p2
!OCL SERIAL
  elemental function p3(z)
    implicit none
    real(rp), intent(in) :: z
    real(rp) :: p3 ! (out)
    !--------------------------------------------------------------------
 
    p3 = 0.976_rp * z + 0.188_rp * z**2 - 1.169_rp * z**3 + 0.755_rp * z**4 - 0.151_rp * z**5
    return
  end function p3
 
end module scale_atm_phy_tb_dgm_common