scale_atm_dyn_dgm_driver_trcadv3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / DGM driver (tracer advection)
!!
!! @par Description
!!      Driver module for tracer advection based on DGM 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_driver_trcadv3d
  !-----------------------------------------------------------------------------
  !
  !++ 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_tracer, only: qa
 
  use scale_timeint_rk, only: timeint_rk
  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_element_operation_base, only: elementoperationbase3d
 
  use scale_meshfield_base, only: &
    meshfieldbase, meshfield2d, meshfield3d
  use scale_localmeshfield_base, only: &
      localmeshfieldbaselist
 
  use scale_element_modalfilter, only: &
    modalfilter
 
  use scale_model_var_manager, only: &
    modelvarmanager, variableinfo
  use scale_model_mesh_manager, only: modelmesh3d
 
  use scale_atm_dyn_dgm_driver_base, only: &
    atmdyndgmdriver_base3d, atmdyndgmdriver_base3d_init, atmdyndgmdriver_base3d_final
 
  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,                 &
    cptot_vid => auxvar_cptot_id, cvtot_vid => auxvar_cvtot_id, rtot_vid => auxvar_rtot_id, &
    pres_vid => auxvar_pres_id, &
    phytend_num1 => phytend_num
  
  use scale_atm_dyn_dgm_bnd, only: atmdynbnd
 
  use scale_atm_dyn_dgm_trcadvect3d_heve, only: &
    atm_dyn_dgm_trcadvect3d_heve_init,                 &
    atm_dyn_dgm_trcadvect3d_heve_final,                &    
    atm_dyn_dgm_trcadvect3d_heve_calc_fct_coef,        &
    atm_dyn_dgm_trcadvect3d_heve_cal_tend,             &
    atm_dyn_dgm_trcadvect3d_tmar,                      &
    atm_dyn_dgm_trcadvect3d_save_massflux,             &
    atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_advtest
 
  use scale_atm_dyn_dgm_driver_nonhydro3d, only: &
    atmdyndgmdriver_nonhydro3d
    
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  type, extends(atmdyndgmdriver_base3d), public :: atmdyndgmdriver_trcadv3d
    integer :: eqs_typeid
    logical :: only_traceradv_flag
    type(sparsemat) :: faceintmat    
 
    ! Limiter
    logical :: disable_limiter
    type(meshfield3d) :: alphadensm, alphadensp ! coeffcient with stabilizion terms in numerical flux (element boundary data)
 
    ! 
    type(meshfield3d), allocatable :: trcvars3d(:)
    type(modelvarmanager) :: trcvar3d_manager
    integer :: trcvar3d_commid
 
    type(meshfield3d), allocatable :: aux_trcvars3d(:)
    type(modelvarmanager) :: auxtrcvar3d_manager
    integer :: auxtrcvar3d_commid
 
    type(meshfield3d), allocatable :: auxtrc_flux_vars3d(:)
    type(modelvarmanager) :: auxtrc_flux_var3d_manager
    integer :: auxtrc_flux_var3d_commid    
 
    ! element-wise modal filter
    logical :: modalfilter_flag
    type(modalfilter) :: modal_filter_3d
 
    ! boundary_condition
    type(atmdynbnd), pointer :: boundary_cond
 
  contains
    procedure :: init => atmdyndgmdriver_trcadv3d_init
    procedure :: final => atmdyndgmdriver_trcadv3d_final
    procedure :: update => atmdyndgmdriver_trcadv3d_update
  end type atmdyndgmdriver_trcadv3d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public, parameter :: eqs_typeid_trcadv3d_heve             = 1
  
  !-
  integer, public, parameter :: auxtrcvars3d_num         = 1
  integer, public, parameter :: auxtrcvars3d_fctcoef_id  = 1
 
  type(variableinfo), public :: atmos_dyn_auxtrcvars3d_vinfo(auxtrcvars3d_num)
  DATA atmos_dyn_auxtrcvars3d_vinfo / &
    variableinfo( auxtrcvars3d_fctcoef_id, 'TRCADV_FCTCOEF', '',  &
                  '1',  3, 'XYZ',  ''                          )  / 
 
  !-
  integer, public, parameter :: trcvars3d_num         = 3
  integer, public, parameter :: trcvars3d_trcadv_id   = 1
  integer, public, parameter :: trcvars3d_dens_id     = 2
  integer, public, parameter :: trcvars3d_dens0_id    = 3
 
  type(variableinfo), public :: atmos_dyn_trcvars3d_vinfo(trcvars3d_num)
  DATA atmos_dyn_trcvars3d_vinfo / &
    variableinfo( trcvars3d_trcadv_id, 'TRCADV', '',     '1',  3, 'XYZ',  '' ), &
    variableinfo( trcvars3d_dens_id  ,   'DENS', '', 'kg/m3',  3, 'XYZ',  '' ), & 
    variableinfo( trcvars3d_dens0_id  ,  'DENS0', '', 'kg/m3',  3, 'XYZ',  '' ) / 
 
  !-
  integer, public, parameter :: mass_flux_num   = 3
  integer, public, parameter :: massflx_z_id    = 1  
  integer, public, parameter :: massflx_x_id    = 2
  integer, public, parameter :: massflx_y_id    = 3
 
  type(variableinfo), public :: atmos_dyn_mass_flux_vinfo(mass_flux_num)
  DATA atmos_dyn_mass_flux_vinfo / &
    variableinfo( massflx_z_id, 'MASSFLX_Z', 'flux in z-direction',  &
                  'kg/s/m2',  3, 'XYZ',  ''                     ),   &
    variableinfo( massflx_x_id, 'MASSFLX_X', 'flux in x-direction',  &
                  'kg/s/m2',  3, 'XYZ',  ''                     ),   & 
    variableinfo( massflx_y_id, 'MASSFLX_Y', 'flux in y-direction',  &
                  'kg/s/m2',  3, 'XYZ',  ''                     )    / 
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
  private :: setup_modalfilter
 
