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.
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#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
  !
  type, public :: timeint_rk
    real(RP), private :: dt
 
    integer, private :: rk_scheme_id_ex
    integer, private :: rk_scheme_id_im
 
    integer, public :: var_num
    integer, private, allocatable :: size_each_var(:)
    integer, public :: nstage
    integer, public :: var_buf_size
    integer, public :: tend_buf_size
 
    ! For Butcher representation (explicit part)
    real(rp), public, allocatable :: coef_a_ex(:,:)
    real(rp), public, allocatable :: coef_b_ex(:)
    real(rp), public, allocatable :: coef_c_ex(:)
 
    ! 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(:,:)
    real(rp), public, allocatable :: coef_b_im(:)
    real(rp), public, allocatable :: coef_c_im(:)
 
    integer, allocatable :: tend_buf_indmap(:)
    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 :: 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 :: 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 :: 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_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
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  !-------------------
 
contains
 
!----------------
 
  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)
    !----------------------------------------
 
    this%dt = dt
    this%ndim = ndim
    this%var_num  = 1
    allocate( this%size_each_var(ndim) )
    this%size_each_var(:) = 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) )
      if ( this%imex_flag ) then
        allocate( this%tend_buf1D_im(size_each_var(1), var_num, this%tend_buf_size) )
      end if
      allocate( this%var0_1D(size_each_var(1), var_num) )
      allocate( this%varTmp_1D(size_each_var(1), var_num) ) 
    case(2)
      allocate( this%tend_buf2D_ex(size_each_var(1),size_each_var(2), var_num, this%tend_buf_size) )
      if ( this%imex_flag ) then
        allocate( this%tend_buf2D_im(size_each_var(1),size_each_var(2), var_num, this%tend_buf_size) )
      end if
      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) ) 
    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) )
      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) )
      end if
      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) ) 
    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)
 
    return
  end subroutine timeint_rk_init
 
  subroutine timeint_rk_final( this )
    implicit none
    class(timeint_rk), intent(inout) :: this
    !----------------------------------------
    
    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
      deallocate( this%coef_a_im, this%coef_b_im, this%coef_c_im )
!    end if
 
    deallocate( this%size_each_var, this%tend_buf_indmap )    
 
    select case(this%ndim)
    case(1)
      deallocate( this%var0_1d )
      deallocate( this%tend_buf1D_ex )
      if (this%imex_flag) deallocate( this%tend_buf1D_im )
      deallocate( this%varTmp_1d )
    case(2)
      deallocate( this%var0_2d )
      deallocate( this%tend_buf2D_ex )
      if (this%imex_flag) deallocate( this%tend_buf2D_im )
      deallocate( this%varTmp_2d )
    case(3)
      deallocate( this%var0_3d )
      deallocate( this%tend_buf3D_ex )
      if (this%imex_flag) deallocate( this%tend_buf3D_im )
      deallocate( this%varTmp_3d )
    end select
 
    return
  end subroutine timeint_rk_final
 
  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
 
  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
  !----------------
 
 
  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 rk_advance_low_storage1d( this, nowstage, q, varid, is, ie )
    else
      call rk_advance_general1d( this, nowstage, q, varid, is, ie  )
    end if
 
    return
  end subroutine timeint_rk_advance1d
 
  subroutine timeint_rk_advance_trcvar1d( this, nowstage, q, varID, is, ie , &
      DENS, DENS0, 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) :: DENS (:)
    real(RP), intent(in) :: DENS0(:)
    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 , &
        dens, dens0, dens_hyd )
    else
      call rk_advance_trcvar_general1d( this, nowstage, q, varid, is, ie , &
        dens, dens0, dens_hyd )
    end if
 
    return
  end subroutine timeint_rk_advance_trcvar1d
 
  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)
    !$omp do
    do i=is, ie
      this%var0_1D(i,varid) = q(i)
    end do
      !$omp end parallel
 
    return
  end subroutine timeint_rk_store_var0_1d
  
  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  )
 
    return
  end subroutine timeint_rk_storeimpl1d
 
  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 rk_advance_low_storage2d( this, nowstage, q, varid, is, ie ,js, je )
    else
      call rk_advance_general2d( this, nowstage, q, varid, is, ie ,js, je  )
    end if
 
    return
  end subroutine timeint_rk_advance2d
 
  subroutine timeint_rk_advance_trcvar2d( this, nowstage, q, varID, is, ie ,js, je , &
      DENS, DENS0, 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) :: DENS (:,:)
    real(RP), intent(in) :: DENS0(:,:)
    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 , &
        dens, dens0, dens_hyd )
    else
      call rk_advance_trcvar_general2d( this, nowstage, q, varid, is, ie ,js, je , &
        dens, dens0, dens_hyd )
    end if
 
    return
  end subroutine timeint_rk_advance_trcvar2d
 
  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)
    !$omp do
    do j=js, je
    do i=is, ie
      this%var0_2D(i,j,varid) = q(i,j)
    end do
    end do
      !$omp end parallel
 
    return
  end subroutine timeint_rk_store_var0_2d
  
  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  )
 
    return
  end subroutine timeint_rk_storeimpl2d
 
  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 rk_advance_low_storage3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke )
    else
      call rk_advance_general3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke  )
    end if
 
    return
  end subroutine timeint_rk_advance3d
 
  subroutine timeint_rk_advance_trcvar3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke , &
      DENS, DENS0, 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) :: DENS (:,:,:)
    real(RP), intent(in) :: DENS0(:,:,:)
    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 , &
        dens, dens0, dens_hyd )
    else
      call rk_advance_trcvar_general3d( this, nowstage, q, varid, is, ie ,js, je ,ks, ke , &
        dens, dens0, dens_hyd )
    end if
 
    return
  end subroutine timeint_rk_advance_trcvar3d
 
  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)
    !$omp do collapse(2)
    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
      !$omp end parallel
 
    return
  end subroutine timeint_rk_store_var0_3d
  
  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  )
 
    return
  end subroutine timeint_rk_storeimpl3d
 
