scale_atmos_dyn_tstep_short_fvm_hevi.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"

#ifdef PROFILE_FAPP
#define PROFILE_START(name) call fapp_start(name, 1, 1)
#define PROFILE_STOP(name)  call fapp_stop (name, 1, 1)
#elif defined(PROFILE_FINEPA)
#define PROFILE_START(name) call start_collection(name)
#define PROFILE_STOP(name)  call stop_collection (name)
#else
#define PROFILE_START(name)
#define PROFILE_STOP(name)
#endif

module scale_atmos_dyn_tstep_short_fvm_hevi
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_index
  use scale_tracer
#if defined DEBUG || defined QUICKDEBUG
  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_short_fvm_hevi_regist
  public :: atmos_dyn_tstep_short_fvm_hevi_setup
  public :: atmos_dyn_tstep_short_fvm_hevi

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
#if 1
#define F2H(k,p,idx) (CDZ(k+p-1)*GSQRT(k+p-1,i,j,idx)/(CDZ(k)*GSQRT(k,i,j,idx)+CDZ(k+1)*GSQRT(k+1,i,j,idx)))
#else
# define F2H(k,p,idx) 0.5_RP
#endif

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,  private, parameter :: NB = 1
  integer,  private, parameter :: VA_FVM_HEVI = 0
  integer                      :: IFS_OFF
  integer                      :: JFS_OFF

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_tstep_short_fvm_hevi_regist( &
       ATMOS_DYN_TYPE, &
       VA_out,         &
       VAR_NAME,       &
       VAR_DESC,       &
       VAR_UNIT        )
    use scale_process, only: &
       prc_mpistop
    implicit none

    character(len=*),       intent(in)  :: atmos_dyn_type
    integer,                intent(out) :: va_out
    character(len=H_SHORT), intent(out) :: var_name(:)
    character(len=H_MID),   intent(out) :: var_desc(:)
    character(len=H_SHORT), intent(out) :: var_unit(:)
    !---------------------------------------------------------------------------

    if ( atmos_dyn_type /= 'FVM-HEVI' .AND. atmos_dyn_type /= 'HEVI' ) then
       write(*,*) 'xxx ATMOS_DYN_TYPE is not FVM-HEVI. Check!'
       call prc_mpistop
    endif

    va_out      = va_fvm_hevi
    var_name(:) = ""
    var_desc(:) = ""
    var_unit(:) = ""

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Register additional prognostic variables (HEVI)'
    if ( va_out < 1 ) then
       if( io_l ) write(io_fid_log,*) '*** => nothing.'
    endif

    return
  end subroutine atmos_dyn_tstep_short_fvm_hevi_regist

  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_tstep_short_fvm_hevi_setup
    implicit none
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*) '*** HEVI Setup'

    return
  end subroutine atmos_dyn_tstep_short_fvm_hevi_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_tstep_short_fvm_hevi( &
       DENS_RK, MOMZ_RK, MOMX_RK, MOMY_RK, RHOT_RK, &
       PROG_RK,                                     &
       mflx_hi, tflx_hi,                            &
       DENS0,   MOMZ0,   MOMX0,   MOMY0,   RHOT0,   &
       DENS,    MOMZ,    MOMX,    MOMY,    RHOT,    &
       DENS_t,  MOMZ_t,  MOMX_t,  MOMY_t,  RHOT_t,  &
       PROG0, PROG,                                 &
       DPRES0, RT2P, CORIOLI,                       &
       num_diff, wdamp_coef, divdmp_coef, DDIV,     &
       FLAG_FCT_MOMENTUM, FLAG_FCT_T,               &
       FLAG_FCT_ALONG_STREAM,                       &
       CDZ, FDZ, FDX, FDY,                          &
       RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY,          &
       PHI, GSQRT, J13G, J23G, J33G, MAPF,          &
       REF_dens, REF_rhot,                          &
       BND_W, BND_E, BND_S, BND_N,                  &
       dtrk, last                                   )
    use scale_grid_index
    use scale_const, only: &
#ifdef dry
       rdry   => const_rdry,  &
       cvdry  => const_cvdry, &
       cpdry  => const_cpdry, &
#endif
       eps    => const_eps,   &
       grav   => const_grav,  &
       p00    => const_pre00
    use scale_atmos_dyn_common, only: &
       atmos_dyn_fct
    use scale_atmos_dyn_fvm_flux, only: &
       atmos_dyn_fvm_flux_valuew_z, &
       atmos_dyn_fvm_fluxz_xyz,     &
       atmos_dyn_fvm_fluxx_xyz,     &
       atmos_dyn_fvm_fluxy_xyz,     &
       atmos_dyn_fvm_fluxz_xyw,     &
       atmos_dyn_fvm_fluxj13_xyw,   &
       atmos_dyn_fvm_fluxj23_xyw,   &
       atmos_dyn_fvm_fluxx_xyw,     &
       atmos_dyn_fvm_fluxy_xyw,     &
       atmos_dyn_fvm_fluxz_uyz,     &
       atmos_dyn_fvm_fluxj13_uyz,   &
       atmos_dyn_fvm_fluxj23_uyz,   &
       atmos_dyn_fvm_fluxx_uyz,     &
       atmos_dyn_fvm_fluxy_uyz,     &
       atmos_dyn_fvm_fluxz_xvz,     &
       atmos_dyn_fvm_fluxj13_xvz,   &
       atmos_dyn_fvm_fluxj23_xvz,   &
       atmos_dyn_fvm_fluxx_xvz,     &
       atmos_dyn_fvm_fluxy_xvz
    use scale_gridtrans, only: &
       i_xyz, &
       i_xyw, &
       i_uyw, &
       i_xvw, &
       i_uyz, &
       i_xvz, &
       i_uvz, &
       i_xy , &
       i_uy , &
       i_xv , &
       i_uv
#ifdef HIST_TEND
    use scale_history, only: &
       hist_in
#endif
    implicit none

    real(RP), intent(out) :: dens_rk(ka,ia,ja)   ! prognostic variables
    real(RP), intent(out) :: momz_rk(ka,ia,ja)   !
    real(RP), intent(out) :: momx_rk(ka,ia,ja)   !
    real(RP), intent(out) :: momy_rk(ka,ia,ja)   !
    real(RP), intent(out) :: rhot_rk(ka,ia,ja)   !

    real(RP), intent(out) :: prog_rk(ka,ia,ja,va)  !

    real(RP), intent(inout) :: mflx_hi(ka,ia,ja,3) ! rho * vel(x,y,z)
    real(RP), intent(out)   :: tflx_hi(ka,ia,ja,3) ! rho * theta * vel(x,y,z)

    real(RP), intent(in),target :: dens0(ka,ia,ja) ! prognostic variables
    real(RP), intent(in),target :: momz0(ka,ia,ja) ! at previous dynamical time step
    real(RP), intent(in),target :: momx0(ka,ia,ja) !
    real(RP), intent(in),target :: momy0(ka,ia,ja) !
    real(RP), intent(in),target :: rhot0(ka,ia,ja) !

    real(RP), intent(in) :: dens(ka,ia,ja)   ! prognostic variables
    real(RP), intent(in) :: momz(ka,ia,ja)   ! at previous RK step
    real(RP), intent(in) :: momx(ka,ia,ja)   !
    real(RP), intent(in) :: momy(ka,ia,ja)   !
    real(RP), intent(in) :: rhot(ka,ia,ja)   !

    real(RP), intent(in) :: dens_t(ka,ia,ja)
    real(RP), intent(in) :: momz_t(ka,ia,ja)
    real(RP), intent(in) :: momx_t(ka,ia,ja)
    real(RP), intent(in) :: momy_t(ka,ia,ja)
    real(RP), intent(in) :: rhot_t(ka,ia,ja)

    real(RP), intent(in)  :: prog0(ka,ia,ja,va)
    real(RP), intent(in)  :: prog (ka,ia,ja,va)

    real(RP), intent(in) :: dpres0(ka,ia,ja)
    real(RP), intent(in) :: rt2p(ka,ia,ja)
    real(RP), intent(in) :: corioli(1,ia,ja)
    real(RP), intent(in) :: num_diff(ka,ia,ja,5,3)
    real(RP), intent(in) :: wdamp_coef(ka)
    real(RP), intent(in) :: divdmp_coef
    real(RP), intent(in) :: ddiv(ka,ia,ja)

    logical,  intent(in) :: flag_fct_momentum
    logical,  intent(in) :: flag_fct_t
    logical,  intent(in) :: flag_fct_along_stream

    real(RP), intent(in) :: cdz(ka)
    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)  :: ref_dens(ka,ia,ja)
    real(RP), intent(in)  :: ref_rhot(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)  :: dtrk
    logical,  intent(in)  :: last


    ! diagnostic variables (work space)
    real(RP) :: pott(ka,ia,ja) ! potential temperature [K]
    real(RP) :: dpres(ka,ia,ja) ! pressure deviation from reference pressure

    real(RP) :: qflx_hi (ka,ia,ja,3)
    real(RP) :: qflx_j13(ka,ia,ja)
    real(RP) :: qflx_j23(ka,ia,ja)

    real(RP) :: advch ! horizontal advection
    real(RP) :: advcv ! vertical advection
    real(RP) :: wdmp  ! Raileight damping
    real(RP) :: div   ! divergence damping
    real(RP) :: pg    ! pressure gradient force
    real(RP) :: cf    ! colioris force
