scale_timeint_rk.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Warning: This file was generated from common/scale_timeint_rk.F90.erb.
!          Do not edit this file.
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!> Module common / Runge-Kutta scheme 
!!
!! @par Description
!!      Driver module to provide various Runge-Kutta schemes. 
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_timeint_rk
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_prc
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type, procedures
  !
 
  !> Derived type to provide RK scheme   
  type, public :: timeint_rk
    real(RP), private :: dt               !< Timestep
 
    integer, private :: rk_scheme_id_ex   !< ID of explicit RK scheme
    integer, private :: rk_scheme_id_im   !< ID of implicit RK scheme
 
    integer, private, allocatable :: size_each_var(:)  !< Dimension information for each variable
    integer, public :: nstage                          !< Number of stage in RK scheme
    integer, public :: tend_buf_size                   !< Buffer size to store tendencies
 
    ! For Butcher representation (explicit part)
    real(rp), public, allocatable :: coef_a_ex(:,:)    !< Coefficients A in the Butcher table (explicit part)
    real(rp), public, allocatable :: coef_b_ex(:)      !< Coefficients b in the Butcher table (explicit part)
    real(rp), public, allocatable :: coef_c_ex(:)      !< Coefficients c in the Butcher table (explicit part)
 
    ! For Shu-Osher representation (explicit part)
    real(rp), public, allocatable :: coef_sig_ex(:,:)
    real(rp), public, allocatable :: coef_gam_ex(:,:)
 
    ! For Butcher representation (implicit part)
    real(rp), public, allocatable :: coef_a_im(:,:)    !< Coefficients A in the Butcher table (implicit part)
    real(rp), public, allocatable :: coef_b_im(:)      !< Coefficients b in the Butcher table (implicit part)
    real(rp), public, allocatable :: coef_c_im(:)      !< Coefficients c in the Butcher table (implicit part)
 
    integer, allocatable :: tend_buf_indmap(:)
    integer :: var_num
    real(rp), public, allocatable :: var_buf1d_ex(:,:,:)
    real(rp), public, allocatable :: tend_buf1d_ex(:,:,:)
    real(rp), public, allocatable :: tend_buf1d_im(:,:,:)
    real(rp), public, allocatable :: var0_1d(:,:)
    real(rp), private, allocatable :: vartmp_1d(:,:)    
    real(rp), public, allocatable :: var_buf2d_ex(:,:,:,:)
    real(rp), public, allocatable :: tend_buf2d_ex(:,:,:,:)
    real(rp), public, allocatable :: tend_buf2d_im(:,:,:,:)
    real(rp), public, allocatable :: var0_2d(:,:,:)
    real(rp), private, allocatable :: vartmp_2d(:,:,:)    
    real(rp), public, allocatable :: var_buf3d_ex(:,:,:,:,:)
    real(rp), public, allocatable :: tend_buf3d_ex(:,:,:,:,:)
    real(rp), public, allocatable :: tend_buf3d_im(:,:,:,:,:)
    real(rp), public, allocatable :: var0_3d(:,:,:,:)
    real(rp), private, allocatable :: vartmp_3d(:,:,:,:)    
 
    logical, private :: low_storage_flag
    logical, public :: imex_flag
    integer, public :: ndim
  contains
    procedure, public :: init => timeint_rk_init
    procedure, public :: final => timeint_rk_final
    procedure, public :: get_implicit_diagfac => timeint_rk_get_implicit_diagfac 
    procedure, public :: get_deltime => timeint_rk_get_deltime
    procedure, public :: advance1d => timeint_rk_advance1d
    procedure, public :: advance1d_varlist => timeint_rk_advance1d_varlist
    procedure, public :: advance_trcvar_1d => timeint_rk_advance_trcvar1d
    procedure, public :: storevar0_1d => timeint_rk_store_var0_1d
    procedure, public :: storeimplicit1d => timeint_rk_storeimpl1d
    procedure, public :: advance2d => timeint_rk_advance2d
    procedure, public :: advance2d_varlist => timeint_rk_advance2d_varlist
    procedure, public :: advance_trcvar_2d => timeint_rk_advance_trcvar2d
    procedure, public :: storevar0_2d => timeint_rk_store_var0_2d
    procedure, public :: storeimplicit2d => timeint_rk_storeimpl2d
    procedure, public :: advance3d => timeint_rk_advance3d
    procedure, public :: advance3d_varlist => timeint_rk_advance3d_varlist
    procedure, public :: advance_trcvar_3d => timeint_rk_advance_trcvar3d
    procedure, public :: storevar0_3d => timeint_rk_store_var0_3d
    procedure, public :: storeimplicit3d => timeint_rk_storeimpl3d
    generic, public :: advance => advance1d, advance2d, advance3d
    generic, public :: advance_varlist => advance1d_varlist, advance2d_varlist, advance3d_varlist
    generic, public :: advance_trcvar => advance_trcvar_1d, advance_trcvar_2d, advance_trcvar_3d
    generic, public :: storevar0 => storevar0_1d, storevar0_2d, storevar0_3d
    generic, public :: storeimplicit => storeimplicit1d, storeimplicit2d, storeimplicit3d
  end type timeint_rk
 
  !> Derived type to store pointer to variable data when using varlist in RK advance procedures
  type, public :: timeint_rk_var
    real(rp), pointer :: var1d(:)
    real(rp), pointer :: var2d(:,:)
    real(rp), pointer :: var3d(:,:,:)
  end type timeint_rk_var
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  !-------------------
 
contains
 
!----------------
 
!> Initialize a object to provide RK scheme 
!!
!! @param rk_scheme_name Name of RK scheme
!! @param dt Timestep
!! @param ndim Number of spatial Dimension
!! @param var_num Number of variables
!! @param size_each_var Array to store size of each dimension
  subroutine timeint_rk_init( this,                  &
    rk_scheme_name, dt, var_num, ndim, size_each_var )
 
    use scale_timeint_rk_butcher_tab, only: &
      timeint_rk_butcher_tab_get_info,      &
      timeint_rk_butcher_tab_get      
    implicit none
    class(timeint_rk), intent(inout) :: this
    character(*), intent(in) :: rk_scheme_name
    real(RP), intent(in) :: dt
    integer, intent(in) :: var_num
    integer, intent(in) :: ndim
    integer, intent(in) :: size_each_var(ndim)
    !----------------------------------------
 
    !$acc enter data create(this)
 
    this%dt = dt
    this%ndim = ndim
    this%var_num = var_num
    allocate( this%size_each_var(ndim) )
    this%size_each_var(:) = size_each_var(:)
    !$acc enter data copyin( this%size_each_var )
    !$acc enter data attach( this%size_each_var )
 
    call timeint_rk_butcher_tab_get_info( rk_scheme_name, & ! (in)
      this%nstage, this%tend_buf_size,                    & ! (out)
      this%low_storage_flag, this%imex_flag               ) ! (out)
    
    allocate ( this%coef_a_ex(this%nstage,this%nstage), this%coef_b_ex(this%nstage), this%coef_c_ex(this%nstage) )
    allocate ( this%coef_sig_ex(this%nstage+1,this%nstage), this%coef_gam_ex(this%nstage+1,this%nstage) )
!    if (this%imex_flag) then
      allocate ( this%coef_a_im(this%nstage,this%nstage), this%coef_b_im(this%nstage), this%coef_c_im(this%nstage) )
!    end if
 
    select case(this%ndim)
    case(1)
      allocate( this%tend_buf1D_ex(size_each_var(1), var_num, this%tend_buf_size) )
      !$acc enter data create( this%tend_buf1D_ex )
      if ( this%imex_flag ) then
        allocate( this%tend_buf1D_im(size_each_var(1), var_num, this%tend_buf_size) )
      else
        allocate( this%tend_buf1D_im(1,1,1) ) ! dummy allocation to avoid error in the case of non-IMEX scheme
      end if
      !$acc enter data create( this%tend_buf1D_im )
      allocate( this%var0_1D(size_each_var(1), var_num) )
      allocate( this%varTmp_1D(size_each_var(1), var_num) ) 
      !$acc enter data create( this%var0_1D, this%varTmp_1D )
    case(2)
      allocate( this%tend_buf2D_ex(size_each_var(1),size_each_var(2), var_num, this%tend_buf_size) )
      !$acc enter data create( this%tend_buf2D_ex )
      if ( this%imex_flag ) then
        allocate( this%tend_buf2D_im(size_each_var(1),size_each_var(2), var_num, this%tend_buf_size) )
      else
        allocate( this%tend_buf2D_im(1,1,1,1) ) ! dummy allocation to avoid error in the case of non-IMEX scheme
      end if
      !$acc enter data create( this%tend_buf2D_im )
      allocate( this%var0_2D(size_each_var(1),size_each_var(2), var_num) )
      allocate( this%varTmp_2D(size_each_var(1),size_each_var(2), var_num) ) 
      !$acc enter data create( this%var0_2D, this%varTmp_2D )
    case(3)
      allocate( this%tend_buf3D_ex(size_each_var(1),size_each_var(2),size_each_var(3), var_num, this%tend_buf_size) )
      !$acc enter data create( this%tend_buf3D_ex )
      if ( this%imex_flag ) then
        allocate( this%tend_buf3D_im(size_each_var(1),size_each_var(2),size_each_var(3), var_num, this%tend_buf_size) )
      else
        allocate( this%tend_buf3D_im(1,1,1,1,1) ) ! dummy allocation to avoid error in the case of non-IMEX scheme
      end if
      !$acc enter data create( this%tend_buf3D_im )
      allocate( this%var0_3D(size_each_var(1),size_each_var(2),size_each_var(3), var_num) )
      allocate( this%varTmp_3D(size_each_var(1),size_each_var(2),size_each_var(3), var_num) ) 
      !$acc enter data create( this%var0_3D, this%varTmp_3D )
    end select
    allocate( this%tend_buf_indmap(this%nstage) )
 
    call timeint_rk_butcher_tab_get( &
      rk_scheme_name, this%nstage, this%imex_flag,    & ! (in)
      this%coef_a_ex, this%coef_b_ex, this%coef_c_ex, & ! (out)
      this%coef_sig_ex, this%coef_gam_ex,             & ! (out)
      this%coef_a_im, this%coef_b_im, this%coef_c_im, & ! (out)
      this%tend_buf_indmap                            ) ! (out)
 
    !$acc enter data copyin( this%coef_a_ex, this%coef_b_ex, this%coef_c_ex, &
    !$acc this%coef_sig_ex, this%coef_gam_ex,                                &
    !$acc this%coef_a_im, this%coef_b_im, this%coef_c_im,                    &
    !$acc this%tend_buf_indmap )      
 
    !$acc enter data attach( this%coef_a_ex, this%coef_b_ex, this%coef_c_ex, &
    !$acc this%coef_sig_ex, this%coef_gam_ex,                                &
    !$acc this%coef_a_im, this%coef_b_im, this%coef_c_im,                    &
    !$acc this%tend_buf_indmap )      
 
    !$acc update device( this%dt, this%rk_scheme_id_ex, this%rk_scheme_id_im, this%nstage, this%tend_buf_size, this%var_num, this%low_storage_flag, this%imex_flag, this%ndim )
    return
  end subroutine timeint_rk_init
 
