scale_atmos_hydrometeor.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_hydrometeor
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_hydrometeor_setup
  public :: atmos_hydrometeor_regist
  public :: atmos_hydrometeor_lhv
  public :: atmos_hydrometeor_lhs
  public :: atmos_hydrometeor_lhf
  public :: atmos_hydrometeor_entr
  public :: atmos_hydrometeor_entr2temp
 
  interface atmos_hydrometeor_regist
     module procedure atmos_hydrometeor_regist
  end interface atmos_hydrometeor_regist
 
  interface atmos_hydrometeor_lhv
     module procedure atmos_hydrometeor_lhv_0d
     module procedure atmos_hydrometeor_lhv_1d
     module procedure atmos_hydrometeor_lhv_2d
     module procedure atmos_hydrometeor_lhv_3d
  end interface atmos_hydrometeor_lhv
 
  interface atmos_hydrometeor_lhs
     module procedure atmos_hydrometeor_lhs_0d
     module procedure atmos_hydrometeor_lhs_1d
     module procedure atmos_hydrometeor_lhs_2d
     module procedure atmos_hydrometeor_lhs_3d
  end interface atmos_hydrometeor_lhs
 
  interface atmos_hydrometeor_lhf
     module procedure atmos_hydrometeor_lhf_0d
     module procedure atmos_hydrometeor_lhf_1d
     module procedure atmos_hydrometeor_lhf_2d
     module procedure atmos_hydrometeor_lhf_3d
  end interface atmos_hydrometeor_lhf
 
  interface atmos_hydrometeor_entr
     module procedure atmos_hydrometeor_entr_0d
     module procedure atmos_hydrometeor_entr_2d
     module procedure atmos_hydrometeor_entr_3d
  end interface atmos_hydrometeor_entr
 
  interface atmos_hydrometeor_entr2temp
     module procedure atmos_hydrometeor_entr2temp_0d
  end interface atmos_hydrometeor_entr2temp
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  logical, public            :: atmos_hydrometeor_ice_phase
 
 
  integer, public            :: i_qv = -1
 
  integer, public, parameter :: n_hyd = 6
 
  integer, public, parameter :: i_hc =  1
  integer, public, parameter :: i_hr =  2
  integer, public, parameter :: i_hi =  3
  integer, public, parameter :: i_hs =  4
  integer, public, parameter :: i_hg =  5
  integer, public, parameter :: i_hh =  6
 
  character(len=H_SHORT), public, parameter :: hyd_name(n_hyd) = &
       (/ "QC", "QR", "QI", "QS", "QG", "QH" /)
  character(len=H_MID),   public, parameter :: hyd_desc(n_hyd) = &
       (/ "cloud    ", "rain     ", "ice water", "snow     ", "graupel  ", "hail     " /)
  character(len=H_SHORT), public, parameter :: num_name(n_hyd) = &
       (/ "NC", "NR", "NI", "NS", "NG", "NH" /)
 
  real(rp), public :: hyd_dens(n_hyd)
 
  logical, public :: atmos_hydrometeor_dry = .true.
 
!  integer, public            :: I_NC = -1
!  integer, public            :: I_NR = -1
!  integer, public            :: I_NI = -1
!  integer, public            :: I_NS = -1
!  integer, public            :: I_NG = -1
!  integer, public            :: I_NH = -1
 