#ifdef HIST_TEND
    real(RP) :: advch_t(ka,ia,ja,5)
    real(RP) :: advcv_t(ka,ia,ja,5)
    real(RP) :: wdmp_t(ka,ia,ja)
    real(RP) :: ddiv_t(ka,ia,ja,3)
    real(RP) :: pg_t(ka,ia,ja,3)
    real(RP) :: cf_t(ka,ia,ja,2)
    logical  :: lhist
#endif

    ! for implicit solver
    real(RP) :: a(ka,ia,ja)
    real(RP) :: b
    real(RP) :: sr(ka,ia,ja)
    real(RP) :: sw(ka,ia,ja)
    real(RP) :: st(ka,ia,ja)
    real(RP) :: pt(ka,ia,ja)
    real(RP) :: c(kmax-1,ia,ja)

    real(RP) :: f1(ka,ia,ja)
    real(RP) :: f2(ka,ia,ja)
    real(RP) :: f3(ka,ia,ja)

    integer :: iis, iie, jjs, jje
    integer :: k, i, j
    integer :: iss, iee

#ifdef DEBUG
    pott(:,:,:) = undef
    dpres(:,:,:) = undef

    pt(:,:,:) = undef

    qflx_hi(:,:,:,:) = undef
    qflx_j13(:,:,:)   = undef
    qflx_j23(:,:,:)   = undef
#endif

#if defined DEBUG || defined QUICKDEBUG
    dens_rk(   1:ks-1,:,:)   = undef
    dens_rk(ke+1:ka  ,:,:)   = undef
    momz_rk(   1:ks-1,:,:)   = undef
    momz_rk(ke+1:ka  ,:,:)   = undef
    momx_rk(   1:ks-1,:,:)   = undef
    momx_rk(ke+1:ka  ,:,:)   = undef
    momy_rk(   1:ks-1,:,:)   = undef
    momy_rk(ke+1:ka  ,:,:)   = undef
    rhot_rk(   1:ks-1,:,:)   = undef
    rhot_rk(ke+1:ka  ,:,:)   = undef
    prog_rk(   1:ks-1,:,:,:) = undef
    prog_rk(ke+1:ka  ,:,:,:) = undef
#endif

#ifdef HIST_TEND
    lhist = last
#endif

#ifdef PROFILE_FIPP
    call fipp_start()
#endif

    ifs_off = 1
    jfs_off = 1
    if ( bnd_w ) ifs_off = 0
    if ( bnd_s ) jfs_off = 0


    do jjs = js, je, jblock
    jje = jjs+jblock-1
    do iis = is, ie, iblock
    iie = iis+iblock-1

       profile_start("hevi_pres")
       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DPRES0,RT2P,RHOT,REF_rhot,DPRES,DENS,PHI)
       do j = jjs, jje+1
       do i = iis, iie+1
          do k = ks, ke
#ifdef DEBUG
             call check( __line__, dpres0(k,i,j) )
             call check( __line__, rt2p(k,i,j) )
             call check( __line__, rhot(k,i,j) )
             call check( __line__, ref_rhot(k,i,j) )
#endif
             dpres(k,i,j) = dpres0(k,i,j) + rt2p(k,i,j) * ( rhot(k,i,j) - ref_rhot(k,i,j) )
          enddo
          dpres(ks-1,i,j) = dpres0(ks-1,i,j) - dens(ks,i,j) * ( phi(ks-1,i,j) - phi(ks+1,i,j) )
          dpres(ke+1,i,j) = dpres0(ke+1,i,j) - dens(ke,i,j) * ( phi(ke+1,i,j) - phi(ke-1,i,j) )
       enddo
       enddo
       profile_stop("hevi_pres")

       !##### continuity equation #####

       profile_start("hevi_mflx_z")
       ! at (x, y, w)
       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JJS,JJE,IIS,IIE,KS,KE,GSQRT,I_XYW,MOMX,MOMY,J13G,J23G,mflx_hi,MAPF,I_XY,num_diff)
       do j = jjs, jje
       do i = iis-1, iie
          mflx_hi(ks-1,i,j,zdir) = 0.0_rp
          do k = ks, ke-1
#ifdef DEBUG
             call check( __line__, momx(k+1,i  ,j) )
             call check( __line__, momx(k+1,i-1,j) )
             call check( __line__, momx(k  ,i  ,j) )
             call check( __line__, momx(k  ,i+1,j) )
             call check( __line__, momy(k+1,i,j) )
             call check( __line__, momy(k+1,i,j-1) )
             call check( __line__, momy(k  ,i,j) )
             call check( __line__, momy(k  ,i,j-1) )
#endif
             mflx_hi(k,i,j,zdir) = j13g(k,i,j,i_xyw) * 0.25_rp / mapf(i,j,2,i_xy) &
                                 * ( momx(k+1,i,j)+momx(k+1,i-1,j) &
                                   + momx(k  ,i,j)+momx(k  ,i-1,j) ) &
                                 + j23g(k,i,j,i_xyw) * 0.25_rp / mapf(i,j,1,i_xy) &
                                 * ( momy(k+1,i,j)+momy(k+1,i,j-1) &
                                   + momy(k  ,i,j)+momy(k  ,i,j-1) ) &
                                 + gsqrt(k,i,j,i_xyw) / ( mapf(i,j,1,i_xy)*mapf(i,j,2,i_xy) ) * num_diff(k,i,j,i_dens,zdir)
          enddo
          mflx_hi(ke  ,i,j,zdir) = 0.0_rp
       enddo
       enddo
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif
       profile_stop("hevi_mflx_z")

       profile_start("hevi_mflx_x")
       iss = max(iis-1,is-ifs_off)
       iee = min(iie,ieh)
       ! at (u, y, z)
       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JJS,JJE,iss,iee,KS,KE,GSQRT,I_UYZ,MOMX,num_diff,mflx_hi,MAPF,I_UY)
       do j = jjs, jje
       do i = iss, iee
       do k = ks, ke
#ifdef DEBUG
          call check( __line__, gsqrt(k,i,j,i_uyz) )
          call check( __line__, momx(k,i,j) )
          call check( __line__, num_diff(k,i,j,i_dens,xdir) )
