scale_atmos_dyn_common.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_dyn_common
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  use scale_index
  use scale_tracer
 
#ifdef DEBUG
  use scale_debug, only: &
     check
  use scale_const, only: &
     undef  => const_undef, &
     iundef => const_undef2
#endif
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_dyn_wdamp_setup
  public :: atmos_dyn_fill_halo
  public :: atmos_dyn_copy_boundary
  public :: atmos_dyn_copy_boundary_tracer
  public :: atmos_dyn_divergence
  public :: atmos_dyn_prep_pres_linearization
 
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
 
contains
 
  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_wdamp_setup( &
       wdamp_coef,              &
       wdamp_tau, wdamp_height, &
       FZ                       )
    use scale_const, only: &
       pi  => const_pi, &
       eps => const_eps
    implicit none
 
    real(rp), intent(out) :: wdamp_coef(ka)
    real(rp), intent(in)  :: wdamp_tau
    real(rp), intent(in)  :: wdamp_height
    real(rp), intent(in)  :: fz(0:ka)
 
    real(rp) :: alpha, sw
 
    integer :: k
    !---------------------------------------------------------------------------
 
    !$acc data copyout(wdamp_coef) copyin(FZ)
 
    if ( wdamp_height < 0.0_rp ) then
       !$acc kernels
       wdamp_coef(:) = 0.0_rp
       !$acc end kernels
    elseif( fz(ke)-wdamp_height < eps ) then
       !$acc kernels
       wdamp_coef(:) = 0.0_rp
       !$acc end kernels
    else
       alpha = 1.0_rp / wdamp_tau
 
       !$acc kernels
       do k = ks, ke
          sw = 0.5_rp + sign( 0.5_rp, fz(k)-wdamp_height )
 
          wdamp_coef(k) = alpha * sw &
                        * 0.5_rp * ( 1.0_rp - cos( pi * (fz(k)-wdamp_height) / (fz(ke)-wdamp_height)) )
       enddo
       !$acc end kernels
       !$acc kernels
       wdamp_coef(   1:ks-1) = wdamp_coef(ks)
       !$acc end kernels
       !$acc kernels
       wdamp_coef(ke+1:ka  ) = wdamp_coef(ke)
       !$acc end kernels
 
       !$acc update host(wdamp_coef)
       log_newline
       log_info("ATMOS_DYN_wdamp_setup",*)                          'Setup Rayleigh damping coefficient'
       log_info_cont('(1x,A)')                   '|=== Rayleigh Damping Coef ===|'
       log_info_cont('(1x,A)')                   '|     k     zh[m]    coef[/s] |'
       do k = ka, ke+1, -1
       log_info_cont('(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' |'
       enddo
       k = ke
       log_info_cont('(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' | KE = TOA'
       do k = ke-1, ks, -1
       log_info_cont('(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' |'
       enddo
       k = ks-1
       log_info_cont('(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' | KS-1 = surface'
       do k = ks-2, 1, -1
       log_info_cont('(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' |'
       enddo
       k = 0
       log_info_cont('(1x,A,I5,F10.2,12x,A)')    '| ',k, fz(k),               ' |'
       log_info_cont('(1x,A)')                   '|=============================|'
    endif
 
    !$acc end data
 
    return
  end subroutine atmos_dyn_wdamp_setup
 
  !-----------------------------------------------------------------------------
 
  subroutine atmos_dyn_fill_halo( var,             &
      fill_constval, lateral_halo, top_bottom_halo )
   implicit none
 
   real(rp), intent(inout) :: var(ka,ia,ja)
   real(rp), intent(in) :: fill_constval
   logical, intent(in) :: lateral_halo
   logical, intent(in) :: top_bottom_halo
 
   integer :: i, j, k
   !----------------------------
 
   !$acc data copy(var)
 
   if (lateral_halo) then
      !$acc kernels
!OCL XFILL
      do j = 1, ja
      do i = 1, isb-1
      do k = 1, ka
         var(k,i,j) = fill_constval
      enddo
      enddo
      do i = ieb+1, ia
      do k = 1, ka
         var(k,i,j) = fill_constval
      enddo
      enddo
      enddo
      !$acc end kernels
      !$acc kernels
!OCL XFILL
      do j = 1, jsb-1
      do i = 1, ia
      do k = 1, ka
         var(k,i,j) = fill_constval
      enddo
      enddo
      enddo
      !$acc end kernels
      !$acc kernels
!OCL XFILL
      do j = jeb+1, ja
      do i = 1, ia
      do k = 1, ka
         var(k,i,j) = fill_constval
      enddo
      enddo
      enddo
      !$acc end kernels
   end if
 
   if (top_bottom_halo) then
      !$acc kernels
!OCL XFILL
      do j = js, je
      do i = is, ie
         var(   1:ks-1,i,j) = fill_constval
         var(ke+1:ka  ,i,j) = fill_constval
      enddo
      enddo
      !$acc end kernels
   end if
 
   !$acc end data
 
   return
  end subroutine atmos_dyn_fill_halo
 
  subroutine atmos_dyn_copy_boundary( &
       DENS,  MOMZ,  MOMX,  MOMY,  RHOT,  PROG, &
       DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, PROG0, &
       BND_W, BND_E, BND_S, BND_N, TwoD )
    implicit none
    real(rp), intent(inout) :: dens (ka,ia,ja)
    real(rp), intent(inout) :: momz (ka,ia,ja)
    real(rp), intent(inout) :: momx (ka,ia,ja)
    real(rp), intent(inout) :: momy (ka,ia,ja)
    real(rp), intent(inout) :: rhot (ka,ia,ja)
    real(rp), intent(inout) :: prog (ka,ia,ja,va)
    real(rp), intent(in)    :: dens0(ka,ia,ja)
    real(rp), intent(in)    :: momz0(ka,ia,ja)
    real(rp), intent(in)    :: momx0(ka,ia,ja)
    real(rp), intent(in)    :: momy0(ka,ia,ja)
    real(rp), intent(in)    :: rhot0(ka,ia,ja)
    real(rp), intent(in)    :: prog0(ka,ia,ja,va)
    logical,  intent(in)    :: bnd_w
    logical,  intent(in)    :: bnd_e
    logical,  intent(in)    :: bnd_s
    logical,  intent(in)    :: bnd_n
    logical,  intent(in)    :: twod
 
    integer :: k, i, j, iv
 
    !$acc data copy(DENS, MOMZ, MOMX, MOMY, RHOT, PROG) &
    !$acc      copyin(DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, PROG0)
 
    if ( bnd_w .and. (.not. twod) ) then
       !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,iv) &
       !$omp shared(JA,IS,KS,KE,DENS,DENS0,MOMZ,MOMZ0,MOMX,MOMX0,MOMY,MOMY0,RHOT,RHOT0,VA,PROG,PROG0)
       !$acc kernels
!OCL XFILL
       do j = 1, ja
       do i = 1, is-1
       do k = ks, ke
          dens(k,i,j) = dens0(k,i,j)
          momz(k,i,j) = momz0(k,i,j)
          momx(k,i,j) = momx0(k,i,j)
          momy(k,i,j) = momy0(k,i,j)
          rhot(k,i,j) = rhot0(k,i,j)
          !$acc loop seq
          do iv = 1, va
             prog(k,i,j,iv) = prog0(k,i,j,iv)
          end do
       enddo
       enddo
       enddo
       !$acc end kernels
    end if
    if ( bnd_e .and. (.not. twod) ) then
       !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,iv) &
       !$omp shared(JA,IE,IA,KS,KE,DENS,DENS0,MOMZ,MOMZ0,MOMX,MOMX0,MOMY,MOMY0,RHOT,RHOT0,VA,PROG,PROG0)
       !$acc kernels
!OCL XFILL
       do j = 1, ja
       do i = ie+1, ia
       do k = ks, ke
          dens(k,i,j) = dens0(k,i,j)
          momz(k,i,j) = momz0(k,i,j)
          momx(k,i,j) = momx0(k,i,j)
          momy(k,i,j) = momy0(k,i,j)
          rhot(k,i,j) = rhot0(k,i,j)
          !$acc loop seq
          do iv = 1, va
             prog(k,i,j,iv) = prog0(k,i,j,iv)
          end do
       enddo
       enddo
       enddo
       !$acc end kernels
       !$omp parallel do private(j,k) OMP_SCHEDULE_ collapse(2)
       !$acc kernels
!OCL XFILL
       do j = 1, ja
       do k = ks, ke
          momx(k,ie,j) = momx0(k,ie,j)
       enddo
       enddo
       !$acc end kernels
    end if
    if ( bnd_s ) then
       !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,iv) &
       !$omp shared(JS,IA,KS,KE,DENS,DENS0,MOMZ,MOMZ0,MOMX,MOMX0,MOMY,MOMY0,RHOT,RHOT0,VA,PROG,PROG0)
       !$acc kernels
!OCL XFILL
       do j = 1, js-1
       do i = 1, ia
       do k = ks, ke
          dens(k,i,j) = dens0(k,i,j)
          momz(k,i,j) = momz0(k,i,j)
          momx(k,i,j) = momx0(k,i,j)
          momy(k,i,j) = momy0(k,i,j)
          rhot(k,i,j) = rhot0(k,i,j)
          !$acc loop seq
          do iv = 1, va
             prog(k,i,j,iv) = prog0(k,i,j,iv)
          end do
       enddo
       enddo
       enddo
       !$acc end kernels
    end if
    if ( bnd_n ) then
       !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp private(i,iv) &
       !$omp shared(JA,JE,IA,KS,KE,DENS,DENS0,MOMZ,MOMZ0,MOMX,MOMX0,MOMY,MOMY0,RHOT,RHOT0,VA,PROG,PROG0)
       !$acc kernels
!OCL XFILL
       do j = je+1, ja
       do i = 1, ia
       do k = ks, ke
          dens(k,i,j) = dens0(k,i,j)
          momz(k,i,j) = momz0(k,i,j)
          momx(k,i,j) = momx0(k,i,j)
          momy(k,i,j) = momy0(k,i,j)
          rhot(k,i,j) = rhot0(k,i,j)
          !$acc loop seq
          do iv = 1, va
             prog(k,i,j,iv) = prog0(k,i,j,iv)
          end do
       enddo
       enddo
       enddo
       !$acc end kernels
       !$omp parallel do private(i,k) OMP_SCHEDULE_ collapse(2)
       !$acc kernels
!OCL XFILL
       do i = 1, ia
       do k = ks, ke
          momy(k,i,je) = momy0(k,i,je)
       enddo
       enddo
       !$acc end kernels
    end if
 
    !$acc end data
 
    return
  end subroutine atmos_dyn_copy_boundary
 
  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_copy_boundary_tracer( &
       QTRC, QTRC0, &
       BND_W, BND_E, BND_S, BND_N, TwoD )
    implicit none
    real(rp), intent(inout) :: qtrc (ka,ia,ja)
    real(rp), intent(in)    :: qtrc0(ka,ia,ja)
    logical,  intent(in)    :: bnd_w
    logical,  intent(in)    :: bnd_e
    logical,  intent(in)    :: bnd_s
    logical,  intent(in)    :: bnd_n
    logical,  intent(in)    :: twod
 
    integer :: k, i, j
 
    !$acc data copy(QTRC) copyin(QTRC0)
 
    if ( bnd_w .and. (.not. twod) ) then
       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JA,IS,KS,KE,QTRC,QTRC0)
       !$acc kernels
!OCL XFILL
       do j = 1, ja
       do i = 1, is-1
       do k = ks, ke
          qtrc(k,i,j) = qtrc0(k,i,j)
       enddo
       enddo
       enddo
       !$acc end kernels
    end if
    if ( bnd_e .and. (.not. twod) ) then
       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JA,IE,IA,KS,KE,QTRC,QTRC0)
       !$acc kernels
