scale_atm_phy_mp_lscond.F90

Go to the documentation of this file.

!> module FElib / Atmosphere / Physics cloud microphysics / Large-scale condensation
!!
!! @par Description
!!     A module to provide a large-scale condensation scheme
!!
!! @par Reference
!!
!! @author Yuta Kawai, Team SCALE
!!
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_mp_lscond
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_prc, only: prc_abort 
  use scale_const, only: &
    rdry => const_rdry,    &
    rvap => const_rvap,    &
    cpdry => const_cpdry,  & 
    cvdry => const_cvdry,  & 
    pres00 => const_pre00, &
    grav => const_grav
  use scale_atmos_hydrometeor, only: &
    cp_vapor, &
    cp_water, &
    cp_ice,   &
    cv_vapor, &
    cv_water, &
    cv_ice,   &
    lhv
  
  use scale_element_base, only: &
    elementbase1d, elementbase3d
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
  public :: atmos_phy_mp_lscond_setup
  public :: atmos_phy_mp_lscond_adjustment
  public :: atmos_phy_mp_lscond_precipitation
  public :: atmos_phy_mp_lscond_precipitation_momentum
 
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !
  integer, private, parameter :: QA_MP  = 2
 
  integer,                parameter, public :: atmos_phy_mp_lscond_ntracers = qa_mp
  integer,                parameter, public :: atmos_phy_mp_lscond_nwaters = 1
  integer,                parameter, public :: atmos_phy_mp_lscond_nices = 0
  character(len=H_SHORT), parameter, public :: atmos_phy_mp_lscond_tracer_names(qa_mp) = (/ &
       'QV', &
       'QC'  /)
  character(len=H_MID)  , parameter, public :: atmos_phy_mp_lscond_tracer_descriptions(qa_mp) = (/ &
       'Ratio of Water Vapor mass to total mass (Specific humidity)', &
       'Ratio of Cloud Water mass to total mass                    '  /)
  character(len=H_SHORT), parameter, public :: atmos_phy_mp_lscond_tracer_units(qa_mp) = (/ &
       'kg/kg',  &
       'kg/kg'   /)
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
 
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,  private, parameter   :: i_qv = 1
  integer,  private, parameter   :: i_qc = 2
 
  integer,  private, parameter   :: i_hyd_qc =  1
 
  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
 
  !-----------------------------------------------------------------------------
  !> Setup a module for a large-scale condensation scheme
  !!
  subroutine atmos_phy_mp_lscond_setup
    use scale_prc, only: &
       prc_abort
    implicit none
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_MP_lscond_setup",*) 'Setup'
    log_info("ATMOS_PHY_MP_lscond_setup",*) 'large-scale condensation'
 
    if( couple_aerosol ) then
       log_error("ATMOS_PHY_MP_lscond_setup",*) 'MP_aerosol_couple should be .false. for large-scale condensation type MP!'
       call prc_abort
    endif
 
    return
  end subroutine atmos_phy_mp_lscond_setup
 
  !> Calculate a state after the saturation process
  !!
