scale_atmos_phy_mp_kessler.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_phy_mp_kessler
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_mp_kessler_setup
  public :: atmos_phy_mp_kessler_adjustment
  public :: atmos_phy_mp_kessler_terminal_velocity
  public :: atmos_phy_mp_kessler_effective_radius
  public :: atmos_phy_mp_kessler_cloud_fraction
  public :: atmos_phy_mp_kessler_qtrc2qhyd
  public :: atmos_phy_mp_kessler_qhyd2qtrc
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, private, parameter :: QA_MP  = 3
 
  integer,                parameter, public :: atmos_phy_mp_kessler_ntracers = qa_mp
  integer,                parameter, public :: atmos_phy_mp_kessler_nwaters = 2
  integer,                parameter, public :: atmos_phy_mp_kessler_nices = 0
  character(len=H_SHORT), parameter, public :: atmos_phy_mp_kessler_tracer_names(qa_mp) = (/ &
       'QV', &
       'QC', &
       'QR'  /)
  character(len=H_MID)  , parameter, public :: atmos_phy_mp_kessler_tracer_descriptions(qa_mp) = (/ &
       'Ratio of Water Vapor mass to total mass (Specific humidity)', &
       'Ratio of Cloud Water mass to total mass                    ', &
       'Ratio of Rain Water mass to total mass                     '/)
  character(len=H_SHORT), parameter, public :: atmos_phy_mp_kessler_tracer_units(qa_mp) = (/ &
       'kg/kg',  &
       'kg/kg',  &
       'kg/kg'   /)
 
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: mp_kessler
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,  private, parameter   :: i_qv = 1
  integer,  private, parameter   :: i_qc = 2
  integer,  private, parameter   :: i_qr = 3
 
  integer,  private, parameter   :: i_hyd_qc =  1
  integer,  private, parameter   :: i_hyd_qr =  2
 
  logical,  private              :: flag_liquid = .true.     ! warm rain
  logical,  private              :: couple_aerosol = .false. ! Consider CCN effect ?
 
  real(rp), private, parameter   :: re_qc =  8.e-6_rp        ! effective radius for cloud water
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_kessler_setup
    use scale_prc, only: &
       prc_abort
    implicit none
 
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_MP_kessler_setup",*) 'Setup'
    log_info("ATMOS_PHY_MP_kessler_setup",*) 'KESSLER-type 1-moment bulk 3 category'
 
    if( couple_aerosol ) then
       log_error("ATMOS_PHY_MP_kessler_setup",*) 'MP_aerosol_couple should be .false. for KESSLER type MP!'
       call prc_abort
    endif
 
    return
  end subroutine atmos_phy_mp_kessler_setup
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_kessler_adjustment( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       DENS, PRES,               &
       dt,                       &
       TEMP, QTRC, CPtot, CVtot, &
       RHOE_t, EVAPORATE         )
    use scale_const, only: &
       dwatr => const_dwatr, &
       pi    => const_pi
    use scale_file_history, only: &
       file_history_in
    use scale_atmos_phy_mp_common, only: &
       mp_saturation_adjustment => atmos_phy_mp_saturation_adjustment
    implicit none
 
    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)       ! density [kg/m3]
    real(rp), intent(in)    :: pres     (ka,ia,ja)       ! pressure [Pa]
    real(dp), intent(in)    :: dt                        ! time interval of microphysics [s]
    real(rp), intent(inout) :: temp     (ka,ia,ja)       ! temperature [K]
    real(rp), intent(inout) :: qtrc     (ka,ia,ja,qa_mp) ! tracer mass concentration [kg/kg]
    real(rp), intent(inout) :: cptot    (ka,ia,ja)       ! total specific heat capacity at constant pressure [J/kg/K]
    real(rp), intent(inout) :: cvtot    (ka,ia,ja)       ! total specific heat capacity at constant volume [J/kg/K]
    real(rp), intent(out)   :: rhoe_t   (ka,ia,ja)       ! tendency of rhoe [J/m3/s]
    real(rp), intent(out)   :: evaporate(ka,ia,ja)       ! number of evaporated cloud [#/m3/s]
 
    real(rp) :: qtrc_dummy(ka,ia,ja)
    real(rp) :: rhoe_d_sat(ka,ia,ja)
    real(rp) :: qc_t_sat  (ka,ia,ja)
 
    integer  :: k, i, j
    !---------------------------------------------------------------------------
 
    log_progress(*) 'atmosphere / physics / microphysics / Kessler'
 
    !##### MP Main #####
    call mp_kessler( &
         ka, ks, ke, ia, is, ie, ja, js, je, &
         dens(:,:,:), pres(:,:,:),   & ! [IN]
         dt,                         & ! [IN]
         temp(:,:,:), qtrc(:,:,:,:), & ! [INOUT]
         cptot(:,:,:), cvtot(:,:,:), & ! [INOUT]
         rhoe_t(:,:,:)               ) ! [OUT]
 
    ! save value before saturation adjustment
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qc_t_sat(k,i,j) = qtrc(k,i,j,i_qc)
    enddo
    enddo
    enddo
 
!OCL XFILL
    qtrc_dummy(:,:,:) = -1.0_rp
 
    call mp_saturation_adjustment( &
         ka, ks, ke, ia, is, ie, ja, js, je,  &
         dens(:,:,:),                         & ! [IN]
         flag_liquid,                         & ! [IN]
         temp(:,:,:),                         & ! [INOUT]
         qtrc(:,:,:,i_qv),                    & ! [INOUT]
         qtrc(:,:,:,i_qc), qtrc_dummy(:,:,:), & ! [INOUT]
         cptot(:,:,:), cvtot(:,:,:),          & ! [INOUT]
         rhoe_d_sat(:,:,:)                    ) ! [OUT]
 
    do j = js, je
    do i = is, ie
    do k = ks, ke
       rhoe_t(k,i,j) = rhoe_t(k,i,j) + rhoe_d_sat(k,i,j) / dt
    enddo
    enddo
    enddo
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qc_t_sat(k,i,j) = ( qtrc(k,i,j,i_qc) - qc_t_sat(k,i,j) ) / dt
    enddo
    enddo
    enddo
 
    call file_history_in( qc_t_sat(:,:,:), 'Pcsat', 'QC production term by satadjust', 'kg/kg/s' )
 
    do j = js, je
    do i = is, ie
    do k = ks, ke
       evaporate(k,i,j) = max( -qc_t_sat(k,i,j), 0.0_rp ) & ! if negative, condensation
                        * dens(k,i,j) / (4.0_rp/3.0_rp*pi*dwatr*re_qc**3) ! mass -> number (assuming constant particle radius as re_qc)
    enddo
    enddo
    enddo
 
    !##### END MP Main #####
 
    return
  end subroutine atmos_phy_mp_kessler_adjustment
 
  !-----------------------------------------------------------------------------
  subroutine mp_kessler( &
         KA, KS, KE, IA, IS, IE, JA, JS, JE, &
         DENS0, PRES0,   &
         dt,             &
         TEMP0, QTRC0,   &
         CPtot0, CVtot0, &
         RHOE_t          )
    use scale_const, only: &
       eps => const_eps
    use scale_atmos_hydrometeor, only: &
       lhv, &
       cp_vapor, &
       cp_water, &
       cv_vapor, &
       cv_water
    use scale_atmos_saturation, only: &
       saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq
    implicit none
    integer,  intent(in) :: ka, ks, ke
    integer,  intent(in) :: ia, is, ie
    integer,  intent(in) :: ja, js, je
 
    real(rp), intent(in) :: dens0(ka,ia,ja)           ! density [kg/m3]
    real(rp), intent(in) :: pres0(ka,ia,ja)           ! pressure [Pa]
    real(dp), intent(in) :: dt                        ! time interval of microphysics [s]
 
    real(rp), intent(inout) :: temp0 (ka,ia,ja)       ! temperature [K]
    real(rp), intent(inout) :: qtrc0 (ka,ia,ja,qa_mp) ! tracer mass concentration [kg/kg]
    real(rp), intent(inout) :: cptot0(ka,ia,ja)       ! total specific heat capacity at constant pressure [J/kg/K]
    real(rp), intent(inout) :: cvtot0(ka,ia,ja)       ! total specific heat capacity at constant volume [J/kg/K]
 
    real(rp), intent(out) :: rhoe_t(ka,ia,ja)         ! tendency of rhoe [J/m3/s]
 
    real(rp) :: dens
    real(rp) :: temp
    real(rp) :: pres
    real(rp) :: cptot, cvtot
    real(rp) :: qv, qc, qr
    real(rp) :: qv_t, qc_t, qr_t
    real(rp) :: e_t, cp_t, cv_t
    real(rp) :: qsatl ! saturated water vapor for liquid water [kg/kg]
    real(rp) :: sliq  ! saturated ratio S for liquid water (0-1)
 
    ! tendency
    real(rp) :: dq_evap ! tendency q (evaporation)
    real(rp) :: dq_auto ! tendency q (autoconversion)
    real(rp) :: dq_accr ! tendency q (accretion)
 
    real(rp) :: vent_factor
    real(rp) :: rdt
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    call prof_rapstart('MP_kessler', 3)
 
    rdt = 1.0_rp / dt
 
    !$omp parallel do default(none) &
    !$omp shared(KS,KE,IS,IE,JS,JE, &
    !$omp        DENS0,TEMP0,PRES0,QTRC0,CPtot0,CVtot0,dt, &
    !$omp        RHOE_t, &
    !$omp        EPS,LHV,CP_VAPOR,CP_WATER,CV_VAPOR,CV_WATER,rdt) &
    !$omp private(k,i,j,dens,temp,pres,cptot,cvtot,qv,qc,qr,qv_t,qc_t,qr_t,e_t,cp_t,cv_t, &
    !$omp         QSATL,Sliq, &
    !$omp         dq_evap,dq_auto,dq_accr,vent_factor) &
    !$omp OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
    do k = ks, ke
 
       ! store to work
       dens = dens0(k,i,j)
       temp = temp0(k,i,j)
       pres = pres0(k,i,j)
       qv   = qtrc0(k,i,j,i_qv)
       qc   = qtrc0(k,i,j,i_qc)
       qr   = qtrc0(k,i,j,i_qr)
 
       call saturation_dens2qsat_liq( &
            temp, dens, & ! [IN]
            qsatl       ) ! [OUT]
 
       sliq = qv / max( qsatl, eps)
 
       ! Auto-conversion (QC->QR)
       dq_auto = 1.e-3_rp * max( qc-1.e-3_rp, 0.0_rp )
 
       ! Accretion (QC->QR)
       dq_accr = 2.2_rp * qc * qr**0.875_rp
 
       ! Evaporation (QR->QV)
       vent_factor = 1.6_rp + 124.9_rp * ( dens*qr )**0.2046_rp
       dq_evap = ( 1.0_rp-min(sliq,1.0_rp) ) / dens * vent_factor  &
               * ( dens*qr )**0.525_rp / ( 5.4e5_rp + 2.55e8_rp / ( pres*qsatl ) )
 
       ! limiter
       qc_t = (-dq_auto-dq_accr         )
       qc_t = max( qc_t, -qc*rdt )
 
       qr_t = ( dq_auto+dq_accr-dq_evap )
       qr_t = max( qr_t, -qr*rdt )
 
       qv_t = - ( qc_t + qr_t )
       qv_t = max( qv_t, -qv*rdt)
 
       ! mass & energy update
       qtrc0(k,i,j,i_qv) = qtrc0(k,i,j,i_qv) + qv_t * dt
       qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) + qc_t * dt
       qtrc0(k,i,j,i_qr) = qtrc0(k,i,j,i_qr) + qr_t * dt
 
       e_t = -lhv * qv_t ! internal energy change
       cp_t = cp_vapor * qv_t &
            + cp_water * ( qc_t + qr_t )
       cptot = cptot0(k,i,j) + cp_t * dt
 
       cv_t = cv_vapor * qv_t &
            + cv_water * ( qc_t + qr_t )
       cvtot = cvtot0(k,i,j) + cv_t * dt
 
       rhoe_t(k,i,j) = dens * e_t
 
       temp0(k,i,j) = ( temp * cvtot0(k,i,j) + e_t * dt ) / cvtot
       cptot0(k,i,j) = cptot
       cvtot0(k,i,j) = cvtot
 
    enddo
    enddo
    enddo
 
    call prof_rapend  ('MP_kessler', 3)
 
    return
  end subroutine mp_kessler
 
  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine atmos_phy_mp_kessler_terminal_velocity( &
       KA, KS, KE,            &
       DENS0, RHOQ0,          &
       REFSTATE_dens_profile, &
       vterm                  )
    implicit none
 
    integer,  intent(in)  :: ka, ks, ke
    real(rp), intent(in)  :: dens0                (ka)         ! density [kg/m3]
    real(rp), intent(in)  :: rhoq0                (ka,qa_mp-1) ! density of each hydrometeor [kg/m3]
    real(rp), intent(in)  :: refstate_dens_profile(ka)         ! reference state density profile [kg/m3]
    real(rp), intent(out) :: vterm                (ka,qa_mp-1) ! terminal velocity of each hydrometeor [m/s]
 
    real(rp) :: qr
    real(rp) :: zerosw
 
    integer  :: k
    !---------------------------------------------------------------------------
 
    do k = ks, ke
       vterm(k,i_hyd_qc) = 0.0_rp
    enddo
 
    do k = ks, ke
       qr     = rhoq0(k,i_hyd_qr) / dens0(k)
       zerosw = 0.5_rp - sign(0.5_rp, qr - 1.e-12_rp )
 
       vterm(k,i_hyd_qr) = - 36.34_rp * ( dens0(k) * ( qr + zerosw ) )**0.1364_rp &
                         * refstate_dens_profile(ks) / refstate_dens_profile(k) * ( 1.0_rp - zerosw )
    enddo
 
    return
  end subroutine atmos_phy_mp_kessler_terminal_velocity
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_kessler_cloud_fraction( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       QTRC,           &
       mask_criterion, &
       cldfrac         )
    implicit none
    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je
 
    real(rp), intent(in)  :: qtrc(ka,ia,ja,qa_mp-1) ! hydrometeor mass concentration [kg/kg]
    real(rp), intent(in)  :: mask_criterion         ! criterion of hydrometeor [kg/kg]
 
    real(rp), intent(out) :: cldfrac(ka,ia,ja)      ! cloud fraction (0 or 1)
 
    real(rp) :: qhydro
    integer  :: k, i, j
    !---------------------------------------------------------------------------
 
    do j  = js, je
    do i  = is, ie
    do k  = ks, ke
       qhydro = qtrc(k,i,j,i_hyd_qc) + qtrc(k,i,j,i_hyd_qr)
       cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
    enddo
    enddo
    enddo
 
    return
  end subroutine atmos_phy_mp_kessler_cloud_fraction
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_kessler_effective_radius( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       DENS0, TEMP0, QTRC0, &
       Re                   )
    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc, &
       i_hr
    implicit none
    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je
 
    real(rp), intent(in)  :: dens0(ka,ia,ja)         ! density [kg/m3]
    real(rp), intent(in)  :: temp0(ka,ia,ja)         ! temperature [K]
    real(rp), intent(in)  :: qtrc0(ka,ia,ja,qa_mp-1) ! hydrometeor mass concentration [kg/kg]
 
    real(rp), intent(out) :: re   (ka,ia,ja,n_hyd)   ! effective radius [cm]
 
    real(rp), parameter :: um2cm = 100.0_rp
    !---------------------------------------------------------------------------
 