!OCL XFILL
       do j = 1, ja
       do i = ie+1, ia
       do k = ks, ke
          qtrc(k,i,j) = qtrc0(k,i,j)
       enddo
       enddo
       enddo
       !$acc end kernels
    end if
    if ( bnd_s ) then
       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JS,IA,KS,KE,QTRC,QTRC0)
       !$acc kernels
!OCL XFILL
       do j = 1, js-1
       do i = 1, ia
       do k = ks, ke
          qtrc(k,i,j) = qtrc0(k,i,j)
       enddo
       enddo
       enddo
       !$acc end kernels
    end if
    if ( bnd_n ) then
       !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JA,JE,IA,KS,KE,QTRC,QTRC0)
       !$acc kernels
!OCL XFILL
       do j = je+1, ja
       do i = 1, ia
       do k = ks, ke
          qtrc(k,i,j) = qtrc0(k,i,j)
       enddo
       enddo
       enddo
       !$acc end kernels
    end if
 
    !$acc end data
 
    return
  end subroutine atmos_dyn_copy_boundary_tracer
 
  !-----------------------------------------------------------------------------
 
  subroutine atmos_dyn_divergence( &
       DDIV, &
       MOMZ, MOMX, MOMY, &
       GSQRT, J13G, J23G, J33G, MAPF, &
       TwoD, &
       RCDZ, RCDX, RCDY, RFDZ, FDZ )
    implicit none
    real(rp), intent(out) :: ddiv(ka,ia,ja)
    real(rp), intent(in)  :: momz(ka,ia,ja)
    real(rp), intent(in)  :: momx(ka,ia,ja)
    real(rp), intent(in)  :: momy(ka,ia,ja)
    real(rp), intent(in)  :: gsqrt(ka,ia,ja,7)
    real(rp), intent(in)  :: j13g(ka,ia,ja,7)
    real(rp), intent(in)  :: j23g(ka,ia,ja,7)
    real(rp), intent(in)  :: j33g
    real(rp), intent(in)  :: mapf(ia,ja,2,7)
    logical,  intent(in)  :: twod
    real(rp), intent(in)  :: rcdz(ka)
    real(rp), intent(in)  :: rcdx(ia)
    real(rp), intent(in)  :: rcdy(ja)
    real(rp), intent(in)  :: rfdz(ka-1)
    real(rp), intent(in)  :: fdz(ka-1)
 
    integer :: k, i, j
 
    call prof_rapstart("DYN_divercence", 2)
 
    !$acc data copyout(DDIV) copyin(MOMZ, MOMX, MOMY, GSQRT, J13G, J23G, MAPF, RCDZ, RCDX, RCDY, RFDZ, FDZ)
 
    ! 3D divergence
 
    if ( twod ) then
       !$omp parallel do private(j,k) OMP_SCHEDULE_
       !$acc kernels
       do j = js, je+1
       do k = ks-1, ke+1
          ddiv(k,is,j) = j33g * ( momz(k,is,j) - momz(k-1,is,j) ) * rcdz(k) &
                              + ( ( momy(k+1,is,j) + momy(k+1,is,j-1) ) * j23g(k+1,is,j,i_xyw) &
                                - ( momy(k-1,is,j) + momy(k-1,is,j-1) ) * j23g(k-1,is,j,i_xyw) ) / ( fdz(k)+fdz(k-1) ) &
                      + mapf(is,j,2,i_xy) &
                      * ( momy(k,is,j  ) * gsqrt(k,is,j  ,i_xvz) / mapf(is,j  ,1,i_xv) &
                        - momy(k,is,j-1) * gsqrt(k,is,j-1,i_xvz) / mapf(is,j-1,1,i_xv) ) * rcdy(j)
       enddo
       enddo
       !$acc end kernels
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif
       !$omp parallel do private(j) OMP_SCHEDULE_
       !$acc kernels
       do j = js, je+1
          ddiv(ks,is,j) = j33g * ( momz(ks,is,j) ) * rcdz(ks) &
                         + ( ( momy(ks+1,is,j) + momy(ks+1,is,j-1) ) * j23g(ks+1,is,j,i_xyw) &
                           - ( momy(ks  ,is,j) + momy(ks  ,is,j-1) ) * j23g(ks  ,is,j,i_xyw) ) * rfdz(ks) &
                       + mapf(is,j,2,i_xy) &
                       * ( momy(ks,is,j  ) * gsqrt(ks,is,j  ,i_xvz) / mapf(is,j  ,1,i_xv) &
                         - momy(ks,is,j-1) * gsqrt(ks,is,j-1,i_xvz) / mapf(is,j-1,1,i_xv) ) * rcdy(j)
          ddiv(ke,is,j) = j33g * ( - momz(ke-1,is,j  ) ) * rcdz(ke) &
                       + ( ( momy(ke  ,is,j) + momy(ke  ,is,j-1) ) * j23g(ke  ,is,j,i_xyw) &
                         - ( momy(ke-1,is,j) + momy(ke-1,is,j-1) ) * j23g(ke-1,is,j,i_xyw) ) * rfdz(ke) &
                       + mapf(is,j,2,i_xy) &
                       * ( momy(ke,is,j  ) * gsqrt(ke,is,j  ,i_xvz) / mapf(is,j  ,1,i_xv) &
                         - momy(ke,is,j-1) * gsqrt(ke,is,j-1,i_xvz) / mapf(is,j-1,1,i_xv) ) * rcdy(j)
       enddo
       !$acc end kernels
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif
    else
       !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
       !$acc kernels
       do j = js, je+1
       do i = is, ie+1
       do k = ks-1, ke+1
          ddiv(k,i,j) = j33g * ( momz(k,i,j) - momz(k-1,i  ,j  ) ) * rcdz(k) &
                      + ( ( momx(k+1,i,j) + momx(k+1,i-1,j  ) ) * j13g(k+1,i,j,i_xyw) &
                        - ( momx(k-1,i,j) + momx(k-1,i-1,j  ) ) * j13g(k-1,i,j,i_xyw) &
                        + ( momy(k+1,i,j) + momy(k+1,i  ,j-1) ) * j23g(k+1,i,j,i_xyw) &
                        - ( momy(k-1,i,j) + momy(k-1,i  ,j-1) ) * j23g(k-1,i,j,i_xyw) ) / ( fdz(k)+fdz(k-1) ) &
                      + mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) &
                      * ( ( momx(k,i  ,j  ) * gsqrt(k,i  ,j  ,i_uyz) / mapf(i  ,j  ,2,i_uy) &
                          - momx(k,i-1,j  ) * gsqrt(k,i-1,j  ,i_uyz) / mapf(i-1,j  ,2,i_uy) ) * rcdx(i) &
                        + ( momy(k,i  ,j  ) * gsqrt(k,i  ,j  ,i_xvz) / mapf(i  ,j  ,1,i_xv) &
                          - momy(k,i,  j-1) * gsqrt(k,i  ,j-1,i_xvz) / mapf(i  ,j-1,1,i_xv) ) * rcdy(j) )
       enddo
       enddo
       enddo
       !$acc end kernels
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif
       !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
       !$acc kernels
       do j = js, je+1
       do i = is, ie+1
          ddiv(ks,i,j) = j33g * ( momz(ks,i,j) ) * rcdz(ks) &
                       + ( ( momx(ks+1,i,j) + momx(ks+1,i-1,j  ) ) * j13g(ks+1,i,j,i_xyw) &
                         - ( momx(ks-1,i,j) + momx(ks  ,i-1,j  ) ) * j13g(ks  ,i,j,i_xyw) &
                         + ( momy(ks+1,i,j) + momy(ks+1,i  ,j-1) ) * j23g(ks+1,i,j,i_xyw) &
                         - ( momy(ks  ,i,j) + momy(ks  ,i  ,j-1) ) * j23g(ks  ,i,j,i_xyw) ) * rfdz(ks) &
                       + mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) &
                       * ( ( momx(ks,i  ,j  ) * gsqrt(ks,i  ,j  ,i_uyz) / mapf(i  ,j  ,2,i_uy) &
                           - momx(ks,i-1,j  ) * gsqrt(ks,i-1,j  ,i_uyz) / mapf(i-1,j  ,2,i_uy) ) * rcdx(i) &
                         + ( momy(ks,i  ,j  ) * gsqrt(ks,i  ,j  ,i_xvz) / mapf(i  ,j  ,1,i_xv) &
                           - momy(ks,i,  j-1) * gsqrt(ks,i  ,j-1,i_xvz) / mapf(i  ,j-1,1,i_xv) ) * rcdy(j) )
          ddiv(ke,i,j) = j33g * ( - momz(ke-1,i  ,j  ) ) * rcdz(ke) &
                       + ( ( momx(ke  ,i,j) + momx(ke  ,i-1,j  ) ) * j13g(ke  ,i,j,i_xyw) &
                         - ( momx(ke-1,i,j) + momx(ke-1,i-1,j  ) ) * j13g(ke-1,i,j,i_xyw) &
                         + ( momy(ke  ,i,j) + momy(ke  ,i  ,j-1) ) * j23g(ke  ,i,j,i_xyw) &
                         - ( momy(ke-1,i,j) + momy(ke-1,i  ,j-1) ) * j23g(ke-1,i,j,i_xyw) ) * rfdz(ke) &
                       + mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) &
                       * ( ( momx(ke,i  ,j  ) * gsqrt(ke,i  ,j  ,i_uyz) / mapf(i  ,j  ,2,i_uy) &
                           - momx(ke,i-1,j  ) * gsqrt(ke,i-1,j  ,i_uyz) / mapf(i-1,j  ,2,i_uy) ) * rcdx(i) &
                         + ( momy(ke,i  ,j  ) * gsqrt(ke,i  ,j  ,i_xvz) / mapf(i  ,j  ,1,i_xv) &
                           - momy(ke,i,  j-1) * gsqrt(ke,i  ,j-1,i_xvz) / mapf(i  ,j-1,1,i_xv) ) * rcdy(j) )
       enddo
       enddo
       !$acc end kernels
