scale_atmos_phy_bl_common.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_phy_bl_common
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
 
#if defined DEBUG || defined QUICKDEBUG
  use scale_debug, only: &
     check
#endif
  use scale_const, only: &
     undef  => const_undef, &
     iundef => const_undef2
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_bl_tendency_tracer
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_bl_tendency_tracer( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       DENS, QTRC, SFLX_Q, &
       Kh, MASS,           &
       CZ, FZ, F2H, DDT,   &
       TRACER_NAME,        &
       RHOQ_t              )
    use scale_matrix, only: &
       matrix_solver_tridiagonal
    use scale_file_history, only: &
       file_history_in
    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je
 
    real(rp),         intent(in) :: dens  (ka,ia,ja)
    real(rp),         intent(in) :: qtrc  (ka,ia,ja)
    real(rp),         intent(in) :: sflx_q(   ia,ja)
    real(rp),         intent(in) :: kh    (ka,ia,ja)
    real(rp),         intent(in) :: mass
    real(rp),         intent(in) :: cz (  ka,ia,ja)
    real(rp),         intent(in) :: fz (0:ka,ia,ja)
    real(rp),         intent(in) :: f2h(ka,2,ia,ja)
    real(dp),         intent(in) :: ddt
    character(len=*), intent(in) :: tracer_name
 
    real(rp), intent(out) :: rhoq_t(ka,ia,ja)
 
    real(rp) :: qtrc_n(ka)
    real(rp) :: rho   (ka)
    real(rp) :: rhokh (ka)
    real(rp) :: a(ka)
    real(rp) :: b(ka)
    real(rp) :: c(ka)
    real(rp) :: d(ka)
    real(rp) :: rho_h
    real(rp) :: ap
    real(rp) :: sf_t
 
    real(rp) :: flx(ka,ia,ja)
 
    real(rp) :: cdz(ka)
    real(rp) :: fdz(ka)
 
    real(rp) :: dt
 
    integer :: ke_pbl
    integer :: k, i, j
 
#ifdef _OPENACC
    real(rp) :: work(ks:ke,4)
#endif
    dt = real( ddt, kind=rp )
 
    !$acc data create(flx)
 
!OCL INDEPENDENT
    !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
    !$omp shared(KA,KS,KE,IS,IE,JS,JE) &
    !$omp shared(RHOQ_t,DENS,QTRC,SFLX_Q,Kh,MASS,CZ,FZ,F2H,DT,flx) &
    !$omp private(QTRC_n,RHO,RHOKh,rho_h,a,b,c,d,ap,sf_t,CDZ,FDZ) &
    !$omp private(KE_PBL,k,i,j)
    !$acc kernels
    !$acc loop independent collapse(2) &
    !$acc private(QTRC_n,RHO,RHOKh,rho_h,a,b,c,d,ap,sf_t,CDZ,FDZ, &
    !$acc         KE_PBL, &
    !$acc         work)
    do j = js, je
    do i = is, ie
 
       ke_pbl = ke-1
       !$acc loop seq
       do k = ke-2, ks+1, -1
          if ( kh(k,i,j) > 0.0_rp ) then
             ke_pbl = k + 1
          else
             exit
          end if
       end do
 
       do k = ks, ke_pbl
          cdz(k) = fz(k  ,i,j) - fz(k-1,i,j)
          fdz(k) = cz(k+1,i,j) - cz(k  ,i,j)
       end do
 
       sf_t = sflx_q(i,j) / cdz(ks)
 
       rho(ks) = dens(ks,i,j) + dt * sf_t * mass
       do k = ks+1, ke_pbl
          rho(k) = dens(k,i,j)
       end do
 
       ! dens * coefficient at the half level
       do k = ks, ke_pbl-1
          rho_h = f2h(k,1,i,j) * dens(k+1,i,j) + f2h(k,2,i,j) * dens(k,i,j)
          rhokh(k) = rho_h * kh(k,i,j)
       end do
 
       d(ks) = ( qtrc(ks,i,j) * dens(ks,i,j) + dt * sf_t ) / rho(ks)
       do k = ks+1, ke_pbl
          d(k) = qtrc(k,i,j)
       end do
 
       c(ks) = 0.0_rp
       do k = ks, ke_pbl-1
          ap = - dt * rhokh(k) / fdz(k)
          a(k) = ap / ( rho(k) * cdz(k) )
#ifdef _OPENACC
          if ( k==ks ) then
             b(k) = - a(k) + 1.0_rp
          else
             b(k) = - a(k) + dt * rhokh(k-1) / ( fdz(k-1) * rho(k) * cdz(k) ) + 1.0_rp
          end if
#else
          b(k) = - a(k) - c(k) + 1.0_rp
#endif
          c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
       end do
       a(ke_pbl) = 0.0_rp
       b(ke_pbl) = - c(ke_pbl) + 1.0_rp
 
       call matrix_solver_tridiagonal( &
               ka, ks, ke_pbl, &
#ifdef _OPENACC
               work(:,:), & ! (wrok)
#endif
               a(:), b(:), c(:), d(:), & ! (in)
               qtrc_n(:)               ) ! (out)
 
       rhoq_t(ks,i,j) = ( qtrc_n(ks) * rho(ks) - qtrc(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
       do k = ks+1, ke_pbl
          rhoq_t(k,i,j) = ( qtrc_n(k) - qtrc(k,i,j) ) * rho(k) / dt
       end do
       do k = ke_pbl+1, ke
          rhoq_t(k,i,j) = 0.0_rp
       end do
 
       do k = ks, ke_pbl-1
          flx(k,i,j) = - rhokh(k) * ( qtrc_n(k+1) - qtrc_n(k) ) / fdz(k)
       end do
 
    end do
    end do
    !$acc end kernels
 
    call file_history_in(flx(:,:,:), 'ZFLX_'//trim(tracer_name)//'_BL', 'Z FLUX of DENS * '//trim(tracer_name)//' (PBL)', 'kg/m2/s', fill_halo=.true.)
 
    !$acc end data
 
    return
  end subroutine atmos_phy_bl_tendency_tracer
 
end module scale_atmos_phy_bl_common