!  integer, public            :: I_QC = -1
!  integer, public            :: I_QR = -1
!  integer, public            :: I_QI = -1
!  integer, public            :: I_QS = -1
!  integer, public            :: I_QG = -1
!  integer, public            :: I_QH = -1
 
  ! hydrometeor (water + ice)
  integer, public            :: qha =  0
  integer, public            :: qhs = -1
  integer, public            :: qhe = -2
  ! water
  integer, public            :: qla =  0
  integer, public            :: qls = -1
  integer, public            :: qle = -2
  ! ice
  integer, public            :: qia =  0
  integer, public            :: qis = -1
  integer, public            :: qie = -2
 
  real(rp), public           :: lhv
  real(rp), public           :: lhs
  real(rp), public           :: lhf
 
  real(rp), public           :: cv_vapor
  real(rp), public           :: cp_vapor
  real(rp), public           :: cv_water
  real(rp), public           :: cp_water
  real(rp), public           :: cv_ice
  real(rp), public           :: cp_ice
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(rp), private :: thermodyn_emask = 1.0_rp
 
  logical,  private :: initialized = .false.
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_setup
    use scale_const, only: &
       const_setup, &
       cpvap          => const_cpvap, &
       cvvap          => const_cvvap, &
       cl             => const_cl,    &
       ci             => const_ci,    &
       lhv00          => const_lhv00, &
       lhs00          => const_lhs00, &
       lhf00          => const_lhf00, &
       lhv0           => const_lhv0,  &
       lhs0           => const_lhs0,  &
       lhf0           => const_lhf0,  &
       dwatr          => const_dwatr, &
       dice           => const_dice,  &
       thermodyn_type => const_thermodyn_type
    use scale_prc, only: &
       prc_abort
    implicit none
    !---------------------------------------------------------------------------
 
    if ( initialized ) return
    initialized = .true.
 
    call const_setup
 
    log_newline
    log_info("ATMOS_HYDROMETEOR_setup",*) 'Setup'
 
    if ( thermodyn_type == 'EXACT' ) then
 
       cv_vapor = cvvap
       cp_vapor = cpvap
       cv_water = cl
       cp_water = cv_water
       cv_ice   = ci
       cp_ice   = cv_ice
 
       lhv = lhv00
       lhs = lhs00
       lhf = lhf00
       thermodyn_emask = 1.0_rp
 
    elseif( thermodyn_type == 'SIMPLE' ) then
 
       cv_vapor = cvvap
       cp_vapor = cpvap
       cv_water = cvvap
       cp_water = cv_water
       cv_ice   = cvvap
       cp_ice   = cv_ice
 
       lhv = lhv0
       lhs = lhs0
       lhf = lhf0
       thermodyn_emask = 0.0_rp
 
    else
       log_error("ATMOS_HYDROMETEOR_setup",*) 'Not appropriate ATMOS_THERMODYN_ENERGY_TYPE. Check!', trim(thermodyn_type)
       call prc_abort
    endif
 
    hyd_dens(:) = (/ dwatr, & ! HC
                     dwatr, & ! HR
                     dice,  & ! HI
                     dice,  & ! HS
                     dice,  & ! HG
                     dice  /) ! HH
 
    return
  end subroutine atmos_hydrometeor_setup
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_regist( &
       NL,   &
       NI,   &
       NAME, &
       DESC, &
       UNIT, &
       Q0,   &
       ADVC  )
    use scale_prc, only: &
      prc_abort
    use scale_tracer, only: &
      tracer_regist
    use scale_const, only: &
      rvap => const_rvap
    implicit none
    integer,          intent(in)  :: nl
    integer,          intent(in)  :: ni
    character(len=*), intent(in)  :: name(1+nl+ni)
    character(len=*), intent(in)  :: desc(1+nl+ni)
    character(len=*), intent(in)  :: unit(1+nl+ni)
 
    integer,          intent(out) :: q0
 
    logical,          intent(in), optional :: advc(1+nl+ni)
 
    real(rp) :: cv   (1+nl+ni)
    real(rp) :: cp   (1+nl+ni)
    real(rp) :: r    (1+nl+ni)
    real(rp) :: ei0  (1+nl+ni)
    logical  :: mass (1+nl+ni)
    logical  :: advc_(1+nl+ni)
 
    integer  :: nq
    integer  :: n
    !---------------------------------------------------------------------------
 
    if ( .not. atmos_hydrometeor_dry ) then
       log_error("ATMOS_HYDROMETEOR_regist",*) 'tracer for hydrometeor is already registerd'
       call prc_abort
    endif
 
    atmos_hydrometeor_dry = .false.
 
    nq = 0
 
    ! vapor
    nq = nq + 1
    cv(nq) = cv_vapor
    cp(nq) = cp_vapor
    r(nq) = rvap
    ei0(nq) = lhv
 
    ! liquid
    do n = 1, nl
       nq = nq + 1
       cv(nq) = cv_water
       cp(nq) = cp_water
       r(nq) = 0.0_rp
       ei0(nq) = 0.0_rp
    end do
 
    ! ice
    do n = 1, ni
       nq = nq + 1
       cv(nq) = cv_ice
       cp(nq) = cp_ice
       r(nq) = 0.0_rp
       ei0(nq) =  - lhf
    end do
 
    ! NQ = 1 + NL + NI,   vapor + liqid + ice
 
    if ( present(advc) ) then
       advc_(:) = advc(:)
    else
       advc_(:) = .true.
    endif
 
    do n = 1, nq
       mass(n) = .true.
    end do
 
    call tracer_regist( q0,             & ! [OUT]
                        nq,             & ! [IN]
                        name,           & ! [IN]
                        desc,           & ! [IN]
                        unit,           & ! [IN]
                        cv, cp, r, ei0, & ! [IN], optional
                        advc_, mass     ) ! [IN], optional
 
    i_qv = q0
 
    if ( nq > 1 ) then
       qhs = q0 + 1
       qha = nl + ni
       qhe = qhs + qha - 1
    endif
 
    if ( nl > 0 ) then
       qls = q0 + 1
       qla = nl
       qle = qls + nl - 1
    endif
 
    if ( ni > 0 ) then
       qis = qle + 1
       qia = ni
       qie = qis + ni - 1
    endif
 
    atmos_hydrometeor_ice_phase = qia > 0
 
    return
  end subroutine atmos_hydrometeor_regist
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhv_0d( &
       temp, &
       lhv   )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhv0  => const_lhv0
    implicit none
    real(RP), intent(in)  :: temp
    real(RP), intent(out) :: lhv
    !---------------------------------------------------------------------------
 
    lhv = lhv0 + ( cp_vapor - cp_water ) * ( temp - tem00 ) * thermodyn_emask
 
    return
  end subroutine atmos_hydrometeor_lhv_0d
 
  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine atmos_hydrometeor_lhv_1d( &
       KA, KS, KE, &
       temp, &
       lhv   )
    implicit none
    integer,  intent(in)  :: KA, KS, KE
 
    real(RP), intent(in)  :: temp(KA)
 
    real(RP), intent(out) :: lhv (KA)
 
    integer :: k
    !---------------------------------------------------------------------------
 
    do k = ks, ke
       call atmos_hydrometeor_lhv_0d( temp(k), lhv(k) )
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhv_1d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhv_2d( &
       IA, IS, IE, JA, JS, JE, &
       temp, &
       lhv   )
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in)  :: temp(IA,JA)
    real(RP), intent(out) :: lhv (IA,JA)
 
    integer :: i, j
    !---------------------------------------------------------------------------
 
    !$omp parallel do OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
       call atmos_hydrometeor_lhv_0d( temp(i,j), lhv(i,j) )
    enddo
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhv_2d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhv_3d( &
       KA, KS, KE, &
       IA, IS, IE, &
       JA, JS, JE, &
       temp, &
       lhv   )
    implicit none
    integer,  intent(in)  :: KA, KS, KE
    integer,  intent(in)  :: IA, IS, IE
    integer,  intent(in)  :: JA, JS, JE
 
    real(RP), intent(in)  :: temp(KA,IA,JA)
    real(RP), intent(out) :: lhv (KA,IA,JA)
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    do j = js, je
    do i = is, ie
    do k = ks, ke
       call atmos_hydrometeor_lhv_0d( temp(k,i,j), lhv(k,i,j) )
    enddo
    enddo
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhv_3d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhs_0d( &
       temp, &
       lhs   )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhs0  => const_lhs0
    implicit none
    real(RP), intent(in)  :: temp
    real(RP), intent(out) :: lhs
    !---------------------------------------------------------------------------
 
    lhs = lhs0 + ( cp_vapor - cp_ice ) * ( temp - tem00 ) * thermodyn_emask
 
    return
  end subroutine atmos_hydrometeor_lhs_0d
 
  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine atmos_hydrometeor_lhs_1d( &
       KA, KS, KE, &
       temp, &
       lhs   )
    implicit none
    integer, intent(in) :: KA, KS, KE
 
    real(RP), intent(in)  :: temp(KA)
 
    real(RP), intent(out) :: lhs (KA)
 
    integer :: k
    !---------------------------------------------------------------------------
 
    do k = ks, ke
       call atmos_hydrometeor_lhs( temp(k), lhs(k) )
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhs_1d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhs_2d( &
       IA, IS, IE, JA, JS, JE, &
       temp, &
       lhs   )
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in)  :: temp(IA,JA)
    real(RP), intent(out) :: lhs (IA,JA)
 
    integer :: i, j
    !---------------------------------------------------------------------------
 
    !$omp parallel do OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
       call atmos_hydrometeor_lhs( temp(i,j), lhs(i,j) )
    enddo
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhs_2d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhs_3d( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       temp, &
       lhs   )
    implicit none
 
    integer,  intent(in)  :: KA, KS, KE
    integer,  intent(in)  :: IA, IS, IE
    integer,  intent(in)  :: JA, JS, JE
 
    real(RP), intent(in)  :: temp(KA,IA,JA)
    real(RP), intent(out) :: lhs (KA,IA,JA)
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    !$omp parallel do OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       call atmos_hydrometeor_lhs( temp(k,i,j), lhs(k,i,j) )
    enddo
    enddo
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhs_3d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhf_0d( &
       temp, &
       lhf   )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhf0  => const_lhf0
    implicit none
    real(RP), intent(in)  :: temp
    real(RP), intent(out) :: lhf
    !---------------------------------------------------------------------------
 
    lhf = lhf0 + ( cp_water - cp_ice ) * ( temp - tem00 ) * thermodyn_emask
 
    return
  end subroutine atmos_hydrometeor_lhf_0d
 
  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine atmos_hydrometeor_lhf_1d( &
       KA, KS, KE, &
       temp, &
       lhf   )
    implicit none
    integer, intent(in) :: KA, KS, KE
 
    real(RP), intent(in)  :: temp(KA)
    real(RP), intent(out) :: lhf (KA)
 
    integer :: k
    !---------------------------------------------------------------------------
 
    do k = ks, ke
       call atmos_hydrometeor_lhf( temp(k), lhf(k) )
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhf_1d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhf_2d( &
       IA, IS, IE, JA, JS, JE, &
       temp, &
       lhf   )
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in)  :: temp(IA,JA)
    real(RP), intent(out) :: lhf (IA,JA)
 
    integer :: i, j
    !---------------------------------------------------------------------------
 
    !$omp parallel do OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
       call atmos_hydrometeor_lhf( temp(i,j), lhf(i,j) )
    enddo
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhf_2d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhf_3d( &
       KA, KS, KE, &
       IA, IS, IE, &
       JA, JS, JE, &
       temp, &
       lhf   )
    implicit none
    integer,  intent(in)  :: KA, KS, KE
    integer,  intent(in)  :: IA, IS, IE
    integer,  intent(in)  :: JA, JS, JE
 
    real(RP), intent(in)  :: temp(KA,IA,JA)
    real(RP), intent(out) :: lhf (KA,IA,JA)
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    !$omp parallel do OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       call atmos_hydrometeor_lhf( temp(k,i,j), lhf(k,i,j) )
    enddo
    enddo
    enddo
 
    return
  end subroutine atmos_hydrometeor_lhf_3d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_entr_0d( &
       TEMP, PRES,   &
       QV, QI, Qdry, &
       Rtot, CPtot,  &
       entr          )
    use scale_const, only: &
       pre00 => const_pre00, &
       tem00 => const_tem00, &
       rdry  => const_rdry,  &
       rvap  => const_rvap,  &
       lhv0  => const_lhv0,  &
       lhf0  => const_lhf0,  &
       psat0 => const_psat0
    implicit none
 
    real(RP), intent(in) :: TEMP
    real(RP), intent(in) :: PRES
    real(RP), intent(in) :: QV
    real(RP), intent(in) :: QI
    real(RP), intent(in) :: Qdry
    real(RP), intent(in) :: Rtot
    real(RP), intent(in) :: CPtot
 
    real(RP), intent(out) :: entr
 
    real(RP) :: pres_dry, pres_vap
    !---------------------------------------------------------------------------
 
    pres_dry = pres * qdry * rdry / rtot
    pres_vap = pres * qv   * rvap / rtot
 
    entr = cptot * log( temp / tem00 ) &
         - qdry * rdry * log( pres_dry / pre00 ) &
         - qv   * rvap * log( pres_vap / psat0 ) &
         + ( qv * lhv0 - qi * lhf0 ) / tem00
 
    return
  end subroutine atmos_hydrometeor_entr_0d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_entr_2d( &
       IA, IS, IE, JA, JS, JE, &
       TEMP, PRES,   &
       QV, QI, Qdry, &
       Rtot, CPtot,  &
       entr          )
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in) :: TEMP (IA,JA)
    real(RP), intent(in) :: PRES (IA,JA)
    real(RP), intent(in) :: QV   (IA,JA)
    real(RP), intent(in) :: QI   (IA,JA)
    real(RP), intent(in) :: Qdry (IA,JA)
    real(RP), intent(in) :: Rtot (IA,JA)
    real(RP), intent(in) :: CPtot(IA,JA)
 
    real(RP), intent(out) :: entr(IA,JA)
 
    integer  :: i, j
    !---------------------------------------------------------------------------
 
    !$omp parallel do OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
       call atmos_hydrometeor_entr_0d( &
            temp(i,j), pres(i,j),        & ! [IN]
            qv(i,j), qi(i,j), qdry(i,j), & ! [IN]
            rtot(i,j), cptot(i,j),       & ! [IN]
            entr(i,j)                    ) ! [OUT]
    enddo
    enddo
 
    return
  end subroutine atmos_hydrometeor_entr_2d
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_entr_3d( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       TEMP, PRES,   &
       QV, QI, Qdry, &
       Rtot, CPtot,  &
       entr          )
    implicit none
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in) :: TEMP (KA,IA,JA)
    real(RP), intent(in) :: PRES (KA,IA,JA)
    real(RP), intent(in) :: QV   (KA,IA,JA)
    real(RP), intent(in) :: QI   (KA,IA,JA)
    real(RP), intent(in) :: Qdry (KA,IA,JA)
    real(RP), intent(in) :: Rtot (KA,IA,JA)
    real(RP), intent(in) :: CPtot(KA,IA,JA)
 
    real(RP), intent(out) :: entr(KA,IA,JA)
 
    integer  :: k, i, j
    !---------------------------------------------------------------------------
 
    !$omp parallel do OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       call atmos_hydrometeor_entr_0d( &
            temp(k,i,j), pres(k,i,j),          & ! [IN]
            qv(k,i,j), qi(k,i,j), qdry(k,i,j), & ! [IN]
            rtot(k,i,j), cptot(k,i,j),         & ! [IN]
            entr(k,i,j)                        ) ! [OUT]
    enddo
    enddo
    enddo
 
    return
  end subroutine atmos_hydrometeor_entr_3d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_entr2temp_0d( &
       entr, pres, &
       qv, qi, qdry, &
       Rtot, CPtot,    &
       temp            )
    use scale_const, only: &
       pre00 => const_pre00, &
       tem00 => const_tem00, &
       rdry  => const_rdry,  &
       rvap  => const_rvap,  &
       lhv0  => const_lhv0,  &
       lhf0  => const_lhf0,  &
       psat0 => const_psat0
    implicit none
 
    real(RP), intent(in) :: entr
    real(RP), intent(in) :: pres
    real(RP), intent(in) :: qv
    real(RP), intent(in) :: qi
    real(RP), intent(in) :: qdry
    real(RP), intent(in) :: Rtot
    real(RP), intent(in) :: CPtot
 
    real(RP), intent(out) :: temp
 
    real(RP) :: pres_dry, pres_vap
    !---------------------------------------------------------------------------
 
    pres_dry = pres * qdry * rdry / rtot
    pres_vap = pres * qv   * rvap / rtot
 
    temp = tem00 &
         * exp( ( entr &
                + qdry * rdry * log( pres_dry / pre00 ) &
                + qv   * rvap * log( pres_vap / psat0 ) &
              - ( qv * lhv0 - qi * lhf0 ) / tem00 ) / cptot )
    return
  end subroutine atmos_hydrometeor_entr2temp_0d
 
end module scale_atmos_hydrometeor