#ifdef DEBUG
       k = iundef; i = iundef; j = iundef
#endif
    end if
 
    !$acc end data
 
    call prof_rapend  ("DYN_divercence", 2)
 
    return
  end subroutine atmos_dyn_divergence
 
 
   !------------------------------------------------------------------------
   ! prepare thermodynamical data
   !   specific heat
   !   pressure data ( linearization )
   !
   ! pres = P0 * ( R * rhot / P0 )**(CP/CV)
   ! d pres / d rhot ~ CP*R/CV * ( R * rhot / P0 )**(R/CV)
   !                 = CP*R/CV * ( pres / P0 )**(R/CP)
   !                 = CP*R/CV * temp / pott
   !                 = CP/CV * pres / rhot
   ! pres ~ P0 * ( R * rhot0 / P0 ) ** (CP/CV) + CV*R/CP * ( pres / P0 )**(R/CP) * rhot'
   !------------------------------------------------------------------------
 
  subroutine atmos_dyn_prep_pres_linearization( &
   DPRES, RT2P, REF_rhot,                            & ! (out)
   rhot, qtrc, ref_pres, aq_r, aq_cv, aq_cp, aq_mass ) ! (in)
   use scale_const, only: &
      p0     => const_pre00, &
      rdry   => const_rdry,  &
      cvdry  => const_cvdry, &
      cpdry  => const_cpdry    
   implicit none
 
   real(rp), intent(out) :: dpres(ka,ia,ja)
   real(rp), intent(out) :: rt2p(ka,ia,ja)
   real(rp), intent(out) :: ref_rhot(ka,ia,ja)
   real(rp), intent(in) :: rhot(ka,ia,ja)
   real(rp), intent(in) :: qtrc(ka,ia,ja,qa)
   real(rp), intent(in) :: ref_pres(ka,ia,ja)
   real(rp), intent(in) :: aq_r(qa)
   real(rp), intent(in) :: aq_cv(qa)
   real(rp), intent(in) :: aq_cp(qa)
   real(rp), intent(in) :: aq_mass(qa)
 
   integer :: i, j, k
   integer :: iq    
   real(rp) :: qdry (ka) ! dry air
   real(rp) :: rtot (ka) ! total R
   real(rp) :: cvtot(ka) ! total CV
   real(rp) :: cptot(ka) ! total CP
   real(rp) :: pres      ! pressure
 