!-------------------
 
 
 
!OCL SERIAL
  subroutine rk_advance_low_storage1d( this, nowstage, q, varID, is, ie  )
    use scale_const, only: &
      eps => const_eps
    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 
 
    integer :: i
    real(RP) :: sig_ss
    real(RP) :: gam_ss
    real(RP) :: one_Minus_sig_ss
    real(RP) :: sig_Ns
    real(RP) :: gam_Ns
 
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_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)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i)
      !$omp do
      do i=is, ie
        q(i) =  this%varTmp_1d(i,varid)                                             &
                + sig_ss * q(i) + gam_ss * this%tend_buf1D_ex(i,varid,1)
      end do
      !$omp end parallel
      call prof_rapend( 'rk_advance_low_storage1D', 3)
 
      return
    end if
      
    !$omp parallel private(i)
    if (nowstage == 1) then
      !$omp do
      do i=is, ie
        this%var0_1D(i,varid)   = q(i)
        this%varTmp_1D(i,varid) = 0.0_rp
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do
      do i=is, ie
        this%varTmp_1d(i,varid) = this%varTmp_1d(i,varid)      &
            + sig_ns * q(i) + gam_ns * this%tend_buf1D_ex(i,varid,1)
      end do
    end if
 
    !$omp do
    do i=is, ie
      q(i) = one_minus_sig_ss * this%var0_1d(i,varid)                &
              + sig_ss * q(i) + gam_ss * this%tend_buf1D_ex(i,varid,1)
    end do
 
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_low_storage1D', 3)
 
    return
  end subroutine rk_advance_low_storage1d
 
  !OCL SERIAL
  subroutine rk_advance_trcvar_low_storage1d( this, nowstage, q, varID, is, ie , &
    DDENS, DDENS0, DENS_hyd )
    use scale_const, only: &
      eps => const_eps
    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(:)
 
    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)
      !$omp do
      do i=is, ie
        dens_ssm1 = dens_hyd(i) + ddens0(i) + c_ssm1 * ( ddens(i) - ddens0(i) )
        q(i) =  ( this%varTmp_1d(i,varid)                                              &
                + sig_ss * q(i) * dens_ssm1 + gam_ss * this%tend_buf1D_ex(i,varid,1) ) &
                / ( dens_hyd(i) + ddens(i) )
      end do
      !$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)
    if (nowstage == 1) then
      !$omp do
      do i=is, ie
        this%var0_1D(i,varid)   = q(i) * ( dens_hyd(i) + ddens0(i) )
        this%varTmp_1D(i,varid) = 0.0_rp
      end do
    end if
 
    if ( abs(sig_ns) > eps .or. abs(gam_ns) > eps ) then
      !$omp do
      do i=is, ie
        dens_ssm1 = dens_hyd(i) + ddens0(i) + c_ssm1 * ( ddens(i) - ddens0(i) )
        this%varTmp_1d(i,varid) = this%varTmp_1d(i,varid)                 &
            + sig_ns * q(i) * dens_ssm1 + gam_ns * this%tend_buf1D_ex(i,varid,1)
      end do
    end if
 
    !$omp do
    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 * this%var0_1d(i,varid)                                &
              + sig_ss * q(i) * dens_ssm1 + gam_ss * this%tend_buf1D_ex(i,varid,1) )   &
              / dens_ss
    end do
 
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_trcvar_low_storage1D', 3)
 
    return
  end subroutine rk_advance_trcvar_low_storage1d
  
 
 
