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

#ifdef HIST_TEND
  real(RP), private, allocatable :: advch_t(:,:,:,:)
  real(RP), private, allocatable :: advcv_t(:,:,:,:)
  real(RP), private, allocatable :: ddiv_t (:,:,:,:)
  real(RP), private, allocatable :: pg_t   (:,:,:,:)
  real(RP), private, allocatable :: cf_t   (:,:,:,:)
#endif
  !-----------------------------------------------------------------------------
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( io_l ) write(io_fid_log,*) '*** HEVI Register'

    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(:) = ""

    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'
#ifdef HEVI_BICGSTAB
    if ( io_l ) write(io_fid_log,*) '*** USING Bi-CGSTAB'
#elif defined(HEVI_LAPACK)
    if ( io_l ) write(io_fid_log,*) '*** USING LAPACK'
#else
    if ( io_l ) write(io_fid_log,*) '*** USING DIRECT'
#endif

#ifdef HIST_TEND
    allocate( advch_t(ka,ia,ja,5) )
    allocate( advcv_t(ka,ia,ja,5) )
    allocate( ddiv_t(ka,ia,ja,3) )
    allocate( pg_t(ka,ia,ja,3) )
    allocate( cf_t(ka,ia,ja,2) )

    advch_t = 0.0_rp
    advcv_t = 0.0_rp
    ddiv_t  = 0.0_rp
    pg_t    = 0.0_rp
    cf_t    = 0.0_rp
#endif

    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, 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, dt                                     )
    use scale_grid_index
    use scale_const, only: &
#ifdef dry
       rdry   => const_rdry,  &
       cvdry  => const_cvdry, &
       cpdry  => const_cpdry, &
#endif
       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) :: 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
    real(RP), intent(in)  :: dt


    ! 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) :: 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) :: 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

#ifdef HIST_TEND
    lhist = dt .eq. dtrk
#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 private(i,j,k) OMP_SCHEDULE_ collapse(2)
       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 private(i,j,k) OMP_SCHEDULE_ collapse(2)
       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 private(i,j,k) OMP_SCHEDULE_ collapse(2)
       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 private(i,j,k) OMP_SCHEDULE_ collapse(2)
       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 private(i,j,k,advch) OMP_SCHEDULE_ collapse(2)
       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
          advch = - ( ( mflx_hi(k,i,j,zdir)-mflx_hi(k-1,i  ,j,  zdir) ) * rcdz(k) &
                    + ( 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) / gsqrt(k,i,j,i_xyz)
          sr(k,i,j) =  advch + dens_t(k,i,j)
#ifdef HIST_TEND
          if ( lhist ) advch_t(k,i,j,i_dens) = advch
#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 private(i,j,k,advcv,advch,cf,div) OMP_SCHEDULE_ collapse(2)
       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) ) * rfdz(k)
          advch = - ( ( qflx_j13(k,i,j)      - qflx_j13(k-1,i  ,j  ) &
                      + qflx_j23(k,i,j)      - qflx_j23(k-1,i  ,j  )    ) * rfdz(k) &
                    + ( 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)
          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) &
                    + 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)
             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 private(i,j,k) OMP_SCHEDULE_ collapse(2)
       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 private(i,j,k,advch) OMP_SCHEDULE_ collapse(2)
       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
          advch = - ( ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i  ,j  ,zdir) ) * rcdz(k) &
                    + ( 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) / gsqrt(k,i,j,i_xyz)
          st(k,i,j) = advch + rhot_t(k,i,j)
#ifdef HIST_TEND
          if ( lhist ) then
             advch_t(k,i,j,i_rhot) = advch
          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)
       !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
       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

#ifdef HEVI_BICGSTAB
          call solve_bicgstab( &
               c(:,i,j),         & ! (inout)
               f1(:,i,j), f2(:,i,j), f3(:,i,j) ) ! (in)
#elif defined(HEVI_LAPACK)
          call solve_lapack( &
               c(:,i,j),        & ! (inout)
#ifdef DEBUG
               i, j,      & ! (in)
#endif
               f1(:,i,j), f2(:,i,j), f3(:,i,j) ) ! (in)