!OCL XFILL
    re(:,:,:,i_hc) =   8.e-6_rp * um2cm
!OCL XFILL
    re(:,:,:,i_hr) = 100.e-6_rp * um2cm
!OCL XFILL
    re(:,:,:,i_hr+1:) = 0.0_rp
 
    return
  end subroutine atmos_phy_mp_kessler_effective_radius
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_kessler_qtrc2qhyd( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       QTRC, &
       Qe    )
    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc, &
       i_hr
    implicit none
    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je
 
    real(rp), intent(in)  :: qtrc(ka,ia,ja,qa_mp-1) ! hydrometeor mass concentration [kg/kg]
 
    real(rp), intent(out) :: qe(ka,ia,ja,n_hyd)   ! mass ratio of each hydrometeor [kg/kg]
    !---------------------------------------------------------------------------
 
!OCL XFILL
    qe(:,:,:,i_hc) = qtrc(:,:,:,i_hyd_qc)
!OCL XFILL
    qe(:,:,:,i_hr) = qtrc(:,:,:,i_hyd_qr)
!OCL XFILL
    qe(:,:,:,i_hr+1:) = 0.0_rp
 
    return
  end subroutine atmos_phy_mp_kessler_qtrc2qhyd
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_kessler_qhyd2qtrc( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       Qe,  &
       QTRC )
    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc, &
       i_hr
    implicit none
    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je
 
    real(rp), intent(in) :: qe(ka,ia,ja,n_hyd)   ! mass ratio of each hydrometeor [kg/kg]
 
    real(rp), intent(out) :: qtrc(ka,ia,ja,qa_mp-1) ! hydrometeor mass concentration [kg/kg]
 
    !---------------------------------------------------------------------------
 
!OCL XFILL
    qtrc(:,:,:,i_hyd_qc) = qe(:,:,:,i_hc)
!OCL XFILL
    qtrc(:,:,:,i_hyd_qr) = qe(:,:,:,i_hr)
 
    ! ignore ice water
 
    return
  end subroutine atmos_phy_mp_kessler_qhyd2qtrc
 
end module scale_atmos_phy_mp_kessler