scale_atmos_dyn_tstep_large_fvm_heve.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_dyn_tstep_large_fvm_heve
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  use scale_index
  use scale_tracer

#ifdef DEBUG
  use scale_debug, only: &
     check
  use scale_const, only: &
     undef  => const_undef, &
     iundef => const_undef2
#endif
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_dyn_tstep_large_fvm_heve_setup
  public :: atmos_dyn_tstep_large_fvm_heve

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  ! tendency
  real(RP), private, allocatable :: dens_t(:,:,:)
  real(RP), private, allocatable :: momz_t(:,:,:)
  real(RP), private, allocatable :: momx_t(:,:,:)
  real(RP), private, allocatable :: momy_t(:,:,:)
  real(RP), private, allocatable :: rhot_t(:,:,:)
  real(RP), private, allocatable, target :: rhoq_t(:,:,:,:)
  real(RP), private, allocatable, target :: zero(:,:,:)
  real(RP), private, pointer     :: rhoq_tn(:,:,:)

  real(RP), private, allocatable :: mflx_hi(:,:,:,:)        ! rho * vel(x,y,z) @ (u,v,w)-face high order


  ! for communication
  integer :: i_comm_dens = 1
  integer :: i_comm_momz = 2
  integer :: i_comm_momx = 3
  integer :: i_comm_momy = 4
  integer :: i_comm_rhot = 5
  integer, allocatable :: i_comm_prog(:)

  integer :: i_comm_dens_t = 1
  integer :: i_comm_momz_t = 2
  integer :: i_comm_momx_t = 3
  integer :: i_comm_momy_t = 4
  integer :: i_comm_rhot_t = 5

  integer, allocatable :: i_comm_rhoq_t(:)
  integer, allocatable :: i_comm_qtrc(:)

  integer :: i_comm_mflx_z = 1
  integer :: i_comm_mflx_x = 2
  integer :: i_comm_mflx_y = 3

  !-----------------------------------------------------------------------------
contains

  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_tstep_large_fvm_heve_setup( &
       DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, PROG, &
       mflx_hi )
    use scale_prc, only: &
       prc_abort
    use scale_prc_cartesc, only: &
       prc_has_e, &
       prc_has_w, &
       prc_has_n, &
       prc_has_s
    use scale_const, only: &
       ohm => const_ohm, &
       undef => const_undef
    use scale_comm_cartesc, only: &
       comm_vars8_init
    implicit none

    ! MPI_RECV_INIT requires intent(inout)
    real(RP),               intent(inout) :: DENS(ka,ia,ja)
    real(RP),               intent(inout) :: MOMZ(ka,ia,ja)
    real(RP),               intent(inout) :: MOMX(ka,ia,ja)
    real(RP),               intent(inout) :: MOMY(ka,ia,ja)
    real(RP),               intent(inout) :: RHOT(ka,ia,ja)
    real(RP),               intent(inout) :: QTRC(ka,ia,ja,qa)
    real(RP),               intent(inout) :: PROG(ka,ia,ja,va)
    real(RP),               intent(inout) :: mflx_hi(ka,ia,ja,3)

    integer :: iv, iq
    !---------------------------------------------------------------------------

    allocate( dens_t(ka,ia,ja) )
    allocate( momz_t(ka,ia,ja) )
    allocate( momx_t(ka,ia,ja) )
    allocate( momy_t(ka,ia,ja) )
    allocate( rhot_t(ka,ia,ja) )
    allocate( rhoq_t(ka,ia,ja,qa) )

    allocate( i_comm_prog(max(va,1)) )
    allocate( i_comm_qtrc(qa) )
    allocate( i_comm_rhoq_t(qa) )

    call comm_vars8_init( 'DENS', dens, i_comm_dens )
    call comm_vars8_init( 'MOMZ', momz, i_comm_momz )
    call comm_vars8_init( 'MOMX', momx, i_comm_momx )
    call comm_vars8_init( 'MOMY', momy, i_comm_momy )
    call comm_vars8_init( 'RHOT', rhot, i_comm_rhot )
    do iv = 1, va
       i_comm_prog(iv) = 5 + iv
       call comm_vars8_init( 'PROG', prog(:,:,:,iv), i_comm_prog(iv) )
    end do

    call comm_vars8_init( 'DENS_t', dens_t, i_comm_dens_t )
    call comm_vars8_init( 'MOMZ_t', momz_t, i_comm_momz_t )
    call comm_vars8_init( 'MOMX_t', momx_t, i_comm_momx_t )
    call comm_vars8_init( 'MOMY_t', momy_t, i_comm_momy_t )
    call comm_vars8_init( 'RHOT_t', rhot_t, i_comm_rhot_t )

    do iq = 1, qa
       i_comm_rhoq_t(iq) = 5 + va + iq
       i_comm_qtrc(iq) = 5 + va + iq

       call comm_vars8_init( 'RHOQ_t', rhoq_t(:,:,:,iq), i_comm_rhoq_t(iq) )
       call comm_vars8_init( 'QTRC',   qtrc(:,:,:,iq), i_comm_qtrc(iq) )
    end do

    call comm_vars8_init( 'mflx_Z', mflx_hi(:,:,:,zdir), i_comm_mflx_z )
    call comm_vars8_init( 'mflx_X', mflx_hi(:,:,:,xdir), i_comm_mflx_x )
    call comm_vars8_init( 'mflx_Y', mflx_hi(:,:,:,ydir), i_comm_mflx_y )

    allocate( zero(ka,ia,ja) )
    zero(:,:,:) = 0.0_rp

    mflx_hi(:,:,:,:) = undef

    return
  end subroutine atmos_dyn_tstep_large_fvm_heve_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_tstep_large_fvm_heve( &
       DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, PROG,             &
       DENS_av, MOMZ_av, MOMX_av, MOMY_av, RHOT_av, QTRC_av, &
       mflx_hi, tflx_hi,                                     &
       num_diff, num_diff_q,                                 &
       QTRC0,                                                &
       DENS_tp, MOMZ_tp, MOMX_tp, MOMY_tp, RHOT_tp, RHOQ_tp, &
       CORIOLI,                                              &
       CDZ, CDX, CDY, FDZ, FDX, FDY,                         &
       RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY,                   &
       PHI, GSQRT,                                           &
       J13G, J23G, J33G, MAPF,                               &
       AQ_R, AQ_CV, AQ_CP, AQ_MASS,                          &
       REF_dens, REF_pott, REF_qv, REF_pres,                 &
       BND_W, BND_E, BND_S, BND_N,                           &
       ND_COEF, ND_COEF_Q, ND_ORDER, ND_SFC_FACT, ND_USE_RS, &
       BND_QA, BND_SMOOTHER_FACT,                            &
       DAMP_DENS,       DAMP_VELZ,       DAMP_VELX,          &
       DAMP_VELY,       DAMP_POTT,       DAMP_QTRC,          &
       DAMP_alpha_DENS, DAMP_alpha_VELZ, DAMP_alpha_VELX,    &
       DAMP_alpha_VELY, DAMP_alpha_POTT, DAMP_alpha_QTRC,    &
       wdamp_coef,                                           &
       divdmp_coef,                                          &
       FLAG_TRACER_SPLIT_TEND,                               &
       FLAG_FCT_MOMENTUM, FLAG_FCT_T, FLAG_FCT_TRACER,       &
       FLAG_FCT_ALONG_STREAM,                                &
       USE_AVERAGE,                                          &
       I_QV,                                                 &
       DTLS, DTSS, Llast                                     )
    use scale_const, only: &
       eps    => const_eps, &
       p0     => const_pre00, &
       rdry   => const_rdry, &
       cvdry  => const_cvdry, &
       cpdry  => const_cpdry
    use scale_comm_cartesc, only: &
       comm_vars8, &
       comm_wait
    use scale_atmos_dyn_common, only: &
       atmos_dyn_divergence, &
       atmos_dyn_numfilter_coef,   &
       atmos_dyn_numfilter_coef_q, &
       atmos_dyn_fct
    use scale_atmos_dyn_fvm_flux_ud1, only: &
       atmos_dyn_fvm_fluxz_xyz_ud1, &
       atmos_dyn_fvm_fluxx_xyz_ud1, &
       atmos_dyn_fvm_fluxy_xyz_ud1
    use scale_atmos_dyn_fvm_flux, only: &
       atmos_dyn_fvm_fluxz_xyz, &
       atmos_dyn_fvm_fluxx_xyz, &
       atmos_dyn_fvm_fluxy_xyz
    use scale_file_history, only: &