!OCL SERIAL
  subroutine rk_advance_low_storage2d( this, nowstage, q, varID, is, ie ,js, je  )
    use scale_const, only: &
      eps => const_eps
    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 
 
    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
 
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_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)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i,j)
      !$omp do
      do j=js, je
      do i=is, ie
        q(i,j) =  this%varTmp_2d(i,j,varid)                                             &
                + sig_ss * q(i,j) + gam_ss * this%tend_buf2D_ex(i,j,varid,1)
      end do
      end do
      !$omp end parallel
      call prof_rapend( 'rk_advance_low_storage2D', 3)
 
      return
    end if
      
    !$omp parallel private(i,j)
    if (nowstage == 1) then
      !$omp do
      do j=js, je
      do i=is, ie
        this%var0_2D(i,j,varid)   = q(i,j)
        this%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
      do j=js, je
      do i=is, ie
        this%varTmp_2d(i,j,varid) = this%varTmp_2d(i,j,varid)      &
            + sig_ns * q(i,j) + gam_ns * this%tend_buf2D_ex(i,j,varid,1)
      end do
      end do
    end if
 
    !$omp do
    do j=js, je
    do i=is, ie
      q(i,j) = one_minus_sig_ss * this%var0_2d(i,j,varid)                &
              + sig_ss * q(i,j) + gam_ss * this%tend_buf2D_ex(i,j,varid,1)
    end do
    end do
 
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_low_storage2D', 3)
 
    return
  end subroutine rk_advance_low_storage2d
 
  !OCL SERIAL
  subroutine rk_advance_trcvar_low_storage2d( this, nowstage, q, varID, is, ie ,js, je , &
    DDENS, DDENS0, DENS_hyd )
    use scale_const, only: &
      eps => const_eps
    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(:,:)
 
    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)
      !$omp do
      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) =  ( this%varTmp_2d(i,j,varid)                                              &
                + sig_ss * q(i,j) * dens_ssm1 + gam_ss * this%tend_buf2D_ex(i,j,varid,1) ) &
                / ( dens_hyd(i,j) + ddens(i,j) )
      end do
      end do
      !$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)
    if (nowstage == 1) then
      !$omp do
      do j=js, je
      do i=is, ie
        this%var0_2D(i,j,varid)   = q(i,j) * ( dens_hyd(i,j) + ddens0(i,j) )
        this%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
      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) )
        this%varTmp_2d(i,j,varid) = this%varTmp_2d(i,j,varid)                 &
            + sig_ns * q(i,j) * dens_ssm1 + gam_ns * this%tend_buf2D_ex(i,j,varid,1)
      end do
      end do
    end if
 
    !$omp do
    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 * this%var0_2d(i,j,varid)                                &
              + sig_ss * q(i,j) * dens_ssm1 + gam_ss * this%tend_buf2D_ex(i,j,varid,1) )   &
              / dens_ss
    end do
    end do
 
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_trcvar_low_storage2D', 3)
 
    return
  end subroutine rk_advance_trcvar_low_storage2d
  
 
 