contains
!OCL SERIAL  
  subroutine atmdyndgmdriver_trcadv3d_init( this, &
    eqs_type_name, tint_type_name, dtsec,               &
    modal_filter_flag, disable_limiter,                 &
    model_mesh3D, boundary_cond, ONLY_TRACERADV_FLAG    )
 
    use scale_localmeshfield_base, only: local_meshfield_type_nodes_faceval
    implicit none
 
    class(atmdyndgmdriver_trcadv3d), intent(inout) :: this
    character(len=*), intent(in) :: eqs_type_name
    character(len=*), intent(in) :: tint_type_name
    real(dp), intent(in) :: dtsec
    logical, intent(in) :: modal_filter_flag
    logical, intent(in) :: disable_limiter
    class(modelmesh3d), intent(inout), target :: model_mesh3d
    class(atmdynbnd), intent(in), target :: boundary_cond
    logical, intent(in) :: only_traceradv_flag
 
    class(meshbase3d), pointer :: mesh3d
    class(hexahedralelement), pointer :: refelem3d
    class(elementbase), pointer :: refelem
 
    integer :: iv
    integer :: ldomid
    !-----------------------------------------------------------------------------
 
    mesh3d => model_mesh3d%ptr_mesh
 
    call atmdyndgmdriver_base3d_init( this, &
      1,                                    &
      tint_type_name, dtsec,                &
      mesh3d                                )
 
    select case(eqs_type_name)
    !-- HEVE ------------------
    case("TRCADV3D_HEVE")
    case default
        log_error("AtmDynDGMDriver_trcadv3d_Init",*) 'Invalid EQS_TYPE in namelist PARAM_ATMOS_DYN. Check!'
        call prc_abort  
    end select
    
    this%ONLY_TRACERADV_FLAG = only_traceradv_flag
    this%disable_limiter = disable_limiter
 
    !- Initialize a module for tracer equations
    call atm_dyn_dgm_trcadvect3d_heve_init( mesh3d, this%FaceIntMat ) 
 
 
    !- Initialize 3D auxiliary variables for tracer advection
 
    call this%TRCVAR3D_manager%Init()
    allocate( this%TRCVARS3D(trcvars3d_num) )
 
    do iv = 1, trcvars3d_num
      call this%TRCVAR3D_manager%Regist( &
        atmos_dyn_trcvars3d_vinfo(iv), mesh3d,     & ! (in) 
        this%TRCVARS3D(iv),                        & ! (inout)
        .false., fill_zero=.true.                  ) ! (in)
    end do
 
    iv = trcvars3d_trcadv_id
    call model_mesh3d%Create_communicator( &
      1, 0, 0,                         & ! (in) 
      this%TRCVAR3D_manager,           & ! (in)
      this%TRCVARS3D(iv:iv),           & ! (in)
      this%TRCVAR3D_commid             ) ! (out)
      
    !- Initialize variables to store time-averaged 3D mass flux
 
    call this%AUXTRC_FLUX_VAR3D_manager%Init()
    allocate( this%AUXTRC_FLUX_VARS3D(mass_flux_num) )
    
    do iv = 1, mass_flux_num
      call this%AUXTRC_FLUX_VAR3D_manager%Regist(  &
        atmos_dyn_mass_flux_vinfo(iv), mesh3d,     & ! (in) 
        this%AUXTRC_FLUX_VARS3D(iv),               & ! (inout)
        .false., fill_zero=.true.                  ) ! (in)
    end do
 
    call model_mesh3d%Create_communicator( &
      1, 1, 0,                           & ! (in) 
      this%AUXTRC_FLUX_VAR3D_manager,    & ! (in)
      this%AUXTRC_FLUX_VARS3D(:),        & ! (in)
      this%AUXTRC_FLUX_VAR3D_commid      ) ! (out)
 
    !- Initialize 3D auxiliary variables for preserving nonnegativity in tracer advection
 
    call this%AUXTRCVAR3D_manager%Init()
    allocate( this%AUX_TRCVARS3D(auxtrcvars3d_num) )
 
    do iv = 1, auxtrcvars3d_num
      call this%AUXTRCVAR3D_manager%Regist( &
        atmos_dyn_auxtrcvars3d_vinfo(iv), mesh3d,  & ! (in) 
        this%AUX_TRCVARS3D(iv),                    & ! (inout)
        .false., fill_zero=.true.                  ) ! (in)
    end do
 
    call model_mesh3d%Create_communicator( &
      1, 0, 0,                         & ! (in) 
      this%AUXTRCVAR3D_manager,        & ! (in)
      this%AUX_TRCVARS3D(:),           & ! (in)
      this%AUXTRCVAR3D_commid          ) ! (out)
 
    !-
      
    call this%alphaDensM%Init( "alphaDensM", "kg/m3.m/s", mesh3d, local_meshfield_type_nodes_faceval )
    call this%alphaDensP%Init( "alphaDensP", "kg/m3.m/s", mesh3d, local_meshfield_type_nodes_faceval )
    !$omp parallel do
    do ldomid=1, mesh3d%LOCAL_MESH_NUM
      this%alphaDensM%local(ldomid)%face_val(:,:) = 0.0_rp
      this%alphaDensP%local(ldomid)%face_val(:,:) = 0.0_rp
    end do
  
    !- initialize an object to manage boundary conditions
    this%boundary_cond => boundary_cond
 
 
    !- initialize an object to manage modal filter
    this%MODALFILTER_FLAG = modal_filter_flag
    if (this%MODALFILTER_FLAG) then
      refelem => mesh3d%refElem
      select type(refelem)
      class is (hexahedralelement)
        refelem3d => refelem
      end select
 
      call setup_modalfilter( this, refelem3d, model_mesh3d%element3D_operation )
    end if
 
    return
  end subroutine atmdyndgmdriver_trcadv3d_init
 
