scale_atmos_dyn.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_atmos_dyn
  !-----------------------------------------------------------------------------
  !
  !++ 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_setup
  public :: atmos_dyn

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public, allocatable :: corioli   (:,:)       ! coriolis term

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  ! for communication
  integer,  private              :: i_comm_dens = 1
  integer,  private              :: i_comm_momz = 2
  integer,  private              :: i_comm_momx = 3
  integer,  private              :: i_comm_momy = 4
  integer,  private              :: i_comm_rhot = 5
  integer,  private, allocatable :: i_comm_prog(:)
  integer,  private, allocatable :: i_comm_qtrc(:)

  logical,  private              :: bnd_w
  logical,  private              :: bnd_e
  logical,  private              :: bnd_s
  logical,  private              :: bnd_n

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

  real(RP), private, allocatable :: num_diff  (:,:,:,:,:)
  real(RP), private, allocatable :: num_diff_q(:,:,:,:)
  real(RP), private, allocatable :: wdamp_coef(:)         ! coefficient for Rayleigh damping of w

  logical,  private              :: dyn_none

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_setup( &
       DYN_Tinteg_Short_TYPE,        &
       DYN_Tinteg_Tracer_TYPE,       &
       DYN_Tinteg_Large_TYPE,        &
       DYN_Tstep_Tracer_TYPE,        &
       DYN_Tstep_Large_TYPE,         &
       DYN_FVM_FLUX_TYPE,            &
       DYN_FVM_FLUX_TYPE_TRACER,     &
       DENS, MOMZ, MOMX, MOMY, RHOT, &
       QTRC, PROG,                   &
       CDZ, CDX, CDY,                &
       FDZ, FDX, FDY,                &
       wdamp_tau,                    &
       wdamp_height,                 &
       FZ,                           &
       enable_coriolis,              &
       lat,                          &
       none                          )
    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
    use scale_atmos_dyn_common, only: &
       atmos_dyn_filter_setup, &
       atmos_dyn_wdamp_setup
    use scale_atmos_dyn_tinteg_short, only: &
       atmos_dyn_tinteg_short_setup
    use scale_atmos_dyn_tinteg_tracer, only: &
       atmos_dyn_tinteg_tracer_setup
    use scale_atmos_dyn_tinteg_large, only: &
       atmos_dyn_tinteg_large_setup
    use scale_atmos_dyn_tstep_short, only: &
       atmos_dyn_tstep_short_setup
    use scale_atmos_dyn_tstep_tracer, only: &
       atmos_dyn_tstep_tracer_setup
    use scale_atmos_dyn_tstep_large, only: &
       atmos_dyn_tstep_large_setup
    use scale_atmos_dyn_fvm_flux, only: &
       atmos_dyn_fvm_flux_setup
    implicit none

    character(len=*),  intent(in)    :: dyn_tinteg_short_type
    character(len=*),  intent(in)    :: dyn_tinteg_tracer_type
    character(len=*),  intent(in)    :: dyn_tinteg_large_type
    character(len=*),  intent(in)    :: dyn_tstep_tracer_type
    character(len=*),  intent(in)    :: dyn_tstep_large_type
    character(len=*),  intent(in)    :: dyn_fvm_flux_type
    character(len=*),  intent(in)    :: dyn_fvm_flux_type_tracer
    real(RP),          intent(inout) :: dens(ka,ia,ja)    ! MPI_RECV_INIT requires intent(inout)
    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(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)    :: wdamp_tau
    real(RP),          intent(in)    :: wdamp_height
    real(RP),          intent(in)    :: fz(0:ka)
    logical,           intent(in)    :: enable_coriolis
    real(RP),          intent(in)    :: lat(ia,ja)
    logical, optional, intent(in)    :: none

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

    dyn_none = .false.

    if ( present(none) ) then
       dyn_none = none
    endif

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

    if ( .NOT. dyn_none ) then
       bnd_w = .NOT. prc_has_w
       bnd_e = .NOT. prc_has_e
       bnd_s = .NOT. prc_has_s
       bnd_n = .NOT. prc_has_n

       allocate( corioli(ia,ja)        )
       allocate( mflx_hi(ka,ia,ja,3)   )
       allocate( num_diff(ka,ia,ja,5,3) )
       allocate( num_diff_q(ka,ia,ja,3)   )
       allocate( wdamp_coef(ka)           )
       mflx_hi(:,:,:,:) = undef

       call atmos_dyn_fvm_flux_setup     ( dyn_fvm_flux_type,            & ! [IN]
                                           dyn_fvm_flux_type_tracer      ) ! [IN]

       call atmos_dyn_tstep_short_setup

       call atmos_dyn_tstep_tracer_setup ( dyn_tstep_tracer_type         ) ! [IN]

       call atmos_dyn_tstep_large_setup  ( dyn_tstep_large_type,         & ! [IN]
                                           dens, momz, momx, momy, rhot, & ! [INOUT]
                                           qtrc, prog,                   & ! [INOUT]
                                           mflx_hi                       ) ! [INOUT]

       call atmos_dyn_tinteg_short_setup ( dyn_tinteg_short_type  ) ! [IN]

       call atmos_dyn_tinteg_tracer_setup( dyn_tinteg_tracer_type ) ! [IN]

       call atmos_dyn_tinteg_large_setup ( dyn_tinteg_large_type  ) ! [IN]

       ! numerical diffusion
       call atmos_dyn_filter_setup( num_diff, num_diff_q,        & ! [INOUT]
                                    cdz, cdx, cdy, fdz, fdx, fdy ) ! [IN]

       ! numerical diffusion
       call atmos_dyn_wdamp_setup( wdamp_coef(:),           & ! [INOUT]
                                   wdamp_tau, wdamp_height, & ! [IN]
                                   fz(:)                    ) ! [IN]

       ! coriolis parameter
       if ( enable_coriolis ) then
          corioli(:,:) = 2.0_rp * ohm * sin( lat(:,:) )
       else
          corioli(:,:) = 0.0_rp
       endif

    else

       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) )
       enddo

       do iq = 1, qa
          i_comm_qtrc(iq) = 5 + va + iq
          call comm_vars8_init( 'QTRC', qtrc(:,:,:,iq), i_comm_qtrc(iq) )
       enddo

    endif

    return
  end subroutine atmos_dyn_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_dyn( &
       DENS,    MOMZ,    MOMX,    MOMY,    RHOT,    QTRC,    &
       PROG,                                                 &
       DENS_av, MOMZ_av, MOMX_av, MOMY_av, RHOT_av, QTRC_av, &
       DENS_tp, MOMZ_tp, MOMX_tp, MOMY_tp, RHOT_tp, RHOQ_tp, &
       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,                 &
       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,    &
       divdmp_coef,                                          &
       FLAG_FCT_MOMENTUM, FLAG_FCT_T, FLAG_FCT_TRACER,       &
       FLAG_FCT_ALONG_STREAM,                                &
       USE_AVERAGE,                                          &
       DTSEC, DTSEC_DYN                                      )
    use scale_comm, only: &
       comm_vars8, &
       comm_wait
    use scale_atmos_boundary, only: &
       bnd_qa
    use scale_atmos_dyn_tinteg_large, only: &
       atmos_dyn_tinteg_large
    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(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)    :: 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)
    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)    :: 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)    :: dtsec
    real(DP), intent(in)    :: dtsec_dyn

    ! for time integartion
    real(RP) :: dens00 (ka,ia,ja)   ! saved density before update

    ! For tracer advection
    real(RP) :: tflx_hi(ka,ia,ja,3) ! rho * theta * vel(x,y,z) @ (u,v,w)-face high order
    real(RP) :: dt

    integer  :: i, j, k, iq
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*) '*** Atmos dynamics step'

    dt = real(dtsec, kind=rp)

    if ( dyn_none ) then

       call prof_rapstart("DYN_Tinteg", 2)

       do j = js, je
       do i = is, ie
       do k = ks, ke
          dens00(k,i,j) = dens(k,i,j) ! save previous value

          dens(k,i,j) = dens(k,i,j) + dens_tp(k,i,j) * dt
       enddo
       enddo
       enddo

       do j = js, je
       do i = is, ie
       do k = ks, ke
          momz(k,i,j) = momz(k,i,j) + momz_tp(k,i,j) * dt
       enddo
       enddo
       enddo

       do j = js, je
       do i = is, ie
       do k = ks, ke
          momx(k,i,j) = momx(k,i,j) + momx_tp(k,i,j) * dt
       enddo
       enddo
       enddo

       do j = js, je
       do i = is, ie
       do k = ks, ke
          momy(k,i,j) = momy(k,i,j) + momy_tp(k,i,j) * dt
       enddo
       enddo
       enddo

       do j = js, je
       do i = is, ie
       do k = ks, ke
          rhot(k,i,j) = rhot(k,i,j) + rhot_tp(k,i,j) * dt
       enddo
       enddo
       enddo

       do iq = 1, qa
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,iq) = max( 0.0_rp, &
                                qtrc(k,i,j,iq) * dens00(k,i,j) + rhoq_tp(k,i,j,iq) * dt &
                              ) / dens(k,i,j)
       enddo
       enddo
       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 iq = 1, qa
          call comm_vars8( qtrc(:,:,:,iq), i_comm_qtrc(iq) )
       enddo

       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 iq = 1, qa
          call comm_wait ( qtrc(:,:,:,iq), i_comm_qtrc(iq), .false. )
       enddo

       if ( use_average ) then
          dens_av(:,:,:)   = dens(:,:,:)
          momz_av(:,:,:)   = momz(:,:,:)
          momx_av(:,:,:)   = momx(:,:,:)
          momy_av(:,:,:)   = momy(:,:,:)
          rhot_av(:,:,:)   = rhot(:,:,:)
          qtrc_av(:,:,:,:) = qtrc(:,:,:,:)
       endif

       call prof_rapend("DYN_Tinteg", 2)

       return
    endif ! if DYN_NONE == .true.

    call prof_rapstart("DYN_Tinteg", 2)

    call atmos_dyn_tinteg_large( dens,    momz,    momx,    momy,    rhot,    qtrc,    & ! [INOUT]
                                 prog,                                                 & ! [INOUT]
                                 dens_av, momz_av, momx_av, momy_av, rhot_av, qtrc_av, & ! [INOUT]
                                 mflx_hi, tflx_hi,                                     & ! [OUT]
                                 num_diff, num_diff_q,                                 & ! [OUT;WORK]
                                 dens_tp, momz_tp, momx_tp, momy_tp, rhot_tp, rhoq_tp, & ! [IN]
                                 corioli,                                              & ! [IN]
                                 cdz, cdx, cdy, fdz, fdx, fdy,                         & ! [IN]
                                 rcdz, rcdx, rcdy, rfdz, rfdx, rfdy,                   & ! [IN]
                                 phi, gsqrt,                                           & ! [IN]
                                 j13g, j23g, j33g, mapf,                               & ! [IN]
                                 aq_r, aq_cv, aq_cp, aq_mass,                          & ! [IN]
                                 ref_dens, ref_pott, ref_qv, ref_pres,                 & ! [IN]
                                 bnd_w, bnd_e, bnd_s, bnd_n,                           & ! [IN]
                                 nd_coef, nd_coef_q, nd_order, nd_sfc_fact, nd_use_rs, & ! [IN]
                                 damp_dens,       damp_velz,       damp_velx,          & ! [IN]
                                 damp_vely,       damp_pott,       damp_qtrc,          & ! [IN]
                                 damp_alpha_dens, damp_alpha_velz, damp_alpha_velx,    & ! [IN]
                                 damp_alpha_vely, damp_alpha_pott, damp_alpha_qtrc,    & ! [IN]
                                 wdamp_coef,                                           & ! [IN]
                                 divdmp_coef,                                          & ! [IN]
                                 flag_fct_momentum, flag_fct_t, flag_fct_tracer,       & ! [IN]
                                 flag_fct_along_stream,                                & ! [IN]
                                 use_average,                                          & ! [IN]
                                 dtsec, dtsec_dyn                                      ) ! [IN]

    call prof_rapend  ("DYN_Tinteg", 2)

#ifdef CHECK_MASS
    call check_mass( dens, damp_dens,           & ! [IN]
                     mflx_hi, tflx_hi,          & ! [IN]
                     gsqrt, mapf,               & ! [IN]
                     rcdx, rcdy,                & ! [IN]
                     dt,                        & ! [IN]
                     bnd_w, bnd_e, bnd_s, bnd_n ) ! [IN]
#endif

    return
  end subroutine atmos_dyn

#ifdef CHECK_MASS
  !-----------------------------------------------------------------------------
  subroutine check_mass( &
       DENS, DAMP_DENS,           &
       mflx_hi, tflx_hi,          &
       GSQRT, MAPF,               &
       RCDX, RCDY,                &
       dt,                        &
       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
    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

       enddo
       enddo
    endif
    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
       enddo
       enddo
    endif
    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
       enddo
       enddo
    endif
    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
       enddo
       enddo
    endif

    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)
    enddo
    enddo
    enddo

    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,ES24.17)') 'total mflx_lb:', 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,ES24.17)') 'total mass   :', 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,ES24.17)') 'total mass2  :', 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