!OCL SERIAL
  subroutine rk_advance_low_storage3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke  )
    use scale_const, only: &
      eps => const_eps
    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 
 
    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
 
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_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)
 
    if ( nowstage == this%nstage ) then
      !$omp parallel private(i,j,k)
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) =  this%varTmp_3d(i,j,k,varid)                                             &
                + sig_ss * q(i,j,k) + gam_ss * this%tend_buf3D_ex(i,j,k,varid,1)
      end do
      end do
      end do
      !$omp end parallel
      call prof_rapend( 'rk_advance_low_storage3D', 3)
 
      return
    end if
      
    !$omp parallel private(i,j,k)
    if (nowstage == 1) then
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        this%var0_3D(i,j,k,varid)   = q(i,j,k)
        this%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)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        this%varTmp_3d(i,j,k,varid) = this%varTmp_3d(i,j,k,varid)      &
            + sig_ns * q(i,j,k) + gam_ns * this%tend_buf3D_ex(i,j,k,varid,1)
      end do
      end do
      end do
    end if
 
    !$omp do collapse(2)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      q(i,j,k) = one_minus_sig_ss * this%var0_3d(i,j,k,varid)                &
              + sig_ss * q(i,j,k) + gam_ss * this%tend_buf3D_ex(i,j,k,varid,1)
    end do
    end do
    end do
 
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_low_storage3D', 3)
 
    return
  end subroutine rk_advance_low_storage3d
 
  !OCL SERIAL
  subroutine rk_advance_trcvar_low_storage3d( this, nowstage, q, varID, is, ie ,js, je ,ks, ke , &
    DDENS, DDENS0, DENS_hyd )
    use scale_const, only: &
      eps => const_eps
    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(:,:,:)
 
    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)
      !$omp do collapse(2)
      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) =  ( this%varTmp_3d(i,j,k,varid)                                              &
                + sig_ss * q(i,j,k) * dens_ssm1 + gam_ss * this%tend_buf3D_ex(i,j,k,varid,1) ) &
                / ( dens_hyd(i,j,k) + ddens(i,j,k) )
      end do
      end do
      end do
      !$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)
    if (nowstage == 1) then
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        this%var0_3D(i,j,k,varid)   = q(i,j,k) * ( dens_hyd(i,j,k) + ddens0(i,j,k) )
        this%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)
      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) )
        this%varTmp_3d(i,j,k,varid) = this%varTmp_3d(i,j,k,varid)                 &
            + sig_ns * q(i,j,k) * dens_ssm1 + gam_ns * this%tend_buf3D_ex(i,j,k,varid,1)
      end do
      end do
      end do
    end if
 
    !$omp do collapse(2)
    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 * this%var0_3d(i,j,k,varid)                                &
              + sig_ss * q(i,j,k) * dens_ssm1 + gam_ss * this%tend_buf3D_ex(i,j,k,varid,1) )   &
              / dens_ss
    end do
    end do
    end do
 
    !$omp end parallel
 
    call prof_rapend( 'rk_advance_trcvar_low_storage3D', 3)
 
    return
  end subroutine rk_advance_trcvar_low_storage3d
  
 
 