!OCL SERIAL
  subroutine atmdyndgmdriver_trcadv3d_update( this, &
    TRC_VARS, PROG_VARS, AUX_VARS, PHYTENDS,             &
    element_operation,                                   &
    Dx, Dy, Dz, Sx, Sy, Sz, Lift, mesh3D,                &
    dyn_driver                                           )
 
    use scale_tracer, only: &
      tracer_advc, tracer_name
    use scale_atm_dyn_dgm_modalfilter, only: &
      atm_dyn_dgm_tracer_modalfilter_apply
    implicit none
 
    class(atmdyndgmdriver_trcadv3d), intent(inout) :: this
    class(meshbase3d), intent(in), target :: mesh3d
    class(modelvarmanager), intent(inout) :: trc_vars
    class(modelvarmanager), intent(inout) :: prog_vars
    class(modelvarmanager), intent(inout) :: aux_vars
    class(modelvarmanager), intent(inout) :: phytends
    class(elementoperationbase3d), intent(in) :: element_operation
    type(sparsemat), intent(in) :: dx, dy, dz
    type(sparsemat), intent(in) :: sx, sy, sz
    type(sparsemat), intent(in) :: lift
    class(atmdyndgmdriver_nonhydro3d), intent(inout) :: dyn_driver
 
    integer :: rkstage
    integer :: tintbuf_ind
    real(rp) :: dt
    real(rp) :: dttmp_trc
 
    class(localmesh3d), pointer :: lcmesh3d
    integer :: n
    integer :: ke
 
    class(meshfield3d), pointer :: ddens, momx, momy, momz, therm
    class(meshfield3d), pointer :: dens_hyd, pres_hyd, rtot, cvtot, cptot, pres
    class(meshfield3d), pointer :: qtrc, rhoq_tp
 
    class(meshfield3d), pointer :: qtrc_tmp, ddens_trc, ddens0_trc
    class(meshfield3d), pointer :: mflx_x_tavg, mflx_y_tavg, mflx_z_tavg
 
    integer :: iq    
    !-----------------------------------------------------------------------------
    
    !-- prepairation
 
    call prof_rapstart( 'ATM_DYN_trc_update_pre', 2)
    call prog_vars%Get3D(dens_vid , ddens)
    call prog_vars%Get3D(therm_vid, therm)
 
    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( rtot_vid, rtot )
    call aux_vars%Get3D( cvtot_vid, cvtot )
    call aux_vars%Get3D( cptot_vid, cptot )
 
    call this%TRCVAR3D_manager%Get3D( trcvars3d_trcadv_id, qtrc_tmp )
    call this%TRCVAR3D_manager%Get3D(  trcvars3d_dens_id,  ddens_trc )
    call this%TRCVAR3D_manager%Get3D( trcvars3d_dens0_id, ddens0_trc )
 
    call this%AUXTRC_FLUX_VAR3D_manager%Get3D( massflx_z_id, mflx_z_tavg )
    call this%AUXTRC_FLUX_VAR3D_manager%Get3D( massflx_x_id, mflx_x_tavg )
    call this%AUXTRC_FLUX_VAR3D_manager%Get3D( massflx_y_id, mflx_y_tavg )
 
    if ( this%ONLY_TRACERADV_FLAG ) then
      call prog_vars%Get3D(momz_vid , momz)
      call prog_vars%Get3D(momx_vid , momx)
      call prog_vars%Get3D(momy_vid , momy)
 
      do n=1, mesh3d%LOCAL_MESH_NUM
        lcmesh3d => mesh3d%lcmesh_list(n)
 
        !$omp parallel do
        do ke=lcmesh3d%NeS, lcmesh3d%NeE
          ddens0_trc %local(n)%val(:,ke) = ddens%local(n)%val(:,ke) 
          ddens_trc  %local(n)%val(:,ke) = ddens%local(n)%val(:,ke)
          mflx_z_tavg%local(n)%val(:,ke) = momz%local(n)%val(:,ke)
          mflx_x_tavg%local(n)%val(:,ke) = momx%local(n)%val(:,ke)
          mflx_y_tavg%local(n)%val(:,ke) = momy%local(n)%val(:,ke)
        end do
        call atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_advtest( &
          this%alphaDensM%local(n)%face_val, this%alphaDensP%local(n)%face_val,            & ! (inout)
          ddens%local(n)%val, momx%local(n)%val, momy%local(n)%val, momz%local(n)%val,     & ! (in)
          dens_hyd%local(n)%val,                                                           & ! (in)
          lcmesh3d%Gsqrt,                                                                  & ! (in)
          lcmesh3d%normal_fn(:,:,1), lcmesh3d%normal_fn(:,:,2), lcmesh3d%normal_fn(:,:,3), & ! (in)
          lcmesh3d%VMapM, lcmesh3d%VMapP, lcmesh3d, lcmesh3d%refElem3D                     ) ! (in)
      end do
    end if
 
    !* Exchange halo data of mass flux
 
    call prof_rapstart( 'ATM_DYN_exchange_mflx', 3)
    call this%AUXTRC_FLUX_VAR3D_manager%MeshFieldComm_Exchange()
    call prof_rapend( 'ATM_DYN_exchange_mflx', 3)
 
    call prof_rapstart( 'ATM_DYN_applyBC_mflux', 3)
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh3d => mesh3d%lcmesh_list(n)
 
      call this%boundary_cond%ApplyBC_PROGVARS_lc( n, & ! (in)
        ddens_trc%local(n)%val(:,:), mflx_x_tavg%local(n)%val, mflx_y_tavg%local(n)%val, mflx_z_tavg%local(n)%val, therm%local(n)%val, & ! (inout)
        dens_hyd%local(n)%val, pres_hyd%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, lcmesh3d, lcmesh3d%refElem3D, lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D   ) ! (in)
    end do
    call prof_rapend( 'ATM_DYN_applyBC_mflux', 3)
 
    call prof_rapend( 'ATM_DYN_trc_update_pre', 2)
 
    !-- Tracer advection ------------------------------------------------
 
    call prof_rapstart( 'ATM_DYN_trc_update', 2)
 
    do iq = 1, qa
    
      call trc_vars%Get3D( iq, qtrc )
      call phytends%Get3D( phytend_num1 + iq, rhoq_tp )
 
      do n=1, mesh3d%LOCAL_MESH_NUM
        lcmesh3d => mesh3d%lcmesh_list(n)
 
        !$omp parallel do
        do ke=lcmesh3d%NeS, lcmesh3d%NeE
          qtrc_tmp%local(n)%val(:,ke) = qtrc%local(n)%val(:,ke)
        end do            
      end do
 
      do rkstage=1, this%tint(1)%nstage
        if ( tracer_advc(iq) ) then
          call prof_rapstart( 'ATM_DYN_exchange_trc', 3)
          call this%TRCVAR3D_manager%MeshFieldComm_Exchange()
          call prof_rapend( 'ATM_DYN_exchange_trc', 3)
 
 
          do n=1, mesh3d%LOCAL_MESH_NUM
            lcmesh3d => mesh3d%lcmesh_list(n)
 
            dt = this%tint(n)%Get_deltime()
            dttmp_trc = dt * this%tint(n)%coef_gam_ex(rkstage+1,rkstage) &
                         / this%tint(n)%coef_sig_ex(rkstage+1,rkstage)
 
            call atm_dyn_dgm_trcadvect3d_heve_calc_fct_coef( &
              this%AUX_TRCVARS3D(auxtrcvars3d_fctcoef_id)%local(n)%val,                     & ! (out)
              qtrc_tmp%local(n)%val,                                                        & ! (in)
              mflx_x_tavg%local(n)%val, mflx_y_tavg%local(n)%val, mflx_z_tavg%local(n)%val, & ! (in) 
              rhoq_tp%local(n)%val,                                                         & ! (in)
              this%alphaDensM%local(n)%face_val, this%alphaDensP%local(n)%face_val,         & ! (in)
              dens_hyd%local(n)%val, ddens_trc%local(n)%val, ddens0_trc%local(n)%val,       & ! (in)
              this%tint(n)%coef_c_ex(rkstage), dttmp_trc,                                   & ! (in) 
              this%FaceIntMat,                                                              & ! (in)
              lcmesh3d, lcmesh3d%refElem3D, lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D, & ! (in)
              this%disable_limiter                                                          ) ! (in)
          end do
 
          call prof_rapstart( 'ATM_DYN_exchange_trc', 3)
          call this%AUXTRCVAR3D_manager%MeshFieldComm_Exchange()
          call prof_rapend( 'ATM_DYN_exchange_trc', 3)
        end if
 
        do n=1, mesh3d%LOCAL_MESH_NUM
          lcmesh3d => mesh3d%lcmesh_list(n)
 
          tintbuf_ind = this%tint(n)%tend_buf_indmap(rkstage)
 
          call prof_rapstart( 'ATM_DYN_update_caltend_ex_trc', 3)  
          if ( tracer_advc(iq) ) then 
            call atm_dyn_dgm_trcadvect3d_heve_cal_tend( &        
              this%tint(n)%tend_buf2D_ex(:,:,1,tintbuf_ind),                                & ! (out)
              qtrc_tmp%local(n)%val,                                                        & ! (in)
              mflx_x_tavg%local(n)%val, mflx_y_tavg%local(n)%val, mflx_z_tavg%local(n)%val, & ! (in) 
              this%alphaDensM%local(n)%face_val, this%alphaDensP%local(n)%face_val,         & ! (in)
              this%AUX_TRCVARS3D(auxtrcvars3d_fctcoef_id)%local(n)%val,                     & ! (out)
              rhoq_tp%local(n)%val,                                                         & ! (in)
              element_operation, this%FaceIntMat,                                           & ! (in)
              lcmesh3d, lcmesh3d%refElem3D, lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D  ) ! (in)
          else
            !$omp parallel do
            do ke=lcmesh3d%NeS, lcmesh3d%NeE
              this%tint(n)%tend_buf2D_ex(:,ke,1,tintbuf_ind) = rhoq_tp%local(n)%val(:,ke)
            end do
          end if
          call prof_rapend( 'ATM_DYN_update_caltend_ex_trc', 3)
 
          call prof_rapstart( 'ATM_DYN_update_advance_trc', 3)                
          call this%tint(n)%Advance_trcvar( &
            rkstage, qtrc_tmp%local(n)%val, 1,                                     &
            1, lcmesh3d%refElem%Np, lcmesh3d%NeS, lcmesh3d%NeE,                    &
            ddens_trc%local(n)%val, ddens0_trc%local(n)%val, dens_hyd%local(n)%val ) 
          call prof_rapend( 'ATM_DYN_update_advance_trc', 3)
 
          if ( rkstage == this%tint(1)%nstage .and. this%MODALFILTER_FLAG ) then
            call prof_rapstart( 'ATM_DYN_update_qtrc_modalfilter', 3)
            call atm_dyn_dgm_tracer_modalfilter_apply( &
              qtrc_tmp%local(n)%val,                                             & ! (inout)
              dens_hyd%local(n)%val, ddens_trc%local(n)%val, ddens%local(n)%val, & ! (in)
              lcmesh3d, lcmesh3d%refElem3D, element_operation                    ) ! (in)
            call prof_rapend( 'ATM_DYN_update_qtrc_modalfilter', 3)
          end if
 
          if ( tracer_advc(iq)                   &
            .and. rkstage == this%tint(1)%nstage &
            .and. this%ONLY_TRACERADV_FLAG       &            
            .and. ( .not. this%disable_limiter ) ) then
 
            call prof_rapstart( 'ATM_DYN_update_trc_TMAR', 3)             
            call atm_dyn_dgm_trcadvect3d_tmar( &
              qtrc_tmp%local(n)%val,                                                       & ! (inout)
              dens_hyd%local(n)%val, ddens_trc%local(n)%val,                               & ! (in)
              lcmesh3d, lcmesh3d%refElem3D, lcmesh3d%lcmesh2D, lcmesh3d%lcmesh2D%refElem2D ) ! (in)
            call prof_rapend( 'ATM_DYN_update_trc_TMAR', 3)  
          end if
 
        end do
      end do ! end for RK loop
 
      do n=1, mesh3d%LOCAL_MESH_NUM
        lcmesh3d => mesh3d%lcmesh_list(n)
        !$omp parallel do
        do ke=lcmesh3d%NeS, lcmesh3d%NeE
          qtrc%local(n)%val(:,ke) = ( dens_hyd%local(n)%val(:,ke) + ddens_trc%local(n)%val(:,ke) ) &
                                  / ( dens_hyd%local(n)%val(:,ke) + ddens%local(n)%val(:,ke) )     &
                                  * qtrc_tmp%local(n)%val(:,ke)
        end do            
      end do
 
    end do ! end do for iq
 
    ! Update pressure
    call update_pressure_specific_heat( pres, rtot, cvtot, cptot, & ! (inout)
      trc_vars, prog_vars, aux_vars, dyn_driver, mesh3d           ) ! (in)
 
    ! Negative fixer
    if (         .not. this%disable_limiter       & 
         .and. ( .not. this%ONLY_TRACERADV_FLAG ) ) then
 
      call prof_rapstart( 'ATM_DYN_trc_negative_fixer', 2)          
      call fix_negative_val( &
        trc_vars, ddens, therm, pres, cvtot, cptot, rtot,               & ! (inout)
        dens_hyd, pres_hyd, dt, mesh3d, dyn_driver%Is_THERMVAR_RHOT()   ) ! (in)
      call prof_rapend( 'ATM_DYN_trc_negative_fixer', 2)
    end if
 
    call prof_rapend( 'ATM_DYN_trc_update', 2)
 
    return
  end subroutine atmdyndgmdriver_trcadv3d_update
 