#endif
          mflx_hi(k,i,j,xdir) = gsqrt(k,i,j,i_uyz) / mapf(i,j,2,i_uy) &
               * ( momx(k,i,j) + num_diff(k,i,j,i_dens,xdir) )
       enddo
       enddo
       enddo
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif
       profile_stop("hevi_mflx_x")

       ! at (x, v, z)
       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JJS,JS,JFS_OFF,JJE,JEH,IIS,IIE,KS,KE,GSQRT,I_XVZ,MOMY,num_diff,mflx_hi,MAPF,I_XV)
       do j = max(jjs-1,js-jfs_off), min(jje,jeh)
       do i = iis, iie
       do k = ks, ke
#ifdef DEBUG
          call check( __line__, gsqrt(k,i,j,i_xvz) )
          call check( __line__, momy(k,i,j) )
          call check( __line__, num_diff(k,i,j,i_dens,ydir) )
#endif
          mflx_hi(k,i,j,ydir) = gsqrt(k,i,j,i_xvz) / mapf(i,j,1,i_xv) &
               * ( momy(k,i,j) + num_diff(k,i,j,i_dens,ydir) )
       enddo
       enddo
       enddo
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif

       !--- update density
       profile_start("hevi_sr")
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,advcv,advch) &
#ifdef HIST_TEND
       !$omp shared(advcv_t,advch_t,lhist) &
#endif
       !$omp shared(JJS,JJE,IIS,IIE,KS,KE) &
       !$omp shared(DENS0,Sr,mflx_hi,DENS_t) &
       !$omp shared(RCDZ,RCDX,RCDY,MAPF,GSQRT,I_XY,I_XYZ)
       do j = jjs, jje
       do i = iis, iie
       do k = ks, ke
#ifdef DEBUG
          call check( __line__, dens0(k,i,j) )
          call check( __line__, mflx_hi(k  ,i  ,j  ,xdir) )
          call check( __line__, mflx_hi(k  ,i-1,j  ,xdir) )
          call check( __line__, mflx_hi(k  ,i  ,j  ,ydir) )
          call check( __line__, mflx_hi(k  ,i  ,j-1,ydir) )
          call check( __line__, dens_t(k,i,j) )
#endif
          advcv = -   ( mflx_hi(k,i,j,zdir)-mflx_hi(k-1,i  ,j,  zdir) ) * rcdz(k)
          advch = - ( ( mflx_hi(k,i,j,xdir)-mflx_hi(k  ,i-1,j,  xdir) ) * rcdx(i) &
                    + ( mflx_hi(k,i,j,ydir)-mflx_hi(k  ,i,  j-1,ydir) ) * rcdy(j) ) &
                * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)
          sr(k,i,j) =  ( advcv + advch ) / gsqrt(k,i,j,i_xyz) + dens_t(k,i,j)
#ifdef HIST_TEND
          if ( lhist ) then
             advcv_t(k,i,j,i_dens) = advcv / gsqrt(k,i,j,i_xyz)
             advch_t(k,i,j,i_dens) = advch / gsqrt(k,i,j,i_xyz)
          end if
#endif
       enddo
       enddo
       enddo
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif
       profile_stop("hevi_sr")

       !##### momentum equation (z) #####

       ! at (x, y, z)
       ! not that z-index is added by -1
       profile_start("hevi_momz_qflxhi_z")
       call atmos_dyn_fvm_fluxz_xyw( qflx_hi(:,:,:,zdir), & ! (out)
            momz, momz, dens, & ! (in)
            gsqrt(:,:,:,i_xyz), j33g, & ! (in)
            num_diff(:,:,:,i_momz,zdir), & ! (in)
            cdz, fdz, dtrk, &
            iis, iie, jjs, jje ) ! (in)
       profile_stop("hevi_momz_qflxhi_z")

       profile_start("hevi_momz_qflxj")
       call atmos_dyn_fvm_fluxj13_xyw( qflx_j13, & ! (out)
            momx, momz, dens, & ! (in)
            gsqrt(:,:,:,i_xyz), j13g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
            cdz, &
            iis, iie, jjs, jje ) ! (in)
       call atmos_dyn_fvm_fluxj23_xyw( qflx_j23, & ! (out)
            momy, momz, dens, & ! (in)
            gsqrt(:,:,:,i_xyz), j23g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
            cdz, &
            iis, iie, jjs, jje ) ! (in)
       profile_stop("hevi_momz_qflxj")

       ! at (u, y, w)
       profile_start("hevi_momz_qflxhi_x")
       call atmos_dyn_fvm_fluxx_xyw( qflx_hi(:,:,:,xdir), & ! (out)
            momx, momz, dens, & ! (in)
            gsqrt(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
            num_diff(:,:,:,i_momz,xdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)
       profile_stop("hevi_momz_qflxhi_x")

       ! at (x, v, w)
       profile_start("hevi_momz_qflxhi_y")
       call atmos_dyn_fvm_fluxy_xyw( qflx_hi(:,:,:,ydir), & ! (out)
            momy, momz, dens, & ! (in)
            gsqrt(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
            num_diff(:,:,:,i_momz,ydir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)
       profile_stop("hevi_momz_qflxhi_y")

       !--- update momentum(z)
       profile_start("hevi_sw")
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,advcv,advch,cf,wdmp,div) &
#ifdef HIST_TEND
       !$omp shared(lhist,advcv_t,advch_t,wdmp_t,ddiv_t) &
#endif
       !$omp shared(JJS,JJE,IIS,IIE,KS,KE) &
       !$omp shared(qflx_hi,qflx_J13,qflx_J23,DDIV,MOMZ0,MOMZ_t,Sw) &
       !$omp shared(RFDZ,RCDX,RCDY,FDZ,dtrk,wdamp_coef,divdmp_coef) &
       !$omp shared(MAPF,GSQRT,I_XY,I_XYW)
       do j = jjs, jje
       do i = iis, iie
       do k = ks, ke-1
#ifdef DEBUG
          call check( __line__, qflx_hi(k  ,i  ,j  ,zdir) )
          call check( __line__, qflx_hi(k-1,i  ,j  ,zdir) )
          call check( __line__, qflx_j13(k  ,i  ,j)        )
          call check( __line__, qflx_j13(k-1,i  ,j)        )
          call check( __line__, qflx_j23(k  ,i  ,j)        )
          call check( __line__, qflx_j23(k-1,i  ,j)        )
          call check( __line__, qflx_hi(k  ,i  ,j  ,xdir) )
          call check( __line__, qflx_hi(k  ,i-1,j  ,xdir) )
          call check( __line__, qflx_hi(k  ,i  ,j  ,ydir) )
          call check( __line__, qflx_hi(k  ,i  ,j-1,ydir) )
          call check( __line__, ddiv(k  ,i,j) )
          call check( __line__, ddiv(k+1,i,j) )
          call check( __line__, momz0(k,i,j) )
          call check( __line__, momz_t(k,i,j) )
#endif
          advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i  ,j  ,zdir) &
                    + qflx_j13(k,i,j)      - qflx_j13(k-1,i  ,j  )      &
                    + qflx_j23(k,i,j)      - qflx_j23(k-1,i  ,j  )      ) * rfdz(k)
          advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k,i-1,j  ,xdir) ) * rcdx(i) &
                    + ( qflx_hi(k,i,j,ydir) - qflx_hi(k,i  ,j-1,ydir) ) * rcdy(j) ) &
                * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)
          wdmp = - wdamp_coef(k) * momz0(k,i,j)
          div = divdmp_coef / dtrk * ( ddiv(k+1,i,j)-ddiv(k,i,j) ) * fdz(k) ! divergence damping
          sw(k,i,j) = ( advcv + advch ) / gsqrt(k,i,j,i_xyw) &
                    + wdmp + div + momz_t(k,i,j)
