module cldwat2m_macro 3,7
!--------------------------------------------------- !
! Purpose : CAM Interface for Cloud Macrophysics !
! Author : Sungsu Park !
! Description : Park et al. 2010. !
! For questions, contact Sungsu Park !
! e-mail : sungsup@ucar.edu !
! phone : 303-497-1375 !
!--------------------------------------------------- !
use shr_kind_mod
, only: r8=>shr_kind_r8
use spmd_utils, only: masterproc
use ppgrid, only: pcols, pver, pverp
use abortutils, only: endrun
use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice, &
vqsatd2_water, vqsatd2_water_single, polysvp
use cam_history, only: addfld, add_default, phys_decomp, outfld
use cam_logfile, only: iulog
implicit none
private
save
public :: ini_macro, mmacro_pcond
! -------------- !
! Set Parameters !
! -------------- !
! -------------------------- !
! Parameters for Ice Stratus !
! -------------------------- !
integer, parameter :: iceopt = 5 ! Option for ice cloud closure. 5 : Modified Slingo Formula. For other options, modify each 2 subroutine.
real(r8), public, parameter :: rhmini = 0.80_r8 ! Minimum rh for ice cloud fraction > 0.
real(r8), public, parameter :: rhmaxi = 1.1_r8 ! rhi at which ice cloud fraction = 1.
real(r8), parameter :: qist_min = 1.e-7_r8 ! Minimum in-stratus ice IWC constraint [ kg/kg ]
real(r8), parameter :: qist_max = 5.e-3_r8 ! Maximum in-stratus ice IWC constraint [ kg/kg ]
! ----------------------------- !
! Parameters for Liquid Stratus !
! ----------------------------- !
logical, parameter :: CAMstfrac = .false. ! If .true. (.false.), use Slingo (triangular PDF-based) liquid stratus fraction
logical, parameter :: freeze_dry = .false. ! If .true., use 'freeze dry' in liquid stratus fraction formula
real(r8), parameter :: qlst_min = 2.e-5_r8 ! Minimum in-stratus LWC constraint [ kg/kg ]
real(r8), parameter :: qlst_max = 3.e-3_r8 ! Maximum in-stratus LWC constraint [ kg/kg ]
real(r8), parameter :: cc = 0.1_r8 ! For newly formed/dissipated in-stratus CWC ( 0 <= cc <= 1 )
real(r8), parameter :: premib = 700.e2_r8 ! Bottom height for mid-level liquid stratus fraction
real(r8), parameter :: premit = 400.e2_r8 ! Top height for mid-level liquid stratus fraction
integer, parameter :: niter = 2 ! For iterative computation of QQ with 'ramda' below.
real(r8), parameter :: ramda = 0.5_r8 ! Explicit : ramda = 0, Implicit : ramda = 1 ( 0<= ramda <= 1 )
real(r8), private :: rhminl ! Critical RH for low-level liquid stratus clouds
real(r8), private :: rhminh ! Critical RH for high-level liquid stratus clouds
contains
! -------------- !
! Initialization !
! -------------- !
subroutine ini_macro 1,2
!--------------------------------------------------------------------- !
! !
! Purpose: Initialize constants for the liquid stratiform macrophysics !
! called from stratiform.F90. !
! Author: Sungsu Park, Dec.01.2009. !
! !
!--------------------------------------------------------------------- !
use cloud_fraction, only: cldfrc_getparams
call cldfrc_getparams(rhminl_out = rhminl , rhminh_out = rhminh)
if( masterproc ) then
write(iulog,*) 'ini_macro: tuning parameters : rhminl = ', rhminl, &
'rhminh = ', rhminh, &
'premit = ', premit, &
'premib = ', premib
end if
return
end subroutine ini_macro
! ------------------------------ !
! Stratiform Liquid Macrophysics !
! ------------------------------ !
! In the version, 'macro --> micro --> advective forcing --> macro...'
! A_...: only 'advective forcing' without 'microphysical tendency'
! C_...: only 'microphysical tendency'
! D_...: only 'detrainment of cumulus condensate'
! So, 'A' and 'C' are exclusive.
subroutine mmacro_pcond( lchnk , ncol , dt , p , dp , & 1,15
T0 , qv0 , ql0 , qi0 , nl0 , ni0 , &
A_T , A_qv , A_ql , A_qi , A_nl , A_ni , &
C_T , C_qv , C_ql , C_qi , C_nl , C_ni , C_qlst , &
D_T , D_qv , D_ql , D_qi , D_nl , D_ni , &
a_cud , a_cu0 , landfrac , snowh , &
s_tendout , qv_tendout , ql_tendout , qi_tendout , nl_tendout , ni_tendout , &
qme , qvadj , qladj , qiadj , qllim , qilim , &
cld , al_st_star , ai_st_star , ql_st_star , qi_st_star )
use constituents, only : qmin
use time_manager, only : is_first_step, get_nstep
use phys_debug_util, only : phys_debug_col
implicit none
integer icol
integer, intent(in) :: lchnk ! Chunk number
integer, intent(in) :: ncol ! Number of active columns
! Input-Output variables
real(r8), intent(inout) :: T0(pcols,pver) ! Temperature [K]
real(r8), intent(inout) :: qv0(pcols,pver) ! Grid-mean water vapor specific humidity [kg/kg]
real(r8), intent(inout) :: ql0(pcols,pver) ! Grid-mean liquid water content [kg/kg]
real(r8), intent(inout) :: qi0(pcols,pver) ! Grid-mean ice water content [kg/kg]
real(r8), intent(inout) :: nl0(pcols,pver) ! Grid-mean number concentration of cloud liquid droplet [#/kg]
real(r8), intent(inout) :: ni0(pcols,pver) ! Grid-mean number concentration of cloud ice droplet [#/kg]
! Input variables
real(r8), intent(in) :: dt ! Model integration time step [s]
real(r8), intent(in) :: p(pcols,pver) ! Pressure at the layer mid-point [Pa]
real(r8), intent(in) :: dp(pcols,pver) ! Pressure thickness [Pa] > 0
real(r8), intent(in) :: A_T(pcols,pver) ! Non-microphysical advective external forcing of T [K/s]
real(r8), intent(in) :: A_qv(pcols,pver) ! Non-microphysical advective external forcing of qv [kg/kg/s]
real(r8), intent(in) :: A_ql(pcols,pver) ! Non-microphysical advective external forcing of ql [kg/kg/s]
real(r8), intent(in) :: A_qi(pcols,pver) ! Non-microphysical advective external forcing of qi [kg/kg/s]
real(r8), intent(in) :: A_nl(pcols,pver) ! Non-microphysical advective external forcing of nl [#/kg/s]
real(r8), intent(in) :: A_ni(pcols,pver) ! Non-microphysical advective external forcing of ni [#/kg/s]
real(r8), intent(in) :: C_T(pcols,pver) ! Microphysical advective external forcing of T [K/s]
real(r8), intent(in) :: C_qv(pcols,pver) ! Microphysical advective external forcing of qv [kg/kg/s]
real(r8), intent(in) :: C_ql(pcols,pver) ! Microphysical advective external forcing of ql [kg/kg/s]
real(r8), intent(in) :: C_qi(pcols,pver) ! Microphysical advective external forcing of qi [kg/kg/s]
real(r8), intent(in) :: C_nl(pcols,pver) ! Microphysical advective external forcing of nl [#/kg/s]
real(r8), intent(in) :: C_ni(pcols,pver) ! Microphysical advective external forcing of ni [#/kg/s]
real(r8), intent(in) :: C_qlst(pcols,pver) ! Microphysical advective external forcing of ql within liquid stratus [kg/kg/s]
real(r8), intent(in) :: D_T(pcols,pver) ! Cumulus detrainment external forcing of T [K/s]
real(r8), intent(in) :: D_qv(pcols,pver) ! Cumulus detrainment external forcing of qv [kg/kg/s]
real(r8), intent(in) :: D_ql(pcols,pver) ! Cumulus detrainment external forcing of ql [kg/kg/s]
real(r8), intent(in) :: D_qi(pcols,pver) ! Cumulus detrainment external forcing of qi [kg/kg/s]
real(r8), intent(in) :: D_nl(pcols,pver) ! Cumulus detrainment external forcing of nl [#/kg/s]
real(r8), intent(in) :: D_ni(pcols,pver) ! Cumulus detrainment external forcing of qi [#/kg/s]
real(r8), intent(in) :: a_cud(pcols,pver) ! Old cumulus fraction before update
real(r8), intent(in) :: a_cu0(pcols,pver) ! New cumulus fraction after update
real(r8), intent(in) :: landfrac(pcols) ! Land fraction
real(r8), intent(in) :: snowh(pcols) ! Snow depth (liquid water equivalent)
! Output variables
real(r8), intent(out) :: s_tendout(pcols,pver) ! Net tendency of grid-mean s from 'Micro+Macro' processes [J/kg/s]
real(r8), intent(out) :: qv_tendout(pcols,pver) ! Net tendency of grid-mean qv from 'Micro+Macro' processes [kg/kg/s]
real(r8), intent(out) :: ql_tendout(pcols,pver) ! Net tendency of grid-mean ql from 'Micro+Macro' processes [kg/kg/s]
real(r8), intent(out) :: qi_tendout(pcols,pver) ! Net tendency of grid-mean qi from 'Micro+Macro' processes [kg/kg/s]
real(r8), intent(out) :: nl_tendout(pcols,pver) ! Net tendency of grid-mean nl from 'Micro+Macro' processes [#/kg/s]
real(r8), intent(out) :: ni_tendout(pcols,pver) ! Net tendency of grid-mean ni from 'Micro+Macro' processes [#/kg/s]
real(r8), intent(out) :: qme (pcols,pver) ! Net condensation rate [kg/kg/s]
real(r8), intent(out) :: qvadj(pcols,pver) ! adjustment tendency from "positive_moisture" call (vapor)
real(r8), intent(out) :: qladj(pcols,pver) ! adjustment tendency from "positive_moisture" call (liquid)
real(r8), intent(out) :: qiadj(pcols,pver) ! adjustment tendency from "positive_moisture" call (ice)
real(r8), intent(out) :: qllim(pcols,pver) ! tendency from "instratus_condensate" call (liquid)
real(r8), intent(out) :: qilim(pcols,pver) ! tendency from "instratus_condensate" call (ice)
real(r8), intent(out) :: cld(pcols,pver) ! Net cloud fraction ( 0 <= cld <= 1 )
real(r8), intent(out) :: al_st_star(pcols,pver) ! Physical liquid stratus fraction
real(r8), intent(out) :: ai_st_star(pcols,pver) ! Physical ice stratus fraction
real(r8), intent(out) :: ql_st_star(pcols,pver) ! In-stratus LWC [kg/kg]
real(r8), intent(out) :: qi_st_star(pcols,pver) ! In-stratus IWC [kg/kg]
! --------------- !
! Local variables !
! --------------- !
integer :: i, j, k, iter, ii, jj ! Loop indexes
! Thermodynamic state variables
real(r8) T(pcols,pver) ! Temperature of equilibrium reference state from which 'Micro & Macro' are computed [K]
real(r8) T1(pcols,pver) ! Temperature after 'fice_force' on T01
real(r8) T_0(pcols,pver) ! Temperature after 'instratus_condensate' on T1
real(r8) T_05(pcols,pver) ! Temperature after 'advection' on T_0
real(r8) T_prime0(pcols,pver) ! Temperature after 'Macrophysics (QQ)' on T_05star
real(r8) T_dprime(pcols,pver) ! Temperature after 'fice_force' on T_prime
real(r8) T_star(pcols,pver) ! Temperature after 'instratus_condensate' on T_dprime
real(r8) qv(pcols,pver) ! Grid-mean qv of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
real(r8) qv1(pcols,pver) ! Grid-mean qv after 'fice_force' on qv01
real(r8) qv_0(pcols,pver) ! Grid-mean qv after 'instratus_condensate' on qv1
real(r8) qv_05(pcols,pver) ! Grid-mean qv after 'advection' on qv_0
real(r8) qv_prime0(pcols,pver) ! Grid-mean qv after 'Macrophysics (QQ)' on qv_05star
real(r8) qv_dprime(pcols,pver) ! Grid-mean qv after 'fice_force' on qv_prime
real(r8) qv_star(pcols,pver) ! Grid-mean qv after 'instratus_condensate' on qv_dprime
real(r8) ql(pcols,pver) ! Grid-mean ql of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
real(r8) ql1(pcols,pver) ! Grid-mean ql after 'fice_force' on ql01
real(r8) ql_0(pcols,pver) ! Grid-mean ql after 'instratus_condensate' on ql1
real(r8) ql_05(pcols,pver) ! Grid-mean ql after 'advection' on ql_0
real(r8) ql_prime0(pcols,pver) ! Grid-mean ql after 'Macrophysics (QQ)' on ql_05star
real(r8) ql_dprime(pcols,pver) ! Grid-mean ql after 'fice_force' on ql_prime
real(r8) ql_star(pcols,pver) ! Grid-mean ql after 'instratus_condensate' on ql_dprime
real(r8) qi(pcols,pver) ! Grid-mean qi of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
real(r8) qi1(pcols,pver) ! Grid-mean qi after 'fice_force' on qi01
real(r8) qi_0(pcols,pver) ! Grid-mean qi after 'instratus_condensate' on qi1
real(r8) qi_05(pcols,pver) ! Grid-mean qi after 'advection' on qi_0
real(r8) qi_prime0(pcols,pver) ! Grid-mean qi after 'Macrophysics (QQ)' on qi_05star
real(r8) qi_dprime(pcols,pver) ! Grid-mean qi after 'fice_force' on qi_prime
real(r8) qi_star(pcols,pver) ! Grid-mean qi after 'instratus_condensate' on qi_dprime
real(r8) nl(pcols,pver) ! Grid-mean nl of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
real(r8) nl1(pcols,pver) ! Grid-mean nl after 'fice_force' on nl01
real(r8) nl_0(pcols,pver) ! Grid-mean nl after 'instratus_condensate' on nl1
real(r8) nl_05(pcols,pver) ! Grid-mean nl after 'advection' on nl_0
real(r8) nl_prime0(pcols,pver) ! Grid-mean nl after 'Macrophysics (QQ)' on nl_05star
real(r8) nl_dprime(pcols,pver) ! Grid-mean nl after 'fice_force' on nl_prime
real(r8) nl_star(pcols,pver) ! Grid-mean nl after 'instratus_condensate' on nl_dprime
real(r8) ni(pcols,pver) ! Grid-mean ni of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
real(r8) ni1(pcols,pver) ! Grid-mean ni after 'fice_force' on ni01
real(r8) ni_0(pcols,pver) ! Grid-mean ni after 'instratus_condensate' on ni1
real(r8) ni_05(pcols,pver) ! Grid-mean ni after 'advection' on ni_0
real(r8) ni_prime0(pcols,pver) ! Grid-mean ni after 'Macrophysics (QQ)' on ni_05star
real(r8) ni_dprime(pcols,pver) ! Grid-mean ni after 'fice_force' on ni_prime
real(r8) ni_star(pcols,pver) ! Grid-mean ni after 'instratus_condensate' on ni_dprime
real(r8) a_st(pcols,pver) ! Stratus fraction of equilibrium reference state
real(r8) a_st_0(pcols,pver) ! Stratus fraction at '_0' state
real(r8) a_st_star(pcols,pver) ! Stratus fraction at '_star' state
real(r8) al_st(pcols,pver) ! Liquid stratus fraction of equilibrium reference state
real(r8) al_st_0(pcols,pver) ! Liquid stratus fraction at '_0' state
real(r8) al_st_nc(pcols,pver) ! Non-physical liquid stratus fraction in the non-cumulus pixels
real(r8) ai_st(pcols,pver) ! Ice stratus fraction of equilibrium reference state
real(r8) ai_st_0(pcols,pver) ! Ice stratus fraction at '_0' state
real(r8) ai_st_nc(pcols,pver) ! Non-physical ice stratus fraction in the non-cumulus pixels
real(r8) ql_st(pcols,pver) ! In-stratus LWC of equilibrium reference state [kg/kg]
real(r8) ql_st_0(pcols,pver) ! In-stratus LWC at '_0' state
real(r8) qi_st(pcols,pver) ! In-stratus IWC of equilibrium reference state [kg/kg]
real(r8) qi_st_0(pcols,pver) ! In-stratus IWC at '_0' state
! Cumulus properties
real(r8) dacudt(pcols,pver)
real(r8) a_cu(pcols,pver)
! Adjustment tendency in association with 'positive_moisture'
real(r8) Tten_pwi1(pcols,pver) ! Pre-process T tendency of input equilibrium state [K/s]
real(r8) qvten_pwi1(pcols,pver) ! Pre-process qv tendency of input equilibrium state [kg/kg/s]
real(r8) qlten_pwi1(pcols,pver) ! Pre-process ql tendency of input equilibrium state [kg/kg/s]
real(r8) qiten_pwi1(pcols,pver) ! Pre-process qi tendency of input equilibrium state [kg/kg/s]
real(r8) nlten_pwi1(pcols,pver) ! Pre-process nl tendency of input equilibrium state [#/kg/s]
real(r8) niten_pwi1(pcols,pver) ! Pre-process ni tendency of input equilibrium state [#/kg/s]
real(r8) Tten_pwi2(pcols,pver) ! Post-process T tendency of provisional equilibrium state [K/s]
real(r8) qvten_pwi2(pcols,pver) ! Post-process qv tendency of provisional equilibrium state [kg/kg/s]
real(r8) qlten_pwi2(pcols,pver) ! Post-process ql tendency of provisional equilibrium state [kg/kg/s]
real(r8) qiten_pwi2(pcols,pver) ! Post-process qi tendency of provisional equilibrium state [kg/kg/s]
real(r8) nlten_pwi2(pcols,pver) ! Post-process nl tendency of provisoonal equilibrium state [#/kg/s]
real(r8) niten_pwi2(pcols,pver) ! Post-process ni tendency of provisional equilibrium state [#/kg/s]
real(r8) A_T_adj(pcols,pver) ! After applying external advective forcing [K/s]
real(r8) A_qv_adj(pcols,pver) ! After applying external advective forcing [kg/kg/s]
real(r8) A_ql_adj(pcols,pver) ! After applying external advective forcing [kg/kg/s]
real(r8) A_qi_adj(pcols,pver) ! After applying external advective forcing [kg/kg/s]
real(r8) A_nl_adj(pcols,pver) ! After applying external advective forcing [#/kg/s]
real(r8) A_ni_adj(pcols,pver) ! After applying external advective forcing [#/kg/s]
! Adjustment tendency in association with 'instratus_condensate'
real(r8) QQw1(pcols,pver) ! Effective adjustive condensation into water due to 'instratus_condensate' [kg/kg/s]
real(r8) QQi1(pcols,pver) ! Effective adjustive condensation into ice due to 'instratus_condensate' [kg/kg/s]
real(r8) QQw2(pcols,pver) ! Effective adjustive condensation into water due to 'instratus_condensate' [kg/kg/s]
real(r8) QQi2(pcols,pver) ! Effective adjustive condensation into ice due to 'instratus_condensate' [kg/kg/s]
real(r8) QQnl1(pcols,pver) ! Tendency of nl associated with QQw1 only when QQw1<0 (net evaporation) [#/kg/s]
real(r8) QQni1(pcols,pver) ! Tendency of ni associated with QQi1 only when QQw1<0 (net evaporation) [#/kg/s]
real(r8) QQnl2(pcols,pver) ! Tendency of nl associated with QQw2 only when QQw2<0 (net evaporation) [#/kg/s]
real(r8) QQni2(pcols,pver) ! Tendency of ni associated with QQi2 only when QQw2<0 (net evaporation) [#/kg/s]
! Macrophysical process tendency variables
real(r8) QQ(pcols,pver) ! Net condensation rate into water+ice [kg/kg/s]
real(r8) QQw(pcols,pver) ! Net condensation rate into water [kg/kg/s]
real(r8) QQi(pcols,pver) ! Net condensation rate into ice [kg/kg/s]
real(r8) QQnl(pcols,pver) ! Tendency of nl associated with QQw both for condensation and evaporation [#/kg/s]
real(r8) QQni(pcols,pver) ! Tendency of ni associated with QQi both for condensation and evaporation [#/kg/s]
real(r8) ACnl(pcols,pver) ! Cloud liquid droplet (nl) activation tendency [#/kg/s]
real(r8) ACni(pcols,pver) ! Cloud ice droplet (ni) activation tendency [#/kg/s]
real(r8) QQw_prev(pcols,pver)
real(r8) QQi_prev(pcols,pver)
real(r8) QQnl_prev(pcols,pver)
real(r8) QQni_prev(pcols,pver)
real(r8) QQw_prog(pcols,pver)
real(r8) QQi_prog(pcols,pver)
real(r8) QQnl_prog(pcols,pver)
real(r8) QQni_prog(pcols,pver)
real(r8) QQ_final(pcols,pver)
real(r8) QQw_final(pcols,pver)
real(r8) QQi_final(pcols,pver)
real(r8) QQn_final(pcols,pver)
real(r8) QQnl_final(pcols,pver)
real(r8) QQni_final(pcols,pver)
real(r8) QQ_all(pcols,pver) ! QQw_all + QQi_all
real(r8) QQw_all(pcols,pver) ! QQw_final + QQw1 + QQw2 + qlten_pwi1 + qlten_pwi2 + A_ql_adj [kg/kg/s]
real(r8) QQi_all(pcols,pver) ! QQi_final + QQi1 + QQi2 + qiten_pwi1 + qiten_pwi2 + A_qi_adj [kg/kg/s]
real(r8) QQn_all(pcols,pver) ! QQnl_all + QQni_all
real(r8) QQnl_all(pcols,pver) ! QQnl_final + QQnl1 + QQnl2 + nlten_pwi1 + nlten_pwi2 + ACnl [#/kg/s]
real(r8) QQni_all(pcols,pver) ! QQni_final + QQni1 + QQni2 + niten_pwi1 + niten_pwi2 + ACni [#/kg/s]
! Coefficient for computing QQ and related processes
real(r8) U(pcols,pver) ! Grid-mean RH
real(r8) U_nc(pcols,pver) ! Mean RH of non-cumulus pixels
real(r8) G_nc(pcols,pver) ! d(U_nc)/d(a_st_nc)
real(r8) F_nc(pcols,pver) ! A function of second parameter for a_st_nc
real(r8) alpha ! = 1/qsat
real(r8) beta ! = (qv/qsat**2)*dqsdT
real(r8) betast ! = alpha*dqsdT
real(r8) gammal ! = alpha + (hlatv/cp)*beta
real(r8) gammai ! = alpha + ((hlatv+hlatf)/cp)*beta
real(r8) gammaQ ! = alpha + (hlatv/cp)*beta
real(r8) deltal ! = 1 + a_st*(hlatv/cp)*(betast/alpha)
real(r8) deltai ! = 1 + a_st*((hlatv+hlatf)/cp)*(betast/alpha)
real(r8) A_Tc ! Advective external forcing of Tc [K/s]
real(r8) A_qt ! Advective external forcing of qt [kg/kg/s]
real(r8) C_Tc ! Microphysical forcing of Tc [K/s]
real(r8) C_qt ! Microphysical forcing of qt [kg/kg/s]
real(r8) dTcdt ! d(Tc)/dt [K/s]
real(r8) dqtdt ! d(qt)/dt [kg/kg/s]
real(r8) dqtstldt ! d(qt_alst)/dt [kg/kg/s]
real(r8) dqidt ! d(qi)/dt [kg/kg/s]
real(r8) dqlstdt ! d(ql_st)/dt [kg/kg/s]
real(r8) dalstdt ! d(al_st)/dt [1/s]
real(r8) dastdt ! d(a_st)/dt [1/s]
real(r8) anic ! Fractional area of non-cumulus and non-ice stratus fraction
real(r8) GG ! G_nc(i,k)/anic
real(r8) aa(2,2)
real(r8) bb(2,1)
real(r8) zeros(pcols,pver)
real(r8) qmin1(pcols,pver)
real(r8) qmin2(pcols,pver)
real(r8) qmin3(pcols,pver)
real(r8) esat_a(pcols) ! Saturation water vapor pressure [Pa]
real(r8) qsat_a(pcols) ! Saturation water vapor specific humidity [kg/kg]
real(r8) dqsdT_a(pcols) ! dqsat/dT [kg/kg/K]
real(r8) Twb_aw(pcols,pver) ! Wet-bulb temperature [K]
real(r8) qvwb_aw(pcols,pver) ! Wet-bulb water vapor specific humidity [kg/kg]
real(r8) esat
real(r8) qsat
real(r8) dqsdT
real(r8) esat_b(pcols)
real(r8) qsat_b(pcols)
real(r8) dqsdT_b(pcols)
real(r8) QQmax,QQmin,QQwmin,QQimin ! For limiting QQ
real(r8) cone ! Number close to but smaller than 1
real(r8) qsmall ! Smallest mixing ratio considered in the macrophysics
cone = 0.999_r8
qsmall = 1.e-18_r8
zeros(:ncol,:) = 0._r8
! ------------------------------------ !
! Global initialization of main output !
! ------------------------------------ !
s_tendout(:ncol,:) = 0._r8
qv_tendout(:ncol,:) = 0._r8
ql_tendout(:ncol,:) = 0._r8
qi_tendout(:ncol,:) = 0._r8
nl_tendout(:ncol,:) = 0._r8
ni_tendout(:ncol,:) = 0._r8
qme(:ncol,:) = 0._r8
cld(:ncol,:) = 0._r8
al_st_star(:ncol,:) = 0._r8
ai_st_star(:ncol,:) = 0._r8
ql_st_star(:ncol,:) = 0._r8
qi_st_star(:ncol,:) = 0._r8
! --------------------------------------- !
! Initialization of internal 2D variables !
! --------------------------------------- !
T(:ncol,:) = 0._r8
T1(:ncol,:) = 0._r8
T_0(:ncol,:) = 0._r8
T_05(:ncol,:) = 0._r8
T_prime0(:ncol,:) = 0._r8
T_dprime(:ncol,:) = 0._r8
T_star(:ncol,:) = 0._r8
qv(:ncol,:) = 0._r8
qv1(:ncol,:) = 0._r8
qv_0(:ncol,:) = 0._r8
qv_05(:ncol,:) = 0._r8
qv_prime0(:ncol,:) = 0._r8
qv_dprime(:ncol,:) = 0._r8
qv_star(:ncol,:) = 0._r8
ql(:ncol,:) = 0._r8
ql1(:ncol,:) = 0._r8
ql_0(:ncol,:) = 0._r8
ql_05(:ncol,:) = 0._r8
ql_prime0(:ncol,:) = 0._r8
ql_dprime(:ncol,:) = 0._r8
ql_star(:ncol,:) = 0._r8
qi(:ncol,:) = 0._r8
qi1(:ncol,:) = 0._r8
qi_0(:ncol,:) = 0._r8
qi_05(:ncol,:) = 0._r8
qi_prime0(:ncol,:) = 0._r8
qi_dprime(:ncol,:) = 0._r8
qi_star(:ncol,:) = 0._r8
nl(:ncol,:) = 0._r8
nl1(:ncol,:) = 0._r8
nl_0(:ncol,:) = 0._r8
nl_05(:ncol,:) = 0._r8
nl_prime0(:ncol,:) = 0._r8
nl_dprime(:ncol,:) = 0._r8
nl_star(:ncol,:) = 0._r8
ni(:ncol,:) = 0._r8
ni1(:ncol,:) = 0._r8
ni_0(:ncol,:) = 0._r8
ni_05(:ncol,:) = 0._r8
ni_prime0(:ncol,:) = 0._r8
ni_dprime(:ncol,:) = 0._r8
ni_star(:ncol,:) = 0._r8
a_st(:ncol,:) = 0._r8
a_st_0(:ncol,:) = 0._r8
a_st_star(:ncol,:) = 0._r8
al_st(:ncol,:) = 0._r8
al_st_0(:ncol,:) = 0._r8
al_st_nc(:ncol,:) = 0._r8
ai_st(:ncol,:) = 0._r8
ai_st_0(:ncol,:) = 0._r8
ai_st_nc(:ncol,:) = 0._r8
ql_st(:ncol,:) = 0._r8
ql_st_0(:ncol,:) = 0._r8
qi_st(:ncol,:) = 0._r8
qi_st_0(:ncol,:) = 0._r8
! Cumulus properties
dacudt(:ncol,:) = 0._r8
a_cu(:ncol,:) = 0._r8
! Adjustment tendency in association with 'positive_moisture'
Tten_pwi1(:ncol,:) = 0._r8
qvten_pwi1(:ncol,:) = 0._r8
qlten_pwi1(:ncol,:) = 0._r8
qiten_pwi1(:ncol,:) = 0._r8
nlten_pwi1(:ncol,:) = 0._r8
niten_pwi1(:ncol,:) = 0._r8
Tten_pwi2(:ncol,:) = 0._r8
qvten_pwi2(:ncol,:) = 0._r8
qlten_pwi2(:ncol,:) = 0._r8
qiten_pwi2(:ncol,:) = 0._r8
nlten_pwi2(:ncol,:) = 0._r8
niten_pwi2(:ncol,:) = 0._r8
A_T_adj(:ncol,:) = 0._r8
A_qv_adj(:ncol,:) = 0._r8
A_ql_adj(:ncol,:) = 0._r8
A_qi_adj(:ncol,:) = 0._r8
A_nl_adj(:ncol,:) = 0._r8
A_ni_adj(:ncol,:) = 0._r8
qvadj (:ncol,:) = 0._r8
qladj (:ncol,:) = 0._r8
qiadj (:ncol,:) = 0._r8
! Adjustment tendency in association with 'instratus_condensate'
QQw1(:ncol,:) = 0._r8
QQi1(:ncol,:) = 0._r8
QQw2(:ncol,:) = 0._r8
QQi2(:ncol,:) = 0._r8
QQnl1(:ncol,:) = 0._r8
QQni1(:ncol,:) = 0._r8
QQnl2(:ncol,:) = 0._r8
QQni2(:ncol,:) = 0._r8
QQnl(:ncol,:) = 0._r8
QQni(:ncol,:) = 0._r8
! Macrophysical process tendency variables
QQ(:ncol,:) = 0._r8
QQw(:ncol,:) = 0._r8
QQi(:ncol,:) = 0._r8
QQnl(:ncol,:) = 0._r8
QQni(:ncol,:) = 0._r8
ACnl(:ncol,:) = 0._r8
ACni(:ncol,:) = 0._r8
QQw_prev(:ncol,:) = 0._r8
QQi_prev(:ncol,:) = 0._r8
QQnl_prev(:ncol,:) = 0._r8
QQni_prev(:ncol,:) = 0._r8
QQw_prog(:ncol,:) = 0._r8
QQi_prog(:ncol,:) = 0._r8
QQnl_prog(:ncol,:) = 0._r8
QQni_prog(:ncol,:) = 0._r8
QQ_final(:ncol,:) = 0._r8
QQw_final(:ncol,:) = 0._r8
QQi_final(:ncol,:) = 0._r8
QQn_final(:ncol,:) = 0._r8
QQnl_final(:ncol,:) = 0._r8
QQni_final(:ncol,:) = 0._r8
QQ_all(:ncol,:) = 0._r8
QQw_all(:ncol,:) = 0._r8
QQi_all(:ncol,:) = 0._r8
QQn_all(:ncol,:) = 0._r8
QQnl_all(:ncol,:) = 0._r8
QQni_all(:ncol,:) = 0._r8
! Coefficient for computing QQ and related processes
U(:ncol,:) = 0._r8
U_nc(:ncol,:) = 0._r8
G_nc(:ncol,:) = 0._r8
F_nc(:ncol,:) = 0._r8
! Other
qmin1(:ncol,:) = 0._r8
qmin2(:ncol,:) = 0._r8
qmin3(:ncol,:) = 0._r8
esat_b(:ncol) = 0._r8
qsat_b(:ncol) = 0._r8
dqsdT_b(:ncol) = 0._r8
esat_a(:ncol) = 0._r8
qsat_a(:ncol) = 0._r8
dqsdT_a(:ncol) = 0._r8
Twb_aw(:ncol,:pver) = 0._r8
qvwb_aw(:ncol,:pver) = 0._r8
! ---------------- !
! Main computation !
! ---------------- !
! ---------------------------------- !
! Compute cumulus-related properties !
! ---------------------------------- !
dacudt(:ncol,:) = (a_cu0(:ncol,:) - a_cud(:ncol,:))/dt
! ---------------------------------------------------------------------- !
! Check if input non-cumulus pixels satisfie a non-negative constraint. !
! If not, force all water vapor substances to be positive in all layers. !
! We should use 'old' cumulus properties for this routine. !
! ---------------------------------------------------------------------- !
T1(:ncol,:) = T0(:ncol,:)
qv1(:ncol,:) = qv0(:ncol,:)
ql1(:ncol,:) = ql0(:ncol,:)
qi1(:ncol,:) = qi0(:ncol,:)
nl1(:ncol,:) = nl0(:ncol,:)
ni1(:ncol,:) = ni0(:ncol,:)
qmin1(:ncol,:) = qmin(1)
qmin2(:ncol,:) = qmin(2)
qmin3(:ncol,:) = qmin(3)
call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, &
qv1, ql1, qi1, T1, qvten_pwi1, qlten_pwi1, qiten_pwi1, Tten_pwi1 )
do k = 1, pver
do i = 1, ncol
if( ql1(i,k) .lt. qsmall ) then
nlten_pwi1(i,k) = -nl1(i,k)/dt
nl1(i,k) = 0._r8
endif
if( qi1(i,k) .lt. qsmall ) then
niten_pwi1(i,k) = -ni1(i,k)/dt
ni1(i,k) = 0._r8
endif
enddo
enddo
! ------------------------------------------------------------- !
! Impose 'in-stratus condensate amount constraint' !
! such that it is bounded by two limiting values. !
! This should also use 'old' cumulus properties since it is !
! before applying external forcings. !
! Below 'QQw1,QQi1' are effective adjustive condensation !
! Although this process also involves freezing of cloud !
! liquid into ice, they can be and only can be expressed !
! in terms of effective condensation. !
! ------------------------------------------------------------- !
do k = 1, pver
call instratus_condensate( lchnk, ncol, k, &
p(:,k), T1(:,k), qv1(:,k), ql1(:,k), qi1(:,k), &
a_cud(:,k), zeros(:,k), zeros(:,k), &
zeros(:,k), zeros(:,k), zeros(:,k), &
landfrac, snowh, &
T_0(:,k), qv_0(:,k), ql_0(:,k), qi_0(:,k), &
al_st_0(:,k), ai_st_0(:,k), ql_st_0(:,k), qi_st_0(:,k) )
a_st_0(:ncol,k) = max(al_st_0(:ncol,k),ai_st_0(:ncol,k))
QQw1(:ncol,k) = (ql_0(:ncol,k) - ql1(:ncol,k))/dt
QQi1(:ncol,k) = (qi_0(:ncol,k) - qi1(:ncol,k))/dt
! -------------------------------------------------- !
! Reduce droplet concentration if evaporation occurs !
! Set a limit such that negative state not happens. !
! -------------------------------------------------- !
do i = 1, ncol
if( QQw1(i,k) .le. 0._r8 ) then
if( ql1(i,k) .gt. qsmall ) then
QQnl1(i,k) = QQw1(i,k)*nl1(i,k)/ql1(i,k)
QQnl1(i,k) = min(0._r8,cone*max(QQnl1(i,k),-nl1(i,k)/dt))
else
QQnl1(i,k) = 0._r8
endif
endif
if( QQi1(i,k) .le. 0._r8 ) then
if( qi1(i,k) .gt. qsmall ) then
QQni1(i,k) = QQi1(i,k)*ni1(i,k)/qi1(i,k)
QQni1(i,k) = min(0._r8,cone*max(QQni1(i,k),-ni1(i,k)/dt))
else
QQni1(i,k) = 0._r8
endif
endif
enddo
enddo
nl_0(:ncol,:) = max(0._r8,nl1(:ncol,:)+QQnl1(:ncol,:)*dt)
ni_0(:ncol,:) = max(0._r8,ni1(:ncol,:)+QQni1(:ncol,:)*dt)
! ----------------------------------------------------------------------------- !
! Check if non-cumulus pixels of '_05' state satisfies non-negative constraint. !
! If not, force all water substances of '_05' state to be positive by imposing !
! adjustive advection. We should use 'new' cumulus properties for this routine. !
! ----------------------------------------------------------------------------- !
T_05(:ncol,:) = T_0(:ncol,:) + ( A_T(:ncol,:) + C_T(:ncol,:) ) * dt
qv_05(:ncol,:) = qv_0(:ncol,:) + ( A_qv(:ncol,:) + C_qv(:ncol,:) ) * dt
ql_05(:ncol,:) = ql_0(:ncol,:) + ( A_ql(:ncol,:) + C_ql(:ncol,:) ) * dt
qi_05(:ncol,:) = qi_0(:ncol,:) + ( A_qi(:ncol,:) + C_qi(:ncol,:) ) * dt
nl_05(:ncol,:) = max(0._r8, nl_0(:ncol,:) + ( A_nl(:ncol,:) + C_nl(:ncol,:) ) * dt )
ni_05(:ncol,:) = max(0._r8, ni_0(:ncol,:) + ( A_ni(:ncol,:) + C_ni(:ncol,:) ) * dt )
call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, &
qv_05, ql_05, qi_05, T_05, A_qv_adj, A_ql_adj, A_qi_adj, A_T_adj )
! -------------------------------------------------------------- !
! Define reference state at the first iteration. This will be !
! continuously updated within the iteration loop below. !
! While equlibrium state properties are already output from the !
! 'instratus_condensate', they will be re-computed within the !
! each iteration process. At the first iteration, they will !
! produce exactly identical results. Note that except at the !
! very first iteration iter = 1, we must use updated cumulus !
! properties at all the other iteration processes. Even at the !
! first iteration, we should use updated cumulus properties !
! when computing limiters for (Q,P,E). !
! -------------------------------------------------------------- !
! -------------------------------------------------------------- !
! Define variables at the reference state of the first iteration !
! -------------------------------------------------------------- !
T(:ncol,:) = T_0(:ncol,:)
qv(:ncol,:) = qv_0(:ncol,:)
ql(:ncol,:) = ql_0(:ncol,:)
qi(:ncol,:) = qi_0(:ncol,:)
al_st(:ncol,:) = al_st_0(:ncol,:)
ai_st(:ncol,:) = ai_st_0(:ncol,:)
a_st(:ncol,:) = a_st_0(:ncol,:)
ql_st(:ncol,:) = ql_st_0(:ncol,:)
qi_st(:ncol,:) = qi_st_0(:ncol,:)
nl(:ncol,:) = nl_0(:ncol,:)
ni(:ncol,:) = ni_0(:ncol,:)
! -------------------------- !
! Main iterative computation !
! -------------------------- !
do iter = 1, niter
! ------------------------------------------ !
! Initialize array within the iteration loop !
! ------------------------------------------ !
QQ(:,:) = 0._r8
QQw(:,:) = 0._r8
QQi(:,:) = 0._r8
QQnl(:,:) = 0._r8
QQni(:,:) = 0._r8
QQw2(:,:) = 0._r8
QQi2(:,:) = 0._r8
QQnl2(:,:) = 0._r8
QQni2(:,:) = 0._r8
nlten_pwi2(:,:) = 0._r8
niten_pwi2(:,:) = 0._r8
ACnl(:,:) = 0._r8
ACni(:,:) = 0._r8
aa(:,:) = 0._r8
bb(:,:) = 0._r8
call findsp_water(lchnk,ncol,qv_05(:,:),T_05(:,:),p(:,:),Twb_aw(:,:),qvwb_aw(:,:))
do k = 1, pver
call vqsatd2_water(T_05(1:ncol,k),p(1:ncol,k),esat_a(1:ncol),qsat_a(1:ncol),dqsdT_a(1:ncol),ncol)
call vqsatd2_water(T(1:ncol,k),p(1:ncol,k),esat_b(1:ncol),qsat_b(1:ncol),dqsdT_b(1:ncol),ncol)
if( iter .eq. 1 ) then
a_cu(:ncol,k) = a_cud(:ncol,k)
else
a_cu(:ncol,k) = a_cu0(:ncol,k)
endif
do i = 1, ncol
U(i,k) = qv(i,k)/qsat_b(i)
U_nc(i,k) = U(i,k)
enddo
if( CAMstfrac ) then
call astG_RHU(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol)
else
call astG_PDF(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol)
endif
call aist_vector(qv(:,k),T(:,k),p(:,k),qi(:,k),landfrac(:),snowh(:),ai_st_nc(:,k),ncol)
ai_st(:ncol,k) = (1._r8-a_cu(:ncol,k))*ai_st_nc(:ncol,k)
al_st(:ncol,k) = (1._r8-a_cu(:ncol,k))*al_st_nc(:ncol,k)
a_st(:ncol,k) = max(al_st(:ncol,k),ai_st(:ncol,k))
do i = 1, ncol
! -------------------------------------------------------- !
! Compute basic thermodynamic coefficients for computing Q !
! -------------------------------------------------------- !
alpha = 1._r8/qsat_b(i)
beta = dqsdt_b(i)*(qv(i,k)/qsat_b(i)**2)
betast = alpha*dqsdT_b(i)
gammal = alpha + (hlatv/cp)*beta
gammai = alpha + ((hlatv+hlatf)/cp)*beta
gammaQ = alpha + (hlatv/cp)*beta
deltal = 1._r8 + a_st(i,k)*(hlatv/cp)*(betast/alpha)
deltai = 1._r8 + a_st(i,k)*((hlatv+hlatf)/cp)*(betast/alpha)
A_Tc = A_T(i,k)+A_T_adj(i,k)-(hlatv/cp)*(A_ql(i,k)+A_ql_adj(i,k))-((hlatv+hlatf)/cp)*(A_qi(i,k)+A_qi_adj(i,k))
A_qt = A_qv(i,k) + A_qv_adj(i,k) + A_ql(i,k) + A_ql_adj(i,k) + A_qi(i,k) + A_qi_adj(i,k)
C_Tc = C_T(i,k) - (hlatv/cp)*C_ql(i,k) - ((hlatv+hlatf)/cp)*C_qi(i,k)
C_qt = C_qv(i,k) + C_ql(i,k) + C_qi(i,k)
dTcdt = A_Tc + C_Tc
dqtdt = A_qt + C_qt
! dqtstldt = A_qt + C_ql(i,k)/max(1.e-2_r8,al_st(i,k)) ! Original
! dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_ql(i,k)/max(1.e-2_r8,al_st(i,k)) ! New 1 on Dec.30.2009.
dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_qlst(i,k) ! New 2 on Dec.30.2009.
! dqtstldt = A_qt + C_qt ! Original Conservative treatment
! dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_qt - C_qi(i,k) ! New Conservative treatment on Dec.30.2009
dqidt = A_qi(i,k) + A_qi_adj(i,k) + C_qi(i,k)
anic = max(1.e-8_r8,(1._r8-a_cu(i,k)))
GG = G_nc(i,k)/anic
aa(1,1) = gammal*al_st(i,k)
aa(1,2) = GG + gammal*cc*ql_st(i,k)
aa(2,1) = alpha + (hlatv/cp)*betast*al_st(i,k)
aa(2,2) = (hlatv/cp)*betast*cc*ql_st(i,k)
bb(1,1) = alpha*dqtdt - beta*dTcdt - gammai*dqidt - GG*al_st_nc(i,k)*dacudt(i,k) + F_nc(i,k)
bb(2,1) = alpha*dqtstldt - betast*(dTcdt + ((hlatv+hlatf)/cp)*dqidt)
call gaussj(aa(1:2,1:2),2,2,bb(1:2,1),1,1)
dqlstdt = bb(1,1)
dalstdt = bb(2,1)
QQ(i,k) = al_st(i,k)*dqlstdt + cc*ql_st(i,k)*dalstdt - ( A_ql(i,k) + A_ql_adj(i,k) + C_ql(i,k) )
! ------------------------------------------------------------ !
! Limiter for QQ !
! Here, 'fice' should be from the reference equilibrium state !
! since QQ itself is computed from the reference state. !
! From the assumption used for derivation of QQ(i), it must be !
! that QQw(i) = QQ(i)*(1._r8-fice(i)), QQi(i) = QQ(i)*fice(i) !
! ------------------------------------------------------------ !
if( QQ(i,k) .ge. 0._r8 ) then
QQmax = (qv_05(i,k) - qmin(1))/dt ! For ghost cumulus & semi-ghost ice stratus
QQmax = max(0._r8,QQmax)
QQ(i,k) = min(QQ(i,k),QQmax)
QQw(i,k) = QQ(i,k)
QQi(i,k) = 0._r8
else
QQmin = 0._r8
if( qv_05(i,k) .lt. qsat_a(i) ) QQmin = min(0._r8,cone*(qv_05(i,k)-qvwb_aw(i,k))/dt)
QQ(i,k) = max(QQ(i,k),QQmin)
QQw(i,k) = QQ(i,k)
QQi(i,k) = 0._r8
QQwmin = min(0._r8,-cone*ql_05(i,k)/dt)
QQimin = min(0._r8,-cone*qi_05(i,k)/dt)
QQw(i,k) = min(0._r8,max(QQw(i,k),QQwmin))
QQi(i,k) = min(0._r8,max(QQi(i,k),QQimin))
endif
! -------------------------------------------------- !
! Reduce droplet concentration if evaporation occurs !
! Note 'QQnl1,QQni1' are computed from the reference !
! equilibrium state but limiter is from 'nl_05'. !
! -------------------------------------------------- !
if( QQw(i,k) .lt. 0._r8 ) then
if( ql_05(i,k) .gt. qsmall ) then
QQnl(i,k) = QQw(i,k)*nl_05(i,k)/ql_05(i,k)
QQnl(i,k) = min(0._r8,cone*max(QQnl(i,k),-nl_05(i,k)/dt))
else
QQnl(i,k) = 0._r8
endif
endif
if( QQi(i,k) .lt. 0._r8 ) then
if( qi_05(i,k) .gt. qsmall ) then
QQni(i,k) = QQi(i,k)*ni_05(i,k)/qi_05(i,k)
QQni(i,k) = min(0._r8,cone*max(QQni(i,k),-ni_05(i,k)/dt))
else
QQni(i,k) = 0._r8
endif
endif
enddo
enddo
! -------------------------------------------------------------------- !
! Until now, we have finished computing all necessary tendencies !
! from the equilibrium input state (T_0). !
! If ramda = 0 : fully explicit scheme !
! ramda = 1 : fully implicit scheme !
! Note that 'ramda = 0.5 with niter = 2' can mimic !
! -------------------------------------------------------------------- !
if( iter .eq. 1 ) then
QQw_prev(:ncol,:) = QQw(:ncol,:)
QQi_prev(:ncol,:) = QQi(:ncol,:)
QQnl_prev(:ncol,:) = QQnl(:ncol,:)
QQni_prev(:ncol,:) = QQni(:ncol,:)
endif
QQw_prog(:ncol,:) = ramda*QQw(:ncol,:) + (1._r8-ramda)*QQw_prev(:ncol,:)
QQi_prog(:ncol,:) = ramda*QQi(:ncol,:) + (1._r8-ramda)*QQi_prev(:ncol,:)
QQnl_prog(:ncol,:) = ramda*QQnl(:ncol,:) + (1._r8-ramda)*QQnl_prev(:ncol,:)
QQni_prog(:ncol,:) = ramda*QQni(:ncol,:) + (1._r8-ramda)*QQni_prev(:ncol,:)
QQw_prev(:ncol,:) = QQw_prog(:ncol,:)
QQi_prev(:ncol,:) = QQi_prog(:ncol,:)
QQnl_prev(:ncol,:) = QQnl_prog(:ncol,:)
QQni_prev(:ncol,:) = QQni_prog(:ncol,:)
! -------------------------------------------------------- !
! Compute final prognostic state on which final diagnostic !
! in-stratus condensate adjustment is applied in the below.!
! Important : I must check whether there are any external !
! advective forcings of 'A_nl(i,k),A_ni(i,k)'. !
! Even they are (i.e., advection of aerosol), !
! actual droplet activation will be performd !
! in microphysics, so it will be completely !
! reasonable to 'A_nl(i,k)=A_ni(i,k)=0'. !
! -------------------------------------------------------- !
do k = 1, pver
do i = 1, ncol
T_prime0(i,k) = T_0(i,k) + dt*( A_T(i,k) + A_T_adj(i,k) + C_T(i,k) + (hlatv*QQw_prog(i,k)+(hlatv+hlatf)*QQi_prog(i,k))/cp )
qv_prime0(i,k) = qv_0(i,k) + dt*( A_qv(i,k) + A_qv_adj(i,k) + C_qv(i,k) - QQw_prog(i,k) - QQi_prog(i,k) )
ql_prime0(i,k) = ql_0(i,k) + dt*( A_ql(i,k) + A_ql_adj(i,k) + C_ql(i,k) + QQw_prog(i,k) )
qi_prime0(i,k) = qi_0(i,k) + dt*( A_qi(i,k) + A_qi_adj(i,k) + C_qi(i,k) + QQi_prog(i,k) )
nl_prime0(i,k) = max(0._r8,nl_0(i,k) + dt*( A_nl(i,k) + C_nl(i,k) + QQnl_prog(i,k) ))
ni_prime0(i,k) = max(0._r8,ni_0(i,k) + dt*( A_ni(i,k) + C_ni(i,k) + QQni_prog(i,k) ))
if( ql_prime0(i,k) .lt. qsmall ) nl_prime0(i,k) = 0._r8
if( qi_prime0(i,k) .lt. qsmall ) ni_prime0(i,k) = 0._r8
enddo
enddo
! -------------------------------------------------- !
! Perform diagnostic 'positive_moisture' constraint. !
! -------------------------------------------------- !
T_dprime(:ncol,:) = T_prime0(:ncol,:)
qv_dprime(:ncol,:) = qv_prime0(:ncol,:)
ql_dprime(:ncol,:) = ql_prime0(:ncol,:)
qi_dprime(:ncol,:) = qi_prime0(:ncol,:)
nl_dprime(:ncol,:) = nl_prime0(:ncol,:)
ni_dprime(:ncol,:) = ni_prime0(:ncol,:)
call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, &
qv_dprime, ql_dprime, qi_dprime, T_dprime, &
qvten_pwi2, qlten_pwi2, qiten_pwi2, Tten_pwi2 )
do k = 1, pver
do i = 1, ncol
if( ql_dprime(i,k) .lt. qsmall ) then
nlten_pwi2(i,k) = -nl_dprime(i,k)/dt
nl_dprime(i,k) = 0._r8
endif
if( qi_dprime(i,k) .lt. qsmall ) then
niten_pwi2(i,k) = -ni_dprime(i,k)/dt
ni_dprime(i,k) = 0._r8
endif
enddo
enddo
! -------------------------------------------------------------- !
! Add tendency associated with detrainment of cumulus condensate !
! This tendency is not used in computing Q !
! Since D_ql,D_qi,D_nl,D_ni > 0, don't need to worry about !
! negative scalar. !
! This tendency is not reflected into Fzs2, which is OK. !
! -------------------------------------------------------------- !
T_dprime(:ncol,:) = T_dprime(:ncol,:) + D_T(:ncol,:) * dt
qv_dprime(:ncol,:) = qv_dprime(:ncol,:) + D_qv(:ncol,:) * dt
ql_dprime(:ncol,:) = ql_dprime(:ncol,:) + D_ql(:ncol,:) * dt
qi_dprime(:ncol,:) = qi_dprime(:ncol,:) + D_qi(:ncol,:) * dt
nl_dprime(:ncol,:) = nl_dprime(:ncol,:) + D_nl(:ncol,:) * dt
ni_dprime(:ncol,:) = ni_dprime(:ncol,:) + D_ni(:ncol,:) * dt
! ---------------------------------------------------------- !
! Impose diagnostic upper and lower limits on the in-stratus !
! condensate amount. This produces a final equilibrium state !
! at the end of each iterative process. !
! ---------------------------------------------------------- !
do k = 1, pver
call instratus_condensate( lchnk , ncol , k , p(:,k) , &
T_dprime(:,k) , qv_dprime(:,k) , ql_dprime(:,k) , qi_dprime(:,k), &
a_cu0(:,k) , zeros(:,k) , zeros(:,k) , &
zeros(:,k) , zeros(:,k) , zeros(:,k) , &
landfrac , snowh , &
T_star(:,k) , qv_star(:,k) , ql_star(:,k) , qi_star(:,k) , &
al_st_star(:,k), ai_st_star(:,k), ql_st_star(:,k), qi_st_star(:,k) )
a_st_star(:ncol,k) = max(al_st_star(:ncol,k),ai_st_star(:ncol,k))
QQw2(:ncol,k) = (ql_star(:ncol,k) - ql_dprime(:ncol,k))/dt
QQi2(:ncol,k) = (qi_star(:ncol,k) - qi_dprime(:ncol,k))/dt
! -------------------------------------------------- !
! Reduce droplet concentration if evaporation occurs !
! -------------------------------------------------- !
do i = 1, ncol
if( QQw2(i,k) .le. 0._r8 ) then
if( ql_dprime(i,k) .ge. qsmall ) then
QQnl2(i,k) = QQw2(i,k)*nl_dprime(i,k)/ql_dprime(i,k)
QQnl2(i,k) = min(0._r8,cone*max(QQnl2(i,k),-nl_dprime(i,k)/dt))
else
QQnl2(i,k) = 0._r8
endif
endif
if( QQi2(i,k) .le. 0._r8 ) then
if( qi_dprime(i,k) .gt. qsmall ) then
QQni2(i,k) = QQi2(i,k)*ni_dprime(i,k)/qi_dprime(i,k)
QQni2(i,k) = min(0._r8,cone*max(QQni2(i,k),-ni_dprime(i,k)/dt))
else
QQni2(i,k) = 0._r8
endif
endif
enddo
enddo
nl_star(:ncol,:) = max(0._r8,nl_dprime(:ncol,:)+QQnl2(:ncol,:)*dt)
ni_star(:ncol,:) = max(0._r8,ni_dprime(:ncol,:)+QQni2(:ncol,:)*dt)
! ------------------------------------------ !
! Final adjustment of droplet concentration. !
! Set # to zero if there is no cloud. !
! ------------------------------------------ !
do k = 1, pver
do i = 1, ncol
if( ql_star(i,k) .lt. qsmall ) then
ACnl(i,k) = - nl_star(i,k)/dt
nl_star(i,k) = 0._r8
endif
if( qi_star(i,k) .lt. qsmall ) then
ACni(i,k) = - ni_star(i,k)/dt
ni_star(i,k) = 0._r8
endif
enddo
enddo
! ----------------------------------------------------- !
! Define equilibrium reference state for next iteration !
! ----------------------------------------------------- !
T(:ncol,:) = T_star(:ncol,:)
qv(:ncol,:) = qv_star(:ncol,:)
ql(:ncol,:) = ql_star(:ncol,:)
qi(:ncol,:) = qi_star(:ncol,:)
al_st(:ncol,:) = al_st_star(:ncol,:)
ai_st(:ncol,:) = ai_st_star(:ncol,:)
a_st(:ncol,:) = a_st_star(:ncol,:)
ql_st(:ncol,:) = ql_st_star(:ncol,:)
qi_st(:ncol,:) = qi_st_star(:ncol,:)
nl(:ncol,:) = nl_star(:ncol,:)
ni(:ncol,:) = ni_star(:ncol,:)
enddo ! End of 'iter' prognostic iterative computation
! ------------------------------------------------------------------------ !
! Compute final tendencies of main output variables and diagnostic outputs !
! Note that the very input state [T0,qv0,ql0,qi0] are !
! marched to [T_star,qv_star,ql_star,qi_star] with equilibrium !
! stratus informations of [a_st_star,ql_st_star,qi_st_star] by !
! below final tendencies and [A_T,A_qv,A_ql,A_qi] !
! ------------------------------------------------------------------------ !
! ------------------ !
! Process tendencies !
! ------------------ !
QQw_final(:ncol,:) = QQw_prog(:ncol,:)
QQi_final(:ncol,:) = QQi_prog(:ncol,:)
QQ_final(:ncol,:) = QQw_final(:ncol,:) + QQi_final(:ncol,:)
QQw_all(:ncol,:) = QQw_prog(:ncol,:) + QQw1(:ncol,:) + QQw2(:ncol,:) + qlten_pwi1(:ncol,:) + qlten_pwi2(:ncol,:) + A_ql_adj(:ncol,:)
QQi_all(:ncol,:) = QQi_prog(:ncol,:) + QQi1(:ncol,:) + QQi2(:ncol,:) + qiten_pwi1(:ncol,:) + qiten_pwi2(:ncol,:) + A_qi_adj(:ncol,:)
QQ_all(:ncol,:) = QQw_all(:ncol,:) + QQi_all(:ncol,:)
QQnl_final(:ncol,:) = QQnl_prog(:ncol,:)
QQni_final(:ncol,:) = QQni_prog(:ncol,:)
QQn_final(:ncol,:) = QQnl_final(:ncol,:) + QQni_final(:ncol,:)
QQnl_all(:ncol,:) = QQnl_prog(:ncol,:) + QQnl1(:ncol,:) + QQnl2(:ncol,:) + nlten_pwi1(:ncol,:) + nlten_pwi2(:ncol,:) + ACnl(:ncol,:) + A_nl_adj(:ncol,:)
QQni_all(:ncol,:) = QQni_prog(:ncol,:) + QQni1(:ncol,:) + QQni2(:ncol,:) + niten_pwi1(:ncol,:) + niten_pwi2(:ncol,:) + ACni(:ncol,:) + A_ni_adj(:ncol,:)
QQn_all(:ncol,:) = QQnl_all(:ncol,:) + QQni_all(:ncol,:)
qme(:ncol,:) = QQ_final(:ncol,:)
qvadj(:ncol,:) = qvten_pwi1(:ncol,:) + qvten_pwi2(:ncol,:) + A_qv_adj(:ncol,:)
qladj(:ncol,:) = qlten_pwi1(:ncol,:) + qlten_pwi2(:ncol,:) + A_ql_adj(:ncol,:)
qiadj(:ncol,:) = qiten_pwi1(:ncol,:) + qiten_pwi2(:ncol,:) + A_qi_adj(:ncol,:)
qllim(:ncol,:) = QQw1 (:ncol,:) + QQw2 (:ncol,:)
qilim(:ncol,:) = QQi1 (:ncol,:) + QQi2 (:ncol,:)
! ----------------- !
! Output tendencies !
! ----------------- !
s_tendout(:ncol,:) = cp*( T_star(:ncol,:) - T0(:ncol,:) )/dt - cp*(A_T(:ncol,:)+C_T(:ncol,:))
qv_tendout(:ncol,:) = ( qv_star(:ncol,:) - qv0(:ncol,:) )/dt - (A_qv(:ncol,:)+C_qv(:ncol,:))
ql_tendout(:ncol,:) = ( ql_star(:ncol,:) - ql0(:ncol,:) )/dt - (A_ql(:ncol,:)+C_ql(:ncol,:))
qi_tendout(:ncol,:) = ( qi_star(:ncol,:) - qi0(:ncol,:) )/dt - (A_qi(:ncol,:)+C_qi(:ncol,:))
nl_tendout(:ncol,:) = ( nl_star(:ncol,:) - nl0(:ncol,:) )/dt - (A_nl(:ncol,:)+C_nl(:ncol,:))
ni_tendout(:ncol,:) = ( ni_star(:ncol,:) - ni0(:ncol,:) )/dt - (A_ni(:ncol,:)+C_ni(:ncol,:))
! ------------------ !
! Net cloud fraction !
! ------------------ !
cld(:ncol,:) = a_st_star(:ncol,:) + a_cu0(:ncol,:)
! --------------------------------- !
! Updated grid-mean state variables !
! --------------------------------- !
T0(:ncol,:) = T_star(:ncol,:)
qv0(:ncol,:) = qv_star(:ncol,:)
ql0(:ncol,:) = ql_star(:ncol,:)
qi0(:ncol,:) = qi_star(:ncol,:)
nl0(:ncol,:) = nl_star(:ncol,:)
ni0(:ncol,:) = ni_star(:ncol,:)
return
end subroutine mmacro_pcond
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine instratus_condensate( lchnk, ncol, k, & 2,23
p_in, T0_in, qv0_in, ql0_in, qi0_in, &
a_dc_in, ql_dc_in, qi_dc_in, &
a_sc_in, ql_sc_in, qi_sc_in, &
landfrac, snowh, &
T_out, qv_out, ql_out, qi_out, &
al_st_out, ai_st_out, ql_st_out, qi_st_out )
! ------------------------------------------------------- !
! Diagnostically force in-stratus condensate to be !
! in the range of 'qlst_min < qc_st < qlst_max' !
! whenever stratus exists in the equilibrium state !
! ------------------------------------------------------- !
use time_manager, only: is_first_step, get_nstep
implicit none
integer, intent(in) :: lchnk ! Chunk identifier
integer, intent(in) :: ncol ! Number of atmospheric columns
integer, intent(in) :: k ! Layer index
real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa]
real(r8), intent(in) :: T0_in(pcols) ! Temperature [K]
real(r8), intent(in) :: qv0_in(pcols) ! Grid-mean water vapor [kg/kg]
real(r8), intent(in) :: ql0_in(pcols) ! Grid-mean LWC [kg/kg]
real(r8), intent(in) :: qi0_in(pcols) ! Grid-mean IWC [kg/kg]
real(r8), intent(in) :: a_dc_in(pcols) ! Deep cumulus cloud fraction
real(r8), intent(in) :: ql_dc_in(pcols) ! In-deep cumulus LWC [kg/kg]
real(r8), intent(in) :: qi_dc_in(pcols) ! In-deep cumulus IWC [kg/kg]
real(r8), intent(in) :: a_sc_in(pcols) ! Shallow cumulus cloud fraction
real(r8), intent(in) :: ql_sc_in(pcols) ! In-shallow cumulus LWC [kg/kg]
real(r8), intent(in) :: qi_sc_in(pcols) ! In-shallow cumulus IWC [kg/kg]
real(r8), intent(in) :: landfrac(pcols) ! Land fraction
real(r8), intent(in) :: snowh(pcols) ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: T_out(pcols) ! Temperature [K]
real(r8), intent(out) :: qv_out(pcols) ! Grid-mean water vapor [kg/kg]
real(r8), intent(out) :: ql_out(pcols) ! Grid-mean LWC [kg/kg]
real(r8), intent(out) :: qi_out(pcols) ! Grid-mean IWC [kg/kg]
real(r8), intent(out) :: al_st_out(pcols) ! Liquid stratus fraction
real(r8), intent(out) :: ai_st_out(pcols) ! Ice stratus fraction
real(r8), intent(out) :: ql_st_out(pcols) ! In-stratus LWC [kg/kg]
real(r8), intent(out) :: qi_st_out(pcols) ! In-stratus IWC [kg/kg]
! Local variables
integer i ! Column index
real(r8) p
real(r8) T0
real(r8) qv0
real(r8) ql0
real(r8) qi0
real(r8) a_dc
real(r8) ql_dc
real(r8) qi_dc
real(r8) a_sc
real(r8) ql_sc
real(r8) qi_sc
real(r8) esat0
real(r8) qsat0
real(r8) dqsdT0
real(r8) U0
real(r8) U0_nc
real(r8) G0_nc
real(r8) al0_st_nc
real(r8) al0_st
real(r8) ai0_st_nc
real(r8) ai0_st
real(r8) a0_st
real(r8) ql0_nc
real(r8) qi0_nc
real(r8) qc0_nc
real(r8) ql0_st
real(r8) qi0_st
real(r8) qc0_st
real(r8) T
real(r8) qv
real(r8) ql
real(r8) qi
real(r8) ql_st
real(r8) qi_st
real(r8) esat
real(r8) qsat
real(r8) dqsdT
real(r8) esat_in(pcols)
real(r8) qsat_in(pcols)
real(r8) dqsdT_in(pcols)
real(r8) U0_in(pcols)
real(r8) al0_st_nc_in(pcols)
real(r8) ai0_st_nc_in(pcols)
real(r8) G0_nc_in(pcols)
integer idxmod
real(r8) U
real(r8) U_nc
real(r8) al_st_nc
real(r8) ai_st_nc
real(r8) G_nc
real(r8) a_st
real(r8) al_st
real(r8) ai_st
real(r8) Tmin0
real(r8) Tmax0
real(r8) Tmin
real(r8) Tmax
! ---------------- !
! Main Computation !
! ---------------- !
call vqsatd2_water(T0_in(1:ncol),p_in(1:ncol),esat_in(1:ncol),qsat_in(1:ncol),dqsdT_in(1:ncol),ncol)
U0_in(:ncol) = qv0_in(:ncol)/qsat_in(:ncol)
if( CAMstfrac ) then
call astG_RHU(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol)
else
call astG_PDF(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol)
endif
call aist_vector(qv0_in(:),T0_in(:),p_in(:),qi0_in(:),landfrac(:),snowh(:),ai0_st_nc_in(:),ncol)
do i = 1, ncol
! ---------------------- !
! Define local variables !
! ---------------------- !
p = p_in(i)
T0 = T0_in(i)
qv0 = qv0_in(i)
ql0 = ql0_in(i)
qi0 = qi0_in(i)
a_dc = a_dc_in(i)
ql_dc = ql_dc_in(i)
qi_dc = qi_dc_in(i)
a_sc = a_sc_in(i)
ql_sc = ql_sc_in(i)
qi_sc = qi_sc_in(i)
ql_dc = 0._r8
qi_dc = 0._r8
ql_sc = 0._r8
qi_sc = 0._r8
esat = esat_in(i)
qsat = qsat_in(i)
dqsdT = dqsdT_in(i)
idxmod = 0
! ------------------------------------------------------------ !
! Force the grid-mean RH to be smaller than 1 if oversaturated !
! In order to be compatible with reduced 3x3 QQ, condensation !
! should occur only into the liquid in gridmean_RH. !
! ------------------------------------------------------------ !
if( qv0 .gt. qsat ) then
call gridmean_RH( lchnk, i, k, p, T0, qv0, ql0, qi0, &
a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, &
landfrac(i), snowh(i) )
call vqsatd2_water_single(T0,p,esat0,qsat0,dqsdT0)
U0 = (qv0/qsat0)
U0_nc = U0
if( CAMstfrac ) then
call astG_RHU_single(U0_nc,p,qv0,landfrac(i),snowh(i),al0_st_nc,G0_nc)
else
call astG_PDF_single(U0_nc,p,qv0,landfrac(i),snowh(i),al0_st_nc,G0_nc)
endif
call aist_single(qv0,T0,p,qi0,landfrac(i),snowh(i),ai0_st_nc)
ai0_st = (1._r8-a_dc-a_sc)*ai0_st_nc
al0_st = (1._r8-a_dc-a_sc)*al0_st_nc
a0_st = max(ai0_st,al0_st)
idxmod = 1
else
ai0_st = (1._r8-a_dc-a_sc)*ai0_st_nc_in(i)
al0_st = (1._r8-a_dc-a_sc)*al0_st_nc_in(i)
endif
a0_st = max(ai0_st,al0_st)
! ----------------------- !
! Handling of input state !
! ----------------------- !
ql0_nc = max(0._r8,ql0-a_dc*ql_dc-a_sc*ql_sc)
qi0_nc = max(0._r8,qi0-a_dc*qi_dc-a_sc*qi_sc)
qc0_nc = ql0_nc + qi0_nc
Tmin0 = T0 - (hlatv/cp)*ql0
Tmax0 = T0 + ((hlatv+hlatf)/cp)*qv0
! ------------------------------------------------------------- !
! Do nothing and just exit if generalized in-stratus condensate !
! condition is satisfied. This includes the case I. !
! For 4x4 liquid stratus, a0_st --> al0_st. !
! ------------------------------------------------------------- !
if( ( ql0_nc .ge. qlst_min*al0_st ) .and. ( ql0_nc .le. qlst_max*al0_st ) ) then
! ------------------ !
! This is the case I !
! ------------------ !
T = T0
qv = qv0
ql = ql0
qi = qi0
goto 10
else
! ----------------------------- !
! This is case II : Dense Cloud !
! ----------------------------- !
if( al0_st .eq. 0._r8 .and. ql0_nc .gt. 0._r8 ) then
! ------------------------------------- !
! Compute hypothetical full evaporation !
! ------------------------------------- !
T = Tmin0
qv = qv0 + ql0
call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
U = qv/qsat
U_nc = U
if( CAMstfrac ) then
call astG_RHU_single(U_nc,p,qv,landfrac(i),snowh(i),al_st_nc,G_nc)
else
call astG_PDF_single(U_nc,p,qv,landfrac(i),snowh(i),al_st_nc,G_nc)
endif
al_st = (1._r8-a_dc-a_sc)*al_st_nc
if( al_st .eq. 0._r8 ) then
ql = 0._r8
qi = qi0
idxmod = 1
goto 10
else
! ------------------------------------------- !
! Evaporate until qc_st decreases to qlst_max !
! ------------------------------------------- !
Tmin = Tmin0
Tmax = T0
call instratus_core( lchnk, i, k, p, &
T0, qv0, ql0, 0._r8, &
a_dc, ql_dc, qi_dc, &
a_sc, ql_sc, qi_sc, ai0_st, &
qlst_max, Tmin, Tmax, landfrac(i), snowh(i), &
T, qv, ql, qi )
idxmod = 1
goto 10
endif
! ------------------------------ !
! This is case III : Empty Cloud !
! ------------------------------ !
elseif( al0_st .gt. 0._r8 .and. ql0_nc .eq. 0._r8 ) then
! ------------------------------------------ !
! Condense until qc_st increases to qlst_min !
! ------------------------------------------ !
Tmin = Tmin0
Tmax = Tmax0
call instratus_core( lchnk, i, k, p, &
T0, qv0, ql0, 0._r8, &
a_dc, ql_dc, qi_dc, &
a_sc, ql_sc, qi_sc, ai0_st, &
qlst_min, Tmin, Tmax, landfrac(i), snowh(i), &
T, qv, ql, qi )
idxmod = 1
goto 10
! --------------- !
! This is case IV !
! --------------- !
elseif( al0_st .gt. 0._r8 .and. ql0_nc .gt. 0._r8 ) then
if( ql0_nc .gt. qlst_max*al0_st ) then
! --------------------------------------- !
! Evaporate until qc_st drops to qlst_max !
! --------------------------------------- !
Tmin = Tmin0
Tmax = Tmax0
call instratus_core( lchnk, i, k, p, &
T0, qv0, ql0, 0._r8, &
a_dc, ql_dc, qi_dc, &
a_sc, ql_sc, qi_sc, ai0_st, &
qlst_max, Tmin, Tmax, landfrac(i), snowh(i), &
T, qv, ql, qi )
idxmod = 1
goto 10
elseif( ql0_nc .lt. qlst_min*al0_st ) then
! -------------------------------------------- !
! Condensate until qc_st increases to qlst_min !
! -------------------------------------------- !
Tmin = Tmin0
Tmax = Tmax0
call instratus_core( lchnk, i, k, p, &
T0, qv0, ql0, 0._r8, &
a_dc, ql_dc, qi_dc, &
a_sc, ql_sc, qi_sc, ai0_st, &
qlst_min, Tmin, Tmax, landfrac(i), snowh(i), &
T, qv, ql, qi )
idxmod = 1
goto 10
else
! ------------------------------------------------ !
! This case should not happen. Issue error message !
! ------------------------------------------------ !
write(iulog,*) 'Impossible case1 in instratus_condensate'
call endrun
endif
! ------------------------------------------------ !
! This case should not happen. Issue error message !
! ------------------------------------------------ !
else
write(iulog,*) 'Impossible case2 in instratus_condensate'
write(iulog,*) al0_st, a_sc, a_dc
write(iulog,*) 1000*ql0_nc, 1000*(ql0+qi0)
call endrun
endif
endif
10 continue
! -------------------------------------------------- !
! Force final energy-moisture conserving consistency !
! -------------------------------------------------- !
qi = qi0
if( idxmod .eq. 1 ) then
call aist_single(qv,T,p,qi,landfrac(i),snowh(i),ai_st_nc)
ai_st = (1._r8-a_dc-a_sc)*ai_st_nc
call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
U = (qv/qsat)
U_nc = U
if( CAMstfrac ) then
call astG_RHU_single(U_nc,p,qv,landfrac(i),snowh(i),al_st_nc,G_nc)
else
call astG_PDF_single(U_nc,p,qv,landfrac(i),snowh(i),al_st_nc,G_nc)
endif
al_st = (1._r8-a_dc-a_sc)*al_st_nc
else
ai_st = (1._r8-a_dc-a_sc)*ai0_st_nc_in(i)
al_st = (1._r8-a_dc-a_sc)*al0_st_nc_in(i)
endif
a_st = max(ai_st,al_st)
if( al_st .eq. 0._r8 ) then
ql_st = 0._r8
else
ql_st = ql/al_st
endif
if( ai_st .eq. 0._r8 ) then
qi_st = 0._r8
else
qi_st = qi/ai_st
endif
qi = ai_st*qi_st
ql = al_st*ql_st
T = T0 - (hlatv/cp)*(ql0-ql) - ((hlatv+hlatf)/cp)*(qi0-qi)
qv = qv0 + ql0 - ql + qi0 - qi
! Below two lines should be removed.
! call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
! qv = max(qv,al_st*qsat+(1._r8-al_st)*1.e-12_r8)
! -------------- !
! Send to output !
! -------------- !
T_out(i) = T
qv_out(i) = qv
ql_out(i) = ql
qi_out(i) = qi
al_st_out(i) = al_st
ai_st_out(i) = ai_st
ql_st_out(i) = ql_st
qi_st_out(i) = qi_st
enddo
return
end subroutine instratus_condensate
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine instratus_core( lchnk, icol, k, p, & 4,6
T0, qv0, ql0, qi0, &
a_dc, ql_dc, qi_dc, &
a_sc, ql_sc, qi_sc, ai_st, &
qcst_crit, Tmin, Tmax, landfrac, snowh, &
T, qv, ql, qi )
! ------------------------------------------------------ !
! Subroutine to find saturation equilibrium state using !
! a Newton iteration method, so that 'qc_st = qcst_crit' !
! is satisfied. !
! ------------------------------------------------------ !
use time_manager, only: is_first_step, get_nstep
implicit none
integer, intent(in) :: lchnk ! Chunk identifier
integer, intent(in) :: icol ! Number of atmospheric columns
integer, intent(in) :: k ! Layer index
real(r8), intent(in) :: p ! Pressure [Pa]
real(r8), intent(in) :: T0 ! Temperature [K]
real(r8), intent(in) :: qv0 ! Grid-mean water vapor [kg/kg]
real(r8), intent(in) :: ql0 ! Grid-mean LWC [kg/kg]
real(r8), intent(in) :: qi0 ! Grid-mean IWC [kg/kg]
real(r8), intent(in) :: a_dc ! Deep cumulus cloud fraction
real(r8), intent(in) :: ql_dc ! In-deep cumulus LWC [kg/kg]
real(r8), intent(in) :: qi_dc ! In-deep cumulus IWC [kg/kg]
real(r8), intent(in) :: a_sc ! Shallow cumulus cloud fraction
real(r8), intent(in) :: ql_sc ! In-shallow cumulus LWC [kg/kg]
real(r8), intent(in) :: qi_sc ! In-shallow cumulus IWC [kg/kg]
real(r8), intent(in) :: ai_st ! Ice stratus fraction (fixed)
real(r8), intent(in) :: Tmin ! Minimum temperature system can have [K]
real(r8), intent(in) :: Tmax ! Maximum temperature system can have [K]
real(r8), intent(in) :: qcst_crit ! Critical in-stratus condensate [kg/kg]
real(r8), intent(in) :: landfrac ! Land fraction
real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: T ! Temperature [K]
real(r8), intent(out) :: qv ! Grid-mean water vapor [kg/kg]
real(r8), intent(out) :: ql ! Grid-mean LWC [kg/kg]
real(r8), intent(out) :: qi ! Grid-mean IWC [kg/kg]
! Local variables
integer i ! Iteration index
real(r8) muQ0, muQ
real(r8) ql_nc0, qi_nc0, qc_nc0, qc_nc
real(r8) fice0, fice
real(r8) ficeg0, ficeg
real(r8) esat0, esat
real(r8) qsat0, qsat
real(r8) dqsdT0, dqsdT
real(r8) dqcncdt, dastdt, dUdt
real(r8) alpha, beta
real(r8) U, U_nc
real(r8) al_st_nc, G_nc
real(r8) al_st
! Variables for root-finding algorithm
integer j
real(r8) x1, x2
real(r8) rtsafe
real(r8) df, dx, dxold, f, fh, fl, temp, xh, xl
real(r8), parameter :: xacc = 1.e-3
! ---------------- !
! Main computation !
! ---------------- !
ql_nc0 = max(0._r8,ql0-a_dc*ql_dc-a_sc*ql_sc)
qi_nc0 = max(0._r8,qi0-a_dc*qi_dc-a_sc*qi_sc)
qc_nc0 = max(0._r8,ql0+qi0-a_dc*(ql_dc+qi_dc)-a_sc*(ql_sc+qi_sc))
fice0 = 0._r8
ficeg0 = 0._r8
muQ0 = 1._r8
! ------------ !
! Root finding !
! ------------ !
x1 = Tmin
x2 = Tmax
call funcd_instratus( x1, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
qcst_crit, landfrac, snowh, &
fl, df, qc_nc, fice, al_st )
call funcd_instratus( x2, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
qcst_crit, landfrac, snowh, &
fh, df, qc_nc, fice, al_st )
if((fl > 0._r8 .and. fh > 0._r8) .or. (fl < 0._r8 .and. fh < 0._r8)) then
call funcd_instratus( T0, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
qcst_crit, landfrac, snowh, &
fl, df, qc_nc, fice, al_st )
rtsafe = T0
goto 10
endif
if( fl == 0._r8) then
rtsafe = x1
goto 10
elseif( fh == 0._r8) then
rtsafe = x2
goto 10
elseif( fl < 0._r8) then
xl = x1
xh = x2
else
xh = x1
xl = x2
end if
rtsafe = 0.5_r8*(x1+x2)
dxold = abs(x2-x1)
dx = dxold
call funcd_instratus( rtsafe, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
qcst_crit, landfrac, snowh, &
f, df, qc_nc, fice, al_st )
do j = 1, 20
if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) > 0._r8 .or. abs(2.0_r8*f) > abs(dxold*df) ) then
dxold = dx
dx = 0.5_r8*(xh-xl)
rtsafe = xl + dx
if(xl == rtsafe) goto 10
else
dxold = dx
dx = f/df
temp = rtsafe
rtsafe = rtsafe - dx
if (temp == rtsafe) goto 10
end if
if(abs(dx) < xacc) goto 10
call funcd_instratus( rtsafe, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
qcst_crit, landfrac, snowh, &
f, df, qc_nc, fice, al_st )
if(f < 0._r8) then
xl = rtsafe
else
xh = rtsafe
endif
enddo
10 continue
! ------------------------------------------- !
! Final safety check before sending to output !
! ------------------------------------------- !
qc_nc = max(0._r8,qc_nc)
T = rtsafe
ql = qc_nc*(1._r8-fice) + a_dc*ql_dc + a_sc*ql_sc
qi = qc_nc*fice + a_dc*qi_dc + a_sc*qi_sc
qv = qv0 + ql0 + qi0 - (qc_nc + a_dc*(ql_dc+qi_dc) + a_sc*(ql_sc+qi_sc))
qv = max(qv,1.e-12_r8)
return
end subroutine instratus_core
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine funcd_instratus( T, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, & 5,3
a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
qcst_crit, landfrac, snowh, &
f, fg, qc_nc, fice, al_st )
! --------------------------------------------------- !
! Subroutine to find function value and gradient at T !
! --------------------------------------------------- !
implicit none
real(r8), intent(in) :: T ! Iteration temperature [K]
real(r8), intent(in) :: p ! Pressure [Pa]
real(r8), intent(in) :: T0 ! Initial temperature [K]
real(r8), intent(in) :: qv0 ! Grid-mean water vapor [kg/kg]
real(r8), intent(in) :: ql0 ! Grid-mean LWC [kg/kg]
real(r8), intent(in) :: qi0 ! Grid-mean IWC [kg/kg]
real(r8), intent(in) :: fice0 !
real(r8), intent(in) :: muQ0 !
real(r8), intent(in) :: qc_nc0 !
real(r8), intent(in) :: a_dc ! Deep cumulus cloud fraction
real(r8), intent(in) :: ql_dc ! In-deep cumulus LWC [kg/kg]
real(r8), intent(in) :: qi_dc ! In-deep cumulus IWC [kg/kg]
real(r8), intent(in) :: a_sc ! Shallow cumulus cloud fraction
real(r8), intent(in) :: ql_sc ! In-shallow cumulus LWC [kg/kg]
real(r8), intent(in) :: qi_sc ! In-shallow cumulus IWC [kg/kg]
real(r8), intent(in) :: ai_st ! Ice stratus fraction (fixed)
real(r8), intent(in) :: qcst_crit ! Critical in-stratus condensate [kg/kg]
real(r8), intent(in) :: landfrac ! Land fraction
real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: f ! Value of minimization function at T
real(r8), intent(out) :: fg ! Gradient of minimization function
real(r8), intent(out) :: qc_nc !
real(r8), intent(out) :: al_st !
real(r8), intent(out) :: fice !
! Local variables
real(r8) esat
real(r8) qsat
real(r8) dqsdT
real(r8) dqcncdt
real(r8) alpha
real(r8) beta
real(r8) U
real(r8) U_nc
real(r8) al_st_nc
real(r8) G_nc
real(r8) dUdt
real(r8) dalstdt
real(r8) qv
! ---------------- !
! Main computation !
! ---------------- !
call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
fice = fice0
qc_nc = (cp/hlatv)*(T-T0)+muQ0*qc_nc0
dqcncdt = (cp/hlatv)
qv = (qv0 + ql0 + qi0 - (qc_nc + a_dc*(ql_dc+qi_dc) + a_sc*(ql_sc+qi_sc)))
alpha = (1._r8/qsat)
beta = (qv/qsat**2._r8)*dqsdT
U = (qv/qsat)
U_nc = U
if( CAMstfrac ) then
call astG_RHU_single(U_nc,p,qv,landfrac,snowh,al_st_nc,G_nc)
else
call astG_PDF_single(U_nc,p,qv,landfrac,snowh,al_st_nc,G_nc)
endif
al_st = (1._r8-a_dc-a_sc)*al_st_nc
dUdt = -(alpha*dqcncdt+beta)
dalstdt = (1._r8/G_nc)*dUdt
if( U_nc .eq. 1._r8 ) dalstdt = 0._r8
f = qc_nc - qcst_crit*al_st
fg = dqcncdt - qcst_crit*dalstdt
return
end subroutine funcd_instratus
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine gridmean_RH( lchnk, icol, k, p, T, qv, ql, qi, & 1,3
a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, &
landfrac, snowh )
! ------------------------------------------------------------- !
! Subroutine to force grid-mean RH = 1 when RH > 1 !
! This is condensation process similar to instratus_condensate. !
! During condensation, we assume 'fice' is maintained in this !
! verison for MG not for RK. !
! ------------------------------------------------------------- !
use time_manager, only: is_first_step, get_nstep
implicit none
integer, intent(in) :: lchnk ! Chunk identifier
integer, intent(in) :: icol ! Number of atmospheric columns
integer, intent(in) :: k ! Layer index
real(r8), intent(in) :: p ! Pressure [Pa]
real(r8), intent(inout) :: T ! Temperature [K]
real(r8), intent(inout) :: qv ! Grid-mean water vapor [kg/kg]
real(r8), intent(inout) :: ql ! Grid-mean LWC [kg/kg]
real(r8), intent(inout) :: qi ! Grid-mean IWC [kg/kg]
real(r8), intent(in) :: a_dc ! Deep cumulus cloud fraction
real(r8), intent(in) :: ql_dc ! In-deep cumulus LWC [kg/kg]
real(r8), intent(in) :: qi_dc ! In-deep cumulus IWC [kg/kg]
real(r8), intent(in) :: a_sc ! Shallow cumulus cloud fraction
real(r8), intent(in) :: ql_sc ! In-shallow cumulus LWC [kg/kg]
real(r8), intent(in) :: qi_sc ! In-shallow cumulus IWC [kg/kg]
real(r8), intent(in) :: landfrac ! Land fraction
real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
! Local variables
integer m ! Iteration index
real(r8) ql_nc0, qi_nc0, qc_nc0
real(r8) Tscale
real(r8) Tc, qt, qc, dqcdt, qc_nc
real(r8) esat, qsat, dqsdT
real(r8) al_st_nc, G_nc
real(r8) f, fg
real(r8), parameter :: xacc = 1.e-3
! ---------------- !
! Main computation !
! ---------------- !
ql_nc0 = max(0._r8,ql-a_dc*ql_dc-a_sc*ql_sc)
qi_nc0 = max(0._r8,qi-a_dc*qi_dc-a_sc*qi_sc)
qc_nc0 = max(0._r8,ql+qi-a_dc*(ql_dc+qi_dc)-a_sc*(ql_sc+qi_sc))
Tc = T - (hlatv/cp)*ql
qt = qv + ql
do m = 1, 20
call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
Tscale = hlatv/cp
qc = (T-Tc)/Tscale
dqcdt = 1._r8/Tscale
f = qsat + qc - qt
fg = dqsdT + dqcdt
fg = sign(1._r8,fg)*max(1.e-10_r8,abs(fg))
if( abs(f/fg) .lt. xacc ) then
goto 10
endif
T = T - f/fg
enddo
write(iulog,*) 'Convergence in gridmean_RH is not fully reached'
write(iulog,*) 'Best guess is written out'
10 continue
call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
qv = qsat
ql = qt - qv
T = Tc + (hlatv/cp)*ql
return
end subroutine gridmean_RH
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine positive_moisture( ncol, dt, qvmin, qlmin, qimin, dp, qv, ql, qi, t, qvten, qlten, qiten, tten ) 4
! ------------------------------------------------------------------------------- !
! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, !
! force them to be larger than minimum value by (1) condensating water vapor !
! into liquid or ice, and (2) by transporting water vapor from the very lower !
! layer. '2._r8' is multiplied to the minimum values for safety. !
! Update final state variables and tendencies associated with this correction. !
! If any condensation happens, update (s,t) too. !
! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
! input tendencies. !
! Be careful the order of k : '1': top layer, 'pver' : near-surface layer !
! ------------------------------------------------------------------------------- !
implicit none
integer, intent(in) :: ncol
real(r8), intent(in) :: dt
real(r8), intent(in) :: dp(pcols,pver), qvmin(pcols,pver), qlmin(pcols,pver), qimin(pcols,pver)
real(r8), intent(inout) :: qv(pcols,pver), ql(pcols,pver), qi(pcols,pver), t(pcols,pver)
real(r8), intent(out) :: qvten(pcols,pver), qlten(pcols,pver), qiten(pcols,pver), tten(pcols,pver)
integer i, k
real(r8) dql, dqi, dqv, sum, aa, dum
tten(:ncol,:pver) = 0._r8
qvten(:ncol,:pver) = 0._r8
qlten(:ncol,:pver) = 0._r8
qiten(:ncol,:pver) = 0._r8
do i = 1, ncol
do k = 1, pver
if( qv(i,k) .lt. qvmin(i,k) .or. ql(i,k) .lt. qlmin(i,k) .or. qi(i,k) .lt. qimin(i,k) ) then
goto 10
endif
enddo
goto 11
10 continue
do k = 1, pver ! From the top to the 1st (lowest) layer from the surface
dql = max(0._r8,1._r8*qlmin(i,k)-ql(i,k))
dqi = max(0._r8,1._r8*qimin(i,k)-qi(i,k))
qlten(i,k) = qlten(i,k) + dql/dt
qiten(i,k) = qiten(i,k) + dqi/dt
qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
tten(i,k) = tten(i,k) + (hlatv/cp)*(dql/dt) + ((hlatv+hlatf)/cp)*(dqi/dt)
ql(i,k) = ql(i,k) + dql
qi(i,k) = qi(i,k) + dqi
qv(i,k) = qv(i,k) - dql - dqi
t(i,k) = t(i,k) + (hlatv * dql + (hlatv+hlatf) * dqi)/cp
dqv = max(0._r8,1._r8*qvmin(i,k)-qv(i,k))
qvten(i,k) = qvten(i,k) + dqv/dt
qv(i,k) = qv(i,k) + dqv
if( k .ne. pver ) then
qv(i,k+1) = qv(i,k+1) - dqv*dp(i,k)/dp(i,k+1)
qvten(i,k+1) = qvten(i,k+1) - dqv*dp(i,k)/dp(i,k+1)/dt
endif
qv(i,k) = max(qv(i,k),qvmin(i,k))
ql(i,k) = max(ql(i,k),qlmin(i,k))
qi(i,k) = max(qi(i,k),qimin(i,k))
end do
! Extra moisture used to satisfy 'qv(i,pver)=qvmin' is proportionally
! extracted from all the layers that has 'qv > 2*qvmin'. This fully
! preserves column moisture.
if( dqv .gt. 1.e-20 ) then
sum = 0._r8
do k = 1, pver
if( qv(i,k) .gt. 2._r8*qvmin(i,k) ) sum = sum + qv(i,k)*dp(i,k)
enddo
aa = dqv*dp(i,pver)/max(1.e-20_r8,sum)
if( aa .lt. 0.5_r8 ) then
do k = 1, pver
if( qv(i,k) .gt. 2._r8*qvmin(i,k) ) then
dum = aa*qv(i,k)
qv(i,k) = qv(i,k) - dum
qvten(i,k) = qvten(i,k) - dum/dt
endif
enddo
else
write(iulog,*) 'Full positive_moisture is impossible in Park Macro'
endif
endif
11 continue
enddo
return
end subroutine positive_moisture
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine astG_PDF_single( U, p, qv, landfrac, snowh, a, Ga ) 4
! --------------------------------------------------------- !
! Compute 'stratus fraction(a)' and Gs=(dU/da) from the !
! analytical formulation of triangular PDF. !
! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', !
! so using constant 'dV' assume that width is proportional !
! to the saturation specific humidity. !
! dV ~ 0.1. !
! cldrh : RH of in-stratus( = 1 if no supersaturation) !
! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
! G is discontinuous across U = 1. In fact, it does not !
! matter whether Ga = 1.e10 or 0 at a = 1: I derived that !
! they will produce the same results. !
! --------------------------------------------------------- !
implicit none
real(r8), intent(in) :: U ! Relative humidity
real(r8), intent(in) :: p ! Pressure [Pa]
real(r8), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg]
real(r8), intent(in) :: landfrac ! Land fraction
real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: a ! Stratus fraction
real(r8), intent(out) :: Ga ! dU/da
! Local variables
integer :: i ! Loop indexes
real(r8) dV ! Width of triangular PDF
real(r8) cldrh ! RH of stratus cloud
real(r8) rhmin ! Critical RH
real(r8) rhwght
! Statement functions
logical land
land = nint(landfrac) == 1
! ---------- !
! Parameters !
! ---------- !
cldrh = 1.0_r8
! ---------------- !
! Main computation !
! ---------------- !
if( p .ge. premib ) then
if( land .and. (snowh.le.0.000001_r8) ) then
rhmin = rhminl - 0.10_r8
else
rhmin = rhminl
endif
dV = cldrh - rhmin
if( U .ge. 1._r8 ) then
a = 1._r8
Ga = 1.e10_r8
elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
(1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
elseif( U .le. (cldrh-dV) ) then
a = 0._r8
Ga = 1.e10_r8
endif
if( freeze_dry ) then
a = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
endif
elseif( p .lt. premit ) then
rhmin = rhminh
dV = cldrh - rhmin
if( U .ge. 1._r8 ) then
a = 1._r8
Ga = 1.e10_r8
elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
(1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
elseif( U .le. (cldrh-dV) ) then
a = 0._r8
Ga = 1.e10_r8
endif
else
rhwght = (premib-(max(p,premit)))/(premib-premit)
! if( land .and. (snowh.le.0.000001_r8) ) then
! rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
! else
rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
! endif
dV = cldrh - rhmin
if( U .ge. 1._r8 ) then
a = 1._r8
Ga = 1.e10_r8
elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
(1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
elseif( U .le. (cldrh-dV) ) then
a = 0._r8
Ga = 1.e10_r8
endif
endif
return
end subroutine astG_PDF_single
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine astG_PDF( U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol ) 2
! --------------------------------------------------------- !
! Compute 'stratus fraction(a)' and Gs=(dU/da) from the !
! analytical formulation of triangular PDF. !
! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', !
! so using constant 'dV' assume that width is proportional !
! to the saturation specific humidity. !
! dV ~ 0.1. !
! cldrh : RH of in-stratus( = 1 if no supersaturation) !
! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
! G is discontinuous across U = 1. In fact, it does not !
! matter whether Ga = 1.e10 or 0 at a = 1: I derived that !
! they will produce the same results. !
! --------------------------------------------------------- !
implicit none
integer, intent(in) :: ncol
real(r8), intent(in) :: U_in(pcols) ! Relative humidity
real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa]
real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor specific humidity [kg/kg]
real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction
real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: a_out(pcols) ! Stratus fraction
real(r8), intent(out) :: Ga_out(pcols) ! dU/da
real(r8) :: U ! Relative humidity
real(r8) :: p ! Pressure [Pa]
real(r8) :: qv ! Grid-mean water vapor specific humidity [kg/kg]
real(r8) :: landfrac ! Land fraction
real(r8) :: snowh ! Snow depth (liquid water equivalent)
real(r8) :: a ! Stratus fraction
real(r8) :: Ga ! dU/da
! Local variables
integer :: i ! Loop indexes
real(r8) dV ! Width of triangular PDF
real(r8) cldrh ! RH of stratus cloud
real(r8) rhmin ! Critical RH
real(r8) rhwght
! Statement functions
logical land
land(i) = nint(landfrac_in(i)) == 1
! ---------- !
! Parameters !
! ---------- !
cldrh = 1.0_r8
! ---------------- !
! Main computation !
! ---------------- !
a_out(:) = 0._r8
Ga_out(:) = 0._r8
do i = 1, ncol
U = U_in(i)
p = p_in(i)
qv = qv_in(i)
landfrac = landfrac_in(i)
snowh = snowh_in(i)
if( p .ge. premib ) then
if( land(i) .and. (snowh.le.0.000001_r8) ) then
rhmin = rhminl - 0.10_r8
else
rhmin = rhminl
endif
dV = cldrh - rhmin
if( U .ge. 1._r8 ) then
a = 1._r8
Ga = 1.e10_r8
elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
(1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
elseif( U .le. (cldrh-dV) ) then
a = 0._r8
Ga = 1.e10_r8
endif
if( freeze_dry ) then
a = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
endif
elseif( p .lt. premit ) then
rhmin = rhminh
dV = cldrh - rhmin
if( U .ge. 1._r8 ) then
a = 1._r8
Ga = 1.e10_r8
elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
(1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
elseif( U .le. (cldrh-dV) ) then
a = 0._r8
Ga = 1.e10_r8
endif
else
rhwght = (premib-(max(p,premit)))/(premib-premit)
! if( land(i) .and. (snowh.le.0.000001_r8) ) then
! rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
! else
rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
! endif
dV = cldrh - rhmin
if( U .ge. 1._r8 ) then
a = 1._r8
Ga = 1.e10_r8
elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
(1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
elseif( U .le. (cldrh-dV) ) then
a = 0._r8
Ga = 1.e10_r8
endif
endif
a_out(i) = a
Ga_out(i) = Ga
enddo
return
end subroutine astG_PDF
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine astG_RHU_single( U, p, qv, landfrac, snowh, a, Ga ) 4
! --------------------------------------------------------- !
! Compute 'stratus fraction(a)' and Gs=(dU/da) from the !
! CAM35 cloud fraction formula. !
! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core !
! For the other cases, I should re-define 'rhminl,rhminh' & !
! 'premib,premit'. !
! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
! G is discontinuous across U = 1. !
! --------------------------------------------------------- !
implicit none
real(r8), intent(in) :: U ! Relative humidity
real(r8), intent(in) :: p ! Pressure [Pa]
real(r8), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg]
real(r8), intent(in) :: landfrac ! Land fraction
real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: a ! Stratus fraction
real(r8), intent(out) :: Ga ! dU/da
! Local variables
real(r8) rhmin ! Critical RH
real(r8) rhdif ! Factor for stratus fraction
real(r8) rhwght
! Statement functions
logical land
land = nint(landfrac) == 1
! ---------------- !
! Main computation !
! ---------------- !
if( p .ge. premib ) then
if( land .and. (snowh.le.0.000001_r8) ) then
rhmin = rhminl - 0.10_r8
else
rhmin = rhminl
endif
rhdif = (U-rhmin)/(1.0_r8-rhmin)
a = min(1._r8,(max(rhdif,0.0_r8))**2)
if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
Ga = 1.e20_r8
else
Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
endif
if( freeze_dry ) then
a = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
endif
elseif( p .lt. premit ) then
rhmin = rhminh
rhdif = (U-rhmin)/(1.0_r8-rhmin)
a = min(1._r8,(max(rhdif,0._r8))**2)
if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
Ga = 1.e20_r8
else
Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
endif
else
rhwght = (premib-(max(p,premit)))/(premib-premit)
! if( land .and. (snowh.le.0.000001_r8) ) then
! rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
! else
rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
! endif
rhdif = (U-rhmin)/(1.0_r8-rhmin)
a = min(1._r8,(max(rhdif,0._r8))**2)
if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
Ga = 1.e10_r8
else
Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
endif
endif
return
end subroutine astG_RHU_single
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine astG_RHU( U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol ) 2
! --------------------------------------------------------- !
! Compute 'stratus fraction(a)' and Gs=(dU/da) from the !
! CAM35 cloud fraction formula. !
! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core !
! For the other cases, I should re-define 'rhminl,rhminh' & !
! 'premib,premit'. !
! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
! G is discontinuous across U = 1. !
! --------------------------------------------------------- !
implicit none
integer, intent(in) :: ncol
real(r8), intent(in) :: U_in(pcols) ! Relative humidity
real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa]
real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor specific humidity [kg/kg]
real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction
real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: a_out(pcols) ! Stratus fraction
real(r8), intent(out) :: Ga_out(pcols) ! dU/da
real(r8) :: U ! Relative humidity
real(r8) :: p ! Pressure [Pa]
real(r8) :: qv ! Grid-mean water vapor specific humidity [kg/kg]
real(r8) :: landfrac ! Land fraction
real(r8) :: snowh ! Snow depth (liquid water equivalent)
real(r8) :: a ! Stratus fraction
real(r8) :: Ga ! dU/da
! Local variables
integer i
real(r8) rhmin ! Critical RH
real(r8) rhdif ! Factor for stratus fraction
real(r8) rhwght
! Statement functions
logical land
land(i) = nint(landfrac_in(i)) == 1
! ---------------- !
! Main computation !
! ---------------- !
a_out(:) = 0._r8
Ga_out(:) = 0._r8
do i = 1, ncol
U = U_in(i)
p = p_in(i)
qv = qv_in(i)
landfrac = landfrac_in(i)
snowh = snowh_in(i)
if( p .ge. premib ) then
if( land(i) .and. (snowh.le.0.000001_r8) ) then
rhmin = rhminl - 0.10_r8
else
rhmin = rhminl
endif
rhdif = (U-rhmin)/(1.0_r8-rhmin)
a = min(1._r8,(max(rhdif,0.0_r8))**2)
if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
Ga = 1.e20_r8
else
Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
endif
if( freeze_dry ) then
a = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
endif
elseif( p .lt. premit ) then
rhmin = rhminh
rhdif = (U-rhmin)/(1.0_r8-rhmin)
a = min(1._r8,(max(rhdif,0._r8))**2)
if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
Ga = 1.e20_r8
else
Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
endif
else
rhwght = (premib-(max(p,premit)))/(premib-premit)
! if( land(i) .and. (snowh.le.0.000001_r8) ) then
! rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
! else
rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
! endif
rhdif = (U-rhmin)/(1.0_r8-rhmin)
a = min(1._r8,(max(rhdif,0._r8))**2)
if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
Ga = 1.e10_r8
else
Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
endif
endif
a_out(i) = a
Ga_out(i) = Ga
enddo
return
end subroutine astG_RHU
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine aist_single( qv, T, p, qi, landfrac, snowh, aist ) 2,5
! --------------------------------------------------------- !
! Compute non-physical ice stratus fraction !
! --------------------------------------------------------- !
use physconst, only: rair
use wv_saturation, only: vqsatd2_water_single
implicit none
real(r8), intent(in) :: qv ! Grid-mean water vapor[kg/kg]
real(r8), intent(in) :: T ! Temperature
real(r8), intent(in) :: p ! Pressure [Pa]
real(r8), intent(in) :: qi ! Grid-mean ice water content [kg/kg]
real(r8), intent(in) :: landfrac ! Land fraction
real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
! Local variables
real(r8) rhmin ! Critical RH
real(r8) rhwght
real(r8) a,b,c,as,bs,cs ! Fit parameters
real(r8) Kc ! Constant for ice cloud calc (wood & field)
real(r8) ttmp ! Limited temperature
real(r8) icicval ! Empirical IWC value [ kg/kg ]
real(r8) rho ! Local air density
real(r8) esl ! Liq sat vapor pressure
real(r8) esi ! Ice sat vapor pressure
real(r8) ncf,phi ! Wilson and Ballard parameters
real(r8) esat, qsat, dqsdT
real(r8) rhi ! grid box averaged relative humidity over ice
real(r8) minice ! minimum grid box avg ice for having a 'cloud'
real(r8) mincld ! minimum ice cloud fraction threshold
real(r8) icimr ! in cloud ice mixing ratio
! real(r8) qist_min ! minimum in cloud ice mixing ratio
! real(r8) qist_max ! maximum in cloud ice mixing ratio
real(r8) rhdif ! working variable for slingo scheme
! integer iceopt ! option for ice cloud closure
! 1=wang & sassen 2=schiller (iciwc)
! 3=wood & field, 4=Wilson (based on smith)
! 5=modified slingo (ssat & empyt cloud)
real(r8) icecrit ! Critical RH for ice clouds in Wilson & Ballard closure ( smaller = more ice clouds )
! Statement functions
logical land
land = nint(landfrac) == 1
! --------- !
! Constants !
! --------- !
! iceopt = 5
! Wang and Sassen IWC paramters ( Option.1 )
a = 26.87_r8
b = 0.569_r8
c = 0.002892_r8
! Schiller parameters ( Option.2 )
as = -68.4202_r8
bs = 0.983917_r8
cs = 2.81795_r8
! Wood and Field parameters ( Option.3 )
Kc = 75._r8
! Wilson & Ballard closure ( Option.4. smaller = more ice clouds)
icecrit = 0.93_r8
! Slingo modified (option 5)
minice = 1.e-12_r8
mincld = 1.e-4_r8
! qist_min = 1.e-7_r8
! qist_max = 5.e-3_r8
! ---------------- !
! Main computation !
! ---------------- !
call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
esl = polysvp(T,0)
esi = polysvp(T,1)
if( iceopt.lt.3 ) then
if( iceopt.eq.1 ) then
ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8
icicval = a + b * ttmp + c * ttmp**2._r8
rho = p/(rair*T)
icicval = icicval * 1.e-6_r8 / rho
else
ttmp = max(190._r8,min(T,273.16_r8))
icicval = 10._r8 **(as * bs**ttmp + cs)
icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
endif
aist = max(0._r8,min(qi/icicval,1._r8))
elseif( iceopt.eq.3 ) then
aist = 1._r8 - exp(-Kc*qi/(qsat*(esi/esl)))
aist = max(0._r8,min(aist,1._r8))
elseif( iceopt.eq.4) then
if( p .ge. premib ) then
if( land .and. (snowh.le.0.000001_r8) ) then
rhmin = rhminl - 0.10_r8
else
rhmin = rhminl
endif
elseif( p .lt. premit ) then
rhmin = rhminh
else
rhwght = (premib-(max(p,premit)))/(premib-premit)
! if( land .and. (snowh.le.0.000001_r8) ) then
! rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
! else
rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
! endif
endif
ncf = qi/((1._r8 - icecrit)*qsat)
if( ncf.le.0._r8 ) then
aist = 0._r8
elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then
aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then
phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
aist = (1._r8 - 4._r8 * cos(phi) * cos(phi))
else
aist = 1._r8
endif
aist = max(0._r8,min(aist,1._r8))
elseif (iceopt.eq.5) then
! set rh ice cloud fraction
rhi= (qv+qi)/qsat * (esl/esi)
rhdif= (rhi-rhmini) / (rhmaxi - rhmini)
aist = min(1.0_r8, max(rhdif,0._r8)**2)
! limiter to remove empty cloud and ice with no cloud
! and set icecld fraction to mincld if ice exists
if (qi.lt.minice) then
aist=0._r8
else
aist=max(mincld,aist)
endif
! enforce limits on icimr
if (qi.ge.minice) then
icimr=qi/aist
!minimum
if (icimr.lt.qist_min) then
aist = max(0._r8,min(1._r8,qi/qist_min))
endif
!maximum
if (icimr.gt.qist_max) then
aist = max(0._r8,min(1._r8,qi/qist_max))
endif
endif
endif
! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate
! computed after updating 'qi_st'.
aist = max(0._r8,min(aist,0.999_r8))
return
end subroutine aist_single
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine aist_vector( qv_in, T_in, p_in, qi_in, landfrac_in, snowh_in, aist_out, ncol ) 2,5
! --------------------------------------------------------- !
! Compute non-physical ice stratus fraction !
! --------------------------------------------------------- !
use physconst, only: rair
use wv_saturation, only: vqsatd2_water
implicit none
integer, intent(in) :: ncol
real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor[kg/kg]
real(r8), intent(in) :: T_in(pcols) ! Temperature
real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa]
real(r8), intent(in) :: qi_in(pcols) ! Grid-mean ice water content [kg/kg]
real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction
real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent)
real(r8), intent(out) :: aist_out(pcols) ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
! Local variables
real(r8) qv ! Grid-mean water vapor[kg/kg]
real(r8) T ! Temperature
real(r8) p ! Pressure [Pa]
real(r8) qi ! Grid-mean ice water content [kg/kg]
real(r8) landfrac ! Land fraction
real(r8) snowh ! Snow depth (liquid water equivalent)
real(r8) aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
real(r8) rhmin ! Critical RH
real(r8) rhwght
real(r8) a,b,c,as,bs,cs ! Fit parameters
real(r8) Kc ! Constant for ice cloud calc (wood & field)
real(r8) ttmp ! Limited temperature
real(r8) icicval ! Empirical IWC value [ kg/kg ]
real(r8) rho ! Local air density
real(r8) esl ! Liq sat vapor pressure
real(r8) esi ! Ice sat vapor pressure
real(r8) ncf,phi ! Wilson and Ballard parameters
real(r8) esat, qsat, dqsdT
real(r8) esat_in(pcols)
real(r8) qsat_in(pcols)
real(r8) dqsdT_in(pcols)
real(r8) rhi ! grid box averaged relative humidity over ice
real(r8) minice ! minimum grid box avg ice for having a 'cloud'
real(r8) mincld ! minimum ice cloud fraction threshold
real(r8) icimr ! in cloud ice mixing ratio
! real(r8) qist_min ! minimum in cloud ice mixing ratio
! real(r8) qist_max ! maximum in cloud ice mixing ratio
real(r8) rhdif ! working variable for slingo scheme
integer i
! integer iceopt ! option for ice cloud closure
! 1=wang & sassen 2=schiller (iciwc)
! 3=wood & field, 4=Wilson (based on smith)
! 5=modified slingo (ssat & empyt cloud)
real(r8) icecrit ! Critical RH for ice clouds in Wilson & Ballard closure ( smaller = more ice clouds )
! Statement functions
logical land
land(i) = nint(landfrac_in(i)) == 1
! --------- !
! Constants !
! --------- !
! iceopt = 5
! Wang and Sassen IWC paramters ( Option.1 )
a = 26.87_r8
b = 0.569_r8
c = 0.002892_r8
! Schiller parameters ( Option.2 )
as = -68.4202_r8
bs = 0.983917_r8
cs = 2.81795_r8
! Wood and Field parameters ( Option.3 )
Kc = 75._r8
! Wilson & Ballard closure ( Option.4. smaller = more ice clouds)
icecrit = 0.93_r8
! Slingo modified (option 5)
minice = 1.e-12_r8
mincld = 1.e-4_r8
! qist_min = 1.e-7_r8
! qist_max = 5.e-3_r8
! ---------------- !
! Main computation !
! ---------------- !
aist_out(:) = 0._r8
esat_in(:) = 0._r8
qsat_in(:) = 0._r8
dqsdT_in(:) = 0._r8
call vqsatd2_water(T_in(1:ncol),p_in(1:ncol),esat_in(1:ncol),qsat_in(1:ncol),dqsdT_in(1:ncol),ncol)
do i = 1, ncol
landfrac = landfrac_in(i)
snowh = snowh_in(i)
T = T_in(i)
qv = qv_in(i)
p = p_in(i)
qi = qi_in(i)
qsat = qsat_in(i)
esl = polysvp(T,0)
esi = polysvp(T,1)
if( iceopt.lt.3 ) then
if( iceopt.eq.1 ) then
ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8
icicval = a + b * ttmp + c * ttmp**2._r8
rho = p/(rair*T)
icicval = icicval * 1.e-6_r8 / rho
else
ttmp = max(190._r8,min(T,273.16_r8))
icicval = 10._r8 **(as * bs**ttmp + cs)
icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
endif
aist = max(0._r8,min(qi/icicval,1._r8))
elseif( iceopt.eq.3 ) then
aist = 1._r8 - exp(-Kc*qi/(qsat*(esi/esl)))
aist = max(0._r8,min(aist,1._r8))
elseif( iceopt.eq.4) then
if( p .ge. premib ) then
if( land(i) .and. (snowh.le.0.000001_r8) ) then
rhmin = rhminl - 0.10_r8
else
rhmin = rhminl
endif
elseif( p .lt. premit ) then
rhmin = rhminh
else
rhwght = (premib-(max(p,premit)))/(premib-premit)
! if( land(i) .and. (snowh.le.0.000001_r8) ) then
! rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
! else
rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
! endif
endif
ncf = qi/((1._r8 - icecrit)*qsat)
if( ncf.le.0._r8 ) then
aist = 0._r8
elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then
aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then
phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
aist = (1._r8 - 4._r8 * cos(phi) * cos(phi))
else
aist = 1._r8
endif
aist = max(0._r8,min(aist,1._r8))
elseif (iceopt.eq.5) then
! set rh ice cloud fraction
rhi= (qv+qi)/qsat * (esl/esi)
rhdif= (rhi-rhmini) / (rhmaxi - rhmini)
aist = min(1.0_r8, max(rhdif,0._r8)**2)
! limiter to remove empty cloud and ice with no cloud
! and set icecld fraction to mincld if ice exists
if (qi.lt.minice) then
aist=0._r8
else
aist=max(mincld,aist)
endif
! enforce limits on icimr
if (qi.ge.minice) then
icimr=qi/aist
!minimum
if (icimr.lt.qist_min) then
aist = max(0._r8,min(1._r8,qi/qist_min))
endif
!maximum
if (icimr.gt.qist_max) then
aist = max(0._r8,min(1._r8,qi/qist_max))
endif
endif
endif
! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate
! computed after updating 'qi_st'.
aist = max(0._r8,min(aist,0.999_r8))
aist_out(i) = aist
enddo
return
end subroutine aist_vector
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine findsp_water_onelayer (lchnk, ncol, q, t, p, tsp, qsp),6
!-----------------------------------------------------------------------
!
! Purpose:
! find the wet bulb temperature for a given t and q
! in a longitude height section
! wet bulb temp is the temperature and spec humidity that is
! just saturated and has the same enthalpy
! if q > qs(t) then tsp > t and qsp = qs(tsp) < q
! if q < qs(t) then tsp < t and qsp = qs(tsp) > q
!
! Method:
! a Newton method is used
! first guess uses an algorithm provided by John Petch from the UKMO
! we exclude points where the physical situation is unrealistic
! e.g. where the temperature is outside the range of validity for the
! saturation vapor pressure, or where the water vapor pressure
! exceeds the ambient pressure, or the saturation specific humidity is
! unrealistic
!
! Author: P. Rasch
!
!-----------------------------------------------------------------------
!
! input arguments
!
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
real(r8), intent(in) :: q(pcols) ! water vapor (kg/kg)
real(r8), intent(in) :: t(pcols) ! temperature (K)
real(r8), intent(in) :: p(pcols) ! pressure (Pa)
!
! output arguments
!
real(r8), intent(out) :: tsp(pcols) ! saturation temp (K)
real(r8), intent(out) :: qsp(pcols) ! saturation mixing ratio (kg/kg)
!
! local variables
!
integer i ! work variable
! integer k ! work variable
logical lflg ! work variable
integer iter ! work variable
integer l ! work variable
logical :: error_found
real(r8) omeps ! 1 minus epsilon
real(r8) trinv ! work variable
real(r8) es ! sat. vapor pressure
real(r8) desdt ! change in sat vap pressure wrt temperature
! real(r8) desdp ! change in sat vap pressure wrt pressure
real(r8) dqsdt ! change in sat spec. hum. wrt temperature
real(r8) dgdt ! work variable
real(r8) g ! work variable
real(r8) weight(pcols) ! work variable
real(r8) hlatsb ! (sublimation)
real(r8) hlatvp ! (vaporization)
real(r8) hltalt(pcols) ! lat. heat. of vap.
real(r8) tterm ! work var.
real(r8) qs ! spec. hum. of water vapor
real(r8) tc ! crit temp of transition to ice
real(r8) tt0 ! Freezing temperature. Needed for findsp
! work variables
real(r8) t1, q1, dt, dq
real(r8) dtm, dqm
real(r8) qvd, a1, tmp
real(r8) rair
real(r8) r1b, c1, c2, c3
real(r8) denom
real(r8) dttol
real(r8) dqtol
integer doit(pcols)
real(r8) enin(pcols), enout(pcols)
real(r8) tlim(pcols)
omeps = 1.0_r8 - epsqs
! trinv = 1.0_r8/ttrice ! put just before it is used, in case ttrice = 0.0 !++xl
a1 = 7.5_r8*log(10._r8)
rair = 287.04_r8
c3 = rair*a1/cp
dtm = 0._r8 ! needed for iter=0 blowup with f90 -ei
dqm = 0._r8 ! needed for iter=0 blowup with f90 -ei
dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
! tmin = 173.16 ! the coldest temperature we can deal with
tt0 = 273.15_r8
!
! max number of times to iterate the calculation
iter = 8
!
! Sungsu comment out k loop
! do k = 1,pver
!
! first guess on the wet bulb temperature
!
do i = 1,ncol
#ifdef DEBUG
if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
write (iulog,*) ' '
write (iulog,*) ' level, t, q, p', t(i), q(i), p(i)
endif
#endif
! limit the temperature range to that relevant to the sat vap pres tables
#if ( ! defined WACCM_MOZART )
tlim(i) = min(max(t(i),173._r8),373._r8)
#else
tlim(i) = min(max(t(i),128._r8),373._r8)
#endif
! es = estblf(tlim(i))
es = polysvp(tlim(i),0) !++xl
denom = p(i) - omeps*es
qs = epsqs*es/denom
doit(i) = 0
enout(i) = 1._r8
! make sure a meaningful calculation is possible
if (p(i) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then
!
! Saturation specific humidity
!
qs = min(epsqs*es/denom,1._r8)
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!++xl
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
! tc = tlim(i) - tt0
! lflg = (tc >= -ttrice .and. tc < 0.0_r8)
! weight(i) = min(-tc*trinv,1.0_r8)
! hlatsb = hlatv + weight(i)*hlatf
! hlatvp = hlatv - 2369.0_r8*tc
! if (tlim(i) < tt0) then
! hltalt(i) = hlatsb
! else
! hltalt(i) = hlatvp
! end if
!
! No icephs or water to ice transition
!
hlatvp = hlatv - 2369.0*(tlim(i)-tt0)
hlatsb = hlatv
if (tlim(i) < tt0) then
hltalt(i) = hlatsb
else
hltalt(i) = hlatvp
end if
!--xl
enin(i) = cp*tlim(i) + hltalt(i)*q(i)
! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
tmp = q(i) - qs
c1 = hltalt(i)*c3
c2 = (tlim(i) + 36._r8)**2
r1b = c2/(c2 + c1*qs)
qvd = r1b*tmp
tsp(i) = tlim(i) + ((hltalt(i)/cp)*qvd)
#ifdef DEBUG
if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
write (iulog,*) ' relative humidity ', q(i)/qs
write (iulog,*) ' first guess ', tsp(i)
endif
#endif
! es = estblf(tsp(i))
es = polysvp(tsp(i),0) !++xl
qsp(i) = min(epsqs*es/(p(i) - omeps*es),1._r8)
else
doit(i) = 1
tsp(i) = tlim(i)
qsp(i) = q(i)
enin(i) = 1._r8
endif
end do ! end do i
!
! now iterate on first guess
!
do l = 1, iter
dtm = 0
dqm = 0
do i = 1,ncol
if (doit(i) == 0) then
! es = estblf(tsp(i))
es = polysvp(tsp(i),0) !++xl
!
! Saturation specific humidity
!
qs = min(epsqs*es/(p(i) - omeps*es),1._r8)
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!
!++xl
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
! tc = tsp(i) - tt0
! lflg = (tc >= -ttrice .and. tc < 0.0_r8)
! weight(i) = min(-tc*trinv,1.0_r8)
! hlatsb = hlatv + weight(i)*hlatf
! hlatvp = hlatv - 2369.0_r8*tc
! if (tsp(i) < tt0) then
! hltalt(i) = hlatsb
! else
! hltalt(i) = hlatvp
! end if
! if (lflg) then
! tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5))))
! else
! tterm = 0.0_r8
! end if
! desdt = hltalt(i)*es/(rgasv*tsp(i)*tsp(i)) + tterm*trinv
!
! No icephs or water to ice transition
!
hlatvp = hlatv - 2369.0*(tsp(i)-tt0)
hlatsb = hlatv
if (tsp(i) < tt0) then
hltalt(i) = hlatsb
else
hltalt(i) = hlatvp
end if
desdt = hltalt(i)*es/(rgasv*tsp(i)*tsp(i))
!--xl
dqsdt = (epsqs + omeps*qs)/(p(i) - omeps*es)*desdt
! g = cp*(tlim(i)-tsp(i)) + hltalt(i)*q(i)- hltalt(i)*qsp(i)
g = enin(i) - (cp*tsp(i) + hltalt(i)*qsp(i))
dgdt = -(cp + hltalt(i)*dqsdt)
t1 = tsp(i) - g/dgdt
dt = abs(t1 - tsp(i))/t1
tsp(i) = max(t1,tmin)
! es = estblf(tsp(i))
es = polysvp(tsp(i),0) !++xl
q1 = min(epsqs*es/(p(i) - omeps*es),1._r8)
dq = abs(q1 - qsp(i))/max(q1,1.e-12_r8)
qsp(i) = q1
#ifdef DEBUG
if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
write (iulog,*) ' rel chg lev, iter, t, q ', l, dt, dq, g
endif
#endif
dtm = max(dtm,dt)
dqm = max(dqm,dq)
! if converged at this point, exclude it from more iterations
if (dt < dttol .and. dq < dqtol) then
doit(i) = 2
endif
enout(i) = cp*tsp(i) + hltalt(i)*qsp(i)
! bail out if we are too near the end of temp range
#if ( ! defined WACCM_MOZART )
if (tsp(i) < 174.16_r8) then
#else
if (tsp(i) < 130.16_r8) then
#endif
doit(i) = 4
endif
else
endif
end do ! do i = 1,ncol
if (dtm < dttol .and. dqm < dqtol) then
go to 10
endif
end do ! do l = 1,iter
10 continue
error_found = .false.
if (dtm > dttol .or. dqm > dqtol) then
do i = 1,ncol
if (doit(i) == 0) error_found = .true.
end do
if (error_found) then
do i = 1,ncol
if (doit(i) == 0) then
write (iulog,*) ' findsp not converging at point i ', i
write (iulog,*) ' t, q, p, enin ', t(i), q(i), p(i), enin(i)
write (iulog,*) ' tsp, qsp, enout ', tsp(i), qsp(i), enout(i)
call endrun ('FINDSP')
endif
end do
endif
endif
do i = 1,ncol
if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
error_found = .true.
endif
end do
if (error_found) then
do i = 1,ncol
if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
write (iulog,*) ' the enthalpy is not conserved for point ', &
i, enin(i), enout(i)
write (iulog,*) ' t, q, p, enin ', t(i), q(i), p(i), enin(i)
write (iulog,*) ' tsp, qsp, enout ', tsp(i), qsp(i), enout(i)
call endrun ('FINDSP')
endif
end do
endif
! end do ! level loop (k=1,pver)
return
end subroutine findsp_water_onelayer
!--xl
! ----------------- !
! End of subroutine !
! ----------------- !
subroutine findsp_water (lchnk, ncol, q, t, p, tsp, qsp) 1,6
!-----------------------------------------------------------------------
!
! Purpose:
! find the wet bulb temperature for a given t and q
! in a longitude height section
! wet bulb temp is the temperature and spec humidity that is
! just saturated and has the same enthalpy
! if q > qs(t) then tsp > t and qsp = qs(tsp) < q
! if q < qs(t) then tsp < t and qsp = qs(tsp) > q
!
! Method:
! a Newton method is used
! first guess uses an algorithm provided by John Petch from the UKMO
! we exclude points where the physical situation is unrealistic
! e.g. where the temperature is outside the range of validity for the
! saturation vapor pressure, or where the water vapor pressure
! exceeds the ambient pressure, or the saturation specific humidity is
! unrealistic
!
! Author: P. Rasch
!
!-----------------------------------------------------------------------
!
! input arguments
!
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
real(r8), intent(in) :: q(pcols,pver) ! water vapor (kg/kg)
real(r8), intent(in) :: t(pcols,pver) ! temperature (K)
real(r8), intent(in) :: p(pcols,pver) ! pressure (Pa)
!
! output arguments
!
real(r8), intent(out) :: tsp(pcols,pver) ! saturation temp (K)
real(r8), intent(out) :: qsp(pcols,pver) ! saturation mixing ratio (kg/kg)
!
! local variables
!
integer i ! work variable
integer k ! work variable
logical lflg ! work variable
integer iter ! work variable
integer l ! work variable
logical :: error_found
real(r8) omeps ! 1 minus epsilon
real(r8) trinv ! work variable
real(r8) es ! sat. vapor pressure
real(r8) desdt ! change in sat vap pressure wrt temperature
! real(r8) desdp ! change in sat vap pressure wrt pressure
real(r8) dqsdt ! change in sat spec. hum. wrt temperature
real(r8) dgdt ! work variable
real(r8) g ! work variable
real(r8) weight(pcols) ! work variable
real(r8) hlatsb ! (sublimation)
real(r8) hlatvp ! (vaporization)
real(r8) hltalt(pcols,pver) ! lat. heat. of vap.
real(r8) tterm ! work var.
real(r8) qs ! spec. hum. of water vapor
real(r8) tc ! crit temp of transition to ice
real(r8) tt0 ! Freezing temperature. Needed for findsp
! work variables
real(r8) t1, q1, dt, dq
real(r8) dtm, dqm
real(r8) qvd, a1, tmp
real(r8) rair
real(r8) r1b, c1, c2, c3
real(r8) denom
real(r8) dttol
real(r8) dqtol
integer doit(pcols)
real(r8) enin(pcols), enout(pcols)
real(r8) tlim(pcols)
omeps = 1.0_r8 - epsqs
! trinv = 1.0_r8/ttrice ! put just before it is used, in case ttrice = 0.0 !++xl
a1 = 7.5_r8*log(10._r8)
rair = 287.04_r8
c3 = rair*a1/cp
dtm = 0._r8 ! needed for iter=0 blowup with f90 -ei
dqm = 0._r8 ! needed for iter=0 blowup with f90 -ei
dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
! tmin = 173.16 ! the coldest temperature we can deal with
tt0 = 273.15_r8
!
! max number of times to iterate the calculation
iter = 8
!
do k = 1,pver
!
! first guess on the wet bulb temperature
!
do i = 1,ncol
#ifdef DEBUG
if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
write (iulog,*) ' '
write (iulog,*) ' level, t, q, p', t(i,k), q(i,k), p(i,k)
endif
#endif
! limit the temperature range to that relevant to the sat vap pres tables
#if ( ! defined WACCM_MOZART )
tlim(i) = min(max(t(i,k),173._r8),373._r8)
#else
tlim(i) = min(max(t(i,k),128._r8),373._r8)
#endif
! es = estblf(tlim(i))
es = polysvp(tlim(i),0) !++xl
denom = p(i,k) - omeps*es
qs = epsqs*es/denom
doit(i) = 0
enout(i) = 1._r8
! make sure a meaningful calculation is possible
if (p(i,k) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then
!
! Saturation specific humidity
!
qs = min(epsqs*es/denom,1._r8)
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!++xl
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
! tc = tlim(i) - tt0
! lflg = (tc >= -ttrice .and. tc < 0.0_r8)
! weight(i) = min(-tc*trinv,1.0_r8)
! hlatsb = hlatv + weight(i)*hlatf
! hlatvp = hlatv - 2369.0_r8*tc
! if (tlim(i) < tt0) then
! hltalt(i,k) = hlatsb
! else
! hltalt(i,k) = hlatvp
! end if
!
! No icephs or water to ice transition
!
hlatvp = hlatv - 2369.0*(tlim(i)-tt0)
hlatsb = hlatv
if (tlim(i) < tt0) then
hltalt(i,k) = hlatsb
else
hltalt(i,k) = hlatvp
end if
!--xl
enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k)
! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
tmp = q(i,k) - qs
c1 = hltalt(i,k)*c3
c2 = (tlim(i) + 36._r8)**2
r1b = c2/(c2 + c1*qs)
qvd = r1b*tmp
tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd)
#ifdef DEBUG
if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
write (iulog,*) ' relative humidity ', q(i,k)/qs
write (iulog,*) ' first guess ', tsp(i,k)
endif
#endif
! es = estblf(tsp(i,k))
es = polysvp(tsp(i,k),0) !++xl
qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
else
doit(i) = 1
tsp(i,k) = tlim(i)
qsp(i,k) = q(i,k)
enin(i) = 1._r8
endif
end do ! end do i
!
! now iterate on first guess
!
do l = 1, iter
dtm = 0
dqm = 0
do i = 1,ncol
if (doit(i) == 0) then
! es = estblf(tsp(i,k))
es = polysvp(tsp(i,k),0) !++xl
!
! Saturation specific humidity
!
qs = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!
!++xl
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
! tc = tsp(i,k) - tt0
! lflg = (tc >= -ttrice .and. tc < 0.0_r8)
! weight(i) = min(-tc*trinv,1.0_r8)
! hlatsb = hlatv + weight(i)*hlatf
! hlatvp = hlatv - 2369.0_r8*tc
! if (tsp(i,k) < tt0) then
! hltalt(i,k) = hlatsb
! else
! hltalt(i,k) = hlatvp
! end if
! if (lflg) then
! tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5))))
! else
! tterm = 0.0_r8
! end if
! desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k)) + tterm*trinv
!
! No icephs or water to ice transition
!
hlatvp = hlatv - 2369.0*(tsp(i,k)-tt0)
hlatsb = hlatv
if (tsp(i,k) < tt0) then
hltalt(i,k) = hlatsb
else
hltalt(i,k) = hlatvp
end if
desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k))
!--xl
dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
! g = cp*(tlim(i)-tsp(i,k)) + hltalt(i,k)*q(i,k)- hltalt(i,k)*qsp(i,k)
g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k))
dgdt = -(cp + hltalt(i,k)*dqsdt)
t1 = tsp(i,k) - g/dgdt
dt = abs(t1 - tsp(i,k))/t1
tsp(i,k) = max(t1,tmin)
! es = estblf(tsp(i,k))
es = polysvp(tsp(i,k),0) !++xl
q1 = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8)
qsp(i,k) = q1
#ifdef DEBUG
if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
write (iulog,*) ' rel chg lev, iter, t, q ', l, dt, dq, g
endif
#endif
dtm = max(dtm,dt)
dqm = max(dqm,dq)
! if converged at this point, exclude it from more iterations
if (dt < dttol .and. dq < dqtol) then
doit(i) = 2
endif
enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)
! bail out if we are too near the end of temp range
#if ( ! defined WACCM_MOZART )
if (tsp(i,k) < 174.16_r8) then
#else
if (tsp(i,k) < 130.16_r8) then
#endif
doit(i) = 4
endif
else
endif
end do ! do i = 1,ncol
if (dtm < dttol .and. dqm < dqtol) then
go to 10
endif
end do ! do l = 1,iter
10 continue
error_found = .false.
if (dtm > dttol .or. dqm > dqtol) then
do i = 1,ncol
if (doit(i) == 0) error_found = .true.
end do
if (error_found) then
do i = 1,ncol
if (doit(i) == 0) then
write (iulog,*) ' findsp not converging at point i ', i
write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
call endrun ('FINDSP')
endif
end do
endif
endif
do i = 1,ncol
if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
error_found = .true.
endif
end do
if (error_found) then
do i = 1,ncol
if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
write (iulog,*) ' the enthalpy is not conserved for point ', &
i, enin(i), enout(i)
write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
call endrun ('FINDSP')
endif
end do
endif
end do ! level loop (k=1,pver)
return
end subroutine findsp_water
!--xl
! ----------------- !
! End of subroutine !
! ----------------- !
SUBROUTINE gaussj(a,n,np,b,m,mp) 1,2
INTEGER m,mp,n,np,NMAX
real(r8) a(np,np),b(np,mp)
real(r8) aa(np,np),bb(np,mp)
PARAMETER (NMAX=50)
INTEGER i,icol,irow,j,k,l,ll,ii,jj,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
real(r8) big,dum,pivinv
aa(:,:) = a(:,:)
bb(:,:) = b(:,:)
do 11 j=1,n
ipiv(j)=0
11 continue
do 22 i=1,n
big=0.
do 13 j=1,n
if(ipiv(j).ne.1)then
do 12 k=1,n
if (ipiv(k).eq.0) then
if (abs(a(j,k)).ge.big)then
big=abs(a(j,k))
irow=j
icol=k
endif
else if (ipiv(k).gt.1) then
write(iulog,*) 'singular matrix in gaussj 1'
do ii = 1, np
do jj = 1, np
write(iulog,*) ii, jj, aa(ii,jj), bb(ii,1)
end do
end do
call endrun
endif
12 continue
endif
13 continue
ipiv(icol)=ipiv(icol)+1
if (irow.ne.icol) then
do 14 l=1,n
dum=a(irow,l)
a(irow,l)=a(icol,l)
a(icol,l)=dum
14 continue
do 15 l=1,m
dum=b(irow,l)
b(irow,l)=b(icol,l)
b(icol,l)=dum
15 continue
endif
indxr(i)=irow
indxc(i)=icol
if (a(icol,icol).eq.0.) then
write(iulog,*) 'singular matrix in gaussj 2'
do ii = 1, np
do jj = 1, np
write(iulog,*) ii, jj, aa(ii,jj), bb(ii,1)
end do
end do
call endrun
endif
pivinv=1./a(icol,icol)
a(icol,icol)=1.
do 16 l=1,n
a(icol,l)=a(icol,l)*pivinv
16 continue
do 17 l=1,m
b(icol,l)=b(icol,l)*pivinv
17 continue
do 21 ll=1,n
if(ll.ne.icol)then
dum=a(ll,icol)
a(ll,icol)=0.
do 18 l=1,n
a(ll,l)=a(ll,l)-a(icol,l)*dum
18 continue
do 19 l=1,m
b(ll,l)=b(ll,l)-b(icol,l)*dum
19 continue
endif
21 continue
22 continue
do 24 l=n,1,-1
if(indxr(l).ne.indxc(l))then
do 23 k=1,n
dum=a(k,indxr(l))
a(k,indxr(l))=a(k,indxc(l))
a(k,indxc(l))=dum
23 continue
endif
24 continue
return
end subroutine gaussj
! ----------------- !
! End of subroutine !
! ----------------- !
end module cldwat2m_macro