scale_atm_phy_tb_dgm_driver.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Physics turbulence / Atmosphere / driver 
!!
!! @par Description
!!      Driver module for turbulent model based on DGM 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_tb_dgm_driver
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_tracer, only: qa
 
  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, only: sparsemat
 
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmesh_base, only: localmeshbase
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
  use scale_element_base, only: &
    elementbase, elementbase2d, elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
  use scale_meshfield_base, only: &
    meshfieldbase, meshfield2d, meshfield3d
  
  use scale_model_mesh_manager, only: modelmesh3d
  use scale_model_var_manager, only: &
    modelvarmanager, variableinfo
 
  use scale_atm_dyn_dgm_bnd, only: atmdynbnd
 
  use scale_atm_phy_tb_dgm_dns, only: &
    atm_phy_tb_dgm_dns_init,          &
    atm_phy_tb_dgm_dns_final,         &
    atm_phy_tb_dgm_dns_cal_grad
 
  use scale_atm_phy_tb_dgm_smg, only: &
    atm_phy_tb_dgm_smg_init,          &
    atm_phy_tb_dgm_smg_final,         &
    atm_phy_tb_dgm_smg_cal_grad
 
  use scale_atm_phy_tb_dgm_globalsmg, only: &
    atm_phy_tb_dgm_globalsmg_init,          &
    atm_phy_tb_dgm_globalsmg_final,         &
    atm_phy_tb_dgm_globalsmg_cal_grad,      &
    atm_phy_tb_dgm_globalsmg_cal_tend,      &
    atm_phy_tb_dgm_globalsmg_cal_grad_qtrc, &
    atm_phy_tb_dgm_globalsmg_cal_tend_qtrc
 
  use scale_atm_dyn_dgm_nonhydro3d_common, only: &
    prgvar_num, &
    dens_vid => prgvar_ddens_id, therm_vid => prgvar_therm_id,&
    momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id,   &
    momz_vid => prgvar_momz_id,                               &
    auxvar_num, &
    preshyd_vid => auxvar_preshydro_id, denshyd_vid => auxvar_denshydro_id, &
    pres_vid => auxvar_pres_id, pt_vid => auxvar_pt_id,                     &
    cptot_vid => auxvar_cptot_id, rtot_vid => auxvar_rtot_id
 
  use scale_atm_phy_tb_dgm_common, only: &
    atm_phy_tb_dgm_common_cal_grad_qtrc, &
    atm_phy_tb_dgm_common_cal_tend,      &
    atm_phy_tb_dgm_common_cal_tend_qtrc, &
    atmos_phy_tb_tends_num1, &
    tb_momx_t_vid => atmos_phy_tb_momx_t_id,tb_momy_t_vid => atmos_phy_tb_momy_t_id,  &
    tb_momz_t_vid => atmos_phy_tb_momz_t_id, tb_rhot_t_vid => atmos_phy_tb_rhot_t_id, &
    atmos_phy_tb_aux_num,                                                   &
    t11_vid => atmos_phy_tb_aux_t11_id, t12_vid => atmos_phy_tb_aux_t12_id, t13_vid => atmos_phy_tb_aux_t13_id, &
    t21_vid => atmos_phy_tb_aux_t21_id, t22_vid => atmos_phy_tb_aux_t22_id, t23_vid => atmos_phy_tb_aux_t23_id, &
    t31_vid => atmos_phy_tb_aux_t31_id, t32_vid => atmos_phy_tb_aux_t32_id, t33_vid => atmos_phy_tb_aux_t33_id, &
    df1_vid => atmos_phy_tb_aux_diffflx1_id, df2_vid => atmos_phy_tb_aux_diffflx2_id,                           &
    df3_vid => atmos_phy_tb_aux_diffflx3_id,                                                                    &
    nu_vid => atmos_phy_tb_diag_nu_id, kh_vid => atmos_phy_tb_diag_kh_id,                                       &
    tke_vid => atmos_phy_tb_diag_tke_id
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
  abstract interface    
    subroutine atm_phy_tb_final()
      implicit none
    end subroutine atm_phy_tb_final
  end interface 
 
  abstract interface
    subroutine atm_phy_tb_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)
      import rp
      import localmesh3d
      import localmesh2d
      import elementbase3d
      import elementbase2d
      import sparsemat
      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), T12(elem%Np,lmesh%NeA), T13(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: T21(elem%Np,lmesh%NeA), T22(elem%Np,lmesh%NeA), T23(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: T31(elem%Np,lmesh%NeA), T32(elem%Np,lmesh%NeA), T33(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: DF1(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: DF2(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: DF3(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: TKE(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: Nu(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: Kh(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)  :: DENS_hyd(elem%Np,lmesh%NeA)
      real(RP), intent(in)  :: PRES_hyd(elem%Np,lmesh%NeA)
      real(RP), intent(in)  :: PRES(elem%Np,lmesh%NeA)
      real(RP), intent(in)  :: PT(elem%Np,lmesh%NeA)
      type(SparseMat), intent(in) :: Dx, Dy, Dz
      type(SparseMat), intent(in) :: Sx, Sy, Sz
      type(SparseMat), intent(in) :: Lift
      logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne)
    end subroutine atm_phy_tb_cal_grad
  end interface
 
  abstract interface
    subroutine atm_phy_tb_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)
      import rp
      import localmesh3d
      import localmesh2d
      import elementbase3d
      import elementbase2d
      import sparsemat
      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)
      real(RP), intent(out) :: DFQ2(elem%Np,lmesh%NeA)
      real(RP), intent(out) :: DFQ3(elem%Np,lmesh%NeA)
      real(RP), intent(inout) :: dRdx(elem%Np,lmesh%NeA)
      real(RP), intent(inout) :: dRdy(elem%Np,lmesh%NeA)
      real(RP), intent(inout) :: dRdz(elem%Np,lmesh%NeA)
      real(RP), intent(in)  :: Kh(elem%Np,lmesh%NeA)
      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)        
      type(SparseMat), intent(in) :: Dx, Dy, Dz
      type(SparseMat), intent(in) :: Sx, Sy, Sz
      type(SparseMat), intent(in) :: Lift
      logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne)
      logical, intent(in) :: cal_grad_dens
    end subroutine atm_phy_tb_cal_grad_qtrc
  end interface
 
  abstract interface
    subroutine atm_phy_tb_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)
 
    import rp
    import localmesh3d
    import localmesh2d
    import elementbase3d
    import elementbase2d
    import sparsemat
    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)
    real(RP), intent(out) :: MOMY_t(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMZ_t(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: RHOT_t(elem%Np,lmesh%NeA)
    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)  :: PRES_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: PT_  (elem%Np,lmesh%NeA)
    type(SparseMat), intent(in) :: Dx, Dy, Dz
    type(SparseMat), intent(in) :: Sx, Sy, Sz
    type(SparseMat), intent(in) :: Lift
    logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne)
 
    end subroutine atm_phy_tb_cal_tend
  end interface
 
  abstract interface 
    subroutine atm_phy_tb_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)
 
      import rp
      import localmesh3d
      import localmesh2d
      import elementbase3d
      import elementbase2d
      import sparsemat
      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)
      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)
      type(SparseMat), intent(in) :: Dx, Dy, Dz
      type(SparseMat), intent(in) :: Sx, Sy, Sz
      type(SparseMat), intent(in) :: Lift
      logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne)
    end subroutine atm_phy_tb_cal_tend_qtrc
  end interface 
 
  type, public :: atmphytbdgmdriver
    integer :: tb_typeid !< Index of the turbulent model
 
    !- Auxiliary variables with tracer
 
    type(meshfield3d), allocatable :: auxtrcvars(:)
    type(modelvarmanager) :: auxtrcvars_manager
    integer :: auxtrcvars_commid
 
    type(meshfield3d) :: grad_dens(3)
 
    procedure(atm_phy_tb_cal_grad), pointer, nopass :: tbsolver_cal_grad => null()
    procedure(atm_phy_tb_cal_grad_qtrc), pointer, nopass :: tbsolver_cal_grad_qtrc => null()
    procedure(atm_phy_tb_final), pointer, nopass :: tbsolver_final => null()
    procedure(atm_phy_tb_cal_tend), pointer, nopass :: tbsolver_cal_tend => null()
    procedure(atm_phy_tb_cal_tend_qtrc), pointer, nopass :: tbsolver_cal_tend_qtrc => null()
  contains
    procedure :: init => atmphytbdgmdriver_init
    procedure :: final => atmphytbdgmdriver_final
    procedure :: tendency => atmphytbdgmdriver_tendency
  end type atmphytbdgmdriver
 
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  integer, parameter :: tb_typeid_dns                 = 1
  integer, parameter :: tb_typeid_smagorinsky         = 2
  integer, parameter :: tb_typeid_smagorinsky_global  = 3
 
  integer, public, parameter :: atmos_phy_tb_auxtrc_dfq3_id    = 1
  integer, public, parameter :: atmos_phy_tb_auxtrc_dfq1_id    = 2  
  integer, public, parameter :: atmos_phy_tb_auxtrc_dfq2_id    = 3
  integer, public, parameter :: atmos_phy_tb_auxtrc_scalar_num = 1 
  integer, public, parameter :: atmos_phy_tb_auxtrc_hvec_num   = 1
  integer, public, parameter :: atmos_phy_tb_auxtrc_num        = 3
  type(variableinfo), public :: atmos_phy_tb_auxtrc_vinfo(atmos_phy_tb_auxtrc_num)
  DATA atmos_phy_tb_auxtrc_vinfo / &
    variableinfo( atmos_phy_tb_auxtrc_dfq3_id, 'DFQ3', 'Kh * gradient of QTRC (z)',    &
                  'm2/s.kg/kg.m-1',  3, 'XYZ',  ''                                 ),  &
    variableinfo( atmos_phy_tb_auxtrc_dfq1_id, 'DFQ1', 'Kh * gradient of QTRC (x)',    &
                  'm2/s.kg/kg.m-1',  3, 'XYZ',  ''                                 ),  &
    variableinfo( atmos_phy_tb_auxtrc_dfq2_id, 'DFQ2', 'Kh * gradient of QTRC (y)',    &
                  'm2/s.kg/kg.m-1',  3, 'XYZ',  ''                                 )   /
 