#ifdef HIST_TEND
          if ( lhist ) then
             advcv_t(k,i,j,i_momz) = advcv / gsqrt(k,i,j,i_xyw)
             advch_t(k,i,j,i_momz) = advch / gsqrt(k,i,j,i_xyw)
             wdmp_t(k,i,j) = wdmp
             ddiv_t(k,i,j,1) = div
          endif
#endif
       enddo
       enddo
       enddo
       profile_stop("hevi_sw")
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif


       !##### Thermodynamic Equation #####

       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JJS,JJE,IIS,IIE,KS,KE,JHALO,IHALO,RHOT,DENS,POTT)
       do j = jjs-jhalo, jje+jhalo
       do i = iis-ihalo, iie+ihalo
       do k = ks, ke
#ifdef DEBUG
          call check( __line__, rhot(k,i,j) )
          call check( __line__, dens(k,i,j) )
#endif
          pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
       enddo
       enddo
       enddo
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif

       ! at (x, y, w)
       call atmos_dyn_fvm_fluxz_xyz( tflx_hi(:,:,:,zdir), & ! (out)
            mflx_hi(:,:,:,zdir), pott, gsqrt(:,:,:,i_xyw), & ! (in)
            num_diff(:,:,:,i_rhot,zdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       ! at (u, y, z)
       call atmos_dyn_fvm_fluxx_xyz( tflx_hi(:,:,:,xdir), & ! (out)
            mflx_hi(:,:,:,xdir), pott, gsqrt(:,:,:,i_uyz), & ! (in)
            num_diff(:,:,:,i_rhot,xdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       ! at (x, v, z)
       call atmos_dyn_fvm_fluxy_xyz( tflx_hi(:,:,:,ydir), & ! (out)
            mflx_hi(:,:,:,ydir), pott, gsqrt(:,:,:,i_xvz), & ! (in)
            num_diff(:,:,:,i_rhot,ydir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)


       profile_start("hevi_st")
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,advcv,advch) &
#ifdef HIST_TEND
       !$omp shared(lhist,advcv_t,advch_t) &
#endif
       !$omp shared(JJS,JJE,IIS,IIE,KS,KE) &
       !$omp shared(tflx_hi,RHOT_t,St,RCDZ,RCDX,RCDY) &
       !$omp shared(MAPF,GSQRT,I_XY,I_XYZ)
       do j = jjs, jje
       do i = iis, iie
       do k = ks, ke
#ifdef DEBUG
          call check( __line__, tflx_hi(k  ,i  ,j  ,zdir) )
          call check( __line__, tflx_hi(k-1,i  ,j  ,zdir) )
          call check( __line__, tflx_hi(k  ,i  ,j  ,xdir) )
          call check( __line__, tflx_hi(k  ,i-1,j  ,xdir) )
          call check( __line__, tflx_hi(k  ,i  ,j  ,ydir) )
          call check( __line__, tflx_hi(k  ,i  ,j-1,ydir) )
          call check( __line__, rhot_t(k,i,j) )
#endif
          advcv = -   ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i  ,j  ,zdir) ) * rcdz(k)
          advch = - ( ( tflx_hi(k,i,j,xdir) - tflx_hi(k  ,i-1,j  ,xdir) ) * rcdx(i) &
                    + ( tflx_hi(k,i,j,ydir) - tflx_hi(k  ,i  ,j-1,ydir) ) * rcdy(j) ) &
                  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)
          st(k,i,j) = ( advcv + advch ) / gsqrt(k,i,j,i_xyz) + rhot_t(k,i,j)
#ifdef HIST_TEND
          if ( lhist ) then
             advcv_t(k,i,j,i_rhot) = advcv / gsqrt(k,i,j,i_xyz)
             advch_t(k,i,j,i_rhot) = advch / gsqrt(k,i,j,i_xyz)
          endif
#endif
       enddo
       enddo
       enddo
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif
       profile_stop("hevi_st")

       ! implicit solver

       profile_start("hevi_solver")

!OCL INDEPENDENT
!OCL PREFETCH_SEQUENTIAL(SOFT)
#ifndef __GFORTRAN__
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,B,pg,advcv) &
#ifdef HIST_TEND
       !$omp shared(lhist,pg_t,advcv_t) &
#endif
#ifdef DEBUG
       !$omp shared(RHOT) &
#endif
       !$omp shared(JJS,JJE,IIS,IIE,KS,KE) &
       !$omp shared(mflx_hi,MOMZ_RK,MOMZ0) &
       !$omp shared(DENS_RK,RHOT_RK,DENS0,RHOT0,DENS,MOMZ,PT,POTT,DPRES) &
       !$omp shared(GRAV,dtrk,A,C,F1,F2,F3,REF_dens,Sr,Sw,St,RT2P) &
       !$omp shared(ATMOS_DYN_FVM_flux_valueW_Z) &
       !$omp shared(MAPF,GSQRT,J33G,I_XY,I_XYZ,I_XYW,CDZ,RCDZ,RFDZ)
#else
       !$omp parallel do default(shared) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp private(B,pg,advcv)
#endif
       do j = jjs, jje
       do i = iis, iie

          call atmos_dyn_fvm_flux_valuew_z( pt(:,i,j), & ! (out)
               momz(:,i,j), pott(:,i,j), gsqrt(:,i,j,i_xyz), & ! (in)
               cdz )

          do k = ks, ke
             a(k,i,j) = dtrk**2 * j33g * rcdz(k) * rt2p(k,i,j) * j33g / gsqrt(k,i,j,i_xyz)
          enddo
          b = grav * dtrk**2 * j33g / ( cdz(ks+1) + cdz(ks) )
          f1(ks,i,j) =        - ( pt(ks+1,i,j) * rfdz(ks) *   a(ks+1,i,j)             + b ) / gsqrt(ks,i,j,i_xyw)
          f2(ks,i,j) = 1.0_rp + ( pt(ks  ,i,j) * rfdz(ks) * ( a(ks+1,i,j)+a(ks,i,j) )     ) / gsqrt(ks,i,j,i_xyw)
          do k = ks+1, ke-2
             b = grav * dtrk**2 * j33g / ( cdz(k+1) + cdz(k) )
             f1(k,i,j) =        - ( pt(k+1,i,j) * rfdz(k) *   a(k+1,i,j)            + b ) / gsqrt(k,i,j,i_xyw)
             f2(k,i,j) = 1.0_rp + ( pt(k  ,i,j) * rfdz(k) * ( a(k+1,i,j)+a(k,i,j) )     ) / gsqrt(k,i,j,i_xyw)
             f3(k,i,j) =        - ( pt(k-1,i,j) * rfdz(k) *              a(k,i,j)   - b ) / gsqrt(k,i,j,i_xyw)
          enddo
          b = grav * dtrk**2 * j33g / ( cdz(ke) + cdz(ke-1) )
          f2(ke-1,i,j) = 1.0_rp + ( pt(ke-1,i,j) * rfdz(ke-1) * ( a(ke,i,j)+a(ke-1,i,j) )    ) / gsqrt(ke-1,i,j,i_xyw)
          f3(ke-1,i,j) =        - ( pt(ke-2,i,j) * rfdz(ke-1) *             a(ke-1,i,j)  - b ) / gsqrt(ke-1,i,j,i_xyw)
          do k = ks, ke-1
             pg = - ( dpres(k+1,i,j) + rt2p(k+1,i,j)*dtrk*st(k+1,i,j) &
                    - dpres(k  ,i,j) - rt2p(k  ,i,j)*dtrk*st(k  ,i,j) ) &
                    * rfdz(k) * j33g / gsqrt(k,i,j,i_xyw) &
                  - grav &
                    * ( f2h(k,1,i_xyz) * ( dens(k+1,i,j) - ref_dens(k+1,i,j) + sr(k+1,i,j) * dtrk ) &
                      + f2h(k,2,i_xyz) * ( dens(k  ,i,j) - ref_dens(k  ,i,j) + sr(k  ,i,j) * dtrk ) )
             c(k-ks+1,i,j) = momz(k,i,j) + dtrk * ( pg + sw(k,i,j) )