!OCL SERIAL
  subroutine atmdyndgmdriver_trcadv3d_final( this )
    implicit none
 
    class(atmdyndgmdriver_trcadv3d), intent(inout) :: this
    !-----------------------------------------------------------------------------
 
    call atm_dyn_dgm_trcadvect3d_heve_final()
 
    nullify( this%boundary_cond )
    
    call this%TRCVAR3D_manager%Final()
    deallocate( this%TRCVARS3D )
 
    call this%AUXTRCVAR3D_manager%Final()
    deallocate( this%AUX_TRCVARS3D )
 
    call this%AUXTRC_FLUX_VAR3D_manager%Final()
    deallocate( this%AUXTRC_FLUX_VARS3D )
 
    call this%alphaDensM%Final()
    call this%alphaDensP%Final()
 
    if ( this%MODALFILTER_FLAG ) then
      call this%modal_filter_3d%Final()
    end if
    
    call this%FaceIntMat%Final()
    call atmdyndgmdriver_base3d_final( this )   
 
    return
  end subroutine atmdyndgmdriver_trcadv3d_final
 
!-- private
!OCL SERIAL
  subroutine fix_negative_val( &
    TRC_VARS, DDENS, THERM, PRES, CVtot, CPtot, Rtot, & ! (inout)
    dens_hyd, pres_hyd, dt, mesh3d, is_thermvar_rhot  ) ! (in)
    use scale_atmos_hydrometeor, only: &
      qla, qia
    use scale_atm_phy_mp_dgm_common, only: &
      atm_phy_mp_dgm_common_negative_fixer
    implicit none
    class(modelvarmanager), intent(inout) :: trc_vars
    class(meshfield3d), intent(inout) :: ddens
    class(meshfield3d), intent(inout) :: therm
    class(meshfield3d), intent(inout) :: pres
    class(meshfield3d), intent(inout) :: cvtot
    class(meshfield3d), intent(inout) :: cptot
    class(meshfield3d), intent(inout) :: rtot
    class(meshfield3d), intent(inout) :: dens_hyd
    class(meshfield3d), intent(inout) :: pres_hyd
    real(rp), intent(in) :: dt
    class(meshbase3d), intent(in), target :: mesh3d
    logical, intent(in) :: is_thermvar_rhot
 
    class(localmesh3d), pointer :: lcmesh3d
    integer :: n
 
    integer :: iq
 
    integer :: trcid_list(qa)
    type(localmeshfieldbaselist) :: lc_qtrc(qa)
    !---------------------------------------------------------------
 
    do iq=1, qa 
      trcid_list(iq) = iq
    end do
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      call trc_vars%GetLocalMeshFieldList( trcid_list, n, lc_qtrc )
      lcmesh3d => mesh3d%lcmesh_list(n)
 
      if ( is_thermvar_rhot ) then
        call atm_phy_mp_dgm_common_negative_fixer( &
          lc_qtrc, ddens%local(n)%val, pres%local(n)%val,            & ! (inout)
          cvtot%local(n)%val, cptot%local(n)%val, rtot%local(n)%val, & ! (inout)
          dens_hyd%local(n)%val, pres_hyd%local(n)%val,              & ! (in)
          dt, lcmesh3d, lcmesh3d%refElem3D, qa, qla, qia,            & ! (in)
          drhot=therm%local(n)%val                                   ) ! (inout)
      else
        call atm_phy_mp_dgm_common_negative_fixer( &
          lc_qtrc, ddens%local(n)%val, pres%local(n)%val,            & ! (inout)
          cvtot%local(n)%val, cptot%local(n)%val, rtot%local(n)%val, & ! (inout)
          dens_hyd%local(n)%val, pres_hyd%local(n)%val,              & ! (in)
          dt, lcmesh3d, lcmesh3d%refElem3D, qa, qla, qia             ) ! (in)
      end if     
    end do
 
    return
  end subroutine fix_negative_val
 
  subroutine update_pressure_specific_heat( &
    PRES, Rtot, CVtot, CPtot,                         &
    TRC_VARS, PROG_VARS, AUX_VARS, dyn_driver, mesh3D )
 
    use scale_tracer, only: &
      tracer_mass, tracer_r, tracer_cv, tracer_cp    
    use scale_atmos_thermodyn, only: &
      atmos_thermodyn_specific_heat
    use scale_atm_dyn_dgm_nonhydro3d_common, only: &
      atm_dyn_dgm_nonhydro3d_common_calc_pressure      
    implicit none
 
    class(meshfield3d), intent(inout) :: pres
    class(meshfield3d), intent(inout) :: rtot
    class(meshfield3d), intent(inout) :: cvtot
    class(meshfield3d), intent(inout) :: cptot
    class(modelvarmanager), intent(inout) :: trc_vars
    class(modelvarmanager), intent(inout) :: prog_vars
    class(modelvarmanager), intent(inout) :: aux_vars
    class(atmdyndgmdriver_nonhydro3d), intent(inout) :: dyn_driver
    class(meshbase3d), intent(in), target :: mesh3d
 
    integer :: n
    integer :: ke
    integer :: iq
 
    class(localmesh3d), pointer :: lcmesh3d
    class(elementbase3d), pointer :: elem3d
 
    real(rp), allocatable :: q_tmp(:,:)
    real(rp), allocatable :: qdry(:)
 
    integer :: trcid_list(qa)
    type(localmeshfieldbaselist) :: lc_qtrc(qa)
    !---------------------------------------------------------------
 
    do iq=1, qa 
      trcid_list(iq) = iq
    end do
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      call trc_vars%GetLocalMeshFieldList( trcid_list, n, lc_qtrc )
      
      lcmesh3d => mesh3d%lcmesh_list(n)
      elem3d => lcmesh3d%refElem3D
      allocate( q_tmp(elem3d%Np,qa), qdry(elem3d%Np) )
 
      !$omp parallel do private(ke, iq, q_tmp, qdry)
      do ke = lcmesh3d%NeS, lcmesh3d%NeE
        do iq = 1, qa
          q_tmp(:,iq) = lc_qtrc(iq)%ptr%val(:,ke)
        end do
        call atmos_thermodyn_specific_heat( &
          elem3d%Np, 1, elem3d%Np, qa,                                                      & ! (in)
          q_tmp(:,:), tracer_mass(:), tracer_r(:), tracer_cv(:), tracer_cp(:),              & ! (in)
          qdry, rtot%local(n)%val(:,ke), cvtot%local(n)%val(:,ke), cptot%local(n)%val(:,ke) ) ! (out)
      end do
      deallocate(q_tmp, qdry)      
    end do
 
    call dyn_driver%calc_pressure( pres, prog_vars, aux_vars )
 
    return
  end subroutine update_pressure_specific_heat
 
  !-- Setup modal filter