!> Finalize a object to provide RK scheme 
!!
  subroutine timeint_rk_final( this )
    implicit none
    class(timeint_rk), intent(inout) :: this
    !----------------------------------------
 
    !$acc exit data delete( this%coef_a_ex, this%coef_b_ex, this%coef_c_ex, &
    !$acc this%coef_sig_ex, this%coef_gam_ex )    
    deallocate( this%coef_a_ex, this%coef_b_ex, this%coef_c_ex )
    deallocate( this%coef_sig_ex, this%coef_gam_ex )
 
!    if (this%imex_flag) then
      !$acc exit data delete( this%coef_a_im, this%coef_b_im, this%coef_c_im )
      deallocate( this%coef_a_im, this%coef_b_im, this%coef_c_im )
!    end if
 
    !$acc exit data delete( this%size_each_var, this%tend_buf_indmap )
    deallocate( this%size_each_var, this%tend_buf_indmap )  
 
    select case(this%ndim)
    case(1)
      !$acc exit data delete( this%var0_1d, this%tend_buf1D_ex )
      deallocate( this%var0_1d )
      deallocate( this%tend_buf1D_ex )
!      if (this%imex_flag) then
        !$acc exit data delete( this%tend_buf1D_im )
        deallocate( this%tend_buf1D_im )
!      end if
      !$acc exit data delete( this%varTmp_1d )
      deallocate( this%varTmp_1d )
    case(2)
      !$acc exit data delete( this%var0_2d, this%tend_buf2D_ex )
      deallocate( this%var0_2d )
      deallocate( this%tend_buf2D_ex )
!      if (this%imex_flag) then
        !$acc exit data delete( this%tend_buf2D_im )
        deallocate( this%tend_buf2D_im )
!      end if
      !$acc exit data delete( this%varTmp_2d )
      deallocate( this%varTmp_2d )
    case(3)
      !$acc exit data delete( this%var0_3d, this%tend_buf3D_ex )
      deallocate( this%var0_3d )
      deallocate( this%tend_buf3D_ex )
!      if (this%imex_flag) then
        !$acc exit data delete( this%tend_buf3D_im )
        deallocate( this%tend_buf3D_im )
!      end if
      !$acc exit data delete( this%varTmp_3d )
      deallocate( this%varTmp_3d )
    end select
 
    !$acc exit data delete(this)
    return
  end subroutine timeint_rk_final
 
!> Get timestep
!!
  elemental function timeint_rk_get_deltime( this ) result(dt)
    implicit none
    class(timeint_rk), intent(in) :: this
    real(RP) :: dt
    !-----------------------------------------------------
    
    dt = this%dt
 
    return
  end function timeint_rk_get_deltime
 
!> Get a diagonal component of coefficient matrix at a stage for implicit RK scheme 
!!
!! @param nowstage Current stage
  elemental function timeint_rk_get_implicit_diagfac( this, nowstage ) result(fac)
    implicit none
    class(timeint_rk), intent(in) :: this
    integer, intent(in) :: nowstage
    real(RP) :: fac
    !-----------------------------------------------------
    
    fac = this%coef_a_im(nowstage,nowstage) * this%dt
 
    return
  end function timeint_rk_get_implicit_diagfac
  !----------------
 
 
!> Advance a variable at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_advance1d( this, nowstage, q, varID, is, ie )
    implicit none  
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
    !----------------------------------------  
    
    if (this%low_storage_flag) then
      call prof_rapstart( 'rk_advance_low_storage1D', 3) 
      call rk_advance_low_storage1d( this, nowstage, q, varid, is, ie , this%size_each_var(1), this%var_num, &
            this%var0_1D, this%varTmp_1d, this%tend_buf1D_ex )
      call prof_rapend( 'rk_advance_low_storage1D', 3) 
    else
      if ( this%imex_flag ) then
        call rk_advance_general1d( this, nowstage, q, varid, is, ie , this%size_each_var(1), this%var_num, &
            this%var0_1D, this%varTmp_1d, this%tend_buf1D_ex, this%tend_buf1D_im )
      else
        call rk_advance_general1d( this, nowstage, q, varid, is, ie , this%size_each_var(1), this%var_num, &
            this%var0_1D, this%varTmp_1d, this%tend_buf1D_ex )
      end if
    end if
 
    return
  end subroutine timeint_rk_advance1d
 
!> Advance variables at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param var_list Array of timeint_rk_var object which saves pointer to variable data
!! @param varIDs Indices of the targeting variables
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine timeint_rk_advance1d_varlist( this, nowstage, var_list, varIDs, is, ie )
    implicit none  
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: varIDs(:)
    type(timeint_rk_var), intent(inout) :: var_list(size(varIDs))
    integer, intent(in) :: is, ie 
 
    integer :: i
    integer :: n
    !----------------------------------------  
 
    n = size(varids)
 
    if (this%low_storage_flag) then
      call prof_rapstart( 'rk_advance_varlist_low_storage1D', 3) 
      i = 1
      do while(i <= n)
         if ( i+1 <= n ) then
          call rk_advance_low_storage1d_var2( this, nowstage, var_list(i)%var1D, var_list(i+1)%var1D, varids(i), varids(i+1), is, ie , this%size_each_var(1), this%var_num, &
                 this%var0_1D, this%varTmp_1d, this%tend_buf1D_ex )
          i = i + 2
        else
          call rk_advance_low_storage1d( this, nowstage, var_list(n)%var1D, varids(n), is, ie , this%size_each_var(1), this%var_num, &
                 this%var0_1D, this%varTmp_1d, this%tend_buf1D_ex )
          i = i + 1
        end if
      end do
      !$acc wait(1)
      call prof_rapend( 'rk_advance_varlist_low_storage1D', 3) 
    else
      do i=1, n
        call this%Advance1D( nowstage, var_list(i)%var1D, varids(i), is, ie  )
      end do
    end if
 
    return
  end subroutine timeint_rk_advance1d_varlist
 
!> Advance tracer variable data at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density
  subroutine timeint_rk_advance_trcvar1d( this, nowstage, q, varID, is, ie , &
      DDENS, DDENS0, DENS_hyd )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
    real(RP), intent(in) :: DDENS (:)
    real(RP), intent(in) :: DDENS0(:)
    real(RP), intent(in) :: DENS_hyd(:)
    !----------------------------------------    
 
    if ( this%imex_flag ) then
      if ( prc_ismaster ) write(*,*) "timeint_rk_advence_trcvar: IMEX is not supported. Check!"
      call prc_abort
    end if
 
    if (this%low_storage_flag) then
      call rk_advance_trcvar_low_storage1d( this, nowstage, q, varid, is, ie , this%size_each_var(1), this%var_num, &
        ddens, ddens0, dens_hyd, this%var0_1D, this%varTmp_1d, this%tend_buf1D_ex )
    else
      call rk_advance_trcvar_general1d( this, nowstage, q, varid, is, ie , this%size_each_var(1), this%var_num, &
        ddens, ddens0, dens_hyd, this%var0_1D, this%varTmp_1d, this%tend_buf1D_ex )
    end if
 
    return
  end subroutine timeint_rk_advance_trcvar1d
 