contains
!OCL SERIAL  
  subroutine atmphytbdgmdriver_init( this, &
    tb_type_name, dtsec,                   &
    model_mesh3D )
    use scale_mesh_cubedspheredom3d, only: meshcubedspheredom3d
    implicit none
 
    class(atmphytbdgmdriver), intent(inout) :: this
    character(len=*), intent(in) :: tb_type_name
    real(dp), intent(in) :: dtsec
    class(modelmesh3d), intent(inout), target :: model_mesh3d
 
    class(meshbase3d), pointer :: mesh3d
    class(meshcubedspheredom3d), pointer :: gm_mesh3d
 
    integer :: iv
    logical :: reg_file_hist
                    
    integer :: idim
    !-----------------------------------------------------------------------------
 
    mesh3d => model_mesh3d%ptr_mesh
 
    !--- Set the type of turbulence scheme
    
    select case(tb_type_name)
    case ('DNS')
      this%TB_TYPEID = tb_typeid_dns
 
      call atm_phy_tb_dgm_dns_init( mesh3d )
      this%tbsolver_cal_grad => atm_phy_tb_dgm_dns_cal_grad
      this%tbsolver_cal_grad_qtrc => atm_phy_tb_dgm_common_cal_grad_qtrc
      this%tbsolver_cal_tend => atm_phy_tb_dgm_common_cal_tend
      this%tbsolver_cal_tend_qtrc => atm_phy_tb_dgm_common_cal_tend_qtrc
      this%tbsolver_final => atm_phy_tb_dgm_dns_final
    case ('SMAGORINSKY')
      this%TB_TYPEID = tb_typeid_smagorinsky
 
      call atm_phy_tb_dgm_smg_init( mesh3d )
      this%tbsolver_cal_grad => atm_phy_tb_dgm_smg_cal_grad
      this%tbsolver_cal_grad_qtrc => atm_phy_tb_dgm_common_cal_grad_qtrc
      this%tbsolver_cal_tend => atm_phy_tb_dgm_common_cal_tend
      this%tbsolver_cal_tend_qtrc => atm_phy_tb_dgm_common_cal_tend_qtrc
      this%tbsolver_final => atm_phy_tb_dgm_smg_final
    case ('SMAGORINSKY_GLOBAL')
      this%TB_TYPEID = tb_typeid_smagorinsky_global
 
      select type(mesh3d)
      class is (meshcubedspheredom3d)
        gm_mesh3d => mesh3d
      end select
      call atm_phy_tb_dgm_globalsmg_init( mesh3d, gm_mesh3d%shallow_approx )
      this%tbsolver_cal_grad => atm_phy_tb_dgm_globalsmg_cal_grad
      this%tbsolver_cal_grad_qtrc => atm_phy_tb_dgm_globalsmg_cal_grad_qtrc
      this%tbsolver_cal_tend => atm_phy_tb_dgm_globalsmg_cal_tend
      this%tbsolver_cal_tend_qtrc => atm_phy_tb_dgm_globalsmg_cal_tend_qtrc
      this%tbsolver_final => atm_phy_tb_dgm_globalsmg_final
    case default
      log_error("AtmPhyTbDGMDriver_Init",*) 'Invalid TB_TYPE in namelist PARAM_ATMOS_TB. Check!'
      call prc_abort
    end select
    
    !- Initialize variables
 
    !-
    if ( qa > 0 ) then
      call this%auxtrcvars_manager%Init()
      allocate( this%auxtrcvars(atmos_phy_tb_auxtrc_num) )
 
      reg_file_hist = .true.    
      do iv = 1, atmos_phy_tb_auxtrc_num
        call this%auxtrcvars_manager%Regist( &
          atmos_phy_tb_auxtrc_vinfo(iv), mesh3d,  & ! (in) 
          this%auxtrcvars(iv),                    & ! (inout)
          reg_file_hist, fill_zero=.true.         ) ! (in)
      end do
 
      !- Setup communication
 
      call model_mesh3d%Create_communicator( &
        atmos_phy_tb_auxtrc_scalar_num, atmos_phy_tb_auxtrc_hvec_num, 0, & ! (in) 
        this%auxtrcvars_manager,         & ! (inout)
        this%auxtrcvars(:),              & ! (in)
        this%auxtrcvars_commid           ) ! (out)
 
      do idim=1, 3
        call this%GRAD_DENS(idim)%Init( "Grad_DENS", "kg/m4", mesh3d )      
      end do
    end if
 
    return
  end subroutine atmphytbdgmdriver_init
 