#ifdef HIST_TEND
             if ( lhist ) pg_t(k,i,j,1) = pg
#endif
          enddo

          call solve_direct( &
               c(:,i,j),        & ! (inout)
               f1(:,i,j), f2(:,i,j), f3(:,i,j) ) ! (in)

          do k = ks, ke-1
#ifdef DEBUG_HEVI2HEVE
          ! for debug (change to explicit integration)
             c(k-ks+1,i,j) = momz(k,i,j)
             mflx_hi(k,i,j,zdir) = mflx_hi(k,i,j,zdir) &
                                 + j33g * momz(k,i,j) / ( mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) )
             momz_rk(k,i,j) = momz0(k,i,j) &
                  + dtrk*( &
                  - j33g * ( dpres(k+1,i,j)-dpres(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw) &
                  - grav * ( f2h(k,2,i_xyz)*(dens(k,i,j)-ref_dens(k,i,j))+f2h(k,1,i_xyz)*(dens(k+1,i,j)-ref_dens(k+1,i,j)) ) &
                  + sw(k,i,j) )
#else
             ! z-momentum flux
             mflx_hi(k,i,j,zdir) = mflx_hi(k,i,j,zdir) &
                                 + j33g * c(k-ks+1,i,j) / ( mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) )
             ! z-momentum
             momz_rk(k,i,j) = momz0(k,i,j) &
                            + ( c(k-ks+1,i,j) - momz(k,i,j) )
#endif
          enddo
          momz_rk(ks-1,i,j) = 0.0_rp
          momz_rk(ke  ,i,j) = 0.0_rp

          ! density and rho*theta
          advcv = - c(1,i,j)            * j33g * rcdz(ks) / gsqrt(ks,i,j,i_xyz) ! C(0) = 0
          dens_rk(ks,i,j) = dens0(ks,i,j) + dtrk * ( advcv + sr(ks,i,j) )
#ifdef HIST_TEND
          if ( lhist ) advcv_t(ks,i,j,i_dens) = advcv
#endif
          advcv = - c(1,i,j)*pt(ks,i,j) * j33g * rcdz(ks) / gsqrt(ks,i,j,i_xyz) ! C(0) = 0
          rhot_rk(ks,i,j) = rhot0(ks,i,j) + dtrk * ( advcv + st(ks,i,j) )
#ifdef HIST_TEND
          if ( lhist ) advcv_t(ks,i,j,i_rhot) = advcv
#endif
          do k = ks+1, ke-1
             advcv = - ( c(k-ks+1,i,j)           - c(k-ks,i,j) ) &
                   * j33g * rcdz(k) / gsqrt(k,i,j,i_xyz)
             dens_rk(k,i,j) = dens0(k,i,j) + dtrk * ( advcv + sr(k,i,j) )
#ifdef HIST_TEND
             if ( lhist ) advcv_t(k,i,j,i_dens) = advcv
#endif
             advcv = - ( c(k-ks+1,i,j)*pt(k,i,j) - c(k-ks,i,j)*pt(k-1,i,j) ) &
                   * j33g * rcdz(k) / gsqrt(k,i,j,i_xyz)
             rhot_rk(k,i,j) = rhot0(k,i,j) + dtrk * ( advcv + st(k,i,j) )
#ifdef HIST_TEND
             if ( lhist ) advcv_t(k,i,j,i_rhot) = advcv
#endif
          enddo
          advcv = c(ke-ks,i,j)                * j33g * rcdz(ke) / gsqrt(ke,i,j,i_xyz) ! C(KE-KS+1) = 0
          dens_rk(ke,i,j) = dens0(ke,i,j) + dtrk * ( advcv + sr(ke,i,j) )
#ifdef HIST_TEND
          if ( lhist ) advcv_t(ke,i,j,i_dens) = advcv
#endif
          advcv = c(ke-ks,i,j) * pt(ke-1,i,j) * j33g * rcdz(ke) / gsqrt(ke,i,j,i_xyz) ! C(KE-KS+1) = 0
          rhot_rk(ke,i,j) = rhot0(ke,i,j) + dtrk * ( advcv + st(ke,i,j) )
#ifdef HIST_TEND
          if ( lhist ) advcv_t(ke,i,j,i_rhot) = advcv
#endif

#ifdef DEBUG
          call check_equation( &
               c(:,i,j), &
               dens(:,i,j), momz(:,i,j), rhot(:,i,j), dpres(:,i,j), &
               ref_dens(:,i,j), &
               sr(:,i,j), sw(:,i,j), st(:,i,j), &
               j33g, gsqrt(:,i,j,:), &
               rt2p(:,i,j), &
               dtrk, i, j )