!OCL SERIAL
  subroutine atmos_phy_mp_lscond_adjustment( &
    KA, KS, KE, IA, IS, IE, JA, JS, JE, & ! (in)
    dens, pres,                         & ! (in)
    dt,                                 & ! (in)
    temp, qtrc, cptot, cvtot,           & ! (inout)
    rhoe_t, evaporate                   ) ! (out)
 
    use scale_atmos_saturation, only: &
      atmos_saturation_pres2qsat_liq
    use scale_const, only: &
      cl => const_cl
    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)
    real(rp), intent(in) :: pres(ka,ia,ja)
    real(rp), intent(in) :: dt
    real(rp), intent(inout) :: temp(ka,ia,ja) 
    real(rp), intent(inout) :: qtrc(ka,ia,ja,qa_mp)
    real(rp), intent(inout) :: cptot(ka,ia,ja) 
    real(rp), intent(inout) :: cvtot(ka,ia,ja) 
    real(rp), intent(out) :: rhoe_t(ka,ia,ja)
    real(rp), intent(out) :: evaporate(ka,ia,ja)
 
    real(rp) :: qsat(ka,ia,ja)
    real(rp) :: d_qc(ka)
    real(rp) :: d_cv(ka)
    real(rp) :: d_cp(ka)
    real(rp) :: d_e(ka)
    real(rp) :: cptot_(ka)
    real(rp) :: cvtot_(ka)
    real(rp) :: coef1(ka)
    real(rp) :: coef2
 
    integer :: i, j
 
    real(rp) :: rdt
    !----------------------------------------------
 
    call atmos_saturation_pres2qsat_liq( &
      ka, ks, ke, ia, is, ie, ja, js, je, & ! (in)
      temp, pres,                         & ! (in)
      qsat                                ) ! (out)
 
 
    coef2 = rvap / lhv
    rdt = 1.0_rp / dt
 
    !$omp parallel do private(i,j, d_qc, d_cv, d_cp, d_e, cptot_, cvtot_, coef1) collapse(2)
    do j=js, je
    do i=is, ie
  !    coef1 = LHV**2 / ( CVDry * Rvap )
      coef1(:) = lhv * ( lhv + ( cv_vapor - cv_water ) * temp(:,i,j) ) / ( cvtot(:,i,j) * rvap )
 
      d_qc(:) = max( ( qtrc(:,i,j,i_qv) - qsat(:,i,j) ) / ( 1.0_rp + coef1 / temp(:,i,j)**2 * ( 1.0_rp - coef2 * temp(:,i,j) ) * qsat(:,i,j) ), &
                     0.0_rp )
 
      qtrc(:,i,j,i_qv) = qtrc(:,i,j,i_qv) - d_qc(:)
      qtrc(:,i,j,i_qc) = qtrc(:,i,j,i_qc) + d_qc(:)
 
      d_cp(:) = - cp_vapor * d_qc(:) &
                + cp_water * d_qc(:)
      cptot_(:) = cptot(:,i,j) + d_cp(:)
 
      d_cv(:) = - cv_vapor * d_qc(:) &
                + cv_water * d_qc(:)
      cvtot_(:) = cvtot(:,i,j) + d_cv(:)
 
      d_e(:) = lhv * d_qc(:)
      rhoe_t(:,i,j) = dens(:,i,j) * d_e(:) * rdt
 
 
      temp(:,i,j) = ( temp(:,i,j) * cvtot(:,i,j) + d_e(:) ) / cvtot_(:)
      cptot(:,i,j) = cptot_(:)
      cvtot(:,i,j) = cvtot_(:)
    end do
    end do
 
    return
  end subroutine atmos_phy_mp_lscond_adjustment
 
  !> Calculate a state after the precipitation process
  !!