!> Store variable data at the time level n for the case of IMEX RK scheme
!!
!! @param q Array to store variable data at the time level n 
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_store_var0_1d( this, q, varID, is, ie  )
    implicit none
    class(timeint_rk), intent(inout) :: this
    real(RP), intent(inout) :: q(:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
 
    integer :: i
    !----------------------------------------    
 
    !$omp parallel private(i)
    !$acc parallel present( q, this%var0_1D )
    !$omp do
    !$acc loop
    do i=is, ie
      this%var0_1D(i,varid) = q(i)
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine timeint_rk_store_var0_1d
  
!> Store variable data after the implicit part of current stage in IMEX RK scheme
!!
!! @param nowstage Current stage
!! @param q Array to store variable data after the implicit part of current stage
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_storeimpl1d( this, nowstage, q, varID, is, ie )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
 
    !----------------------------------------    
 
    call rk_storeimpl_general1d( this, nowstage, q, varid, is, ie , this%size_each_var(1), this%var_num, &
      this%var0_1D, this%varTmp_1d, this%tend_buf1D_im )
 
    return
  end subroutine timeint_rk_storeimpl1d
 
!> Advance a variable at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_advance2d( this, nowstage, q, varID, is, ie ,js, je )
    implicit none  
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:,:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
    !----------------------------------------  
    
    if (this%low_storage_flag) then
      call prof_rapstart( 'rk_advance_low_storage2D', 3) 
      call rk_advance_low_storage2d( this, nowstage, q, varid, is, ie ,js, je , this%size_each_var(1),this%size_each_var(2), this%var_num, &
            this%var0_2D, this%varTmp_2d, this%tend_buf2D_ex )
      call prof_rapend( 'rk_advance_low_storage2D', 3) 
    else
      if ( this%imex_flag ) then
        call rk_advance_general2d( this, nowstage, q, varid, is, ie ,js, je , this%size_each_var(1),this%size_each_var(2), this%var_num, &
            this%var0_2D, this%varTmp_2d, this%tend_buf2D_ex, this%tend_buf2D_im )
      else
        call rk_advance_general2d( this, nowstage, q, varid, is, ie ,js, je , this%size_each_var(1),this%size_each_var(2), this%var_num, &
            this%var0_2D, this%varTmp_2d, this%tend_buf2D_ex )
      end if
    end if
 
    return
  end subroutine timeint_rk_advance2d
 
!> Advance variables at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param var_list Array of timeint_rk_var object which saves pointer to variable data
!! @param varIDs Indices of the targeting variables
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine timeint_rk_advance2d_varlist( this, nowstage, var_list, varIDs, is, ie ,js, je )
    implicit none  
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: varIDs(:)
    type(timeint_rk_var), intent(inout) :: var_list(size(varIDs))
    integer, intent(in) :: is, ie ,js, je 
 
    integer :: i
    integer :: n
    !----------------------------------------  
 
    n = size(varids)
 
    if (this%low_storage_flag) then
      call prof_rapstart( 'rk_advance_varlist_low_storage2D', 3) 
      i = 1
      do while(i <= n)
         if ( i+1 <= n ) then
          call rk_advance_low_storage2d_var2( this, nowstage, var_list(i)%var2D, var_list(i+1)%var2D, varids(i), varids(i+1), is, ie ,js, je , this%size_each_var(1),this%size_each_var(2), this%var_num, &
                 this%var0_2D, this%varTmp_2d, this%tend_buf2D_ex )
          i = i + 2
        else
          call rk_advance_low_storage2d( this, nowstage, var_list(n)%var2D, varids(n), is, ie ,js, je , this%size_each_var(1),this%size_each_var(2), this%var_num, &
                 this%var0_2D, this%varTmp_2d, this%tend_buf2D_ex )
          i = i + 1
        end if
      end do
      !$acc wait(1)
      call prof_rapend( 'rk_advance_varlist_low_storage2D', 3) 
    else
      do i=1, n
        call this%Advance2D( nowstage, var_list(i)%var2D, varids(i), is, ie ,js, je  )
      end do
    end if
 
    return
  end subroutine timeint_rk_advance2d_varlist
 
!> Advance tracer variable data at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density
  subroutine timeint_rk_advance_trcvar2d( this, nowstage, q, varID, is, ie ,js, je , &
      DDENS, DDENS0, DENS_hyd )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:,:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
    real(RP), intent(in) :: DDENS (:,:)
    real(RP), intent(in) :: DDENS0(:,:)
    real(RP), intent(in) :: DENS_hyd(:,:)
    !----------------------------------------    
 
    if ( this%imex_flag ) then
      if ( prc_ismaster ) write(*,*) "timeint_rk_advence_trcvar: IMEX is not supported. Check!"
      call prc_abort
    end if
 
    if (this%low_storage_flag) then
      call rk_advance_trcvar_low_storage2d( this, nowstage, q, varid, is, ie ,js, je , this%size_each_var(1),this%size_each_var(2), this%var_num, &
        ddens, ddens0, dens_hyd, this%var0_2D, this%varTmp_2d, this%tend_buf2D_ex )
    else
      call rk_advance_trcvar_general2d( this, nowstage, q, varid, is, ie ,js, je , this%size_each_var(1),this%size_each_var(2), this%var_num, &
        ddens, ddens0, dens_hyd, this%var0_2D, this%varTmp_2d, this%tend_buf2D_ex )
    end if
 
    return
  end subroutine timeint_rk_advance_trcvar2d
 
!> Store variable data at the time level n for the case of IMEX RK scheme
!!
!! @param q Array to store variable data at the time level n 
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_store_var0_2d( this, q, varID, is, ie ,js, je  )
    implicit none
    class(timeint_rk), intent(inout) :: this
    real(RP), intent(inout) :: q(:,:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
 
    integer :: i,j
    !----------------------------------------    
 
    !$omp parallel private(i,j)
    !$acc parallel present( q, this%var0_2D )
    !$omp do
    !$acc loop collapse(2)
    do j=js, je
    do i=is, ie
      this%var0_2D(i,j,varid) = q(i,j)
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine timeint_rk_store_var0_2d
  
!> Store variable data after the implicit part of current stage in IMEX RK scheme
!!
!! @param nowstage Current stage
!! @param q Array to store variable data after the implicit part of current stage
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_storeimpl2d( this, nowstage, q, varID, is, ie ,js, je )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:,:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
 
    !----------------------------------------    
 
    call rk_storeimpl_general2d( this, nowstage, q, varid, is, ie ,js, je , this%size_each_var(1),this%size_each_var(2), this%var_num, &
      this%var0_2D, this%varTmp_2d, this%tend_buf2D_im )
 
    return
  end subroutine timeint_rk_storeimpl2d
 
!> Advance a variable at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_advance3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke )
    implicit none  
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:,:,:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
    !----------------------------------------  
    
    if (this%low_storage_flag) then
      call prof_rapstart( 'rk_advance_low_storage3D', 3) 
      call rk_advance_low_storage3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke , this%size_each_var(1),this%size_each_var(2),this%size_each_var(3), this%var_num, &
            this%var0_3D, this%varTmp_3d, this%tend_buf3D_ex )
      call prof_rapend( 'rk_advance_low_storage3D', 3) 
    else
      if ( this%imex_flag ) then
        call rk_advance_general3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke , this%size_each_var(1),this%size_each_var(2),this%size_each_var(3), this%var_num, &
            this%var0_3D, this%varTmp_3d, this%tend_buf3D_ex, this%tend_buf3D_im )
      else
        call rk_advance_general3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke , this%size_each_var(1),this%size_each_var(2),this%size_each_var(3), this%var_num, &
            this%var0_3D, this%varTmp_3d, this%tend_buf3D_ex )
      end if
    end if
 
    return
  end subroutine timeint_rk_advance3d
 
!> Advance variables at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param var_list Array of timeint_rk_var object which saves pointer to variable data
!! @param varIDs Indices of the targeting variables
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine timeint_rk_advance3d_varlist( this, nowstage, var_list, varIDs, is, ie ,js, je ,ks, ke )
    implicit none  
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: varIDs(:)
    type(timeint_rk_var), intent(inout) :: var_list(size(varIDs))
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
 
    integer :: i
    integer :: n
    !----------------------------------------  
 
    n = size(varids)
 
    if (this%low_storage_flag) then
      call prof_rapstart( 'rk_advance_varlist_low_storage3D', 3) 
      i = 1
      do while(i <= n)
         if ( i+1 <= n ) then
          call rk_advance_low_storage3d_var2( this, nowstage, var_list(i)%var3D, var_list(i+1)%var3D, varids(i), varids(i+1), is, ie ,js, je ,ks, ke , this%size_each_var(1),this%size_each_var(2),this%size_each_var(3), this%var_num, &
                 this%var0_3D, this%varTmp_3d, this%tend_buf3D_ex )
          i = i + 2
        else
          call rk_advance_low_storage3d( this, nowstage, var_list(n)%var3D, varids(n), is, ie ,js, je ,ks, ke , this%size_each_var(1),this%size_each_var(2),this%size_each_var(3), this%var_num, &
                 this%var0_3D, this%varTmp_3d, this%tend_buf3D_ex )
          i = i + 1
        end if
      end do
      !$acc wait(1)
      call prof_rapend( 'rk_advance_varlist_low_storage3D', 3) 
    else
      do i=1, n
        call this%Advance3D( nowstage, var_list(i)%var3D, varids(i), is, ie ,js, je ,ks, ke  )
      end do
    end if
 
    return
  end subroutine timeint_rk_advance3d_varlist
 
!> Advance tracer variable data at the current stage with explicit RK part
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density
  subroutine timeint_rk_advance_trcvar3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke , &
      DDENS, DDENS0, DENS_hyd )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:,:,:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
    real(RP), intent(in) :: DDENS (:,:,:)
    real(RP), intent(in) :: DDENS0(:,:,:)
    real(RP), intent(in) :: DENS_hyd(:,:,:)
    !----------------------------------------    
 
    if ( this%imex_flag ) then
      if ( prc_ismaster ) write(*,*) "timeint_rk_advence_trcvar: IMEX is not supported. Check!"
      call prc_abort
    end if
 
    if (this%low_storage_flag) then
      call rk_advance_trcvar_low_storage3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke , this%size_each_var(1),this%size_each_var(2),this%size_each_var(3), this%var_num, &
        ddens, ddens0, dens_hyd, this%var0_3D, this%varTmp_3d, this%tend_buf3D_ex )
    else
      call rk_advance_trcvar_general3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke , this%size_each_var(1),this%size_each_var(2),this%size_each_var(3), this%var_num, &
        ddens, ddens0, dens_hyd, this%var0_3D, this%varTmp_3d, this%tend_buf3D_ex )
    end if
 
    return
  end subroutine timeint_rk_advance_trcvar3d
 
!> Store variable data at the time level n for the case of IMEX RK scheme
!!
!! @param q Array to store variable data at the time level n 
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_store_var0_3d( this, q, varID, is, ie ,js, je ,ks, ke  )
    implicit none
    class(timeint_rk), intent(inout) :: this
    real(RP), intent(inout) :: q(:,:,:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
 
    integer :: i,j,k
    !----------------------------------------    
 
    !$omp parallel private(i,j,k)
    !$acc parallel present( q, this%var0_3D )
    !$omp do collapse(2)
    !$acc loop collapse(3)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      this%var0_3D(i,j,k,varid) = q(i,j,k)
    end do
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine timeint_rk_store_var0_3d
  
!> Store variable data after the implicit part of current stage in IMEX RK scheme
!!
!! @param nowstage Current stage
!! @param q Array to store variable data after the implicit part of current stage
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
  subroutine timeint_rk_storeimpl3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    real(RP), intent(inout) :: q(:,:,:)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
 
    !----------------------------------------    
 
    call rk_storeimpl_general3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke , this%size_each_var(1),this%size_each_var(2),this%size_each_var(3), this%var_num, &
      this%var0_3D, this%varTmp_3d, this%tend_buf3D_im )
 
    return
  end subroutine timeint_rk_storeimpl3d
 