#ifdef DRY
   real(rp) :: cpovcv
#endif
 
   !--------------------------------------
 
#ifdef DRY
   cpovcv = cpdry / cvdry
#endif
 
!OCL XFILL
   !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
   !$omp shared(JA,IA,KS,KE) &
   !$omp shared(P0,Rdry,RHOT,AQ_R,AQ_CV,AQ_CP,QTRC,AQ_MASS,REF_rhot,REF_pres,CPdry,CVdry,QA,RT2P,DPRES) &
#ifdef DRY
   !$omp shared(CPovCV) &
#endif
   !$omp private(i,j,k,iq) &
   !$omp private(PRES,Rtot,CVtot,CPtot,QDRY)
   !$acc kernels copyout(DPRES, RT2P, REF_rhot) copyin(RHOT, QTRC, REF_pres, AQ_R, AQ_CV, AQ_CP, AQ_MASS)
   do j = 1, ja
   !$acc loop private(Rtot, CVtot, CPtot, QDRY)
   do i = 1, ia
      do k = ks, ke
         rtot(k) = 0.0_rp
         cvtot(k) = 0.0_rp
         cptot(k) = 0.0_rp
         qdry(k) = 1.0_rp
      end do
      !$acc loop seq
      do iq = 1, qa
         do k = ks, ke
            rtot(k) = rtot(k) + aq_r(iq) * qtrc(k,i,j,iq)
            cvtot(k) = cvtot(k) + aq_cv(iq) * qtrc(k,i,j,iq)
            cptot(k) = cptot(k) + aq_cp(iq) * qtrc(k,i,j,iq)
            qdry(k) = qdry(k) - qtrc(k,i,j,iq) * aq_mass(iq)
         enddo
      end do
      do k = ks, ke
         rtot(k) = rtot(k) + rdry  * qdry(k)
         cvtot(k) = cvtot(k) + cvdry * qdry(k)
         cptot(k) = cptot(k) + cpdry * qdry(k)
      end do
      do k = ks, ke
        pres = p0 * ( rtot(k) * rhot(k,i,j) / p0 )**( cptot(k) / cvtot(k) )
        rt2p(k,i,j) = cptot(k) / cvtot(k) * pres / rhot(k,i,j)
        dpres(k,i,j) = pres - ref_pres(k,i,j)
        ref_rhot(k,i,j) = rhot(k,i,j)
      end do
      dpres(ks-1,i,j) = dpres(ks+1,i,j) + ( ref_pres(ks+1,i,j) - ref_pres(ks-1,i,j) )
      dpres(ke+1,i,j) = dpres(ke-1,i,j) + ( ref_pres(ke-1,i,j) - ref_pres(ke+1,i,j) )
   end do
   end do
   !$acc end kernels
 
   return
 end subroutine atmos_dyn_prep_pres_linearization
 
end module scale_atmos_dyn_common