!OCL SERIAL
  subroutine atmphytbdgmdriver_tendency( this, &
    TB_TENDS, PROG_VARS, TRC_VARS, AUX_VARS,   &
    AUX_TB_VARS,  DIAG_TB_VARS,                &
    boundary_cond,                             &
    Dx, Dy, Dz, Sx, Sy, Sz, Lift, mesh3D,      &
    gfilter_flag, gfilter, diag_nu_kh,         &
    modal_filter_flag, modal_filter )
 
    use scale_tracer, only: &
      qa, tracer_advc, tracer_name
      
    use scale_model_var_manager, only: modelvarmanager
 
    use scale_meshfield_filter_operation_3d, only: &
      meshfieldfilteroperation3d
    use scale_element_modalfilter, only: modalfilter
    implicit none
 
    class(atmphytbdgmdriver), intent(inout) :: this
    class(meshbase3d), intent(in), target :: mesh3d
    class(modelvarmanager), intent(inout) :: tb_tends
    class(modelvarmanager), intent(inout) :: prog_vars
    class(modelvarmanager), intent(inout) :: trc_vars
    class(modelvarmanager), intent(inout) :: aux_vars
    class(modelvarmanager), intent(inout) :: aux_tb_vars
    class(modelvarmanager), intent(inout) :: diag_tb_vars
    type(atmdynbnd), intent(in) :: boundary_cond
    type(sparsemat), intent(in) :: dx, dy, dz
    type(sparsemat), intent(in) :: sx, sy, sz
    type(sparsemat), intent(in) :: lift
    logical, intent(in) :: gfilter_flag
    class(meshfieldfilteroperation3d), intent(inout) :: gfilter
    type(meshfield3d), intent(inout) :: diag_nu_kh(2)
    logical, intent(in) :: modal_filter_flag
    type(modalfilter), intent(in) :: modal_filter
 
    class(localmesh3d), pointer :: lcmesh3d
    integer :: n
    integer :: ke
 
    class(meshfield3d), pointer :: ddens, momx, momy, momz, therm
    class(meshfield3d), pointer :: pres_hyd, dens_hyd, pres, pt, rtot, cptot
    class(meshfield3d), pointer :: t11, t12, t13, t21, t22, t23, t31, t32, t33
    class(meshfield3d), pointer :: df1, df2, df3
    class(meshfield3d), pointer :: dfq1, dfq2, dfq3
    class(meshfield3d), pointer :: nu, kh, tke
    class(meshfield3d), pointer :: tb_momx_t, tb_momy_t, tb_momz_t, tb_rhot_t
    class(meshfield3d), pointer :: qtrc
    class(meshfield3d), pointer :: tb_rhoq_t
 
    type dyn_bndinfo
      logical, allocatable :: is_bound(:,:)
    end type
    type(dyn_bndinfo), allocatable :: bnd_info(:)
 
    integer :: iq
    logical :: cal_grad_flag
    !-----------------------------------------------------------------------------
    
    !-
    call prof_rapstart( 'ATM_TB_tendency_pre', 2)
    call prog_vars%Get3D(dens_vid , ddens)
    call prog_vars%Get3D(therm_vid, therm)
    call prog_vars%Get3D(momz_vid , momz )
    call prog_vars%Get3D(momx_vid , momx )
    call prog_vars%Get3D(momy_vid , momy )
 
    call aux_vars%Get3D( preshyd_vid, pres_hyd )
    call aux_vars%Get3D( denshyd_vid, dens_hyd )
    call aux_vars%Get3D( pres_vid, pres )
    call aux_vars%Get3D( pt_vid, pt )
    call aux_vars%Get3D( rtot_vid, rtot )
    call aux_vars%Get3D( cptot_vid, cptot )
 
    call aux_tb_vars%Get3D( t11_vid, t11 )
    call aux_tb_vars%Get3D( t12_vid, t12 )
    call aux_tb_vars%Get3D( t13_vid, t13 )
    call aux_tb_vars%Get3D( t21_vid, t21 )
    call aux_tb_vars%Get3D( t22_vid, t22 )
    call aux_tb_vars%Get3D( t23_vid, t23 )
    call aux_tb_vars%Get3D( t31_vid, t31 )
    call aux_tb_vars%Get3D( t32_vid, t32 )
    call aux_tb_vars%Get3D( t33_vid, t33 )
    call aux_tb_vars%Get3D( df1_vid, df1 )
    call aux_tb_vars%Get3D( df2_vid, df2 )
    call aux_tb_vars%Get3D( df3_vid, df3 )
 
    call diag_tb_vars%Get3D( tke_vid, tke )
    call diag_tb_vars%Get3D( nu_vid, nu )
    call diag_tb_vars%Get3D( kh_vid, kh )
 
    call tb_tends%Get3D( tb_momx_t_vid, tb_momx_t )
    call tb_tends%Get3D( tb_momy_t_vid, tb_momy_t )
    call tb_tends%Get3D( tb_momz_t_vid, tb_momz_t )
    call tb_tends%Get3D( tb_rhot_t_vid, tb_rhot_t )
 
    if ( qa > 0 ) then
      call this%auxtrcvars_manager%Get3D( atmos_phy_tb_auxtrc_dfq1_id, dfq1 )
      call this%auxtrcvars_manager%Get3D( atmos_phy_tb_auxtrc_dfq2_id, dfq2 )
      call this%auxtrcvars_manager%Get3D( atmos_phy_tb_auxtrc_dfq3_id, dfq3 )
    end if
 
    allocate( bnd_info(mesh3d%LOCAL_MESH_NUM) )
    call prof_rapend( 'ATM_TB_tendency_pre', 2)
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh3d => mesh3d%lcmesh_list(n)
 
      !- Apply boundary conditions
      call prof_rapstart('ATM_PHY_TB_bnd', 2)
      allocate( bnd_info(n)%is_bound(lcmesh3d%refElem3D%NfpTot,lcmesh3d%Ne) )
      call boundary_cond%Inquire_bound_flag(  bnd_info(n)%is_bound, & ! (out)
        n, lcmesh3d%VMapM, lcmesh3d%VMapP, lcmesh3d%VMapB,          & ! (in)
        lcmesh3d, lcmesh3d%refElem3D                                ) ! (in)
      
      call boundary_cond%ApplyBC_Grad_TBVARS_lc( n, &
        ddens%local(n)%val, momx%local(n)%val, momy%local(n)%val, momz%local(n)%val, pt%local(n)%val, pres%local(n)%val, & ! (inout)
        dens_hyd%local(n)%val, pres_hyd%local(n)%val, rtot%local(n)%val, cptot%local(n)%val,                             & ! (in)
        lcmesh3d%Gsqrt(:,:), lcmesh3d%GsqrtH(:,:), lcmesh3d%GIJ(:,:,1,1), lcmesh3d%GIJ(:,:,1,2), lcmesh3d%GIJ(:,:,2,2),  & ! (in)
        lcmesh3d%GI3(:,:,1), lcmesh3d%GI3(:,:,2),                                                                        & ! (in)
        lcmesh3d%normal_fn(:,:,1), lcmesh3d%normal_fn(:,:,2), lcmesh3d%normal_fn(:,:,3),                                 & ! (in)
        lcmesh3d%vmapM, lcmesh3d%vmapP, lcmesh3d%vmapB,                                                                  & ! (in)
        lcmesh3d, lcmesh3d%refElem3D, lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D                                     ) ! (in)        
      call prof_rapend('ATM_PHY_TB_bnd', 2)
 
      call prof_rapstart('ATM_PHY_TB_cal_grad', 2)
      call this%tbsolver_cal_grad( &
          t11%local(n)%val, t12%local(n)%val, t13%local(n)%val,         & ! (out)
          t21%local(n)%val, t22%local(n)%val, t23%local(n)%val,         & ! (out)
          t31%local(n)%val, t32%local(n)%val, t33%local(n)%val,         & ! (out)
          df1%local(n)%val, df2%local(n)%val, df3%local(n)%val,         & ! (out)
          tke%local(n)%val, nu%local(n)%val, kh%local(n)%val,           & ! (out)
          ddens%local(n)%val, momx%local(n)%val, momy%local(n)%val, momz%local(n)%val, therm%local(n)%val,            & ! (in)
          dens_hyd%local(n)%val, pres_hyd%local(n)%val, pres%local(n)%val, pt%local(n)%val,                           & ! (in)
          dx, dy, dz, sx, sy, sz, lift, lcmesh3d, lcmesh3d%refElem3D, lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D, & ! (in)
          bnd_info(n)%is_bound                                                                                        ) ! (in)
      call prof_rapend('ATM_PHY_TB_cal_grad', 2)
    end do
 
    !* Exchange halo data
    call prof_rapstart('ATM_PHY_TB_exchange_prgv', 2)
    call aux_tb_vars%MeshFieldComm_Exchange()
    call prof_rapend('ATM_PHY_TB_exchange_prgv', 2)
 
    if ( gfilter_flag ) then
      call prof_rapstart('ATM_PHY_TB_filter_NU_KH', 2)
      call gfilter%Apply( diag_nu_kh, mesh3d )
      call prof_rapend('ATM_PHY_TB_filter_NU_KH', 2)
    end if
    if ( modal_filter_flag ) then
      call prof_rapstart('ATM_PHY_TB_filter_NU_KH', 2)
      do n=1, mesh3d%LOCAL_MESH_NUM
        lcmesh3d => mesh3d%lcmesh_list(n)
        !$omp parallel do
        do ke=lcmesh3d%NeS, lcmesh3d%NeE
          nu%local(n)%val(:,ke) = matmul( modal_filter%FilterMat, nu%local(n)%val(:,ke) )
          kh%local(n)%val(:,ke) = matmul( modal_filter%FilterMat, kh%local(n)%val(:,ke) )
        end do
      end do
      call prof_rapend('ATM_PHY_TB_filter_NU_KH', 2)
    end if
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh3d => mesh3d%lcmesh_list(n)
 
      !- Apply boundary conditions for stress tensor
      call prof_rapstart('ATM_PHY_TB_bnd', 2)
      call boundary_cond%ApplyBC_Grad_TBStress_lc( n, &
        t11%local(n)%val, t12%local(n)%val, t13%local(n)%val,                                                            & ! (inout)
        t21%local(n)%val, t22%local(n)%val, t23%local(n)%val,                                                            & ! (inout)
        t31%local(n)%val, t32%local(n)%val, t33%local(n)%val,                                                            & ! (inout)
        df1%local(n)%val, df2%local(n)%val, df3%local(n)%val,                                                            & ! (inout)
        lcmesh3d%Gsqrt(:,:), lcmesh3d%GsqrtH(:,:), lcmesh3d%GIJ(:,:,1,1), lcmesh3d%GIJ(:,:,1,2), lcmesh3d%GIJ(:,:,2,2),  & ! (in)
        lcmesh3d%GI3(:,:,1), lcmesh3d%GI3(:,:,2),                                                                        & ! (in)
        lcmesh3d%normal_fn(:,:,1), lcmesh3d%normal_fn(:,:,2), lcmesh3d%normal_fn(:,:,3),                                 & ! (in)
        lcmesh3d%vmapM, lcmesh3d%vmapP, lcmesh3d%vmapB,                                                                  & ! (in)
        lcmesh3d, lcmesh3d%refElem3D, lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D                                     ) ! (in)        
      call prof_rapend('ATM_PHY_TB_bnd', 2)
 
      call prof_rapstart('ATM_PHY_TB_cal_tend', 2)
      call this%tbsolver_cal_tend( &
        tb_momx_t%local(n)%val, tb_momy_t%local(n)%val, tb_momz_t%local(n)%val, tb_rhot_t%local(n)%val,             & ! (out)
        t11%local(n)%val, t12%local(n)%val, t13%local(n)%val,                                                       & ! (in)
        t21%local(n)%val, t22%local(n)%val, t23%local(n)%val,                                                       & ! (in)
        t31%local(n)%val, t32%local(n)%val, t33%local(n)%val,                                                       & ! (in)
        df1%local(n)%val, df2%local(n)%val, df3%local(n)%val, nu%local(n)%val, kh%local(n)%val,                     & ! (in)
        ddens%local(n)%val, momx%local(n)%val, momy%local(n)%val, momz%local(n)%val, therm%local(n)%val,            & ! (in)
        dens_hyd%local(n)%val, pres_hyd%local(n)%val, pres%local(n)%val, pt%local(n)%val,                           & ! (in)
        dx, dy, dz, sx, sy, sz, lift, lcmesh3d, lcmesh3d%refElem3D, lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D, & ! (in)
        bnd_info(n)%is_bound                                                                                        ) ! (in)
      call prof_rapend('ATM_PHY_TB_cal_tend', 2)
    end do
 
    cal_grad_flag = .true.
    do iq = 1, qa
      if ( .not. tracer_advc(iq) ) cycle
      
      call trc_vars%Get3D( iq, qtrc )
      call tb_tends%Get3D( atmos_phy_tb_tends_num1 + iq, tb_rhoq_t )
      
      do n=1, mesh3d%LOCAL_MESH_NUM
        lcmesh3d => mesh3d%lcmesh_list(n)
 
        call prof_rapstart('ATM_PHY_TB_cal_grad_qtrc', 2)
        call this%tbsolver_cal_grad_qtrc( &
          dfq1%local(n)%val, dfq2%local(n)%val, dfq3%local(n)%val,                       & ! (out)
          this%GRAD_DENS(1)%local(n)%val, this%GRAD_DENS(2)%local(n)%val,                & ! (inout)
          this%GRAD_DENS(3)%local(n)%val,                                                & ! (inout)
          kh%local(n)%val, qtrc%local(n)%val, ddens%local(n)%val, dens_hyd%local(n)%val, & ! (in) 
          dx, dy, dz, sx, sy, sz, lift, lcmesh3d, lcmesh3d%refElem3D,                    & ! (in)
          lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D,                                & ! (in)
          bnd_info(n)%is_bound, cal_grad_flag )                                            ! (in)
        call prof_rapend('ATM_PHY_TB_cal_grad_qtrc', 2)
      end do
      cal_grad_flag = .false.
 
      call this%auxtrcvars_manager%MeshFieldComm_Exchange()
 
      do n=1, mesh3d%LOCAL_MESH_NUM
        lcmesh3d => mesh3d%lcmesh_list(n)
 
        call prof_rapstart('ATM_PHY_TB_cal_tend_qtrc', 2)
        call this%tbsolver_cal_tend_qtrc( tb_rhoq_t%local(n)%val,          & ! (out)
          dfq1%local(n)%val, dfq2%local(n)%val, dfq3%local(n)%val,         & ! (out)
          kh%local(n)%val, ddens%local(n)%val, dens_hyd%local(n)%val,      & ! (in) 
          dx, dy, dz, sx, sy, sz, lift, lcmesh3d, lcmesh3d%refElem3D,      & ! (in)
          lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D,                  & ! (in)
          bnd_info(n)%is_bound )                                             ! (in)
        call prof_rapend('ATM_PHY_TB_cal_tend_qtrc', 2)
      end do        
    end do
 
    do n = 1, mesh3d%LOCAL_MESH_NUM
      deallocate( bnd_info(n)%is_bound )
    end do
  
    return
  end subroutine atmphytbdgmdriver_tendency
 
!OCL SERIAL
  subroutine atmphytbdgmdriver_final( this )
    implicit none
 
    class(atmphytbdgmdriver), intent(inout) :: this
 
    integer :: idim
    !-----------------------------------------------------------------------------
 
    call this%tbsolver_final()
 
    if ( qa > 0 ) then
      do idim=1, 3
        call this%GRAD_DENS(idim)%Final()
      end do
  
      call this%auxtrcvars_manager%Final()
    end if
 
    return
  end subroutine atmphytbdgmdriver_final
 
end module scale_atm_phy_tb_dgm_driver