#endif

       enddo
       enddo
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif

       profile_stop("hevi_solver")


       !##### momentum equation (x) #####

       profile_start("hevi_momx")
       ! at (u, y, w)
       call atmos_dyn_fvm_fluxz_uyz( qflx_hi(:,:,:,zdir), & ! (out)
            momz, momx, dens, & ! (in)
            gsqrt(:,:,:,i_uyw), j33g, & ! (in)
            num_diff(:,:,:,i_momx,zdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)
       call atmos_dyn_fvm_fluxj13_uyz( qflx_j13, & ! (out)
            momx, momx, dens, & ! (in)
            gsqrt(:,:,:,i_uyz), j13g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)
       call atmos_dyn_fvm_fluxj23_uyz( qflx_j23, & ! (out)
            momy, momx, dens, & ! (in)
            gsqrt(:,:,:,i_uyz), j23g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       ! at (x, y, z)
       ! note that x-index is added by -1
       call atmos_dyn_fvm_fluxx_uyz( qflx_hi(:,:,:,xdir), & ! (out)
            momx, momx, dens, & ! (in)
            gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
            num_diff(:,:,:,i_momx,xdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       ! at (u, v, z)
       call atmos_dyn_fvm_fluxy_uyz( qflx_hi(:,:,:,ydir), & ! (out)
            momy, momx, dens, & ! (in)
            gsqrt(:,:,:,i_uvz), mapf(:,:,1,i_uv), & ! (in)
            num_diff(:,:,:,i_momx,ydir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       !--- update momentum(x)
       iee = min(iie,ieh)
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,advch,advcv,pg,cf,div) &
#ifdef HIST_TEND
       !$omp shared(lhist,advch_t,advcv_t,pg_t,cf_t,ddiv_t) &
#endif
       !$omp shared(JJS,JJE,IIS,iee,KS,KE) &
       !$omp shared(MOMX_RK,DPRES,DENS,MOMX,MOMY,DDIV,MOMX0,MOMX_t) &
       !$omp shared(qflx_hi,qflx_J13,qflx_J23) &
       !$omp shared(RCDZ,RCDY,RFDX,CDZ,FDX) &
       !$omp shared(MAPF,GSQRT,J13G,I_XY,I_UY,I_UV,I_XYZ,I_UYW,I_UYZ) &
       !$omp shared(dtrk,CORIOLI,divdmp_coef)
       do j = jjs, jje
       do i = iis, iee
       do k = ks, ke
#ifdef DEBUG
          call check( __line__, qflx_hi(k  ,i  ,j  ,zdir) )
          call check( __line__, qflx_hi(k-1,i  ,j  ,zdir) )
          call check( __line__, qflx_hi(k  ,i  ,j  ,xdir) )
          call check( __line__, qflx_hi(k  ,i-1,j  ,xdir) )
          call check( __line__, qflx_hi(k  ,i  ,j  ,ydir) )
          call check( __line__, qflx_hi(k  ,i  ,j-1,ydir) )
          call check( __line__, dpres(k,i+1,j) )
          call check( __line__, dpres(k,i  ,j) )
          call check( __line__, corioli(1,i  ,j) )
          call check( __line__, corioli(1,i+1,j) )
          call check( __line__, momy(k,i  ,j  ) )
          call check( __line__, momy(k,i+1,j  ) )
          call check( __line__, momy(k,i  ,j-1) )
          call check( __line__, momy(k,i+1,j-1) )
          call check( __line__, ddiv(k,i+1,j) )
          call check( __line__, ddiv(k,i  ,j) )
          call check( __line__, momx0(k,i,j) )
#endif
          advcv = -   ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i  ,j  ,zdir) &
                      + qflx_j13(k,i,j)      - qflx_j13(k-1,i  ,j       ) &
                      + qflx_j23(k,i,j)      - qflx_j23(k-1,i  ,j       ) ) * rcdz(k)
          advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k  ,i-1,j  ,xdir) ) * rfdx(i) &
                    + ( qflx_hi(k,i,j,ydir) - qflx_hi(k  ,i  ,j-1,ydir) ) * rcdy(j) ) &
                * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy)
          pg = ( ( gsqrt(k,i+1,j,i_xyz) * dpres(k,i+1,j) & ! [x,y,z]
                 - gsqrt(k,i  ,j,i_xyz) * dpres(k,i  ,j) & ! [x,y,z]
                 ) * rfdx(i) &
               + ( j13g(k  ,i,j,i_uyw) &
                 * 0.5_rp * ( f2h(k,1,i_uyz) * ( dpres(k+1,i+1,j)+dpres(k+1,i,j) ) &
                            + f2h(k,2,i_uyz) * ( dpres(k  ,i+1,j)+dpres(k  ,i,j) ) ) & ! [x,y,z->u,y,w]
                 - j13g(k-1,i,j,i_uyw) &
                 * 0.5_rp * ( f2h(k,1,i_uyz) * ( dpres(k  ,i+1,j)+dpres(k  ,i,j) ) &
                            + f2h(k,2,i_uyz) * ( dpres(k-1,i+1,j)+dpres(k-1,i,j) ) ) & ! [x,y,z->u,y,w]
                 ) * rcdz(k) ) &
               * mapf(i,j,1,i_uy)
          cf = 0.125_rp * ( corioli(1,i+1,j  )+corioli(1,i,j  ) ) & ! [x,y,z->u,y,z]
             * ( momy(k,i+1,j  )+momy(k,i,j  ) &
               + momy(k,i+1,j-1)+momy(k,i,j-1) ) &  ! [x,v,z->u,y,z]
             + 0.25_rp * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) &
             * ( momy(k,i,j) + momy(k,i,j-1) + momy(k,i+1,j) + momy(k,i+1,j-1) ) &
             * ( ( momy(k,i,j) + momy(k,i,j-1) + momy(k,i+1,j) + momy(k,i+1,j-1) ) * 0.25_rp &
                 * ( 1.0_rp/mapf(i+1,j,2,i_xy) - 1.0_rp/mapf(i,j,2,i_xy) ) * rfdx(i) &
                   - momx(k,i,j) &
                   * ( 1.0_rp/mapf(i,j,1,i_uv) - 1.0_rp/mapf(i,j-1,1,i_uv) ) * rcdy(j) ) &
             * 2.0_rp / ( dens(k,i+1,j) + dens(k,i,j) ) ! metric term
          div = divdmp_coef / dtrk * ( ddiv(k,i+1,j)/mapf(i+1,j,2,i_xy) - ddiv(k,i,j)/mapf(i,j,1,i_xy) ) &
              * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) * fdx(i) ! divergence damping
          momx_rk(k,i,j) = momx0(k,i,j) &
               + dtrk * ( ( advcv + advch - pg ) / gsqrt(k,i,j,i_uyz) + cf + div + momx_t(k,i,j) )
#ifdef HIST_TEND
          if ( lhist ) then
             advcv_t(k,i,j,i_momx) = advcv / gsqrt(k,i,j,i_uyz)
             advch_t(k,i,j,i_momx) = advch / gsqrt(k,i,j,i_uyz)
             pg_t(k,i,j,2) = - pg / gsqrt(k,i,j,i_uyz)
             cf_t(k,i,j,1) = cf
             ddiv_t(k,i,j,2) = div
          endif