#else
          call solve_direct( &
               c(:,i,j),        & ! (inout)
               f1(:,i,j), f2(:,i,j), f3(:,i,j) ) ! (in)
#endif

          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 private(i,j,k,advch,advcv,pg,cf,div) OMP_SCHEDULE_ collapse(2)
       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) ) * rcdz(k)
          advch = - ( ( qflx_j13(k,i,j)      - qflx_j13(k-1,i  ,j  ) &
                      + qflx_j23(k,i,j)      - qflx_j23(k-1,i  ,j  )      ) * rfdz(k) &
                    + ( 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) ) * rcdx(i) &
                   - momx(k,i,j) &
                   * ( 1.0_rp/mapf(i,j,1,i_uv) - 1.0_rp/mapf(i,j-1,1,i_uv) ) * rfdy(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 private(i,j,k,advch,advcv,pg,cf,div) OMP_SCHEDULE_ collapse(2)
       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) ) * rcdz(k)
          advch = - ( ( qflx_j13(k,i,j)      - qflx_j13(k-1,i  ,j  ) &
                      + qflx_j23(k,i,j)      - qflx_j23(k-1,i  ,j  )      ) * rcdz(k) &
                    + ( 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(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

#ifdef HEVI_BICGSTAB
!OCL SERIAL
  subroutine solve_bicgstab( &
       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) :: r0(kmax-1)

    real(RP) :: M(3,kmax-1)
    real(RP) :: p(kmax-1)
    real(RP) :: ap(kmax-1)
    real(RP) :: s(kmax-1)
    real(RP) :: as(kmax-1)
    real(RP) :: al, be, w

    real(RP), pointer :: r(:)
    real(RP), pointer :: rn(:)
    real(RP), pointer :: swap(:)
    real(RP), target :: v0(kmax-1)
    real(RP), target :: v1(kmax-1)
    real(RP) :: r0r
    real(RP) :: norm, error, epsilon

    integer :: k, iter

    epsilon = 0.1_rp**(rp-1)

    m(3,1) = f1(ks)
    m(2,1) = f2(ks)
    do k = ks+1, ke-2
       m(3,k-ks+1) = f1(k)
       m(2,k-ks+1) = f2(k)
       m(1,k-ks+1) = f3(k)
    enddo
    m(2,ke-ks) = f2(ke-1)
    m(1,ke-ks) = f3(ke-1)

    norm = 0.0_rp
    do k = 1, kmax-1
       norm = norm + c(k)**2
    enddo

    r  => v0
    rn => v1

    call mul_matrix( v1, m, c )

    do k = 1, kmax-1
       r(k) = c(k) - v1(k)
       r0(k) = r(k)
       p(k) = r(k)
    enddo

    r0r = r0(1) * r(1)
    do k = 2, kmax-1
       r0r = r0r + r0(k)*r(k)
    enddo
    do iter = 1, kmax-1
       error = 0.0_rp
       do k = 1, kmax-1
          error = error + r(k)**2
       enddo

!       if ( error < epsilon .OR. error / norm < epsilon ) then
       if ( error/norm < epsilon ) then
#ifdef DEBUG
!          write(*,*) "Bi-CGSTAB converged:", iter
#endif
          exit
       endif

       call mul_matrix( ap, m, p )
       al = r0(1) * ap(1)
       do k = 2, kmax-1
          al = al + r0(k)*ap(k)
       enddo
       al = r0r / al ! (r0,r) / (r0,Mp)
       s(:) = r(:) - al*ap(:)
       call mul_matrix( as, m, s )
       be = as(1) * s(1)  ! be is used as just work variable here
       w =  as(1) * as(1)
       do k = 2, kmax-1
          be = be + as(k)*s(k)
          w  = w  + as(k)*as(k)
       enddo
       w = be / w ! (as,s) / (as,as)

       c(:) = c(:) + al*p(:) + w*s(:)
       rn(:) = s(:) - w*as(:)
       be = al/w / r0r
       r0r = r0(1) * rn(1)
       do k = 2, kmax-1
          r0r = r0r + r0(k)*rn(k)
       enddo
       be = be * r0r ! al/w * (r0,rn)/(r0,r)
       p(:) = rn(:) + be * ( p(:) - w*ap(:) )

       swap => rn
       rn => r
       r => swap
    enddo

    if ( iter >= kmax-1 ) then
       write(*,*) 'xxx [atmos_dyn_hevi] Bi-CGSTAB'
       write(*,*) 'xxx not converged', error, norm
       call prc_mpistop
    endif

    return
  end subroutine solve_bicgstab

  subroutine mul_matrix(V, M, C)
    implicit none
    real(RP), intent(out) :: V(kmax-1)
    real(RP), intent(in)  :: M(3, kmax-1)
    real(RP), intent(in)  :: C(kmax-1)

    integer :: k

    ! k = 1
    v(1) = m(3,1)*c(2) + m(2,1)*c(1)
    do k = 2, kmax-2
       v(k) = m(3,k)*c(k+1) + m(2,k)*c(k) + m(1,k)*c(k-1)
    enddo
    ! k = KE-1
    v(kmax-1) = m(2,kmax-1)*c(kmax-1) + m(1,kmax-1)*c(kmax-2)

    return
  end subroutine mul_matrix

#elif defined(HEVI_LAPACK)
!OCL SERIAL
  subroutine solve_lapack( &
       C,   & ! (inout)
#ifdef DEBUG
       i, j, &
#endif
       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)
#ifdef DEBUG
    integer , intent(in)    :: i
    integer , intent(in)    :: j
#endif

    real(RP) :: M(nb*3+1,kmax-1)
    integer  :: IPIV(kmax-1)
    integer  :: INFO

    integer :: k

#ifdef DEBUG
    real(RP) :: M2(kmax-1,kmax-1)
    real(RP) :: C2(kmax-1)
    real(RP) :: sum
#endif


    ! band matrix
    do k = ks+1, ke-2
       ! k-1, +1
       m(nb+1,k-ks+1) = f1(k-1)
       ! k, 0
       m(nb+2,k-ks+1) = f2(k)
       ! k+1, -1
       m(nb+3,k-ks+1) = f3(k+1)
    enddo
    ! KS, 0
    m(nb+2,1    ) = f2(ks)
    ! KS+1, -1
    m(nb+3,1    ) = f3(ks+1)
    ! KE-2, +1
    m(nb+1,ke-ks) = f1(ke-2)
    ! KE-1, 0
    m(nb+2,ke-ks) = f2(ke-1)

#ifdef DEBUG
          m2(:,:) = 0.0_rp
          do k = 1, kmax-1
             if (k>1) m2(k-1,k) = m(nb+3,k-1)
             m2(k,k) = m(nb+2,k)
             if (k<kmax-1) m2(k+1,k) = m(nb+1,k+1)
             c2(k) = c(k)
          enddo
#endif
          if ( rp == dp ) then
             call dgbsv( kmax-1, nb, nb, 1, m, nb*3+1, ipiv, c, kmax-1, info)
          else
             call sgbsv( kmax-1, nb, nb, 1, m, nb*3+1, ipiv, c, kmax-1, info)
          endif
          ! C is (\rho w)^{n+1}
#ifdef DEBUG
          if ( info /= 0 ) then
             write(*,*) "DGBSV was failed", info
             call prc_mpistop
          endif

          do k = 1, kmax-1
             sum = 0.0_rp
             if (k>1) sum = sum + m2(k-1,k)*c(k-1)
             sum = sum + m2(k,k)*c(k)
             if (k<kmax-1) sum = sum + m2(k+1,k)*c(k+1)
             if ( abs(sum-c2(k)) > 1e-10_rp ) then
                write(*,*) "sum is different"
                write(*,*) k+2, i, j, sum, c2(k)
!                write(*,*) M2(k-1:k+1,k), c(k-1:k+1)
                call prc_mpistop
             endif
          enddo
#endif
  end subroutine solve_lapack
#else
!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

#endif

#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