#ifdef hist_tend
       file_history_in, &
#endif
       file_history_set_disable
    use scale_atmos_dyn_tinteg_short, only: &
       atmos_dyn_tinteg_short
    use scale_atmos_dyn_tinteg_tracer, only: &
       atmos_dyn_tinteg_tracer
    implicit none

    real(RP), intent(inout) :: DENS(ka,ia,ja)
    real(RP), intent(inout) :: MOMZ(ka,ia,ja)
    real(RP), intent(inout) :: MOMX(ka,ia,ja)
    real(RP), intent(inout) :: MOMY(ka,ia,ja)
    real(RP), intent(inout) :: RHOT(ka,ia,ja)
    real(RP), intent(inout) :: QTRC(ka,ia,ja,qa)
    real(RP), intent(inout) :: PROG(ka,ia,ja,va)

    real(RP), intent(inout) :: DENS_av(ka,ia,ja)
    real(RP), intent(inout) :: MOMZ_av(ka,ia,ja)
    real(RP), intent(inout) :: MOMX_av(ka,ia,ja)
    real(RP), intent(inout) :: MOMY_av(ka,ia,ja)
    real(RP), intent(inout) :: RHOT_av(ka,ia,ja)
    real(RP), intent(inout) :: QTRC_av(ka,ia,ja,qa)

    real(RP), intent(out)   :: mflx_hi(ka,ia,ja,3)
    real(RP), intent(out)   :: tflx_hi(ka,ia,ja,3)

    real(RP), intent(out)   :: num_diff(ka,ia,ja,5,3)
    real(RP), intent(out)   :: num_diff_q(ka,ia,ja,3)

    real(RP), intent(in)    :: QTRC0(ka,ia,ja,qa)

    real(RP), intent(in)    :: DENS_tp(ka,ia,ja)
    real(RP), intent(in)    :: MOMZ_tp(ka,ia,ja)
    real(RP), intent(in)    :: MOMX_tp(ka,ia,ja)
    real(RP), intent(in)    :: MOMY_tp(ka,ia,ja)
    real(RP), intent(in)    :: RHOT_tp(ka,ia,ja)
    real(RP), intent(in)    :: RHOQ_tp(ka,ia,ja,qa)

    real(RP), intent(in)    :: CORIOLI(ia,ja)

    real(RP), intent(in)    :: CDZ (ka)
    real(RP), intent(in)    :: CDX (ia)
    real(RP), intent(in)    :: CDY (ja)
    real(RP), intent(in)    :: FDZ (ka-1)
    real(RP), intent(in)    :: FDX (ia-1)
    real(RP), intent(in)    :: FDY (ja-1)
    real(RP), intent(in)    :: RCDZ(ka)
    real(RP), intent(in)    :: RCDX(ia)
    real(RP), intent(in)    :: RCDY(ja)
    real(RP), intent(in)    :: RFDZ(ka-1)
    real(RP), intent(in)    :: RFDX(ia-1)
    real(RP), intent(in)    :: RFDY(ja-1)

    real(RP), intent(in)    :: PHI  (ka,ia,ja)
    real(RP), intent(in)    :: GSQRT(ka,ia,ja,7)
    real(RP), intent(in)    :: J13G (ka,ia,ja,7)
    real(RP), intent(in)    :: J23G (ka,ia,ja,7)
    real(RP), intent(in)    :: J33G
    real(RP), intent(in)    :: MAPF (ia,ja,2,4)

    real(RP), intent(in)    :: AQ_R   (qa)
    real(RP), intent(in)    :: AQ_CV  (qa)
    real(RP), intent(in)    :: AQ_CP  (qa)
    real(RP), intent(in)    :: AQ_MASS(qa)

    real(RP), intent(in)    :: REF_dens(ka,ia,ja)
    real(RP), intent(in)    :: REF_pott(ka,ia,ja)
    real(RP), intent(in)    :: REF_qv  (ka,ia,ja)
    real(RP), intent(in)    :: REF_pres(ka,ia,ja)

    logical,  intent(in)    :: BND_W
    logical,  intent(in)    :: BND_E
    logical,  intent(in)    :: BND_S
    logical,  intent(in)    :: BND_N

    real(RP), intent(in)    :: ND_COEF
    real(RP), intent(in)    :: ND_COEF_Q
    integer,  intent(in)    :: ND_ORDER
    real(RP), intent(in)    :: ND_SFC_FACT
    logical,  intent(in)    :: ND_USE_RS

    integer,  intent(in)    :: BND_QA
    real(RP), intent(in)    :: BND_SMOOTHER_FACT

    real(RP), intent(in)    :: DAMP_DENS(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_VELZ(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_VELX(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_VELY(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_POTT(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_QTRC(ka,ia,ja,bnd_qa)

    real(RP), intent(in)    :: DAMP_alpha_DENS(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_alpha_VELZ(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_alpha_VELX(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_alpha_VELY(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_alpha_POTT(ka,ia,ja)
    real(RP), intent(in)    :: DAMP_alpha_QTRC(ka,ia,ja,bnd_qa)

    real(RP), intent(in)    :: wdamp_coef(ka)
    real(RP), intent(in)    :: divdmp_coef

    logical,  intent(in)    :: FLAG_TRACER_SPLIT_TEND
    logical,  intent(in)    :: FLAG_FCT_MOMENTUM
    logical,  intent(in)    :: FLAG_FCT_T
    logical,  intent(in)    :: FLAG_FCT_TRACER
    logical,  intent(in)    :: FLAG_FCT_ALONG_STREAM

    logical,  intent(in)    :: USE_AVERAGE

    integer,  intent(in)    :: I_QV

    real(DP), intent(in)    :: DTLS
    real(DP), intent(in)    :: DTSS
    logical , intent(in)    :: Llast

    ! for time integartion
    real(RP) :: DENS00  (ka,ia,ja) ! saved density before small step loop

    ! diagnostic variables
    real(RP) :: DDIV    (ka,ia,ja) ! 3 dimensional divergence
    real(RP) :: DPRES0  (ka,ia,ja) ! pressure deviation
    real(RP) :: RT2P    (ka,ia,ja) ! factor of RHOT to PRES
    real(RP) :: REF_rhot(ka,ia,ja) ! reference of RHOT

    real(RP) :: QDRY  ! dry air
    real(RP) :: Rtot  ! total R
    real(RP) :: CVtot ! total CV
    real(RP) :: CPtot ! total CP
    real(RP) :: PRES  ! pressure
#ifdef DRY
    real(RP) :: CPovCV
#endif

    real(RP) :: DENS_tq(ka,ia,ja)
    real(RP) :: diff(ka,ia,ja)
    real(RP) :: damp
#ifdef HIST_TEND
    real(RP) :: damp_t(ka,ia,ja)
#endif

    ! For tracer advection
    real(RP) :: mflx_av  (ka,ia,ja,3)  ! rho * vel(x,y,z) @ (u,v,w)-face average

    real(RP) :: dtl
    real(RP) :: dts
    integer  :: nstep

    integer  :: i, j, k, iq, step
    integer  :: iv
    integer  :: n
    !---------------------------------------------------------------------------

    call prof_rapstart("DYN_Large_Preparation", 2)

    dts   = real(DTSS, kind=RP)            ! short time step
    dtl   = real(DTLS, kind=RP)            ! large time step
    nstep = ceiling( ( dtl - eps ) / dts )
    dts   = dtl / nstep                    ! dts is divisor of dtl and smaller or equal to dtss

#ifdef DEBUG
    log_info("ATMOS_DYN_Tstep_large_fvm_heve",*)                         'Dynamics large time step'
    log_info_cont('(1x,A,F0.2,A,F0.2,A,I0)') &
    '-> DT_large, DT_small, DT_large/DT_small : ', dtl, ', ', dts, ', ', nstep

    dens00(:,:,:) = undef

    num_diff(:,:,:,:,:) = undef

    mflx_hi(:,:,:,:) = undef
    tflx_hi(:,:,:,:) = undef
#endif

!OCL XFILL
    dens00(:,:,:) = dens(:,:,:)

    if ( use_average ) then
!OCL XFILL
       dens_av(:,:,:) = 0.0_rp
!OCL XFILL
       momz_av(:,:,:) = 0.0_rp
!OCL XFILL
       momx_av(:,:,:) = 0.0_rp
!OCL XFILL
       momy_av(:,:,:) = 0.0_rp
!OCL XFILL
       rhot_av(:,:,:) = 0.0_rp
    endif

#ifndef DRY
!OCL XFILL
    mflx_av(:,:,:,:) = 0.0_rp

#endif

#ifdef DRY
    cpovcv = cpdry / cvdry
#endif
    !------------------------------------------------------------------------
    ! prepare thermodynamical data
    !   specific heat
    !   pressure data ( linearization )
    !
    ! pres = P0 * ( R * rhot / P0 )**(CP/CV)
    ! d pres / d rhot ~ CP*R/CV * ( R * rhot / P0 )**(R/CV)
    !                 = CP*R/CV * ( pres / P0 )**(R/CP)
    !                 = CP*R/CV * temp / pott
    !                 = CP/CV * pres / rhot
    ! pres ~ P0 * ( R * rhot0 / P0 ) ** (CP/CV) + CV*R/CP * ( pres / P0 )**(R/CP) * rhot'
    !------------------------------------------------------------------------
!OCL XFILL
    !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
    !$omp shared(JA,IA,KS,KE) &
    !$omp shared(P0,Rdry,RHOT,AQ_R,AQ_CV,AQ_CP,QTRC,AQ_MASS,REF_rhot,REF_pres,CPdry,CVdry,QA,RT2P,DPRES0) &
#ifdef DRY
    !$omp shared(CPovCV) &
#endif
    !$omp private(i,j,k,iq) &
    !$omp private(PRES,Rtot,CVtot,CPtot,QDRY)
    do j = 1, ja
    do i = 1, ia
       do k = ks, ke
#ifdef DRY
          pres = p0 * ( rdry * rhot(k,i,j) / p0 )**cpovcv
          rt2p(k,i,j) = cpovcv * pres / rhot(k,i,j)
#else
          rtot  = 0.0_rp
          cvtot = 0.0_rp
          cptot = 0.0_rp
          qdry  = 1.0_rp
          do iq = 1, qa
             rtot  = rtot + aq_r(iq) * qtrc(k,i,j,iq)
             cvtot = cvtot + aq_cv(iq) * qtrc(k,i,j,iq)
             cptot = cptot + aq_cp(iq) * qtrc(k,i,j,iq)
             qdry  = qdry  - qtrc(k,i,j,iq) * aq_mass(iq)
          enddo
          rtot  = rtot  + rdry  * qdry
          cvtot = cvtot + cvdry * qdry
          cptot = cptot + cpdry * qdry
          pres = p0 * ( rtot * rhot(k,i,j) / p0 )**( cptot / cvtot )
          rt2p(k,i,j) = cptot / cvtot * pres / rhot(k,i,j)
#endif
          dpres0(k,i,j) = pres - ref_pres(k,i,j)
          ref_rhot(k,i,j) = rhot(k,i,j)
       end do
       dpres0(ks-1,i,j) = dpres0(ks+1,i,j) + ( ref_pres(ks+1,i,j) - ref_pres(ks-1,i,j) )
       dpres0(ke+1,i,j) = dpres0(ke-1,i,j) + ( ref_pres(ke-1,i,j) - ref_pres(ke+1,i,j) )
    end do
    end do

    call prof_rapend  ("DYN_Large_Preparation", 2)

    !###########################################################################
    ! Update DENS,MONZ,MOMX,MOMY,MOMZ,RHOT
    !###########################################################################

    call prof_rapstart("DYN_Large_Tendency", 2)

#ifdef HIST_TEND
!OCL XFILL
    damp_t(:,:,:) = 0.0_rp
#endif

!OCL XFILL
    dens_tq(:,:,:) = 0.0_rp

    do iq = 1, qa

       if ( iq >= i_qv .and. iq < i_qv+bnd_qa ) then

          !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
!OCL XFILL
          do j = js-1, je+1
          do i = is-1, ie+1
          do k = ks, ke
             diff(k,i,j) = qtrc(k,i,j,iq) - damp_qtrc(k,i,j,iq-i_qv+1)
          enddo
          enddo
          enddo
          !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
          !$omp private(i,j,k,damp) &
#ifdef HIST_TEND
          !$omp shared(damp_t) &
#endif
          !$omp shared(JS,JE,IS,IE,KS,KE,iq) &
          !$omp shared(RHOQ_t,RHOQ_tp,DENS_tq,DAMP_alpha_QTRC,diff,BND_SMOOTHER_FACT,DENS00,TRACER_MASS,I_QV)
!OCL XFILL
          do j = js, je
          do i = is, ie
          do k = ks, ke
             damp = - damp_alpha_qtrc(k,i,j,iq-i_qv+1) &
                  * ( diff(k,i,j) & ! rayleigh damping
                    - ( diff(k,i-1,j) + diff(k,i+1,j) + diff(k,i,j-1) + diff(k,i,j+1) - diff(k,i,j)*4.0_rp ) &
                    * 0.125_rp * bnd_smoother_fact ) ! horizontal smoother
#ifdef HIST_TEND
             damp_t(k,i,j) = damp
#endif
             damp = damp * dens00(k,i,j)
             rhoq_t(k,i,j,iq) = rhoq_tp(k,i,j,iq) + damp
             dens_tq(k,i,j) = dens_tq(k,i,j) + damp * tracer_mass(iq) ! only for mass tracer
          enddo
          enddo
          enddo
#ifdef HIST_TEND
          call file_history_in(rhoq_tp(:,:,:,iq), trim(tracer_name(iq))//'_t_phys', &
               'tendency of '//trim(tracer_name(iq))//' due to physics (w/ HIST_TEND)', 'kg/kg/s' )
          call file_history_in(damp_t,            trim(tracer_name(iq))//'_t_damp', &
               'tendency of '//trim(tracer_name(iq))//' due to damping (w/ HIST_TEND)', 'kg/kg/s' )
#endif
!OCL XFILL
          do j = js, je
          do i = is, ie
             rhoq_t(   1:ks-1,i,j,iq) = 0.0_rp
             rhoq_t(ke+1:ka  ,i,j,iq) = 0.0_rp
          enddo
          enddo

          call comm_vars8( rhoq_t(:,:,:,iq), i_comm_rhoq_t(iq) )
          call comm_wait ( rhoq_t(:,:,:,iq), i_comm_rhoq_t(iq), .false. )

       else

          !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(3)
!OCL XFILL
          do j = 1, ja
          do i = 1, ia
          do k = 1, ka
             rhoq_t(k,i,j,iq) = rhoq_tp(k,i,j,iq)
          enddo
          enddo
          enddo

       end if

    end do

    call prof_rapend  ("DYN_Large_Tendency", 2)

    call prof_rapstart("DYN_Large_Boundary", 2)

    if ( bnd_w ) then
       !$omp parallel do private(j,k) OMP_SCHEDULE_ collapse(2)
       do j = js, je
       do k = ks, ke
          mflx_hi(k,is-1,j,xdir) = gsqrt(k,is-1,j,i_uyz) * momx(k,is-1,j) / mapf(is-1,j,2,i_uy)
       enddo
       enddo
    end if
    if ( bnd_e ) then
       !$omp parallel do private(j,k) OMP_SCHEDULE_ collapse(2)
       do j = js, je
       do k = ks, ke
          mflx_hi(k,ie,j,xdir) = gsqrt(k,ie,j,i_uyz) * momx(k,ie,j) / mapf(ie,j,2,i_uy)
       enddo
       enddo
    end if
    if ( bnd_s ) then
       !$omp parallel do private(i,k) OMP_SCHEDULE_ collapse(2)
       do i = is, ie
       do k = ks, ke
          mflx_hi(k,i,js-1,ydir) = gsqrt(k,i,js-1,i_xvz) * momy(k,i,js-1) / mapf(i,js-1,1,i_xv)
       enddo
       enddo
    end if
    if ( bnd_n ) then
       !$omp parallel do private(i,k) OMP_SCHEDULE_ collapse(2)
       do i = is, ie
       do k = ks, ke
          mflx_hi(k,i,je,ydir) = gsqrt(k,i,je,i_xvz) * momy(k,i,je) / mapf(i,je,1,i_xv)
       enddo
       enddo
    end if

    call prof_rapend  ("DYN_Large_Boundary", 2)


    do step = 1, nstep

       call file_history_set_disable( .not. (llast .AND. step == nstep ) )

       !-----< prepare tendency >-----

       call prof_rapstart("DYN_Large_Tendency", 2)

       !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
!OCL XFILL
       do j = js-1, je+1
       do i = is-1, ie+1
       do k = ks, ke
          diff(k,i,j) = dens(k,i,j) - damp_dens(k,i,j)
       enddo
       enddo
       enddo
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,damp) &
#ifdef HIST_TEND
       !$omp shared(damp_t) &
#endif
       !$omp shared(JS,JE,IS,IE,KS,KE) &
       !$omp shared(DAMP_alpha_DENS,diff,DENS_tq,DENS_t,DENS_tp,BND_SMOOTHER_FACT)
!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ks, ke
          damp = - damp_alpha_dens(k,i,j) &
               * ( diff(k,i,j) & ! rayleigh damping
                 - ( diff(k,i-1,j) + diff(k,i+1,j) + diff(k,i,j-1) + diff(k,i,j+1) - diff(k,i,j)*4.0_rp ) &
                 * 0.125_rp * bnd_smoother_fact ) & ! horizontal smoother
               + dens_tq(k,i,j) ! dencity change due to rayleigh damping for tracers
          dens_t(k,i,j) = dens_tp(k,i,j) & ! tendency from physical step
                        + damp
#ifdef HIST_TEND
          damp_t(k,i,j) = damp
#endif
       enddo
       enddo
       enddo
!OCL XFILL
       do j = js, je
       do i = is, ie
          dens_t(   1:ks-1,i,j) = 0.0_rp
          dens_t(ke+1:ka  ,i,j) = 0.0_rp
       enddo
       enddo
       call comm_vars8( dens_t(:,:,:), i_comm_dens_t )
#ifdef HIST_TEND
       call file_history_in(dens_tp, 'DENS_t_phys', 'tendency of dencity due to physics (w/ HIST_TEND)', 'kg/m3/s' )
       call file_history_in(damp_t,  'DENS_t_damp', 'tendency of dencity due to damping (w/ HIST_TEND)', 'kg/m3/s' )
#endif

       !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
!OCL XFILL
       do j = js-1, je+1
       do i = is-1, ie+1
       do k = ks, ke-1
          diff(k,i,j) = momz(k,i,j) - damp_velz(k,i,j) * ( dens(k,i,j)+dens(k+1,i,j) ) * 0.5_rp
       enddo
       enddo
       enddo
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,damp) &
#ifdef HIST_TEND
       !$omp shared(damp_t) &
#endif
       !$omp shared(JS,JE,IS,IE,KS,KE) &
       !$omp shared(DAMP_alpha_VELZ,diff,BND_SMOOTHER_FACT,MOMZ_t,MOMZ_tp)
!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ks, ke-1
          damp = - damp_alpha_velz(k,i,j) &
               * ( diff(k,i,j) & ! rayleigh damping
                 - ( diff(k,i-1,j) + diff(k,i+1,j) + diff(k,i,j-1) + diff(k,i,j+1) - diff(k,i,j)*4.0_rp ) &
                 * 0.125_rp * bnd_smoother_fact ) ! horizontal smoother

          momz_t(k,i,j) = momz_tp(k,i,j) & ! tendency from physical step
                        + damp
#ifdef HIST_TEND
          damp_t(k,i,j) = damp
#endif
       enddo
       enddo
       enddo
!OCL XFILL
       do j = js, je
       do i = is, ie
          momz_t( 1:ks-1,i,j) = 0.0_rp
          momz_t(ke:ka  ,i,j) = 0.0_rp
       enddo
       enddo
       call comm_vars8( momz_t(:,:,:), i_comm_momz_t )
#ifdef HIST_TEND
       call file_history_in(momz_tp, 'MOMZ_t_phys', 'tendency of momentum z due to physics (w/ HIST_TEND)', 'kg/m2/s2', dim_type='ZHXY' )
       call file_history_in(damp_t,  'MOMZ_t_damp', 'tendency of momentum z due to damping (w/ HIST_TEND)', 'kg/m2/s2', dim_type='ZHXY' )
#endif

       !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
!OCL XFILL
       do j = js-1, je+1
       do i = is-1, ie+1
       do k = ks, ke
          diff(k,i,j) = momx(k,i,j) - damp_velx(k,i,j) * ( dens(k,i,j)+dens(k,i+1,j) ) * 0.5_rp
       enddo
       enddo
       enddo
!OCL XFILL
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,damp) &
#ifdef HIST_TEND
       !$omp shared(damp_t) &
#endif
       !$omp shared(JS,JE,IS,IE,KS,KE) &
       !$omp shared(DAMP_alpha_VELX,diff,BND_SMOOTHER_FACT,MOMX_tp,MOMX_t)
       do j = js, je
       do i = is, ie
       do k = ks, ke
          damp = - damp_alpha_velx(k,i,j) &
               * ( diff(k,i,j) & ! rayleigh damping
                 - ( diff(k,i-1,j) + diff(k,i+1,j) + diff(k,i,j-1) + diff(k,i,j+1) - diff(k,i,j)*4.0_rp ) &
                 * 0.125_rp * bnd_smoother_fact ) ! horizontal smoother
          momx_t(k,i,j) = momx_tp(k,i,j) & ! tendency from physical step
                        + damp
#ifdef HIST_TEND
          damp_t(k,i,j) = damp
#endif
       enddo
       enddo
       enddo
!OCL XFILL
       do j = js, je
       do i = is, ie
          momx_t(   1:ks-1,i,j) = 0.0_rp
          momx_t(ke+1:ka  ,i,j) = 0.0_rp
       enddo
       enddo
       call comm_vars8( momx_t(:,:,:), i_comm_momx_t )
#ifdef HIST_TEND
       call file_history_in(momx_tp, 'MOMX_t_phys', 'tendency of momentum x due to physics (w/ HIST_TEND)', 'kg/m2/s2', dim_type='ZXHY' )
       call file_history_in(damp_t,  'MOMX_t_damp', 'tendency of momentum x due to damping (w/ HIST_TEND)', 'kg/m2/s2', dim_type='ZXHY' )
#endif

       !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
!OCL XFILL
       do j = js-1, je+1
       do i = is-1, ie+1
       do k = ks, ke
          diff(k,i,j) = momy(k,i,j) - damp_vely(k,i,j) * ( dens(k,i,j)+dens(k,i,j+1) ) * 0.5_rp
       enddo
       enddo
       enddo
!OCL XFILL
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,damp) &
#ifdef HIST_TEND
       !$omp shared(damp_t) &
#endif
       !$omp shared(JS,JE,IS,IE,KS,KE) &
       !$omp shared(DAMP_alpha_VELY,diff,BND_SMOOTHER_FACT,MOMY_tp,MOMY_t)
       do j = js, je
       do i = is, ie
       do k = ks, ke
          damp = - damp_alpha_vely(k,i,j) &
               * ( diff(k,i,j) & ! rayleigh damping
                 - ( diff(k,i-1,j) + diff(k,i+1,j) + diff(k,i,j-1) + diff(k,i,j+1) - diff(k,i,j)*4.0_rp ) &
                 * 0.125_rp * bnd_smoother_fact ) ! horizontal smoother
          momy_t(k,i,j) = momy_tp(k,i,j) & ! tendency from physical step
                        + damp
#ifdef HIST_TEND
          damp_t(k,i,j) = damp
#endif
       enddo
       enddo
       enddo
!OCL XFILL
       do j = js, je
       do i = is, ie
          momy_t(   1:ks-1,i,j) = 0.0_rp
          momy_t(ke+1:ka  ,i,j) = 0.0_rp
       enddo
       enddo
       call comm_vars8( momy_t(:,:,:), i_comm_momy_t )
#ifdef HIST_TEND
       call file_history_in(momy_tp, 'MOMY_t_phys', 'tendency of momentum y due to physics (w/ HIST_TEND)', 'kg/m2/s2', dim_type='ZXYH' )
       call file_history_in(damp_t,  'MOMY_t_damp', 'tendency of momentum y due to damping (w/ HIST_TEND)', 'kg/m2/s2', dim_type='ZXYH' )
#endif

       !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
!OCL XFILL
       do j = js-1, je+2
       do i = is-1, ie+2
       do k = ks, ke
          diff(k,i,j) = rhot(k,i,j) - damp_pott(k,i,j) * dens(k,i,j)
       enddo
       enddo
       enddo
!OCL XFILL
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,damp) &
#ifdef HIST_TEND
       !$omp shared(damp_t) &
#endif
       !$omp shared(JS,JE,IS,IE,KS,KE) &
       !$omp shared(DAMP_alpha_POTT,diff,BND_SMOOTHER_FACT,RHOT_t,RHOT_tp)
       do j = js, je
       do i = is, ie
       do k = ks, ke
          damp = - damp_alpha_pott(k,i,j) &
               * ( diff(k,i,j) & ! rayleigh damping
                 - ( diff(k,i-1,j) + diff(k,i+1,j) + diff(k,i,j-1) + diff(k,i,j+1) - diff(k,i,j)*4.0_rp ) &
                 * 0.125_rp * bnd_smoother_fact ) ! horizontal smoother
          rhot_t(k,i,j) = rhot_tp(k,i,j) & ! tendency from physical step
                        + damp
#ifdef HIST_TEND
          damp_t(k,i,j) = damp
#endif
       enddo
       enddo
       enddo
!OCL XFILL
       do j = js, je
       do i = is, ie
          rhot_t(   1:ks-1,i,j) = 0.0_rp
          rhot_t(ke+1:ka  ,i,j) = 0.0_rp
       enddo
       enddo
       call comm_vars8( rhot_t(:,:,:), i_comm_rhot_t )
#ifdef HIST_TEND
       call file_history_in(rhot_tp, 'RHOT_t_phys', 'tendency of rho*theta temperature due to physics (w/ HIST_TEND)', 'K kg/m3/s' )
       call file_history_in(damp_t,  'RHOT_t_damp', 'tendency of rho*theta temperature due to damping (w/ HIST_TEND)', 'K kg/m3/s' )
#endif

       call comm_wait ( dens_t(:,:,:), i_comm_dens_t, .false. )
       call comm_wait ( momz_t(:,:,:), i_comm_momz_t, .false. )
       call comm_wait ( momx_t(:,:,:), i_comm_momx_t, .false. )
       call comm_wait ( momy_t(:,:,:), i_comm_momy_t, .false. )
       call comm_wait ( rhot_t(:,:,:), i_comm_rhot_t, .false. )

       call prof_rapend  ("DYN_Large_Tendency", 2)

       call prof_rapstart("DYN_Large_Numfilter", 2)

       !-----< prepare numerical diffusion coefficient >-----

       if ( nd_coef == 0.0_rp ) then
!OCL XFILL
          num_diff(:,:,:,:,:) = 0.0_rp
       else
          call atmos_dyn_numfilter_coef( num_diff(:,:,:,:,:),                    & ! [OUT]
                                         dens, momz, momx, momy, rhot,           & ! [IN]
                                         cdz, cdx, cdy, fdz, fdx, fdy, dts,      & ! [IN]
                                         ref_dens, ref_pott,                     & ! [IN]
                                         nd_coef, nd_order, nd_sfc_fact, nd_use_rs ) ! [IN]
       endif

       call prof_rapend  ("DYN_Large_Numfilter", 2)

       if ( divdmp_coef > 0.0_rp ) then

          call atmos_dyn_divergence( ddiv, & ! (out)
               momz, momx, momy, & ! (in)
               gsqrt, j13g, j23g, j33g, mapf, & ! (in)
               rcdz, rcdx, rcdy, rfdz, fdz ) ! (in)

       else

!XFILL
          ddiv = 0.0_rp

       end if

       !------------------------------------------------------------------------
       ! Start short time integration
       !------------------------------------------------------------------------

       call prof_rapstart("DYN_Short_Tinteg", 2)

       call atmos_dyn_tinteg_short( dens, momz, momx, momy, rhot, prog,       & ! (inout)
                                    mflx_hi, tflx_hi,                         & ! (inout)
                                    dens_t, momz_t, momx_t, momy_t, rhot_t,   & ! (in)
                                    dpres0, rt2p, corioli,                    & ! (in)
                                    num_diff, wdamp_coef, divdmp_coef, ddiv,  & ! (in)
                                    flag_fct_momentum, flag_fct_t,            & ! (in)
                                    flag_fct_along_stream,                    & ! (in)
                                    cdz, fdz, fdx, fdy,                       & ! (in)
                                    rcdz, rcdx, rcdy, rfdz, rfdx, rfdy,       & ! (in)
                                    phi, gsqrt, j13g, j23g, j33g, mapf,       & ! (in)
                                    ref_dens, ref_rhot,                       & ! (in)
                                    bnd_w, bnd_e, bnd_s, bnd_n,               & ! (in)
                                    dts                                       ) ! (in)

       call prof_rapend  ("DYN_Short_Tinteg", 2)

#ifdef CHECK_MASS
       call check_mass( &
            dens, damp_dens, &
            mflx_hi, tflx_hi, &
            gsqrt, mapf, &
            rcdx, rcdy, &
            dts, step, &
            bnd_w, bnd_e, bnd_s, bnd_n )
#endif

       !$omp parallel do default(none) private(i,j,iv) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JSB,JEB,ISB,IEB,KS,KA,DENS,MOMZ,MOMX,MOMY,RHOT,VA,PROG,KE)
       do j  = jsb, jeb
       do i  = isb, ieb
          dens(   1:ks-1,i,j) = dens(ks,i,j)
          momz(   1:ks-1,i,j) = 0.0_rp
          momx(   1:ks-1,i,j) = momx(ks,i,j)
          momy(   1:ks-1,i,j) = momy(ks,i,j)
          rhot(   1:ks-1,i,j) = rhot(ks,i,j)
          do iv = 1, va
             prog(   1:ks-1,i,j,iv) = prog(ks,i,j,iv)
          end do
          dens(ke+1:ka,  i,j) = dens(ke,i,j)
          momz(ke+1:ka,  i,j) = 0.0_rp
          momx(ke+1:ka,  i,j) = momx(ke,i,j)
          momy(ke+1:ka,  i,j) = momy(ke,i,j)
          rhot(ke+1:ka,  i,j) = rhot(ke,i,j)
          do iv = 1, va
             prog(ke+1:ka,  i,j,iv) = prog(ke,i,j,iv)
          end do
       enddo
       enddo

       call comm_vars8( dens(:,:,:), i_comm_dens )
       call comm_vars8( momz(:,:,:), i_comm_momz )
       call comm_vars8( momx(:,:,:), i_comm_momx )
       call comm_vars8( momy(:,:,:), i_comm_momy )
       call comm_vars8( rhot(:,:,:), i_comm_rhot )
       do iv = 1, va
          call comm_vars8( prog(:,:,:,iv), i_comm_prog(iv) )
       end do
       call comm_wait ( dens(:,:,:), i_comm_dens, .false. )
       call comm_wait ( momz(:,:,:), i_comm_momz, .false. )
       call comm_wait ( momx(:,:,:), i_comm_momx, .false. )
       call comm_wait ( momy(:,:,:), i_comm_momy, .false. )
       call comm_wait ( rhot(:,:,:), i_comm_rhot, .false. )
       do iv = 1, va
          call comm_wait ( prog(:,:,:,iv), i_comm_prog(iv), .false. )
       end do

       if ( use_average ) then
          do j = jsb, jeb
          do i = isb, ieb
          do k = ks, ke
             dens_av(k,i,j) = dens_av(k,i,j) + dens(k,i,j)
             momz_av(k,i,j) = momz_av(k,i,j) + momz(k,i,j)
             momx_av(k,i,j) = momx_av(k,i,j) + momx(k,i,j)
             momy_av(k,i,j) = momy_av(k,i,j) + momy(k,i,j)
             rhot_av(k,i,j) = rhot_av(k,i,j) + rhot(k,i,j)
          end do
          end do
          end do
       endif

#ifndef DRY
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          mflx_av(k,i,j,:) = mflx_av(k,i,j,:) + mflx_hi(k,i,j,:)
       end do
       end do
       end do
#endif

    enddo ! dynamical steps

    if ( use_average ) then
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          dens_av(k,i,j) = dens_av(k,i,j) / nstep
       end do
       end do
       end do
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          momz_av(k,i,j) = momz_av(k,i,j) / nstep
       end do
       end do
       end do
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          momx_av(k,i,j) = momx_av(k,i,j) / nstep
       end do
       end do
       end do
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          momy_av(k,i,j) = momy_av(k,i,j) / nstep
       end do
       end do
       end do
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          rhot_av(k,i,j) = rhot_av(k,i,j) / nstep
       end do
       end do
       end do
    endif

#ifndef DRY
    !###########################################################################
    ! Update Tracers
    !###########################################################################

!OCL XFILL
    do n = 1, 3
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       mflx_hi(k,i,j,n) = mflx_av(k,i,j,n) / nstep
    end do
    end do
    end do
    end do

    call comm_vars8( mflx_hi(:,:,:,zdir), i_comm_mflx_z )
    call comm_vars8( mflx_hi(:,:,:,xdir), i_comm_mflx_x )
    call comm_vars8( mflx_hi(:,:,:,ydir), i_comm_mflx_y )
    call comm_wait ( mflx_hi(:,:,:,zdir), i_comm_mflx_z, .false. )
    call comm_wait ( mflx_hi(:,:,:,xdir), i_comm_mflx_x, .false. )
    call comm_wait ( mflx_hi(:,:,:,ydir), i_comm_mflx_y, .false. )

    if ( use_average ) then
!OCL XFILL
       qtrc_av(:,:,:,:) = 0.0_rp
    endif

    !------------------------------------------------------------------------
    ! Update each tracer
    !------------------------------------------------------------------------

#ifdef SDM
    do iq = 1, i_qv
#else
    do iq = 1, qa
#endif

       if ( tracer_advc(iq) ) then

          call prof_rapstart("DYN_Large_Numfilter", 2)

          if ( nd_coef_q == 0.0_rp ) then
!OCL XFILL
             num_diff_q(:,:,:,:) = 0.0_rp
          else
             call atmos_dyn_numfilter_coef_q( num_diff_q(:,:,:,:),                    & ! [OUT]
                                              dens00, qtrc(:,:,:,iq), iq==i_qv,       & ! [IN]
                                              cdz, cdx, cdy, dtl,                     & ! [IN]
                                              ref_qv, iq,                             & ! [IN]
                                              nd_coef_q, nd_order, nd_sfc_fact, nd_use_rs ) ! [IN]
          endif

          call prof_rapend  ("DYN_Large_Numfilter", 2)

          call prof_rapstart("DYN_Tracer_Tinteg", 2)

          if ( flag_tracer_split_tend ) then
             !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
             !$omp shared(KA,IA,JA,iq,QTRC,RHOQ_t,DENS,dtl)
             do j = 1, ja
             do i = 1, ia
             do k = 1, ka
                qtrc(k,i,j,iq) = qtrc(k,i,j,iq) &
                               + rhoq_t(k,i,j,iq) * dtl / dens(k,i,j)
             end do
             end do
             end do
             rhoq_tn => zero
          else
             rhoq_tn => rhoq_t(:,:,:,iq)
          end if

          call atmos_dyn_tinteg_tracer( &
               qtrc(:,:,:,iq), & ! (inout)
               qtrc0(:,:,:,iq), rhoq_tn, &! (in)
               dens00, dens, & ! (in)
               mflx_hi, num_diff_q, & ! (in)
               gsqrt, mapf(:,:,:,i_xy), & ! (in)
               cdz, rcdz, rcdx, rcdy,             & ! (in)
               bnd_w, bnd_e, bnd_s, bnd_n, & ! (in)
               dtl, & ! (in)
               llast .AND. flag_fct_tracer, & ! (in)
               flag_fct_along_stream         ) ! (in)

          call prof_rapend  ("DYN_Tracer_Tinteg", 2)

       else

          do j = js, je
          do i = is, ie
          do k = ks, ke
             qtrc(k,i,j,iq) = ( qtrc0(k,i,j,iq) * dens00(k,i,j) &
                              + rhoq_t(k,i,j,iq) * dtl          ) / dens(k,i,j)
          end do
          end do
          end do

       end if

       if ( use_average ) then
          do j = jsb, jeb
          do i = isb, ieb
          do k = ks, ke
             qtrc_av(k,i,j,iq) = qtrc(k,i,j,iq)
          end do
          end do
          end do
       endif

       call comm_vars8( qtrc(:,:,:,iq), i_comm_qtrc(iq) )

    enddo ! scalar quantities loop

    do iq = 1, qa
       call comm_wait ( qtrc(:,:,:,iq), i_comm_qtrc(iq), .false. )
    enddo
#endif

    return
  end subroutine atmos_dyn_tstep_large_fvm_heve

#ifdef CHECK_MASS
  subroutine check_mass( &
       DENS, DAMP_DENS, &
       mflx_hi, tflx_hi, &
       GSQRT, MAPF, &
       RCDX, RCDY, &
       dt, step, &
       BND_W, BND_E, BND_S, BND_N &
       )
    use mpi
    use scale_atmos_grid_cartesc_real, only: &
       vol => atmos_grid_cartesc_real_vol
    use scale_comm_cartesc, only: &
       comm_datatype, &
       comm_world
    use scale_file_history, only: &
       file_history_in
    implicit none
    real(RP), intent(in) :: DENS     (ka,ia,ja)
    real(RP), intent(in) :: DAMP_DENS(ka,ia,ja)
    real(RP), intent(in) :: mflx_hi  (ka,ia,ja,3)
    real(RP), intent(in) :: tflx_hi  (ka,ia,ja,3)
    real(RP), intent(in) :: GSQRT    (ka,ia,ja,7)
    real(RP), intent(in) :: MAPF     (   ia,ja,2,7)
    real(RP), intent(in) :: RCDX(ia)
    real(RP), intent(in) :: RCDY(ja)
    real(RP), intent(in) :: dt
    integer,  intent(in) :: step
    logical,  intent(in) :: BND_W
    logical,  intent(in) :: BND_E
    logical,  intent(in) :: BND_S
    logical,  intent(in) :: BND_N

    ! lateral boundary flux
    real(RP) :: mflx_lb_horizontal(ka)
    real(RP) :: allmflx_lb_horizontal(ka)
    real(RP) :: mflx_lb_total
    real(RP) :: mass_total
    real(RP) :: mass_total2
    real(RP) :: allmflx_lb_total
    real(RP) :: allmass_total
    real(RP) :: allmass_total2

    integer :: k, i, j
    integer :: ierr


    call file_history_in(mflx_hi(:,:,:,zdir), 'MFLXZ', 'momentum flux of z-direction (w/ CHECK_MASS)', 'kg/m2/s', dim_type='ZHXY' )
    call file_history_in(mflx_hi(:,:,:,xdir), 'MFLXX', 'momentum flux of x-direction (w/ CHECK_MASS)', 'kg/m2/s', dim_type='ZXHY' )
    call file_history_in(mflx_hi(:,:,:,ydir), 'MFLXY', 'momentum flux of y-direction (w/ CHECK_MASS)', 'kg/m2/s', dim_type='ZXYH' )

    call file_history_in(tflx_hi(:,:,:,zdir), 'TFLXZ', 'potential temperature flux of z-direction (w/ CHECK_MASS)', 'K*kg/m2/s', dim_type='ZHXY' )
    call file_history_in(tflx_hi(:,:,:,xdir), 'TFLXX', 'potential temperature flux of x-direction (w/ CHECK_MASS)', 'K*kg/m2/s', dim_type='ZXHY' )
    call file_history_in(tflx_hi(:,:,:,ydir), 'TFLXY', 'potential temperature flux of y-direction (w/ CHECK_MASS)', 'K*kg/m2/s', dim_type='ZXYH' )

    mflx_lb_total            = 0.0_rp
    mflx_lb_horizontal(:)    = 0.0_rp
    allmflx_lb_horizontal(:) = 0.0_rp

    if ( bnd_w ) then ! for western boundary
       i = is
       do j = js, je
       do k = ks, ke
          mflx_lb_total = mflx_lb_total + mflx_hi(k,i-1,j,xdir) * rcdx(i) * vol(k,i,j) &
               * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) * dt
          mflx_lb_horizontal(k) = mflx_lb_horizontal(k) + mflx_hi(k,i-1,j,xdir) * rcdx(i) * vol(k,i,j) &
               * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) * dt

       end do
       end do
    end if
    if ( bnd_e ) then ! for eastern boundary
       i = ie
       do j = js, je
       do k = ks, ke
          mflx_lb_total = mflx_lb_total - mflx_hi(k,i,j,xdir) * rcdx(i) * vol(k,i,j) &
               * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) * dt
          mflx_lb_horizontal(k) = mflx_lb_horizontal(k) - mflx_hi(k,i,j,xdir) * rcdx(i) * vol(k,i,j) &
               * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) * dt
       end do
       end do
    end if
    if ( bnd_s ) then ! for sourthern boundary
       j = js
       do i = is, ie
       do k = ks, ke
          mflx_lb_total = mflx_lb_total + mflx_hi(k,i,j-1,ydir) * rcdy(j) * vol(k,i,j) &
               * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) * dt
          mflx_lb_horizontal(k) = mflx_lb_horizontal(k) + mflx_hi(k,i,j-1,ydir) * rcdy(j) * vol(k,i,j) &
               * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) * dt
       end do
       end do
    end if
    if ( bnd_n ) then ! for northern boundary
       j = je
       do i = is, ie
       do k = ks, ke
          mflx_lb_total = mflx_lb_total - mflx_hi(k,i,j,ydir) * rcdy(j) * vol(k,i,j) &
               * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) * dt
          mflx_lb_horizontal(k) = mflx_lb_horizontal(k) - mflx_hi(k,i,j,ydir) * rcdy(j) * vol(k,i,j) &
               * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) * dt
       end do
       end do
    end if

    mass_total  = 0.0_rp
    mass_total2 = 0.0_rp

    ! check total mass in the inner region
    do j = js, je
    do i = is, ie
    do k = ks, ke
       mass_total  = mass_total  + dens(k,i,j) * vol(k,i,j)
       mass_total2 = mass_total2 + damp_dens(k,i,j) * vol(k,i,j)
    end do
    end do
    end do

    call mpi_allreduce( mflx_lb_total,        &
                        allmflx_lb_total,     &
                        1,                    &
                        comm_datatype,        &
                        mpi_sum,              &
                        comm_world,           &
                        ierr                  )

    log_info("check_mass",'(A,1x,I1,1x,ES24.17)') 'total mflx_lb:', step, allmflx_lb_total

    call mpi_allreduce( mass_total,           &
                        allmass_total,        &
                        1,                    &
                        comm_datatype,        &
                        mpi_sum,              &
                        comm_world,           &
                        ierr                  )

    log_info("check_mass",'(A,1x,I1,1x,ES24.17)') 'total mass   :', step, allmass_total

    call mpi_allreduce( mass_total2,          &
                        allmass_total2,       &
                        1,                    &
                        comm_datatype,        &
                        mpi_sum,              &
                        comm_world,           &
                        ierr                  )

    log_info("check_mass",'(A,1x,I1,1x,ES24.17)') 'total mass2  :', step, allmass_total2

    call mpi_allreduce( mflx_lb_horizontal(ks:ke),    &
                        allmflx_lb_horizontal(ks:ke), &
                        kmax,                         &
                        comm_datatype,                &
                        mpi_sum,                      &
                        comm_world,                   &
                        ierr                          )

    call file_history_in(allmflx_lb_horizontal(:), 'ALLMOM_lb_hz',                           &
                    'horizontally total momentum flux from lateral boundary (w/ CHECK_MASS)', 'kg/m2/s' )

    return
  end subroutine check_mass
#endif

end module scale_atmos_dyn_tstep_large_fvm_heve