#endif
       enddo
       enddo
       enddo
       profile_stop("hevi_momx")
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif

       !##### momentum equation (y) #####
       profile_start("hevi_momy")
       ! at (x, v, w)
       call atmos_dyn_fvm_fluxz_xvz( qflx_hi(:,:,:,zdir), & ! (out)
            momz, momy, dens, & ! (in)
            gsqrt(:,:,:,i_xvw), j33g, & ! (in)
            num_diff(:,:,:,i_momy,zdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)
       call atmos_dyn_fvm_fluxj13_xvz( qflx_j13, & ! (out)
            momx, momy, dens, & ! (in)
            gsqrt(:,:,:,i_xvz), j13g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)
       call atmos_dyn_fvm_fluxj23_xvz( qflx_j23, & ! (out)
            momy, momy, dens, & ! (in)
            gsqrt(:,:,:,i_xvz), j23g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       ! at (u, v, z)
       call atmos_dyn_fvm_fluxx_xvz( qflx_hi(:,:,:,xdir), & ! (out)
            momx, momy, dens, & ! (in)
            gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_uv), & ! (in)
            num_diff(:,:,:,i_momy,xdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       ! at (x, y, z)
       ! note that y-index is added by -1
       call atmos_dyn_fvm_fluxy_xvz( qflx_hi(:,:,:,ydir), & ! (out)
            momy, momy, dens, & ! (in)
            gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
            num_diff(:,:,:,i_momy,ydir), & ! (in
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       !--- update momentum(y)
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,j,k,advch,advcv,pg,cf,div) &
#ifdef HIST_TEND
       !$omp shared(lhist,advch_t,advcv_t,pg_t,cf_t,ddiv_t) &
#endif
       !$omp shared(JJS,JJE,JEH,IIS,IIE,KS,KE) &
       !$omp shared(MOMY_RK,DPRES,DENS,MOMX,MOMY,DDIV,MOMY0,MOMY_t) &
       !$omp shared(qflx_hi,qflx_J13,qflx_J23) &
       !$omp shared(RCDZ,RCDX,RFDY,CDZ,FDY) &
       !$omp shared(MAPF,GSQRT,J23G,I_XY,I_XV,I_UV,I_XYZ,I_XVW,I_XVZ) &
       !$omp shared(dtrk,CORIOLI,divdmp_coef)
       do j = jjs, min(jje,jeh)
       do i = iis, iie
       do k = ks, ke
#ifdef DEBUG
          call check( __line__, qflx_hi(k  ,i  ,j  ,zdir) )
          call check( __line__, qflx_hi(k-1,i  ,j  ,zdir) )
          call check( __line__, qflx_hi(k  ,i  ,j  ,xdir) )
          call check( __line__, qflx_hi(k  ,i-1,j  ,xdir) )
          call check( __line__, qflx_hi(k  ,i  ,j  ,ydir) )
          call check( __line__, qflx_hi(k  ,i  ,j-1,ydir) )
          call check( __line__, dpres(k,i,j  ) )
          call check( __line__, dpres(k,i,j+1) )
          call check( __line__, corioli(1,i,j  ) )
          call check( __line__, corioli(1,i,j+1) )
          call check( __line__, momx(k,i  ,j  ) )
          call check( __line__, momx(k,i  ,j+1) )
          call check( __line__, momx(k,i-1,j  ) )
          call check( __line__, momx(k,i-1,j+1) )
          call check( __line__, ddiv(k,i,j+1) )
          call check( __line__, ddiv(k,i,j  ) )
          call check( __line__, momy_t(k,i,j) )
          call check( __line__, momy0(k,i,j) )
#endif
          advcv = -   ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i  ,j  ,zdir) &
                      + qflx_j13(k,i,j)      - qflx_j13(k-1,i  ,j  )      &
                      + qflx_j23(k,i,j)      - qflx_j23(k-1,i  ,j  )      ) * rcdz(k)
          advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k  ,i-1,j  ,xdir) ) * rcdx(i) &
                    + ( qflx_hi(k,i,j,ydir) - qflx_hi(k  ,i  ,j-1,ydir) ) * rfdy(j) ) &
                  * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv)
          pg = ( ( gsqrt(k,i,j+1,i_xyz) * dpres(k,i,j+1) & ! [x,y,z]
                 - gsqrt(k,i,j  ,i_xyz) * dpres(k,i,j  ) & ! [x,y,z]
                 ) * rfdy(j) &
               + ( j23g(k  ,i,j,i_xvw) &
                 * 0.5_rp * ( f2h(k  ,1,i_xvz) * ( dpres(k+1,i,j+1)+dpres(k+1,i,j) ) &
                            + f2h(k  ,2,i_xvz) * ( dpres(k  ,i,j+1)+dpres(k  ,i,j) ) ) & ! [x,y,z->x,v,w]
                 - j23g(k-1,i,j,i_xvw) &
                 * 0.5_rp * ( f2h(k-1,1,i_xvz) * ( dpres(k  ,i,j+1)+dpres(k  ,i,j) ) &
                            + f2h(k-1,2,i_xvz) * ( dpres(k-1,i,j+1)+dpres(k-1,i,j) ) ) & ! [x,y,z->x,v,w]
                 ) * rcdz(k) ) &
               * mapf(i,j,2,i_xv)
          cf = - 0.125_rp * ( corioli(1,i  ,j+1)+corioli(1,i  ,j) ) & ! [x,y,z->x,v,z]
                          * ( momx(k,i  ,j+1)+momx(k,i  ,j) &
                            + momx(k,i-1,j+1)+momx(k,i-1,j) ) & ! [u,y,z->x,v,z]
               - 0.25_rp * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) &
               * ( momx(k,i,j) + momx(k,i-1,j) + momx(k,i,j+1) + momx(k,i-1,j+1) ) &
               * ( momy(k,i,j) &
                 * ( 1.0_rp/mapf(i,j,2,i_uv) - 1.0_rp/mapf(i-1,j,2,i_uv) ) * rcdx(i) &
                 - 0.25_rp * ( momx(k,i,j)+momx(k,i-1,j)+momx(k,i,j+1)+momx(k,i-1,j+1) ) &
                 * ( 1.0_rp/mapf(i,j+1,1,i_xy) - 1.0_rp/mapf(i,j,1,i_xy) ) * rfdy(j) ) &
               * 2.0_rp / ( dens(k,i,j+1) + dens(k,i,j) ) ! metoric term
          div = divdmp_coef / dtrk * ( ddiv(k,i,j+1)/mapf(i,j+1,1,i_xy) - ddiv(k,i,j)/mapf(i,j,1,i_xy) ) &
              * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) * fdy(j) ! divergence damping
          momy_rk(k,i,j) = momy0(k,i,j) &
                         + dtrk * ( ( advcv + advch - pg ) / gsqrt(k,i,j,i_xvz) + cf + div + momy_t(k,i,j) )
#ifdef HIST_TEND
          if ( lhist ) then
             advcv_t(k,i,j,i_momy) = advcv / gsqrt(k,i,j,i_xvz)
             advch_t(k,i,j,i_momy) = advch / gsqrt(k,i,j,i_xvz)
             pg_t(k,i,j,3) = - pg / gsqrt(k,i,j,i_xvz)
             cf_t(k,i,j,2) = cf
             ddiv_t(k,i,j,3) = div
          endif
#endif
       enddo
       enddo
       enddo
       profile_stop("hevi_momy")
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif

    enddo
    enddo

#ifdef PROFILE_FIPP
       call fipp_stop()
#endif

#ifdef HIST_TEND
    if ( lhist ) then
       call hist_in(advcv_t(:,:,:,i_dens), 'DENS_t_advcv', 'tendency of density    (vert. advection)',    'kg/m3/s'   )
       call hist_in(advcv_t(:,:,:,i_momz), 'MOMZ_t_advcv', 'tendency of momentum z (vert. advection)',    'kg/m2/s2', zdim='half')
       call hist_in(advcv_t(:,:,:,i_momx), 'MOMX_t_advcv', 'tendency of momentum x (vert. advection)',    'kg/m2/s2', xdim='half')
       call hist_in(advcv_t(:,:,:,i_momy), 'MOMY_t_advcv', 'tendency of momentum y (vert. advection)',    'kg/m2/s2', ydim='half')
       call hist_in(advcv_t(:,:,:,i_rhot), 'RHOT_t_advcv', 'tendency of rho*theta  (vert. advection)',    'K kg/m3/s' )

       call hist_in(advch_t(:,:,:,i_dens), 'DENS_t_advch', 'tendency of density    (horiz. advection)',   'kg/m3/s'   )
       call hist_in(advch_t(:,:,:,i_momz), 'MOMZ_t_advch', 'tendency of momentum z (horiz. advection)',   'kg/m2/s2', zdim='half')
       call hist_in(advch_t(:,:,:,i_momx), 'MOMX_t_advch', 'tendency of momentum x (horiz. advection)',   'kg/m2/s2', xdim='half')
       call hist_in(advch_t(:,:,:,i_momy), 'MOMY_t_advch', 'tendency of momentum y (horiz. advection)',   'kg/m2/s2', ydim='half')
       call hist_in(advch_t(:,:,:,i_rhot), 'RHOT_t_advch', 'tendency of rho*theta  (horiz. advection)',   'K kg/m3/s' )

       call hist_in(pg_t(:,:,:,1),      'MOMZ_t_pg',    'tendency of momentum z (pressure gradient)',  'kg/m2/s2', zdim='half')
       call hist_in(pg_t(:,:,:,2),      'MOMX_t_pg',    'tendency of momentum x (pressure gradient)',  'kg/m2/s2', xdim='half')
       call hist_in(pg_t(:,:,:,3),      'MOMY_t_pg',    'tendency of momentum y (pressure gradient)',  'kg/m2/s2', ydim='half')

       call hist_in(wdmp_t(:,:,:),        'MOMZ_t_wdamp', 'tendency of momentum z (Rayleigh damping)',   'kg/m2/s2', zdim='half')

       call hist_in(ddiv_t(:,:,:,1),      'MOMZ_t_ddiv',  'tendency of momentum z (divergence damping)', 'kg/m2/s2', zdim='half')
       call hist_in(ddiv_t(:,:,:,2),      'MOMX_t_ddiv',  'tendency of momentum x (divergence damping)', 'kg/m2/s2', xdim='half')
       call hist_in(ddiv_t(:,:,:,3),      'MOMY_t_ddiv',  'tendency of momentum y (divergence damping)', 'kg/m2/s2', ydim='half')

       call hist_in(cf_t(:,:,:,1),      'MOMX_t_cf',    'tendency of momentum x (coliolis force)',     'kg/m2/s2', xdim='half')
       call hist_in(cf_t(:,:,:,2),      'MOMY_t_cf',    'tendency of momentum y (coliolis force)',     'kg/m2/s2', ydim='half')
    endif