!-------------------
 
 
!> Advance variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_low_storage1d( this, nowstage, q, varID, is, ie , IA, var_num, &
      var0_1D, varTmp_1d, tend_buf1D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none      
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA
    real(RP), intent(inout) :: q(IA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_1d(IA,var_num)
    real(RP), intent(inout) :: varTmp_1d(IA,var_num)
    real(RP), intent(inout) :: tend_buf1D_ex(IA,var_num,this%tend_buf_size)
 
    integer :: i
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
 
    !----------------------------------------    
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i)
      !$omp do
      !$acc parallel present( q, varTmp_1d, tend_buf1D_ex ) async(1)
      !$acc loop
      do i=is, ie
        q(i) =  vartmp_1d(i,varid)                                             &
                + sig_ss * q(i) + gam_ss * tend_buf1d_ex(i,varid,1)
      end do
      !$acc end parallel
      !$omp end parallel
 
      return
    end if
      
    !$omp parallel private(i)
    !$acc parallel present( q, var0_1D, varTmp_1D, tend_buf1D_ex ) async(1)
    if (nowstage == 1) then
      !$omp do
      !$acc loop
      do i=is, ie
        var0_1d(i,varid)   = q(i)
        vartmp_1d(i,varid) = 0.0_rp
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do
      !$acc loop
      do i=is, ie
        vartmp_1d(i,varid) = vartmp_1d(i,varid)      &
            + sig_ns * q(i) + gam_ns * tend_buf1d_ex(i,varid,1)
      end do
    end if
 
    !$omp do
    !$acc loop
    do i=is, ie
      q(i) = one_minus_sig_ss * var0_1d(i,varid)                &
              + sig_ss * q(i) + gam_ss * tend_buf1d_ex(i,varid,1)
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_advance_low_storage1d
 
!> Advance variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q1 Array to store variable data
!! @param q2 Array to store variable data
!! @param varID1 Index of the targeting variable
!! @param varID2 Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_low_storage1d_var2( this, nowstage, q1, q2, varID1, varID2, is, ie , IA, var_num, &
      var0_1D, varTmp_1d, tend_buf1D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none      
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA
    real(RP), intent(inout) :: q1(IA)
    real(RP), intent(inout) :: q2(IA)
    integer, intent(in) :: varID1
    integer, intent(in) :: varID2
    integer, intent(in) :: is, ie 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_1d(IA,var_num)
    real(RP), intent(inout) :: varTmp_1d(IA,var_num)
    real(RP), intent(inout) :: tend_buf1D_ex(IA,var_num,this%tend_buf_size)
 
    integer :: i
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
 
    !----------------------------------------    
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i)
      !$omp do
      !$acc parallel present( q1, q2, varTmp_1d, tend_buf1D_ex ) async(1)
      !$acc loop
      do i=is, ie
        q1(i) =  vartmp_1d(i,varid1)                                             &
                + sig_ss * q1(i) + gam_ss * tend_buf1d_ex(i,varid1,1)
        q2(i) =  vartmp_1d(i,varid2)                                             &
                + sig_ss * q2(i) + gam_ss * tend_buf1d_ex(i,varid2,1)
      end do
      !$acc end parallel
      !$omp end parallel
 
      return
    end if
      
    !$omp parallel private(i)
    !$acc parallel present( q1, q2, var0_1D, varTmp_1D, tend_buf1D_ex ) async(1)
    if (nowstage == 1) then
      !$omp do
      !$acc loop
      do i=is, ie
        var0_1d(i,varid1)   = q1(i)
        var0_1d(i,varid2)   = q2(i)
 
        vartmp_1d(i,varid1) = 0.0_rp
        vartmp_1d(i,varid2) = 0.0_rp
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do
      !$acc loop
      do i=is, ie
        vartmp_1d(i,varid1) = vartmp_1d(i,varid1)      &
            + sig_ns * q1(i) + gam_ns * tend_buf1d_ex(i,varid1,1)
        vartmp_1d(i,varid2) = vartmp_1d(i,varid2)      &
            + sig_ns * q2(i) + gam_ns * tend_buf1d_ex(i,varid2,1)
      end do
    end if
 
    !$omp do
    !$acc loop
    do i=is, ie
      q1(i) = one_minus_sig_ss * var0_1d(i,varid1)                &
              + sig_ss * q1(i) + gam_ss * tend_buf1d_ex(i,varid1,1)
      q2(i) = one_minus_sig_ss * var0_1d(i,varid2)                &
              + sig_ss * q2(i) + gam_ss * tend_buf1d_ex(i,varid2,1)
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_advance_low_storage1d_var2
 
!> Advance tracer variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density
!OCL SERIAL
  subroutine rk_advance_trcvar_low_storage1d( this, nowstage, q, varID, is, ie , IA, var_num, &
    DDENS, DDENS0, DENS_hyd,                         &
    var0_1D, varTmp_1d, tend_buf1D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA
    real(RP), intent(inout) :: q(IA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
    integer, intent(in) :: var_num
    real(RP), intent(in) :: DDENS (IA)
    real(RP), intent(in) :: DDENS0(IA)
    real(RP), intent(in) :: DENS_hyd(IA)
    real(RP), intent(inout) :: var0_1d(IA,var_num)
    real(RP), intent(inout) :: varTmp_1d(IA,var_num)
    real(RP), intent(inout) :: tend_buf1D_ex(IA,var_num,this%tend_buf_size)
 
    integer :: i
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
    real(RP) :: c_ssm1
    real(RP) :: c_ss
 
    real(RP) :: dens_ssm1
    real(RP) :: dens_ss
   
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_trcvar_low_storage1D', 3) 
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)   
    c_ssm1 = this%coef_c_ex(nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i)
      !$acc parallel present( q, varTmp_1d, tend_buf1D_ex, DDENS, DDENS0, DENS_hyd )
      !$omp do
      !$acc loop
      do i=is, ie
        dens_ssm1 = dens_hyd(i) + ddens0(i) + c_ssm1 * ( ddens(i) - ddens0(i) )
        q(i) =  ( vartmp_1d(i,varid)                                              &
                + sig_ss * q(i) * dens_ssm1 + gam_ss * tend_buf1d_ex(i,varid,1) ) &
                / ( dens_hyd(i) + ddens(i) )
      end do
      !$acc end parallel
      !$omp end parallel
      call prof_rapend( 'rk_advance_trcvar_low_storage1D', 3)
 
      return
    end if
 
    c_ss   = this%coef_c_ex(nowstage+1) 
 
    !$omp parallel private(i,dens_ss,dens_ssm1)
    !$acc parallel present( q, varTmp_1d, tend_buf1D_ex, DDENS, DDENS0, DENS_hyd )
    if (nowstage == 1) then
      !$omp do
      !$acc loop
      do i=is, ie
        var0_1d(i,varid)   = q(i) * ( dens_hyd(i) + ddens0(i) )
        vartmp_1d(i,varid) = 0.0_rp
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do
      !$acc loop
      do i=is, ie
        dens_ssm1 = dens_hyd(i) + ddens0(i) + c_ssm1 * ( ddens(i) - ddens0(i) )
        vartmp_1d(i,varid) = vartmp_1d(i,varid)                      &
            + sig_ns * q(i) * dens_ssm1 + gam_ns * tend_buf1d_ex(i,varid,1)
      end do
    end if
 
    !$omp do
      !$acc loop
    do i=is, ie
      dens_ssm1 = dens_hyd(i) + ddens0(i) + c_ssm1 * ( ddens(i) - ddens0(i) )
      dens_ss   = dens_hyd(i) + ddens0(i) + c_ss   * ( ddens(i) - ddens0(i) )
 
      q(i) = ( one_minus_sig_ss * var0_1d(i,varid)                                &
              + sig_ss * q(i) * dens_ssm1 + gam_ss * tend_buf1d_ex(i,varid,1) )   &
              / dens_ss
    end do
    !$acc end parallel
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_trcvar_low_storage1D', 3)
 
    return
  end subroutine rk_advance_trcvar_low_storage1d
  
 
!> Advance variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_low_storage2d( this, nowstage, q, varID, is, ie ,js, je , IA,JA, var_num, &
      var0_2D, varTmp_2d, tend_buf2D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none      
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA
    real(RP), intent(inout) :: q(IA,JA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_2d(IA,JA,var_num)
    real(RP), intent(inout) :: varTmp_2d(IA,JA,var_num)
    real(RP), intent(inout) :: tend_buf2D_ex(IA,JA,var_num,this%tend_buf_size)
 
    integer :: i,j
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
 
    !----------------------------------------    
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i,j)
      !$omp do
      !$acc parallel present( q, varTmp_2d, tend_buf2D_ex ) async(1)
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        q(i,j) =  vartmp_2d(i,j,varid)                                             &
                + sig_ss * q(i,j) + gam_ss * tend_buf2d_ex(i,j,varid,1)
      end do
      end do
      !$acc end parallel
      !$omp end parallel
 
      return
    end if
      
    !$omp parallel private(i,j)
    !$acc parallel present( q, var0_2D, varTmp_2D, tend_buf2D_ex ) async(1)
    if (nowstage == 1) then
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        var0_2d(i,j,varid)   = q(i,j)
        vartmp_2d(i,j,varid) = 0.0_rp
      end do
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        vartmp_2d(i,j,varid) = vartmp_2d(i,j,varid)      &
            + sig_ns * q(i,j) + gam_ns * tend_buf2d_ex(i,j,varid,1)
      end do
      end do
    end if
 
    !$omp do
    !$acc loop collapse(2)
    do j=js, je
    do i=is, ie
      q(i,j) = one_minus_sig_ss * var0_2d(i,j,varid)                &
              + sig_ss * q(i,j) + gam_ss * tend_buf2d_ex(i,j,varid,1)
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_advance_low_storage2d
 
!> Advance variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q1 Array to store variable data
!! @param q2 Array to store variable data
!! @param varID1 Index of the targeting variable
!! @param varID2 Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_low_storage2d_var2( this, nowstage, q1, q2, varID1, varID2, is, ie ,js, je , IA,JA, var_num, &
      var0_2D, varTmp_2d, tend_buf2D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none      
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA
    real(RP), intent(inout) :: q1(IA,JA)
    real(RP), intent(inout) :: q2(IA,JA)
    integer, intent(in) :: varID1
    integer, intent(in) :: varID2
    integer, intent(in) :: is, ie ,js, je 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_2d(IA,JA,var_num)
    real(RP), intent(inout) :: varTmp_2d(IA,JA,var_num)
    real(RP), intent(inout) :: tend_buf2D_ex(IA,JA,var_num,this%tend_buf_size)
 
    integer :: i,j
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
 
    !----------------------------------------    
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i,j)
      !$omp do
      !$acc parallel present( q1, q2, varTmp_2d, tend_buf2D_ex ) async(1)
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        q1(i,j) =  vartmp_2d(i,j,varid1)                                             &
                + sig_ss * q1(i,j) + gam_ss * tend_buf2d_ex(i,j,varid1,1)
        q2(i,j) =  vartmp_2d(i,j,varid2)                                             &
                + sig_ss * q2(i,j) + gam_ss * tend_buf2d_ex(i,j,varid2,1)
      end do
      end do
      !$acc end parallel
      !$omp end parallel
 
      return
    end if
      
    !$omp parallel private(i,j)
    !$acc parallel present( q1, q2, var0_2D, varTmp_2D, tend_buf2D_ex ) async(1)
    if (nowstage == 1) then
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        var0_2d(i,j,varid1)   = q1(i,j)
        var0_2d(i,j,varid2)   = q2(i,j)
 
        vartmp_2d(i,j,varid1) = 0.0_rp
        vartmp_2d(i,j,varid2) = 0.0_rp
      end do
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        vartmp_2d(i,j,varid1) = vartmp_2d(i,j,varid1)      &
            + sig_ns * q1(i,j) + gam_ns * tend_buf2d_ex(i,j,varid1,1)
        vartmp_2d(i,j,varid2) = vartmp_2d(i,j,varid2)      &
            + sig_ns * q2(i,j) + gam_ns * tend_buf2d_ex(i,j,varid2,1)
      end do
      end do
    end if
 
    !$omp do
    !$acc loop collapse(2)
    do j=js, je
    do i=is, ie
      q1(i,j) = one_minus_sig_ss * var0_2d(i,j,varid1)                &
              + sig_ss * q1(i,j) + gam_ss * tend_buf2d_ex(i,j,varid1,1)
      q2(i,j) = one_minus_sig_ss * var0_2d(i,j,varid2)                &
              + sig_ss * q2(i,j) + gam_ss * tend_buf2d_ex(i,j,varid2,1)
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_advance_low_storage2d_var2
 
!> Advance tracer variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density
!OCL SERIAL
  subroutine rk_advance_trcvar_low_storage2d( this, nowstage, q, varID, is, ie ,js, je , IA,JA, var_num, &
    DDENS, DDENS0, DENS_hyd,                         &
    var0_2D, varTmp_2d, tend_buf2D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA
    real(RP), intent(inout) :: q(IA,JA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
    integer, intent(in) :: var_num
    real(RP), intent(in) :: DDENS (IA,JA)
    real(RP), intent(in) :: DDENS0(IA,JA)
    real(RP), intent(in) :: DENS_hyd(IA,JA)
    real(RP), intent(inout) :: var0_2d(IA,JA,var_num)
    real(RP), intent(inout) :: varTmp_2d(IA,JA,var_num)
    real(RP), intent(inout) :: tend_buf2D_ex(IA,JA,var_num,this%tend_buf_size)
 
    integer :: i,j
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
    real(RP) :: c_ssm1
    real(RP) :: c_ss
 
    real(RP) :: dens_ssm1
    real(RP) :: dens_ss
   
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_trcvar_low_storage2D', 3) 
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)   
    c_ssm1 = this%coef_c_ex(nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i,j)
      !$acc parallel present( q, varTmp_2d, tend_buf2D_ex, DDENS, DDENS0, DENS_hyd )
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        dens_ssm1 = dens_hyd(i,j) + ddens0(i,j) + c_ssm1 * ( ddens(i,j) - ddens0(i,j) )
        q(i,j) =  ( vartmp_2d(i,j,varid)                                              &
                + sig_ss * q(i,j) * dens_ssm1 + gam_ss * tend_buf2d_ex(i,j,varid,1) ) &
                / ( dens_hyd(i,j) + ddens(i,j) )
      end do
      end do
      !$acc end parallel
      !$omp end parallel
      call prof_rapend( 'rk_advance_trcvar_low_storage2D', 3)
 
      return
    end if
 
    c_ss   = this%coef_c_ex(nowstage+1) 
 
    !$omp parallel private(i,j,dens_ss,dens_ssm1)
    !$acc parallel present( q, varTmp_2d, tend_buf2D_ex, DDENS, DDENS0, DENS_hyd )
    if (nowstage == 1) then
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        var0_2d(i,j,varid)   = q(i,j) * ( dens_hyd(i,j) + ddens0(i,j) )
        vartmp_2d(i,j,varid) = 0.0_rp
      end do
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        dens_ssm1 = dens_hyd(i,j) + ddens0(i,j) + c_ssm1 * ( ddens(i,j) - ddens0(i,j) )
        vartmp_2d(i,j,varid) = vartmp_2d(i,j,varid)                      &
            + sig_ns * q(i,j) * dens_ssm1 + gam_ns * tend_buf2d_ex(i,j,varid,1)
      end do
      end do
    end if
 
    !$omp do
      !$acc loop collapse(2)
    do j=js, je
    do i=is, ie
      dens_ssm1 = dens_hyd(i,j) + ddens0(i,j) + c_ssm1 * ( ddens(i,j) - ddens0(i,j) )
      dens_ss   = dens_hyd(i,j) + ddens0(i,j) + c_ss   * ( ddens(i,j) - ddens0(i,j) )
 
      q(i,j) = ( one_minus_sig_ss * var0_2d(i,j,varid)                                &
              + sig_ss * q(i,j) * dens_ssm1 + gam_ss * tend_buf2d_ex(i,j,varid,1) )   &
              / dens_ss
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_trcvar_low_storage2D', 3)
 
    return
  end subroutine rk_advance_trcvar_low_storage2d
  
 
!> Advance variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_low_storage3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke , IA,JA,KA, var_num, &
      var0_3D, varTmp_3d, tend_buf3D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none      
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA,KA
    real(RP), intent(inout) :: q(IA,JA,KA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: varTmp_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: tend_buf3D_ex(IA,JA,KA,var_num,this%tend_buf_size)
 
    integer :: i,j,k
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
 
    !----------------------------------------    
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i,j,k)
      !$omp do collapse(2)
      !$acc parallel present( q, varTmp_3d, tend_buf3D_ex ) async(1)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) =  vartmp_3d(i,j,k,varid)                                             &
                + sig_ss * q(i,j,k) + gam_ss * tend_buf3d_ex(i,j,k,varid,1)
      end do
      end do
      end do
      !$acc end parallel
      !$omp end parallel
 
      return
    end if
      
    !$omp parallel private(i,j,k)
    !$acc parallel present( q, var0_3D, varTmp_3D, tend_buf3D_ex ) async(1)
    if (nowstage == 1) then
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        var0_3d(i,j,k,varid)   = q(i,j,k)
        vartmp_3d(i,j,k,varid) = 0.0_rp
      end do
      end do
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        vartmp_3d(i,j,k,varid) = vartmp_3d(i,j,k,varid)      &
            + sig_ns * q(i,j,k) + gam_ns * tend_buf3d_ex(i,j,k,varid,1)
      end do
      end do
      end do
    end if
 
    !$omp do collapse(2)
    !$acc loop collapse(3)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      q(i,j,k) = one_minus_sig_ss * var0_3d(i,j,k,varid)                &
              + sig_ss * q(i,j,k) + gam_ss * tend_buf3d_ex(i,j,k,varid,1)
    end do
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_advance_low_storage3d
 
!> Advance variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q1 Array to store variable data
!! @param q2 Array to store variable data
!! @param varID1 Index of the targeting variable
!! @param varID2 Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_low_storage3d_var2( this, nowstage, q1, q2, varID1, varID2, is, ie ,js, je ,ks, ke , IA,JA,KA, var_num, &
      var0_3D, varTmp_3d, tend_buf3D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none      
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA,KA
    real(RP), intent(inout) :: q1(IA,JA,KA)
    real(RP), intent(inout) :: q2(IA,JA,KA)
    integer, intent(in) :: varID1
    integer, intent(in) :: varID2
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: varTmp_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: tend_buf3D_ex(IA,JA,KA,var_num,this%tend_buf_size)
 
    integer :: i,j,k
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
 
    !----------------------------------------    
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i,j,k)
      !$omp do collapse(2)
      !$acc parallel present( q1, q2, varTmp_3d, tend_buf3D_ex ) async(1)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q1(i,j,k) =  vartmp_3d(i,j,k,varid1)                                             &
                + sig_ss * q1(i,j,k) + gam_ss * tend_buf3d_ex(i,j,k,varid1,1)
        q2(i,j,k) =  vartmp_3d(i,j,k,varid2)                                             &
                + sig_ss * q2(i,j,k) + gam_ss * tend_buf3d_ex(i,j,k,varid2,1)
      end do
      end do
      end do
      !$acc end parallel
      !$omp end parallel
 
      return
    end if
      
    !$omp parallel private(i,j,k)
    !$acc parallel present( q1, q2, var0_3D, varTmp_3D, tend_buf3D_ex ) async(1)
    if (nowstage == 1) then
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        var0_3d(i,j,k,varid1)   = q1(i,j,k)
        var0_3d(i,j,k,varid2)   = q2(i,j,k)
 
        vartmp_3d(i,j,k,varid1) = 0.0_rp
        vartmp_3d(i,j,k,varid2) = 0.0_rp
      end do
      end do
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        vartmp_3d(i,j,k,varid1) = vartmp_3d(i,j,k,varid1)      &
            + sig_ns * q1(i,j,k) + gam_ns * tend_buf3d_ex(i,j,k,varid1,1)
        vartmp_3d(i,j,k,varid2) = vartmp_3d(i,j,k,varid2)      &
            + sig_ns * q2(i,j,k) + gam_ns * tend_buf3d_ex(i,j,k,varid2,1)
      end do
      end do
      end do
    end if
 
    !$omp do collapse(2)
    !$acc loop collapse(3)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      q1(i,j,k) = one_minus_sig_ss * var0_3d(i,j,k,varid1)                &
              + sig_ss * q1(i,j,k) + gam_ss * tend_buf3d_ex(i,j,k,varid1,1)
      q2(i,j,k) = one_minus_sig_ss * var0_3d(i,j,k,varid2)                &
              + sig_ss * q2(i,j,k) + gam_ss * tend_buf3d_ex(i,j,k,varid2,1)
    end do
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_advance_low_storage3d_var2
 
!> Advance tracer variable data at the current stage with explicit RK part (low-storage case)
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density
!OCL SERIAL
  subroutine rk_advance_trcvar_low_storage3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke , IA,JA,KA, var_num, &
    DDENS, DDENS0, DENS_hyd,                         &
    var0_3D, varTmp_3d, tend_buf3D_ex )
    use scale_const, only: &
      eps => const_eps
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA,KA
    real(RP), intent(inout) :: q(IA,JA,KA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
    integer, intent(in) :: var_num
    real(RP), intent(in) :: DDENS (IA,JA,KA)
    real(RP), intent(in) :: DDENS0(IA,JA,KA)
    real(RP), intent(in) :: DENS_hyd(IA,JA,KA)
    real(RP), intent(inout) :: var0_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: varTmp_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: tend_buf3D_ex(IA,JA,KA,var_num,this%tend_buf_size)
 
    integer :: i,j,k
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
    real(RP) :: c_ssm1
    real(RP) :: c_ss
 
    real(RP) :: dens_ssm1
    real(RP) :: dens_ss
   
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_trcvar_low_storage3D', 3) 
 
    sig_ss = this%coef_sig_ex(nowstage+1,nowstage)
    sig_ns = this%coef_sig_ex(this%nstage+1,nowstage)
    one_minus_sig_ss = 1.0_rp - sig_ss
    gam_ss = this%dt * this%coef_gam_ex(nowstage+1,nowstage)
    gam_ns = this%dt * this%coef_gam_ex(this%nstage+1,nowstage)   
    c_ssm1 = this%coef_c_ex(nowstage)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i,j,k)
      !$acc parallel present( q, varTmp_3d, tend_buf3D_ex, DDENS, DDENS0, DENS_hyd )
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        dens_ssm1 = dens_hyd(i,j,k) + ddens0(i,j,k) + c_ssm1 * ( ddens(i,j,k) - ddens0(i,j,k) )
        q(i,j,k) =  ( vartmp_3d(i,j,k,varid)                                              &
                + sig_ss * q(i,j,k) * dens_ssm1 + gam_ss * tend_buf3d_ex(i,j,k,varid,1) ) &
                / ( dens_hyd(i,j,k) + ddens(i,j,k) )
      end do
      end do
      end do
      !$acc end parallel
      !$omp end parallel
      call prof_rapend( 'rk_advance_trcvar_low_storage3D', 3)
 
      return
    end if
 
    c_ss   = this%coef_c_ex(nowstage+1) 
 
    !$omp parallel private(i,j,k,dens_ss,dens_ssm1)
    !$acc parallel present( q, varTmp_3d, tend_buf3D_ex, DDENS, DDENS0, DENS_hyd )
    if (nowstage == 1) then
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        var0_3d(i,j,k,varid)   = q(i,j,k) * ( dens_hyd(i,j,k) + ddens0(i,j,k) )
        vartmp_3d(i,j,k,varid) = 0.0_rp
      end do
      end do
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        dens_ssm1 = dens_hyd(i,j,k) + ddens0(i,j,k) + c_ssm1 * ( ddens(i,j,k) - ddens0(i,j,k) )
        vartmp_3d(i,j,k,varid) = vartmp_3d(i,j,k,varid)                      &
            + sig_ns * q(i,j,k) * dens_ssm1 + gam_ns * tend_buf3d_ex(i,j,k,varid,1)
      end do
      end do
      end do
    end if
 
    !$omp do collapse(2)
      !$acc loop collapse(3)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      dens_ssm1 = dens_hyd(i,j,k) + ddens0(i,j,k) + c_ssm1 * ( ddens(i,j,k) - ddens0(i,j,k) )
      dens_ss   = dens_hyd(i,j,k) + ddens0(i,j,k) + c_ss   * ( ddens(i,j,k) - ddens0(i,j,k) )
 
      q(i,j,k) = ( one_minus_sig_ss * var0_3d(i,j,k,varid)                                &
              + sig_ss * q(i,j,k) * dens_ssm1 + gam_ss * tend_buf3d_ex(i,j,k,varid,1) )   &
              / dens_ss
    end do
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_trcvar_low_storage3D', 3)
 
    return
  end subroutine rk_advance_trcvar_low_storage3d
  
 
 
!> Advance variable data at the current stage with explicit RK part (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_general1d( this, nowstage, q, varID, is, ie , IA, var_num, &
    var0_1D, varTmp_1d, tend_buf1D_ex, tend_buf1D_im )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA
    real(RP), intent(inout) :: q(IA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_1d(IA,var_num)
    real(RP), intent(inout) :: varTmp_1d(IA,var_num)
    real(RP), intent(inout) :: tend_buf1D_ex(IA,var_num,this%tend_buf_size)
    real(RP), intent(inout), optional :: tend_buf1D_im(IA,var_num,this%tend_buf_size)
 
    integer :: i
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: coef_a_ex_dt
    real(RP) :: coef_a_im_dt
    real(RP) :: coef_b_ex_dt
    real(RP) :: coef_b_im_dt
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_general1D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp parallel do
      !$acc parallel loop collapse(1) present( q, var0_1d, varTmp_1d )
      do i=is, ie
        var0_1d(i,varid) = q(i)
        vartmp_1d(i,varid) = q(i)
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      if ( this%imex_flag ) then
        coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
        coef_b_im_dt = this%coef_b_im(nowstage) * this%dt
        !$omp parallel do
        !$acc parallel loop collapse(1) present( q, varTmp_1d, tend_buf1D_ex, tend_buf1D_im )
        do i=is, ie
          q(i) =  vartmp_1d(i,varid) &
              + coef_b_ex_dt * tend_buf1d_ex(i,varid,tintbuf_ind) & 
              + coef_b_im_dt * tend_buf1d_im(i,varid,tintbuf_ind)   
        end do
      else
        coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
 
        !$omp parallel do
        !$acc parallel loop collapse(1) present( q, varTmp_1d, tend_buf1D_ex )
        do i=is, ie
          q(i) =  vartmp_1d(i,varid)                       &
                  + coef_b_ex_dt * tend_buf1d_ex(i,varid,tintbuf_ind)
        end do
      end if    
      call prof_rapend( 'rk_advance_general1D', 3)
 
      return
    end if 
 
    !$omp parallel private( s, coef_a_ex_dt, coef_a_im_dt, coef_b_ex_dt, coef_b_im_dt ) &
    !$omp private( i )
    !$acc parallel present( q, var0_1D, varTmp_1D, tend_buf1D_ex, tend_buf1D_im )
 
    if ( this%imex_flag ) then
      coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
      coef_b_im_dt = this%coef_b_im(nowstage) * this%dt
 
      !$omp do
      !$acc loop
      do i=is, ie
        q(i) = var0_1d(i,varid)
        vartmp_1d(i,varid) =  vartmp_1d(i,varid)   &
              + coef_b_ex_dt * tend_buf1d_ex(i,varid,tintbuf_ind)  & 
              + coef_b_im_dt * tend_buf1d_im(i,varid,tintbuf_ind)
      end do
    else
      coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
 
      !$omp do
      !$acc loop
      do i=is, ie
        q(i) = var0_1d(i,varid)
        vartmp_1d(i,varid) =  vartmp_1d(i,varid) &
              + coef_b_ex_dt * tend_buf1d_ex(i,varid,tintbuf_ind)
      end do
    end if 
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,nowstage)
      !$omp do
      !$acc loop
      do i=is, ie
        q(i) = var0_1d(i,varid)             &
            +  coef_a_ex_dt * tend_buf1d_ex(i,varid,1)
      end do
    else if ( .not. this%imex_flag ) then      
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
      !$omp do
      !$acc loop
      do i=is, ie
        q(i) = q(i)                            &
            +  coef_a_ex_dt * tend_buf1d_ex(i,varid,s)
      end do
      end do
    else ! IMEX   
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
        coef_a_im_dt = this%dt * this%coef_a_im(nowstage+1,s)            
      !$omp do
      !$acc loop
        do i=is, ie
          q(i) = q(i)                                                 &
            + coef_a_ex_dt * tend_buf1d_ex(i,varid,s) &
            + coef_a_im_dt * tend_buf1d_im(i,varid,s)  
        end do
      end do
    end if
    !$acc end parallel
    !$omp end parallel
    call prof_rapend( 'rk_advance_general1D', 3)
 
    return
  end subroutine rk_advance_general1d
 
!> Advance tracer variable data at the current stage with explicit RK part (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density  
!OCL SERIAL
  subroutine rk_advance_trcvar_general1d( this, nowstage, q, varID, is, ie , IA, var_num, &
    DDENS, DDENS0, DENS_hyd,                         &
    var0_1D, varTmp_1d, tend_buf1D_ex )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA
    real(RP), intent(inout) :: q(IA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
    integer, intent(in) :: var_num
    real(RP), intent(in) :: DDENS (IA)
    real(RP), intent(in) :: DDENS0(IA)
    real(RP), intent(in) :: DENS_hyd(IA)
    real(RP), intent(inout) :: var0_1d(IA,var_num)
    real(RP), intent(inout) :: varTmp_1d(IA,var_num)
    real(RP), intent(inout) :: tend_buf1D_ex(IA,var_num,this%tend_buf_size)
 
    integer :: i
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: dens_
 
    real(RP) :: coef_a_ex
    real(RP) :: coef_a_ex_dt
    real(RP) :: coef_b_ex_dt
    real(RP) :: c_ss
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_trcvar_general1D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( this%nstage == 1 ) then
      !$omp parallel do
      !$acc parallel loop present( q, varTmp_1d, DDENS0, DENS_hyd )
      do i=is, ie
        vartmp_1d(i,varid) = q(i) * ( dens_hyd(i) + ddens0(i) )
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      coef_b_ex_dt = this%dt * this%coef_b_ex(nowstage)
      !$omp parallel do
      !$acc parallel loop present( q, varTmp_1d, DDENS, DENS_hyd )
      do i=is, ie
        q(i) =  ( vartmp_1d(i,varid)                    &
            + coef_b_ex_dt * tend_buf1d_ex(i,varid,tintbuf_ind)  ) &
            / ( dens_hyd(i) + ddens(i) )
      end do
      call prof_rapend( 'rk_advance_trcvar_general1D', 3)
 
      return
    end if 
 
    c_ss = this%coef_c_ex(nowstage+1) 
 
    !$omp parallel private( s, dens_, coef_a_ex, coef_a_ex_dt, coef_b_ex_dt ) &
    !$omp private( i )
    !$acc parallel present( q, var0_1D, varTmp_1D, tend_buf1D_ex, DDENS, DDENS0, DENS_hyd )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      !$acc loop
      do i=is, ie
        var0_1d(i,varid) = q(i) * ( dens_hyd(i) + ddens0(i) )
        vartmp_1d(i,varid) = var0_1d(i,varid)
      end do
    end if 
 
    coef_b_ex_dt = this%dt * this%coef_b_ex(nowstage)
    !$omp do
    !$acc loop
    do i=is, ie
      q(i) = this%var0_1d(i,varid)
      vartmp_1d(i,varid) =  vartmp_1d(i,varid)                                             &
            + coef_b_ex_dt * tend_buf1d_ex(i,varid,tintbuf_ind)
    end do
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      coef_a_ex = this%coef_a_ex(nowstage+1,nowstage)
      coef_b_ex_dt = this%dt * this%coef_a_ex(nowstage+1,nowstage)
      !$omp do
      !$acc loop
      do i=is, ie
        dens_ = ( dens_hyd(i) + ddens0(i) ) &
            + coef_a_ex * ( ddens(i) - ddens0(i) )
        q(i) = ( var0_1d(i,varid)            &
            + coef_b_ex_dt * tend_buf1d_ex(i,varid,1) ) &
            / dens_
      end do
    else if ( .not. this%imex_flag ) then    
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
        !$omp do
        !$acc loop
        do i=is, ie
          q(i) = q(i)                           &
            + coef_a_ex_dt * tend_buf1d_ex(i,varid,s)
        end do
      end do
      !$omp do
      !$acc loop
      do i=is, ie
        dens_ = dens_hyd(i) + ddens0(i) + c_ss * ( ddens(i) - ddens0(i) )
        q(i) = q(i) / dens_
      end do
    end if
    !$acc end parallel
    !$omp end parallel
    call prof_rapend( 'rk_advance_trcvar_generall1D', 3)
 
    return
  end subroutine rk_advance_trcvar_general1d
 
!> Store variable data after the implicit part of current stage in IMEX RK scheme (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data after the implicit part of current stage
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_storeimpl_general1d( this, nowstage, q, varID, is, ie , IA, var_num, &
    var0_1D, varTmp_1d, tend_buf1D_im )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA
    real(RP), intent(inout) :: q(IA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_1d(IA,var_num)
    real(RP), intent(inout) :: varTmp_1d(IA,var_num)
    real(RP), intent(inout) :: tend_buf1D_im(IA,var_num,this%tend_buf_size)
 
    integer :: i
    integer :: tintbuf_ind
    real(RP) :: coef_a_im_dt
    !----------------------------------------    
 
    if (.not. this%imex_flag ) return
    
    tintbuf_ind = this%tend_buf_indmap(nowstage)   
    coef_a_im_dt = this%dt * this%coef_a_im(nowstage,nowstage)
 
    !$omp parallel
    !$acc parallel present( q, var0_1D, varTmp_1D, tend_buf1D_im )
 
    if ( nowstage == 1 ) then
      !$omp do
      !$acc loop
      do i=is, ie
        var0_1d(i,varid) = q(i)
        vartmp_1d(i,varid) = q(i)
      end do
    end if
            
    !$omp do
    !$acc loop
    do i=is, ie
      q(i) = q(i)                                                                              &
            + coef_a_im_dt * tend_buf1d_im(i,varid,tintbuf_ind)  
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_storeimpl_general1d
 
 
!> Advance variable data at the current stage with explicit RK part (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_general2d( this, nowstage, q, varID, is, ie ,js, je , IA,JA, var_num, &
    var0_2D, varTmp_2d, tend_buf2D_ex, tend_buf2D_im )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA
    real(RP), intent(inout) :: q(IA,JA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_2d(IA,JA,var_num)
    real(RP), intent(inout) :: varTmp_2d(IA,JA,var_num)
    real(RP), intent(inout) :: tend_buf2D_ex(IA,JA,var_num,this%tend_buf_size)
    real(RP), intent(inout), optional :: tend_buf2D_im(IA,JA,var_num,this%tend_buf_size)
 
    integer :: i,j
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: coef_a_ex_dt
    real(RP) :: coef_a_im_dt
    real(RP) :: coef_b_ex_dt
    real(RP) :: coef_b_im_dt
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_general2D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp parallel do
      !$acc parallel loop collapse(2) present( q, var0_2d, varTmp_2d )
      do j=js, je
      do i=is, ie
        var0_2d(i,j,varid) = q(i,j)
        vartmp_2d(i,j,varid) = q(i,j)
      end do
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      if ( this%imex_flag ) then
        coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
        coef_b_im_dt = this%coef_b_im(nowstage) * this%dt
        !$omp parallel do
        !$acc parallel loop collapse(2) present( q, varTmp_2d, tend_buf2D_ex, tend_buf2D_im )
        do j=js, je
        do i=is, ie
          q(i,j) =  vartmp_2d(i,j,varid) &
              + coef_b_ex_dt * tend_buf2d_ex(i,j,varid,tintbuf_ind) & 
              + coef_b_im_dt * tend_buf2d_im(i,j,varid,tintbuf_ind)   
        end do
        end do
      else
        coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
 
        !$omp parallel do
        !$acc parallel loop collapse(2) present( q, varTmp_2d, tend_buf2D_ex )
        do j=js, je
        do i=is, ie
          q(i,j) =  vartmp_2d(i,j,varid)                       &
                  + coef_b_ex_dt * tend_buf2d_ex(i,j,varid,tintbuf_ind)
        end do
        end do
      end if    
      call prof_rapend( 'rk_advance_general2D', 3)
 
      return
    end if 
 
    !$omp parallel private( s, coef_a_ex_dt, coef_a_im_dt, coef_b_ex_dt, coef_b_im_dt ) &
    !$omp private( i,j )
    !$acc parallel present( q, var0_2D, varTmp_2D, tend_buf2D_ex, tend_buf2D_im )
 
    if ( this%imex_flag ) then
      coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
      coef_b_im_dt = this%coef_b_im(nowstage) * this%dt
 
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        q(i,j) = var0_2d(i,j,varid)
        vartmp_2d(i,j,varid) =  vartmp_2d(i,j,varid)   &
              + coef_b_ex_dt * tend_buf2d_ex(i,j,varid,tintbuf_ind)  & 
              + coef_b_im_dt * tend_buf2d_im(i,j,varid,tintbuf_ind)
      end do
      end do
    else
      coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
 
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        q(i,j) = var0_2d(i,j,varid)
        vartmp_2d(i,j,varid) =  vartmp_2d(i,j,varid) &
              + coef_b_ex_dt * tend_buf2d_ex(i,j,varid,tintbuf_ind)
      end do
      end do
    end if 
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,nowstage)
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        q(i,j) = var0_2d(i,j,varid)             &
            +  coef_a_ex_dt * tend_buf2d_ex(i,j,varid,1)
      end do
      end do
    else if ( .not. this%imex_flag ) then      
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        q(i,j) = q(i,j)                            &
            +  coef_a_ex_dt * tend_buf2d_ex(i,j,varid,s)
      end do
      end do
      end do
    else ! IMEX   
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
        coef_a_im_dt = this%dt * this%coef_a_im(nowstage+1,s)            
      !$omp do
      !$acc loop collapse(2)
        do j=js, je
        do i=is, ie
          q(i,j) = q(i,j)                                                 &
            + coef_a_ex_dt * tend_buf2d_ex(i,j,varid,s) &
            + coef_a_im_dt * tend_buf2d_im(i,j,varid,s)  
        end do
        end do
      end do
    end if
    !$acc end parallel
    !$omp end parallel
    call prof_rapend( 'rk_advance_general2D', 3)
 
    return
  end subroutine rk_advance_general2d
 
!> Advance tracer variable data at the current stage with explicit RK part (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density  
!OCL SERIAL
  subroutine rk_advance_trcvar_general2d( this, nowstage, q, varID, is, ie ,js, je , IA,JA, var_num, &
    DDENS, DDENS0, DENS_hyd,                         &
    var0_2D, varTmp_2d, tend_buf2D_ex )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA
    real(RP), intent(inout) :: q(IA,JA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
    integer, intent(in) :: var_num
    real(RP), intent(in) :: DDENS (IA,JA)
    real(RP), intent(in) :: DDENS0(IA,JA)
    real(RP), intent(in) :: DENS_hyd(IA,JA)
    real(RP), intent(inout) :: var0_2d(IA,JA,var_num)
    real(RP), intent(inout) :: varTmp_2d(IA,JA,var_num)
    real(RP), intent(inout) :: tend_buf2D_ex(IA,JA,var_num,this%tend_buf_size)
 
    integer :: i,j
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: dens_
 
    real(RP) :: coef_a_ex
    real(RP) :: coef_a_ex_dt
    real(RP) :: coef_b_ex_dt
    real(RP) :: c_ss
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_trcvar_general2D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( this%nstage == 1 ) then
      !$omp parallel do
      !$acc parallel loop collapse(2) present( q, varTmp_2d, DDENS0, DENS_hyd )
      do j=js, je
      do i=is, ie
        vartmp_2d(i,j,varid) = q(i,j) * ( dens_hyd(i,j) + ddens0(i,j) )
      end do
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      coef_b_ex_dt = this%dt * this%coef_b_ex(nowstage)
      !$omp parallel do
      !$acc parallel loop collapse(2) present( q, varTmp_2d, DDENS, DENS_hyd )
      do j=js, je
      do i=is, ie
        q(i,j) =  ( vartmp_2d(i,j,varid)                    &
            + coef_b_ex_dt * tend_buf2d_ex(i,j,varid,tintbuf_ind)  ) &
            / ( dens_hyd(i,j) + ddens(i,j) )
      end do
      end do
      call prof_rapend( 'rk_advance_trcvar_general2D', 3)
 
      return
    end if 
 
    c_ss = this%coef_c_ex(nowstage+1) 
 
    !$omp parallel private( s, dens_, coef_a_ex, coef_a_ex_dt, coef_b_ex_dt ) &
    !$omp private( i,j )
    !$acc parallel present( q, var0_2D, varTmp_2D, tend_buf2D_ex, DDENS, DDENS0, DENS_hyd )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        var0_2d(i,j,varid) = q(i,j) * ( dens_hyd(i,j) + ddens0(i,j) )
        vartmp_2d(i,j,varid) = var0_2d(i,j,varid)
      end do
      end do
    end if 
 
    coef_b_ex_dt = this%dt * this%coef_b_ex(nowstage)
    !$omp do
      !$acc loop collapse(2)
    do j=js, je
    do i=is, ie
      q(i,j) = this%var0_2d(i,j,varid)
      vartmp_2d(i,j,varid) =  vartmp_2d(i,j,varid)                                             &
            + coef_b_ex_dt * tend_buf2d_ex(i,j,varid,tintbuf_ind)
    end do
    end do
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      coef_a_ex = this%coef_a_ex(nowstage+1,nowstage)
      coef_b_ex_dt = this%dt * this%coef_a_ex(nowstage+1,nowstage)
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        dens_ = ( dens_hyd(i,j) + ddens0(i,j) ) &
            + coef_a_ex * ( ddens(i,j) - ddens0(i,j) )
        q(i,j) = ( var0_2d(i,j,varid)            &
            + coef_b_ex_dt * tend_buf2d_ex(i,j,varid,1) ) &
            / dens_
      end do
      end do
    else if ( .not. this%imex_flag ) then    
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
        !$omp do
        !$acc loop collapse(2)
        do j=js, je
        do i=is, ie
          q(i,j) = q(i,j)                           &
            + coef_a_ex_dt * tend_buf2d_ex(i,j,varid,s)
        end do
        end do
      end do
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        dens_ = dens_hyd(i,j) + ddens0(i,j) + c_ss * ( ddens(i,j) - ddens0(i,j) )
        q(i,j) = q(i,j) / dens_
      end do
      end do
    end if
    !$acc end parallel
    !$omp end parallel
    call prof_rapend( 'rk_advance_trcvar_generall2D', 3)
 
    return
  end subroutine rk_advance_trcvar_general2d
 
!> Store variable data after the implicit part of current stage in IMEX RK scheme (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data after the implicit part of current stage
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_storeimpl_general2d( this, nowstage, q, varID, is, ie ,js, je , IA,JA, var_num, &
    var0_2D, varTmp_2d, tend_buf2D_im )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA
    real(RP), intent(inout) :: q(IA,JA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_2d(IA,JA,var_num)
    real(RP), intent(inout) :: varTmp_2d(IA,JA,var_num)
    real(RP), intent(inout) :: tend_buf2D_im(IA,JA,var_num,this%tend_buf_size)
 
    integer :: i,j
    integer :: tintbuf_ind
    real(RP) :: coef_a_im_dt
    !----------------------------------------    
 
    if (.not. this%imex_flag ) return
    
    tintbuf_ind = this%tend_buf_indmap(nowstage)   
    coef_a_im_dt = this%dt * this%coef_a_im(nowstage,nowstage)
 
    !$omp parallel
    !$acc parallel present( q, var0_2D, varTmp_2D, tend_buf2D_im )
 
    if ( nowstage == 1 ) then
      !$omp do
      !$acc loop collapse(2)
      do j=js, je
      do i=is, ie
        var0_2d(i,j,varid) = q(i,j)
        vartmp_2d(i,j,varid) = q(i,j)
      end do
      end do
    end if
            
    !$omp do
    !$acc loop collapse(2)
    do j=js, je
    do i=is, ie
      q(i,j) = q(i,j)                                                                              &
            + coef_a_im_dt * tend_buf2d_im(i,j,varid,tintbuf_ind)  
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_storeimpl_general2d
 
 
!> Advance variable data at the current stage with explicit RK part (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_advance_general3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke , IA,JA,KA, var_num, &
    var0_3D, varTmp_3d, tend_buf3D_ex, tend_buf3D_im )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA,KA
    real(RP), intent(inout) :: q(IA,JA,KA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: varTmp_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: tend_buf3D_ex(IA,JA,KA,var_num,this%tend_buf_size)
    real(RP), intent(inout), optional :: tend_buf3D_im(IA,JA,KA,var_num,this%tend_buf_size)
 
    integer :: i,j,k
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: coef_a_ex_dt
    real(RP) :: coef_a_im_dt
    real(RP) :: coef_b_ex_dt
    real(RP) :: coef_b_im_dt
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_general3D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp parallel do collapse(2)
      !$acc parallel loop collapse(3) present( q, var0_3d, varTmp_3d )
      do k=ks, ke
      do j=js, je
      do i=is, ie
        var0_3d(i,j,k,varid) = q(i,j,k)
        vartmp_3d(i,j,k,varid) = q(i,j,k)
      end do
      end do
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      if ( this%imex_flag ) then
        coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
        coef_b_im_dt = this%coef_b_im(nowstage) * this%dt
        !$omp parallel do collapse(2)
        !$acc parallel loop collapse(3) present( q, varTmp_3d, tend_buf3D_ex, tend_buf3D_im )
        do k=ks, ke
        do j=js, je
        do i=is, ie
          q(i,j,k) =  vartmp_3d(i,j,k,varid) &
              + coef_b_ex_dt * tend_buf3d_ex(i,j,k,varid,tintbuf_ind) & 
              + coef_b_im_dt * tend_buf3d_im(i,j,k,varid,tintbuf_ind)   
        end do
        end do
        end do
      else
        coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
 
        !$omp parallel do collapse(2)
        !$acc parallel loop collapse(3) present( q, varTmp_3d, tend_buf3D_ex )
        do k=ks, ke
        do j=js, je
        do i=is, ie
          q(i,j,k) =  vartmp_3d(i,j,k,varid)                       &
                  + coef_b_ex_dt * tend_buf3d_ex(i,j,k,varid,tintbuf_ind)
        end do
        end do
        end do
      end if    
      call prof_rapend( 'rk_advance_general3D', 3)
 
      return
    end if 
 
    !$omp parallel private( s, coef_a_ex_dt, coef_a_im_dt, coef_b_ex_dt, coef_b_im_dt ) &
    !$omp private( i,j,k )
    !$acc parallel present( q, var0_3D, varTmp_3D, tend_buf3D_ex, tend_buf3D_im )
 
    if ( this%imex_flag ) then
      coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
      coef_b_im_dt = this%coef_b_im(nowstage) * this%dt
 
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = var0_3d(i,j,k,varid)
        vartmp_3d(i,j,k,varid) =  vartmp_3d(i,j,k,varid)   &
              + coef_b_ex_dt * tend_buf3d_ex(i,j,k,varid,tintbuf_ind)  & 
              + coef_b_im_dt * tend_buf3d_im(i,j,k,varid,tintbuf_ind)
      end do
      end do
      end do
    else
      coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt
 
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = var0_3d(i,j,k,varid)
        vartmp_3d(i,j,k,varid) =  vartmp_3d(i,j,k,varid) &
              + coef_b_ex_dt * tend_buf3d_ex(i,j,k,varid,tintbuf_ind)
      end do
      end do
      end do
    end if 
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,nowstage)
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = var0_3d(i,j,k,varid)             &
            +  coef_a_ex_dt * tend_buf3d_ex(i,j,k,varid,1)
      end do
      end do
      end do
    else if ( .not. this%imex_flag ) then      
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = q(i,j,k)                            &
            +  coef_a_ex_dt * tend_buf3d_ex(i,j,k,varid,s)
      end do
      end do
      end do
      end do
    else ! IMEX   
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
        coef_a_im_dt = this%dt * this%coef_a_im(nowstage+1,s)            
      !$omp do collapse(2)
      !$acc loop collapse(3)
        do k=ks, ke
        do j=js, je
        do i=is, ie
          q(i,j,k) = q(i,j,k)                                                 &
            + coef_a_ex_dt * tend_buf3d_ex(i,j,k,varid,s) &
            + coef_a_im_dt * tend_buf3d_im(i,j,k,varid,s)  
        end do
        end do
        end do
      end do
    end if
    !$acc end parallel
    !$omp end parallel
    call prof_rapend( 'rk_advance_general3D', 3)
 
    return
  end subroutine rk_advance_general3d
 
!> Advance tracer variable data at the current stage with explicit RK part (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store tracer variable data
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
!! @param DDENS  Array to store density perturbation data at the time level n+1
!! @param DDENS0  Array to store density perturbation  data at the time level n
!! @param DENS_hyd Array to hydrostatic part of density  
!OCL SERIAL
  subroutine rk_advance_trcvar_general3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke , IA,JA,KA, var_num, &
    DDENS, DDENS0, DENS_hyd,                         &
    var0_3D, varTmp_3d, tend_buf3D_ex )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA,KA
    real(RP), intent(inout) :: q(IA,JA,KA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
    integer, intent(in) :: var_num
    real(RP), intent(in) :: DDENS (IA,JA,KA)
    real(RP), intent(in) :: DDENS0(IA,JA,KA)
    real(RP), intent(in) :: DENS_hyd(IA,JA,KA)
    real(RP), intent(inout) :: var0_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: varTmp_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: tend_buf3D_ex(IA,JA,KA,var_num,this%tend_buf_size)
 
    integer :: i,j,k
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: dens_
 
    real(RP) :: coef_a_ex
    real(RP) :: coef_a_ex_dt
    real(RP) :: coef_b_ex_dt
    real(RP) :: c_ss
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_trcvar_general3D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( this%nstage == 1 ) then
      !$omp parallel do collapse(2)
      !$acc parallel loop collapse(3) present( q, varTmp_3d, DDENS0, DENS_hyd )
      do k=ks, ke
      do j=js, je
      do i=is, ie
        vartmp_3d(i,j,k,varid) = q(i,j,k) * ( dens_hyd(i,j,k) + ddens0(i,j,k) )
      end do
      end do
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      coef_b_ex_dt = this%dt * this%coef_b_ex(nowstage)
      !$omp parallel do collapse(2)
      !$acc parallel loop collapse(3) present( q, varTmp_3d, DDENS, DENS_hyd )
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) =  ( vartmp_3d(i,j,k,varid)                    &
            + coef_b_ex_dt * tend_buf3d_ex(i,j,k,varid,tintbuf_ind)  ) &
            / ( dens_hyd(i,j,k) + ddens(i,j,k) )
      end do
      end do
      end do
      call prof_rapend( 'rk_advance_trcvar_general3D', 3)
 
      return
    end if 
 
    c_ss = this%coef_c_ex(nowstage+1) 
 
    !$omp parallel private( s, dens_, coef_a_ex, coef_a_ex_dt, coef_b_ex_dt ) &
    !$omp private( i,j,k )
    !$acc parallel present( q, var0_3D, varTmp_3D, tend_buf3D_ex, DDENS, DDENS0, DENS_hyd )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        var0_3d(i,j,k,varid) = q(i,j,k) * ( dens_hyd(i,j,k) + ddens0(i,j,k) )
        vartmp_3d(i,j,k,varid) = var0_3d(i,j,k,varid)
      end do
      end do
      end do
    end if 
 
    coef_b_ex_dt = this%dt * this%coef_b_ex(nowstage)
    !$omp do collapse(2)
      !$acc loop collapse(3)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      q(i,j,k) = this%var0_3d(i,j,k,varid)
      vartmp_3d(i,j,k,varid) =  vartmp_3d(i,j,k,varid)                                             &
            + coef_b_ex_dt * tend_buf3d_ex(i,j,k,varid,tintbuf_ind)
    end do
    end do
    end do
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      coef_a_ex = this%coef_a_ex(nowstage+1,nowstage)
      coef_b_ex_dt = this%dt * this%coef_a_ex(nowstage+1,nowstage)
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        dens_ = ( dens_hyd(i,j,k) + ddens0(i,j,k) ) &
            + coef_a_ex * ( ddens(i,j,k) - ddens0(i,j,k) )
        q(i,j,k) = ( var0_3d(i,j,k,varid)            &
            + coef_b_ex_dt * tend_buf3d_ex(i,j,k,varid,1) ) &
            / dens_
      end do
      end do
      end do
    else if ( .not. this%imex_flag ) then    
      do s=1, nowstage
        coef_a_ex_dt = this%dt * this%coef_a_ex(nowstage+1,s)
        !$omp do collapse(2)
        !$acc loop collapse(3)
        do k=ks, ke
        do j=js, je
        do i=is, ie
          q(i,j,k) = q(i,j,k)                           &
            + coef_a_ex_dt * tend_buf3d_ex(i,j,k,varid,s)
        end do
        end do
        end do
      end do
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        dens_ = dens_hyd(i,j,k) + ddens0(i,j,k) + c_ss * ( ddens(i,j,k) - ddens0(i,j,k) )
        q(i,j,k) = q(i,j,k) / dens_
      end do
      end do
      end do
    end if
    !$acc end parallel
    !$omp end parallel
    call prof_rapend( 'rk_advance_trcvar_generall3D', 3)
 
    return
  end subroutine rk_advance_trcvar_general3d
 
!> Store variable data after the implicit part of current stage in IMEX RK scheme (general case)
!!
!! @param nowstage Current stage
!! @param q Array to store variable data after the implicit part of current stage
!! @param varID Index of the targeting variable
!! @param is Index at which begins the loop for the corresponding direction
!! @param ie Index at which finishes the loop for the corresponding direction
!! @param js Index at which begins the loop for the corresponding direction
!! @param je Index at which finishes the loop for the corresponding direction
!! @param ks Index at which begins the loop for the corresponding direction
!! @param ke Index at which finishes the loop for the corresponding direction
!OCL SERIAL
  subroutine rk_storeimpl_general3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke , IA,JA,KA, var_num, &
    var0_3D, varTmp_3d, tend_buf3D_im )
    implicit none
    class(timeint_rk), intent(inout) :: this
    integer, intent(in) :: nowstage
    integer, intent(in) :: IA,JA,KA
    real(RP), intent(inout) :: q(IA,JA,KA)
    integer, intent(in) :: varID
    integer, intent(in) :: is, ie ,js, je ,ks, ke 
    integer, intent(in) :: var_num
    real(RP), intent(inout) :: var0_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: varTmp_3d(IA,JA,KA,var_num)
    real(RP), intent(inout) :: tend_buf3D_im(IA,JA,KA,var_num,this%tend_buf_size)
 
    integer :: i,j,k
    integer :: tintbuf_ind
    real(RP) :: coef_a_im_dt
    !----------------------------------------    
 
    if (.not. this%imex_flag ) return
    
    tintbuf_ind = this%tend_buf_indmap(nowstage)   
    coef_a_im_dt = this%dt * this%coef_a_im(nowstage,nowstage)
 
    !$omp parallel
    !$acc parallel present( q, var0_3D, varTmp_3D, tend_buf3D_im )
 
    if ( nowstage == 1 ) then
      !$omp do collapse(2)
      !$acc loop collapse(3)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        var0_3d(i,j,k,varid) = q(i,j,k)
        vartmp_3d(i,j,k,varid) = q(i,j,k)
      end do
      end do
      end do
    end if
            
    !$omp do collapse(2)
    !$acc loop collapse(3)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      q(i,j,k) = q(i,j,k)                                                                              &
            + coef_a_im_dt * tend_buf3d_im(i,j,k,varid,tintbuf_ind)  
    end do
    end do
    end do
    !$acc end parallel
    !$omp end parallel
 
    return
  end subroutine rk_storeimpl_general3d
 
end module scale_timeint_rk