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)
    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

    ! liquid
    do n = 1, nl
       nq = nq + 1
       cv(nq) = cv_water
       cp(nq) = cp_water
       r(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
    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,  & ! [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