#endif

  end subroutine atmos_dyn_tstep_short_fvm_hevi

!OCL SERIAL
  subroutine solve_direct( &
       C,         & ! (inout)
       f1, f2, f3 ) ! (in)

    use scale_process, only: &
       prc_mpistop
    implicit none
    real(RP), intent(inout) :: C(KMAX-1)
    real(RP), intent(in)    :: F1(KA)
    real(RP), intent(in)    :: F2(KA)
    real(RP), intent(in)    :: F3(KA)

    real(RP) :: e(KMAX-2)
    real(RP) :: f(KMAX-2)

    real(RP) :: rdenom

    integer :: k

    rdenom = 1.0_rp / f2(ks)
    e(1) = - f1(ks) * rdenom
    f(1) = c(1) * rdenom
    do k = 2, kmax-2
       rdenom = 1.0_rp / ( f2(k+ks-1) + f3(k+ks-1) * e(k-1) )
       e(k) = - f1(k+ks-1) * rdenom
       f(k) = ( c(k) - f3(k+ks-1) * f(k-1) ) * rdenom
    enddo

    ! C = \rho w
    c(kmax-1) = ( c(kmax-1) - f3(ke-1) * f(kmax-2) ) &
              / ( f2(ke-1) + f3(ke-1) * e(kmax-2) ) ! C(KMAX-1) = f(KMAX-1)
    do k = kmax-2, 1, -1
       c(k) = e(k) * c(k+1) + f(k)
    enddo

    return
  end subroutine solve_direct

#ifdef DEBUG
  subroutine check_equation( &
       VECT, &
       DENS, MOMZ, RHOT, DPRES, &
       REF_dens, &
       Sr, Sw, St, &
       J33G, G, &
       RT2P, &
       dt, i, j )
    use scale_const, only: &
         eps => const_eps, &
         grav => const_grav
    use scale_process, only: &
         prc_mpistop
    use scale_grid, only: &
         rcdz => grid_rcdz, &
         rfdz => grid_rfdz
    use scale_gridtrans, only: &
       i_xyz, &
       i_xyw
    implicit none
    real(RP), intent(in) :: VECT(KMAX-1)
    real(RP), intent(in) :: DENS(KA)
    real(RP), intent(in) :: MOMZ(KA)
    real(RP), intent(in) :: RHOT(KA)
    real(RP), intent(in) :: DPRES(KA)
    real(RP), intent(in) :: REF_dens(KA)
    real(RP), intent(in) :: Sr(KA)
    real(RP), intent(in) :: Sw(KA)
    real(RP), intent(in) :: St(KA)
    real(RP), intent(in) :: J33G
    real(RP), intent(in) :: G(KA,8)
    real(RP), intent(in) :: RT2P(KA)
    real(RP), intent(in) :: dt
    integer , intent(in) :: i
    integer , intent(in) :: j

    real(RP), parameter :: small = 1e-6_rp

    real(RP) :: MOMZ_N(KA)
    real(RP) :: DENS_N(KA)
    real(RP) :: RHOT_N(KA)
    real(RP) :: DPRES_N(KA)

    real(RP) :: POTT(KA)
    real(RP) :: PT(KA)

    real(RP) :: error, lhs, rhs
    integer :: k


    do k = ks, ke-1
       momz_n(k) = vect(k-ks+1)
    enddo
    momz_n(:ks-1) = 0.0_rp
    momz_n(ke:) = 0.0_rp

    ! density
    do k = ks+1, ke-1
       dens_n(k) = dens(k) &
            + dt * ( - j33g * ( momz_n(k) - momz_n(k-1) ) * rcdz(k) / g(k,i_xyz) + sr(k) )
    enddo
    dens_n(ks) = dens(ks) &
         + dt * ( - j33g * momz_n(ks) * rcdz(ks) / g(ks,i_xyz) + sr(ks) )
    dens_n(ke) = dens(ke) &
         + dt * ( j33g * momz_n(ke-1) * rcdz(ke) / g(ke,i_xyz) + sr(ke) )

    ! rho*theta
    do k = ks, ke
       pott(k) = rhot(k) / dens(k)
    enddo
    do k = ks+1, ke-2
       pt(k) = ( 7.0_rp * ( pott(k+1) + pott(k  ) ) &
                 -        ( pott(k+2) + pott(k-1) ) ) / 12.0_rp
    enddo
    pt(ks-1) = 0.0_rp
    pt(ks  ) = ( pott(ks+1) + pott(ks  ) ) * 0.5_rp
    pt(ke-1) = ( pott(ke  ) + pott(ke-1) ) * 0.5_rp
    pt(ke  ) = 0.0_rp
    do k = ks+1, ke-1
       rhot_n(k) = rhot(k) &
            + dt * ( - j33g * ( momz_n(k)*pt(k) - momz_n(k-1)*pt(k-1) ) * rcdz(k) / g(k,i_xyz) &
                     + st(k) )
    enddo
    rhot_n(ks) = rhot(ks) &
         + dt * ( - j33g * momz_n(ks)*pt(ks) * rcdz(ks) / g(ks,i_xyz) + st(ks) )
    rhot_n(ke) = rhot(ke) &
         + dt * ( j33g * momz_n(ke-1)*pt(ke-1) * rcdz(ke) / g(ke-1,i_xyz) + st(ke) )


    do k = ks, ke
       dpres_n(k) = dpres(k) + rt2p(k) * ( rhot_n(k) - rhot(k) )
    enddo

    do k = ks, ke
       lhs = ( dens_n(k) - dens(k) ) / dt
       rhs = - j33g * ( momz_n(k) - momz_n(k-1) ) * rcdz(k) / g(k,i_xyz) + sr(k)
       if ( abs(lhs) < small ) then
          error = rhs
       else
          error = ( lhs - rhs ) / lhs
       endif
       if ( abs(error) > small ) then
          write(*,*)"HEVI: DENS error", k, i, j, error, lhs, rhs
          write(*,*)eps
          call prc_mpistop
       endif
    enddo

    do k = ks, ke-1
       lhs = ( momz_n(k) - momz(k) ) / dt
       rhs = - j33g * ( dpres_n(k+1) - dpres_n(k) ) * rfdz(k) / g(k,i_xyw) &
             - grav * ( dens_n(k+1) - ref_dens(k+1) + dens_n(k) - ref_dens(k) ) * 0.5_rp &
             + sw(k)
       if ( abs(lhs) < small ) then
          error = rhs
       else
          error = ( lhs - rhs ) / lhs
       endif
       if ( abs(error) > small ) then
          write(*,*)"HEVI: MOMZ error", k, i, j, error, lhs, rhs
          write(*,*) momz_n(k), momz(k), dt
          write(*,*) - j33g * ( dpres(k+1) - dpres(k) ) * rfdz(k) / g(k,i_xyw) &
             - grav * ( dens(k+1) -ref_dens(k+1) + dens(k) -ref_dens(k) ) * 0.5_rp &
             + sw(k)
          call prc_mpistop
       endif
    enddo

    do k = ks, ke
       lhs = ( rhot_n(k) - rhot(k) ) / dt
       rhs = - j33g * ( momz_n(k)*pt(k) - momz_n(k-1)*pt(k-1) ) * rcdz(k) / g(k,i_xyz) + st(k)
       if ( abs(lhs) < small ) then
          error = rhs
       else
          error = ( lhs - rhs ) / lhs
       endif
       if ( abs(error) > small ) then
          write(*,*)"HEVI: RHOT error", k, i, j, error, lhs, rhs
          call prc_mpistop
       endif
    enddo

    return
  end subroutine check_equation
#endif

end module scale_atmos_dyn_tstep_short_fvm_hevi