!OCL SERIAL
  subroutine atmos_phy_mp_lscond_precipitation( &
    DENS, RHOQ, CPtot, CVtot, RHOE,            & ! (inout)
    sflx_rain, sflx_snow, esflx,               & ! (inout)
    temp, dt,                                  & ! (in)
    qha, qla, qia, lcmesh, elem, elem1d        ) ! (in)
    
    use scale_const, only: &
      undef => const_undef8    
    implicit none
 
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: qha                   !< hydrometeor (water + ice)
    real(rp), intent(inout) :: dens (elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(inout) :: rhoq (elem%np,lcmesh%nez,lcmesh%ne2d,qha)
    real(rp), intent(inout) :: cptot(elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(inout) :: cvtot(elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(inout) :: rhoe (elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(inout) :: sflx_rain(elem%nfp_v,lcmesh%ne2da)
    real(rp), intent(inout) :: sflx_snow(elem%nfp_v,lcmesh%ne2da)
    real(rp), intent(inout) :: esflx    (elem%nfp_v,lcmesh%ne2da)
    real(rp), intent(in) :: temp (elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(in) :: dt
    integer, intent(in) :: qla, qia
    class(elementbase1d), intent(in) :: elem1d
 
    real(rp) :: rhocp(elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: rhocv(elem%np,lcmesh%nez,lcmesh%ne2d)
 
    real(rp) :: ddens(elem%np)
    real(rp) :: dinternalen(elem%np)
 
    real(rp) :: vint_weight(elem%nnode_v,elem%nnode_h1d**2)
    real(rp) :: condens_vint_lc(elem%nnode_h1d**2)
    real(rp) :: ien_vint_lc(elem%nnode_h1d**2)
 
    real(rp) :: eflx(elem%np)
    real(rp) :: cp(qha)
    real(rp) :: cv(qha)
 
    integer :: ke2d
    integer :: p2d
 
    integer :: ke_z
    integer :: ke
    integer :: iq
 
    real(rp) :: rdt
    !-------------------------------------------------------
 
    do iq = 1, qha
      if ( iq > qla + qia ) then
        cp(iq) = undef 
        cv(iq) = undef
      else if ( iq > qla ) then ! ice water
        cp(iq) = cp_ice
        cv(iq) = cv_ice
      else                      ! liquid water
        cp(iq) = cp_water
        cv(iq) = cv_water
      end if
    end do
 
    rdt = 1.0_rp / dt
 
    !$omp parallel do collapse(2)
    do ke2d = 1, lcmesh%Ne2D
    do ke_z = 1, lcmesh%NeZ
      rhocp(:,ke_z,ke2d) = cptot(:,ke_z,ke2d) * dens(:,ke_z,ke2d)
      rhocv(:,ke_z,ke2d) = cvtot(:,ke_z,ke2d) * dens(:,ke_z,ke2d)
    end do
    end do
 
    do iq = 1, qha
      !$omp parallel do private(ke2D,ke_z,ke,p2D, &
      !$omp dDENS, vint_weight, condens_vint_lc, dInternalEn, ien_vint_lc)
      do ke2d = 1, lcmesh%Ne2D
      do ke_z = 1, lcmesh%NeZ
        ke = ke2d + (ke_z-1)*lcmesh%Ne2D
 
        ddens(:) = - rhoq(:,ke_z,ke2d,iq)
        rhoq(:,ke_z,ke2d,iq) = 0.0_rp
 
        do p2d=1, elem%Nnode_h1D**2
          vint_weight(:,p2d) = 0.5_rp * elem1d%IntWeight_lgl(:) * ( lcmesh%zlev(elem%Colmask(elem%Nnode_v,p2d),ke) - lcmesh%zlev(elem%Colmask(1,p2d),ke) )
        end do
 
        do p2d=1, elem%Nnode_h1D**2
          condens_vint_lc(p2d) = sum( vint_weight(:,p2d) * ddens(elem%Colmask(:,p2d)) )
        end do
      
        if ( iq > qla ) then ! ice water
            sflx_snow(:,ke2d) = sflx_snow(:,ke2d)  &
                              + condens_vint_lc(:) * rdt
        else                 ! liquid water
            sflx_rain(:,ke2d) = sflx_rain(:,ke2d)  &
                              + condens_vint_lc(:) * rdt
        end if
 
        !--- update density
 
        rhocp(:,ke_z,ke2d) = rhocp(:,ke_z,ke2d) + cp(iq) * ddens(:)
        rhocv(:,ke_z,ke2d) = rhocv(:,ke_z,ke2d) + cv(iq) * ddens(:)
        dens(:,ke_z,ke2d) = dens(:,ke_z,ke2d) + ddens(:)
 
        !--- update internal energy   
 
        dinternalen(:) = cp(iq) * ddens(:) * temp(:,ke_z,ke2d)
 
        do p2d=1, elem%Nnode_h1D**2
          ien_vint_lc(p2d) = sum( vint_weight(:,p2d) * dinternalen(elem%Colmask(:,p2d))  )
        end do
        esflx(:,ke2d) = esflx(:,ke2d) &
                      + ien_vint_lc(:) * rdt
 
        rhoe(:,ke_z,ke2d) = rhoe(:,ke_z,ke2d) + dinternalen(:)
      end do
      end do
    end do
 
    !$omp parallel do collapse(2)
    do ke2d = 1, lcmesh%Ne2D
    do ke_z = 1, lcmesh%NeZ  
      cptot(:,ke_z,ke2d) = rhocp(:,ke_z,ke2d) / dens(:,ke_z,ke2d)
      cvtot(:,ke_z,ke2d) = rhocv(:,ke_z,ke2d) / dens(:,ke_z,ke2d)
    end do
    end do
 
    return
  end subroutine atmos_phy_mp_lscond_precipitation
 
  !> Calculate a tendency of momentum due to the precipitation process
  !!
!OCL SERIAL
  subroutine atmos_phy_mp_lscond_precipitation_momentum( &
    MOMU_t, MOMV_t, MOMZ_t,                & ! (out)
    dens, momu, momv, momz, dens_new,      & ! (in)
    rdt_mp, lcmesh, elem                   ) ! (in)
    implicit none
 
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: momu_t(elem%np,lcmesh%nea)
    real(rp), intent(out) :: momv_t(elem%np,lcmesh%nea)
    real(rp), intent(out) :: momz_t(elem%np,lcmesh%nea)
    real(rp), intent(in) :: dens(elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(in) :: momu(elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(in) :: momv(elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(in) :: momz(elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(in) :: dens_new(elem%np,lcmesh%nez,lcmesh%ne2d)
    real(rp), intent(in) :: rdt_mp
    
    integer :: ke2d
    integer :: ke_z
    integer :: ke
    real(rp) :: coef(elem%np)
    !----------------------------------------------------------
 
    !$omp parallel do collapse(2) private( &
    !$omp ke2D, ke_z, ke, coef )
    do ke2d = 1, lcmesh%Ne2D
    do ke_z = 1, lcmesh%NeZ
      ke = ke2d + (ke_z-1)*lcmesh%Ne2D
      coef(:) = ( dens_new(:,ke_z,ke2d) / dens(:,ke_z,ke2d) - 1.0_rp ) * rdt_mp
 
      momu_t(:,ke) = coef(:) * momu(:,ke_z,ke2d)
      momv_t(:,ke) = coef(:) * momv(:,ke_z,ke2d)
      momz_t(:,ke) = coef(:) * momz(:,ke_z,ke2d)
    end do
    end do
    return
  end subroutine atmos_phy_mp_lscond_precipitation_momentum
  
end module scale_atm_phy_mp_lscond