!OCL SERIAL
  subroutine setup_modalfilter( this, refElem3D, element_operation )
    use scale_element_line, only: lineelement
    implicit none
 
    class(atmdyndgmdriver_trcadv3d), target, intent(inout) :: this
    class(hexahedralelement), target, intent(in) :: refelem3d
    class(elementoperationbase3d), intent(inout) :: element_operation
 
    real(rp) :: mf_etac_h  = 2.0_rp/3.0_rp
    real(rp) :: mf_alpha_h = 36.0_rp
    integer  :: mf_order_h = 16
    real(rp) :: mf_etac_v  = 2.0_rp/3.0_rp
    real(rp) :: mf_alpha_v = 36.0_rp
    integer  :: mf_order_v = 16
 
    namelist /param_atmos_dyn_tracer_modalfilter/ &
      mf_etac_h, mf_alpha_h, mf_order_h,   &
      mf_etac_v, mf_alpha_v, mf_order_v    
 
    integer :: ierr
 
    type(lineelement) :: elemv1d
    !---------------------------------------------------------------
 
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_dyn_tracer_modalfilter,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_DYN_TRCADV3D_setup_modalfilter",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_DYN_TRCADV3D_setup_modalfilter",*) 'Invalid names in namelist PARAM_ATMOS_DYN_TRACER_MODALFILTER. Check!'
      call prc_abort
    endif
 
    log_nml(param_atmos_dyn_tracer_modalfilter)
 
    call this%modal_filter_3d%Init( &
        refelem3d,                           & ! (in)
        mf_etac_h, mf_alpha_h, mf_order_h,   & ! (in)
        mf_etac_v, mf_alpha_v, mf_order_v    ) ! (in)
 
    call element_operation%Setup_ModalFilter_tracer( &
      mf_etac_h, mf_alpha_h, mf_order_h,   & ! (in)
      mf_etac_v, mf_alpha_v, mf_order_v    ) ! (in)
 
    return
  end subroutine setup_modalfilter
 
end module scale_atm_dyn_dgm_driver_trcadv3d