!OCL SERIAL
  subroutine rk_advance_general1d( 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 
 
    integer :: i
    integer :: s
    integer :: tintbuf_ind
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_general1D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( this%nstage == 1 ) then
      !$omp parallel do
      do i=is, ie
        this%varTmp_1d(i,varid) = q(i)
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      if ( this%imex_flag ) then
 
      !$omp parallel do
        do i=is, ie
          q(i) =  this%varTmp_1d(i,varid) + this%dt * ( &
              + this%coef_b_ex(nowstage) * this%tend_buf1D_ex(i,varid,tintbuf_ind)   & 
              + this%coef_b_im(nowstage) * this%tend_buf1D_im(i,varid,tintbuf_ind)   )
        end do
      else
 
        !$omp parallel do
        do i=is, ie
          q(i) =  this%varTmp_1d(i,varid)                                             &
                  + this%dt * this%coef_b_ex(nowstage) * this%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 ) &
    !$omp private( i )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      do i=is, ie
        this%var0_1D(i,varid) = q(i)
        this%varTmp_1D(i,varid) = q(i)
      end do
    end if 
 
    if ( this%imex_flag ) then
      !$omp do
      do i=is, ie
        q(i) = this%var0_1d(i,varid)
        this%varTmp_1d(i,varid) =  this%varTmp_1d(i,varid) + this%dt * ( &
              + this%coef_b_ex(nowstage) * this%tend_buf1D_ex(i,varid,tintbuf_ind)   & 
              + this%coef_b_im(nowstage) * this%tend_buf1D_im(i,varid,tintbuf_ind)   )
      end do
    else
      !$omp do
      do i=is, ie
        q(i) = this%var0_1d(i,varid)
        this%varTmp_1d(i,varid) =  this%varTmp_1d(i,varid)                                             &
                  + this%dt * this%coef_b_ex(nowstage) * this%tend_buf1D_ex(i,varid,tintbuf_ind)
      end do
    end if 
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      do i=is, ie
        q(i) = this%var0_1d(i,varid)                                                &
            +  this%dt * this%coef_a_ex(nowstage+1,nowstage)*this%tend_buf1D_ex(i,varid,1)
      end do
    else if ( .not. this%imex_flag ) then      
      do s=1, nowstage
      !$omp do
      do i=is, ie
        q(i) = q(i)                                  &
            +  this%dt * this%coef_a_ex(nowstage+1,s)*this%tend_buf1D_ex(i,varid,s)
      end do
      end do
    else ! IMEX   
      do s=1, nowstage
      !$omp do
      do i=is, ie
        q(i) = q(i)                                        &
            +  this%dt * ( this%coef_a_ex(nowstage+1,s)*this%tend_buf1D_ex(i,varid,s)  &
                         + this%coef_a_im(nowstage+1,s)*this%tend_buf1D_im(i,varid,s)  )                         
      end do
      end do
    end if
 
    !$omp end parallel
    call prof_rapend( 'rk_advance_general1D', 3)
 
    return
  end subroutine rk_advance_general1d
 
  !OCL SERIAL
  subroutine rk_advance_trcvar_general1d( 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(:)
 
    integer :: i
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: dens_
    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
      do i=is, ie
        this%varTmp_1d(i,varid) = q(i) * ( dens_hyd(i) + ddens0(i) )
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      !$omp parallel do
      do i=is, ie
        q(i) =  ( this%varTmp_1d(i,varid)                                          &
            + this%dt * this%coef_b_ex(nowstage) * this%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_ ) &
    !$omp private( i )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      do i=is, ie
        this%var0_1D(i,varid) = q(i) * ( dens_hyd(i) + ddens0(i) )
        this%varTmp_1D(i,varid) = this%var0_1D(i,varid)
      end do
    end if 
 
    !$omp do
    do i=is, ie
      q(i) = this%var0_1d(i,varid)
      this%varTmp_1d(i,varid) =  this%varTmp_1d(i,varid)                                             &
            + this%dt * this%coef_b_ex(nowstage) * this%tend_buf1D_ex(i,varid,tintbuf_ind)
    end do
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      do i=is, ie
        dens_ = ( dens_hyd(i) + ddens0(i) ) &
            + this%coef_a_ex(nowstage+1,nowstage)*( ddens(i) - ddens0(i) )
        q(i) = ( this%var0_1d(i,varid)                                           &
            + this%dt * this%coef_a_ex(nowstage+1,nowstage)*this%tend_buf1D_ex(i,varid,1) ) &
            / dens_
      end do
    else if ( .not. this%imex_flag ) then      
      do s=1, nowstage
      !$omp do
      do i=is, ie
        q(i) = q(i)                                  &
            + this%dt * this%coef_a_ex(nowstage+1,s)*this%tend_buf1D_ex(i,varid,s)
      end do
      end do
      !$omp do
      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
 
    !$omp end parallel
    call prof_rapend( 'rk_advance_trcvar_generall1D', 3)
 
    return
  end subroutine rk_advance_trcvar_general1d
 
!OCL SERIAL
  subroutine rk_storeimpl_general1d( 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 
 
    integer :: i
    integer :: tintbuf_ind
    !----------------------------------------    
 
    if (.not. this%imex_flag ) return
    
    tintbuf_ind = this%tend_buf_indmap(nowstage)   
 
    !$omp parallel
 
    if ( nowstage == 1 ) then
      !$omp do
      do i=is, ie
        this%var0_1D(i,varid) = q(i)
        this%varTmp_1D(i,varid) = q(i)
      end do
    end if
            
    !$omp do
    do i=is, ie
      q(i) = q(i)                                                                              &
            +  this%dt * this%coef_a_im(nowstage,nowstage)*this%tend_buf1D_im(i,varid,tintbuf_ind)  
    end do
    !$omp end parallel
 
    return
  end subroutine rk_storeimpl_general1d
 
 
!OCL SERIAL
  subroutine rk_advance_general2d( 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 
 
    integer :: i,j
    integer :: s
    integer :: tintbuf_ind
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_general2D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( this%nstage == 1 ) then
      !$omp parallel do
      do j=js, je
      do i=is, ie
        this%varTmp_2d(i,j,varid) = q(i,j)
      end do
      end do
    end if
 
    if ( nowstage ==  this%nstage ) then
      if ( this%imex_flag ) then
 
      !$omp parallel do
        do j=js, je
        do i=is, ie
          q(i,j) =  this%varTmp_2d(i,j,varid) + this%dt * ( &
              + this%coef_b_ex(nowstage) * this%tend_buf2D_ex(i,j,varid,tintbuf_ind)   & 
              + this%coef_b_im(nowstage) * this%tend_buf2D_im(i,j,varid,tintbuf_ind)   )
        end do
        end do
      else
 
        !$omp parallel do
        do j=js, je
        do i=is, ie
          q(i,j) =  this%varTmp_2d(i,j,varid)                                             &
                  + this%dt * this%coef_b_ex(nowstage) * this%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 ) &
    !$omp private( i,j )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      do j=js, je
      do i=is, ie
        this%var0_2D(i,j,varid) = q(i,j)
        this%varTmp_2D(i,j,varid) = q(i,j)
      end do
      end do
    end if 
 
    if ( this%imex_flag ) then
      !$omp do
      do j=js, je
      do i=is, ie
        q(i,j) = this%var0_2d(i,j,varid)
        this%varTmp_2d(i,j,varid) =  this%varTmp_2d(i,j,varid) + this%dt * ( &
              + this%coef_b_ex(nowstage) * this%tend_buf2D_ex(i,j,varid,tintbuf_ind)   & 
              + this%coef_b_im(nowstage) * this%tend_buf2D_im(i,j,varid,tintbuf_ind)   )
      end do
      end do
    else
      !$omp do
      do j=js, je
      do i=is, ie
        q(i,j) = this%var0_2d(i,j,varid)
        this%varTmp_2d(i,j,varid) =  this%varTmp_2d(i,j,varid)                                             &
                  + this%dt * this%coef_b_ex(nowstage) * this%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
      !$omp do
      do j=js, je
      do i=is, ie
        q(i,j) = this%var0_2d(i,j,varid)                                                &
            +  this%dt * this%coef_a_ex(nowstage+1,nowstage)*this%tend_buf2D_ex(i,j,varid,1)
      end do
      end do
    else if ( .not. this%imex_flag ) then      
      do s=1, nowstage
      !$omp do
      do j=js, je
      do i=is, ie
        q(i,j) = q(i,j)                                  &
            +  this%dt * this%coef_a_ex(nowstage+1,s)*this%tend_buf2D_ex(i,j,varid,s)
      end do
      end do
      end do
    else ! IMEX   
      do s=1, nowstage
      !$omp do
      do j=js, je
      do i=is, ie
        q(i,j) = q(i,j)                                        &
            +  this%dt * ( this%coef_a_ex(nowstage+1,s)*this%tend_buf2D_ex(i,j,varid,s)  &
                         + this%coef_a_im(nowstage+1,s)*this%tend_buf2D_im(i,j,varid,s)  )                         
      end do
      end do
      end do
    end if
 
    !$omp end parallel
    call prof_rapend( 'rk_advance_general2D', 3)
 
    return
  end subroutine rk_advance_general2d
 
  !OCL SERIAL
  subroutine rk_advance_trcvar_general2d( 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(:,:)
 
    integer :: i,j
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: dens_
    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
      do j=js, je
      do i=is, ie
        this%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
      !$omp parallel do
      do j=js, je
      do i=is, ie
        q(i,j) =  ( this%varTmp_2d(i,j,varid)                                          &
            + this%dt * this%coef_b_ex(nowstage) * this%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_ ) &
    !$omp private( i,j )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      do j=js, je
      do i=is, ie
        this%var0_2D(i,j,varid) = q(i,j) * ( dens_hyd(i,j) + ddens0(i,j) )
        this%varTmp_2D(i,j,varid) = this%var0_2D(i,j,varid)
      end do
      end do
    end if 
 
    !$omp do
    do j=js, je
    do i=is, ie
      q(i,j) = this%var0_2d(i,j,varid)
      this%varTmp_2d(i,j,varid) =  this%varTmp_2d(i,j,varid)                                             &
            + this%dt * this%coef_b_ex(nowstage) * this%tend_buf2D_ex(i,j,varid,tintbuf_ind)
    end do
    end do
 
    if ( this%tend_buf_size == 1 .and. (.not. this%imex_flag) ) then
      !$omp do
      do j=js, je
      do i=is, ie
        dens_ = ( dens_hyd(i,j) + ddens0(i,j) ) &
            + this%coef_a_ex(nowstage+1,nowstage)*( ddens(i,j) - ddens0(i,j) )
        q(i,j) = ( this%var0_2d(i,j,varid)                                           &
            + this%dt * this%coef_a_ex(nowstage+1,nowstage)*this%tend_buf2D_ex(i,j,varid,1) ) &
            / dens_
      end do
      end do
    else if ( .not. this%imex_flag ) then      
      do s=1, nowstage
      !$omp do
      do j=js, je
      do i=is, ie
        q(i,j) = q(i,j)                                  &
            + this%dt * this%coef_a_ex(nowstage+1,s)*this%tend_buf2D_ex(i,j,varid,s)
      end do
      end do
      end do
      !$omp do
      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
 
    !$omp end parallel
    call prof_rapend( 'rk_advance_trcvar_generall2D', 3)
 
    return
  end subroutine rk_advance_trcvar_general2d
 
!OCL SERIAL
  subroutine rk_storeimpl_general2d( 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 
 
    integer :: i,j
    integer :: tintbuf_ind
    !----------------------------------------    
 
    if (.not. this%imex_flag ) return
    
    tintbuf_ind = this%tend_buf_indmap(nowstage)   
 
    !$omp parallel
 
    if ( nowstage == 1 ) then
      !$omp do
      do j=js, je
      do i=is, ie
        this%var0_2D(i,j,varid) = q(i,j)
        this%varTmp_2D(i,j,varid) = q(i,j)
      end do
      end do
    end if
            
    !$omp do
    do j=js, je
    do i=is, ie
      q(i,j) = q(i,j)                                                                              &
            +  this%dt * this%coef_a_im(nowstage,nowstage)*this%tend_buf2D_im(i,j,varid,tintbuf_ind)  
    end do
    end do
    !$omp end parallel
 
    return
  end subroutine rk_storeimpl_general2d
 
 
!OCL SERIAL
  subroutine rk_advance_general3d( 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 
 
    integer :: i,j,k
    integer :: s
    integer :: tintbuf_ind
    !----------------------------------------    
 
    call prof_rapstart( 'rk_advance_general3D', 3)
 
    tintbuf_ind = this%tend_buf_indmap(nowstage)
 
    if ( this%nstage == 1 ) then
      !$omp parallel do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        this%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
 
      !$omp parallel do collapse(2)
        do k=ks, ke
        do j=js, je
        do i=is, ie
          q(i,j,k) =  this%varTmp_3d(i,j,k,varid) + this%dt * ( &
              + this%coef_b_ex(nowstage) * this%tend_buf3D_ex(i,j,k,varid,tintbuf_ind)   & 
              + this%coef_b_im(nowstage) * this%tend_buf3D_im(i,j,k,varid,tintbuf_ind)   )
        end do
        end do
        end do
      else
 
        !$omp parallel do collapse(2)
        do k=ks, ke
        do j=js, je
        do i=is, ie
          q(i,j,k) =  this%varTmp_3d(i,j,k,varid)                                             &
                  + this%dt * this%coef_b_ex(nowstage) * this%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 ) &
    !$omp private( i,j,k )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        this%var0_3D(i,j,k,varid) = q(i,j,k)
        this%varTmp_3D(i,j,k,varid) = q(i,j,k)
      end do
      end do
      end do
    end if 
 
    if ( this%imex_flag ) then
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = this%var0_3d(i,j,k,varid)
        this%varTmp_3d(i,j,k,varid) =  this%varTmp_3d(i,j,k,varid) + this%dt * ( &
              + this%coef_b_ex(nowstage) * this%tend_buf3D_ex(i,j,k,varid,tintbuf_ind)   & 
              + this%coef_b_im(nowstage) * this%tend_buf3D_im(i,j,k,varid,tintbuf_ind)   )
      end do
      end do
      end do
    else
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = this%var0_3d(i,j,k,varid)
        this%varTmp_3d(i,j,k,varid) =  this%varTmp_3d(i,j,k,varid)                                             &
                  + this%dt * this%coef_b_ex(nowstage) * this%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
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = this%var0_3d(i,j,k,varid)                                                &
            +  this%dt * this%coef_a_ex(nowstage+1,nowstage)*this%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
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = q(i,j,k)                                  &
            +  this%dt * this%coef_a_ex(nowstage+1,s)*this%tend_buf3D_ex(i,j,k,varid,s)
      end do
      end do
      end do
      end do
    else ! IMEX   
      do s=1, nowstage
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = q(i,j,k)                                        &
            +  this%dt * ( this%coef_a_ex(nowstage+1,s)*this%tend_buf3D_ex(i,j,k,varid,s)  &
                         + this%coef_a_im(nowstage+1,s)*this%tend_buf3D_im(i,j,k,varid,s)  )                         
      end do
      end do
      end do
      end do
    end if
 
    !$omp end parallel
    call prof_rapend( 'rk_advance_general3D', 3)
 
    return
  end subroutine rk_advance_general3d
 
  !OCL SERIAL
  subroutine rk_advance_trcvar_general3d( 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(:,:,:)
 
    integer :: i,j,k
    integer :: s
    integer :: tintbuf_ind
 
    real(RP) :: dens_
    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)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        this%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
      !$omp parallel do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) =  ( this%varTmp_3d(i,j,k,varid)                                          &
            + this%dt * this%coef_b_ex(nowstage) * this%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_ ) &
    !$omp private( i,j,k )
 
    if ( nowstage == 1 .and. (.not. this%imex_flag) ) then
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        this%var0_3D(i,j,k,varid) = q(i,j,k) * ( dens_hyd(i,j,k) + ddens0(i,j,k) )
        this%varTmp_3D(i,j,k,varid) = this%var0_3D(i,j,k,varid)
      end do
      end do
      end do
    end if 
 
    !$omp do collapse(2)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      q(i,j,k) = this%var0_3d(i,j,k,varid)
      this%varTmp_3d(i,j,k,varid) =  this%varTmp_3d(i,j,k,varid)                                             &
            + this%dt * this%coef_b_ex(nowstage) * this%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
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        dens_ = ( dens_hyd(i,j,k) + ddens0(i,j,k) ) &
            + this%coef_a_ex(nowstage+1,nowstage)*( ddens(i,j,k) - ddens0(i,j,k) )
        q(i,j,k) = ( this%var0_3d(i,j,k,varid)                                           &
            + this%dt * this%coef_a_ex(nowstage+1,nowstage)*this%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
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        q(i,j,k) = q(i,j,k)                                  &
            + this%dt * this%coef_a_ex(nowstage+1,s)*this%tend_buf3D_ex(i,j,k,varid,s)
      end do
      end do
      end do
      end do
      !$omp do collapse(2)
      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
 
    !$omp end parallel
    call prof_rapend( 'rk_advance_trcvar_generall3D', 3)
 
    return
  end subroutine rk_advance_trcvar_general3d
 
!OCL SERIAL
  subroutine rk_storeimpl_general3d( 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 
 
    integer :: i,j,k
    integer :: tintbuf_ind
    !----------------------------------------    
 
    if (.not. this%imex_flag ) return
    
    tintbuf_ind = this%tend_buf_indmap(nowstage)   
 
    !$omp parallel
 
    if ( nowstage == 1 ) then
      !$omp do collapse(2)
      do k=ks, ke
      do j=js, je
      do i=is, ie
        this%var0_3D(i,j,k,varid) = q(i,j,k)
        this%varTmp_3D(i,j,k,varid) = q(i,j,k)
      end do
      end do
      end do
    end if
            
    !$omp do collapse(2)
    do k=ks, ke
    do j=js, je
    do i=is, ie
      q(i,j,k) = q(i,j,k)                                                                              &
            +  this%dt * this%coef_a_im(nowstage,nowstage)*this%tend_buf3D_im(i,j,k,varid,tintbuf_ind)  
    end do
    end do
    end do
    !$omp end parallel
 
    return
  end subroutine rk_storeimpl_general3d
 
end module scale_timeint_rk