scale_atmos_dyn_tstep_large_fvm_heve.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_atmos_dyn_tstep_large_fvm_heve
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_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 :: RHOQ_t(:,:,:,:)

  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_process, only: &
       prc_mpistop
    use scale_rm_process, 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, 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 )

    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, &
       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_FCT_MOMENTUM, FLAG_FCT_T, FLAG_FCT_TRACER,       &
       FLAG_FCT_ALONG_STREAM,                                &
       USE_AVERAGE,                                          &
       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, only: &
       comm_vars8, &
       comm_wait
    use scale_gridtrans, only: &
       i_xyz, &
       i_xyw, &
       i_uyz, &
       i_xvz, &
       i_xy,  &
       i_uy,  &
       i_xv
    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_atmos_boundary, only: &
       bnd_qa, &
       bnd_smoother_fact => atmos_boundary_smoother_fact
    use scale_history, only: &
#ifdef hist_tend
       hist_in, &
#endif
       hist_switch
    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

    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_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

    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) :: qflx_hi  (ka,ia,ja,3)  ! rho * vel(x,y,z) * phi @ (u,v,w)-face high order
    real(RP) :: qflx_lo  (ka,ia,ja,3)  ! rho * vel(x,y,z) * phi,  monotone flux
    real(RP) :: qflx_anti(ka,ia,ja,3)  ! anti-diffusive flux

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

    integer  :: iis, iie
    integer  :: jjs, jje
    integer  :: i, j, k, iq, step
    integer  :: iv

    real(RP) :: diff_coef
    !---------------------------------------------------------------------------

    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
    if( io_l ) write(io_fid_log,*)                         '*** Dynamics large time step'
    if( io_l ) write(io_fid_log,'(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

    qflx_hi(:,:,:,:) = undef
    qflx_lo(:,:,:,:) = 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, bnd_qa

       !$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) = qtrc(k,i,j,iq) - damp_qtrc(k,i,j,iq)
       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)
!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ks, ke
          damp = - damp_alpha_qtrc(k,i,j,iq) &
               * ( 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 hist_in(rhoq_tp(:,:,:,iq), trim(tracer_name(iq))//'_t_phys', &
                    'tendency of '//trim(tracer_name(iq))//' due to physics', 'kg/kg/s' )
       call hist_in(damp_t,            trim(tracer_name(iq))//'_t_damp', &
                    'tendency of '//trim(tracer_name(iq))//' due to damping', '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. )

    end do

    !$omp parallel do private(i,j,k,iq) OMP_SCHEDULE_ collapse(3)
!OCL XFILL
    do iq = bnd_qa+1, qa
       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 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 hist_switch( 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 hist_in(dens_tp, 'DENS_t_phys', 'tendency of dencity due to physics', 'kg/m3/s' )
       call hist_in(damp_t,  'DENS_t_damp', 'tendency of dencity due to damping', '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 hist_in(momz_tp, 'MOMZ_t_phys', 'tendency of momentum z due to physics',  'kg/m2/s2', zdim='half' )
       call hist_in(damp_t,  'MOMZ_t_damp', 'tendency of momentum z due to damping', 'kg/m2/s2', zdim='half' )
#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 hist_in(momx_tp, 'MOMX_t_phys', 'tendency of momentum x due to physics', 'kg/m2/s2', xdim='half' )
       call hist_in(damp_t,  'MOMX_t_damp', 'tendency of momentum x due to damping', 'kg/m2/s2', xdim='half' )
#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 hist_in(momy_tp, 'MOMY_t_phys', 'tendency of momentum y due to physics', 'kg/m2/s2', ydim='half' )
       call hist_in(damp_t,  'MOMY_t_damp', 'tendency of momentum y due to damping', 'kg/m2/s2', ydim='half' )
#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 hist_in(rhot_tp, 'RHOT_t_phys', 'tendency of rho*theta temperature due to physics', 'K kg/m3/s' )
       call hist_in(damp_t,  'RHOT_t_damp', 'tendency of rho*theta temperature due to damping', '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) = momz(ks,i,j)
          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) = momz(ke,i,j)
          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
          dens_av(:,:,:) = dens_av(:,:,:) + dens(:,:,:)
          momz_av(:,:,:) = momz_av(:,:,:) + momz(:,:,:)
          momx_av(:,:,:) = momx_av(:,:,:) + momx(:,:,:)
          momy_av(:,:,:) = momy_av(:,:,:) + momy(:,:,:)
          rhot_av(:,:,:) = rhot_av(:,:,:) + rhot(:,:,:)
       endif

#ifndef DRY
       mflx_av(:,:,:,:) = mflx_av(:,:,:,:) + mflx_hi(:,:,:,:)
#endif

    enddo ! dynamical steps

    if ( use_average ) then
       dens_av(:,:,:) = dens_av(:,:,:) / nstep
       momz_av(:,:,:) = momz_av(:,:,:) / nstep
       momx_av(:,:,:) = momx_av(:,:,:) / nstep
       momy_av(:,:,:) = momy_av(:,:,:) / nstep
       rhot_av(:,:,:) = rhot_av(:,:,:) / nstep
    endif

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

!OCL XFILL
    mflx_hi(:,:,:,:) = mflx_av(:,:,:,:) / nstep

    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),                 & ! [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)

          call atmos_dyn_tinteg_tracer( &
               qtrc(:,:,:,iq), & ! (inout)
               qtrc0(:,:,:,iq), rhoq_t(:,:,:,iq), &! (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
          qtrc_av(:,:,:,iq) = qtrc(:,:,:,iq)
       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_grid_real, only: &
       vol => real_vol
    use scale_comm, only: &
       comm_datatype, &
       comm_world
    use scale_history, only: &
       hist_in
    use scale_gridtrans, only: &
       i_xyz, &
       i_xy
    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 hist_in(mflx_hi(:,:,:,zdir), 'MFLXZ', 'momentum flux of z-direction', 'kg/m2/s', zdim='half' )
    call hist_in(mflx_hi(:,:,:,xdir), 'MFLXX', 'momentum flux of x-direction', 'kg/m2/s', xdim='half' )
    call hist_in(mflx_hi(:,:,:,ydir), 'MFLXY', 'momentum flux of y-direction', 'kg/m2/s', ydim='half' )

    call hist_in(tflx_hi(:,:,:,zdir), 'TFLXZ', 'potential temperature flux of z-direction', 'K*kg/m2/s', zdim='half' )
    call hist_in(tflx_hi(:,:,:,xdir), 'TFLXX', 'potential temperature flux of x-direction', 'K*kg/m2/s', xdim='half' )
    call hist_in(tflx_hi(:,:,:,ydir), 'TFLXY', 'potential temperature flux of y-direction', 'K*kg/m2/s', ydim='half' )

    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                  )

    if( io_l ) write(io_fid_log,'(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                  )

    if( io_l ) write(io_fid_log,'(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                  )

    if( io_l ) write(io_fid_log,'(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 hist_in(allmflx_lb_horizontal(:), 'ALLMOM_lb_hz',                           &
                    'horizontally total momentum flux from lateral boundary', 'kg/m2/s' )

    return
  end subroutine check_mass
#endif

end module scale_atmos_dyn_tstep_large_fvm_heve