Line data Source code
1 :
2 : module cldwat2m_macro
3 :
4 : !--------------------------------------------------- !
5 : ! Purpose : CAM Interface for Cloud Macrophysics !
6 : ! Author : Sungsu Park !
7 : ! Description : Park et al. 2010. !
8 : ! For questions, contact Sungsu Park !
9 : ! e-mail : sungsup@ucar.edu !
10 : ! phone : 303-497-1375 !
11 : !--------------------------------------------------- !
12 :
13 : use shr_kind_mod, only: r8=>shr_kind_r8
14 : use spmd_utils, only: masterproc
15 : use ppgrid, only: pcols, pver, pverp
16 : use cam_abortutils, only: endrun
17 : use physconst, only: cpair, latvap, latice, rh2o, gravit, rair
18 : use wv_saturation, only: qsat_water, svp_water, svp_ice, qsat_ice
19 : use cam_logfile, only: iulog
20 : use ref_pres, only: top_lev=>trop_cloud_top_lev
21 : use cldfrc2m, only: astG_PDF_single, astG_PDF, astG_RHU_single, &
22 : astG_RHU, aist_single, aist_vector, &
23 : rhmini_const, rhmaxi_const
24 :
25 : implicit none
26 : private
27 : save
28 :
29 : public :: &
30 : ini_macro, &
31 : mmacro_pcond
32 :
33 : ! -------------- !
34 : ! Set Parameters !
35 : ! -------------- !
36 :
37 : ! ------------------------------------------------------------------------------- !
38 : ! Parameter used for selecting generalized critical RH for liquid and ice stratus !
39 : ! ------------------------------------------------------------------------------- !
40 :
41 : integer :: i_rhminl ! This is for liquid stratus fraction.
42 : ! If 0 : Original fixed critical RH from the namelist.
43 : ! If 1 : Add convective detrainment effect on the above '0' option.
44 : ! In this case, 'tau_detw' [s] should be specified below.
45 : ! If 2 : Use fully scale-adaptive method.
46 : ! In this case, 'tau_detw' [s] and 'c_aniso' [no unit] should
47 : ! be specified below.
48 :
49 : integer :: i_rhmini ! This is for ice stratus fraction.
50 : ! If 0 : Original fixed critical RH from the namelist.
51 : ! If 1 : Add convective detrainment effect on the above '0' option.
52 : ! In this case, 'tau_deti' [s] should be specified below.
53 : ! If 2 : Use fully scale-adaptive method.
54 : ! In this case, 'tau_deti' [s] and 'c_aniso' [no unit] should
55 : ! be specified below.
56 : ! Note that 'micro_mg_cam' is using below 'rhmini_const', regardless
57 : ! of 'i_rhmini'. This connection should be built in future.
58 :
59 : real(r8), parameter :: tau_detw =100._r8 ! Dissipation time scale of convective liquid condensate detrained
60 : ! into the clear portion. [hr]. 0.5-3 hr is possible.
61 : real(r8), parameter :: tau_deti = 1._r8 ! Dissipation time scale of convective ice condensate detrained
62 : ! into the clear portion. [hr]. 0.5-3 hr is possible.
63 : real(r8), parameter :: c_aniso = 1._r8 ! Inverse of anisotropic factor of PBL turbulence
64 :
65 : ! ----------------------------- !
66 : ! Parameters for Liquid Stratus !
67 : ! ----------------------------- !
68 :
69 : logical, parameter :: CAMstfrac = .false. ! If .true. (.false.),
70 : ! use Slingo (triangular PDF-based) liquid stratus fraction
71 : real(r8), parameter :: qlst_min = 2.e-5_r8 ! Minimum in-stratus LWC constraint [ kg/kg ]
72 : real(r8), parameter :: qlst_max = 3.e-3_r8 ! Maximum in-stratus LWC constraint [ kg/kg ]
73 : real(r8), parameter :: cc = 0.1_r8 ! For newly formed/dissipated in-stratus CWC ( 0 <= cc <= 1 )
74 : integer, parameter :: niter = 2 ! For iterative computation of QQ with 'ramda' below.
75 : real(r8), parameter :: ramda = 0.5_r8 ! Explicit : ramda = 0, Implicit : ramda = 1 ( 0<= ramda <= 1 )
76 : real(r8), private :: rhminl_const ! Critical RH for low-level liquid stratus clouds
77 : real(r8), private :: rhminl_adj_land_const ! rhminl adjustment for snowfree land
78 : real(r8), private :: rhminh_const ! Critical RH for high-level liquid stratus clouds
79 : real(r8), private :: premit ! Top height for mid-level liquid stratus fraction
80 : real(r8), private :: premib ! Bottom height for mid-level liquid stratus fraction
81 :
82 : real(r8), parameter :: qsmall = 1.e-18_r8 ! Smallest mixing ratio considered in the macrophysics
83 :
84 : contains
85 :
86 : ! -------------- !
87 : ! Initialization !
88 : ! -------------- !
89 :
90 0 : subroutine ini_macro(rhminl_opt_in, rhmini_opt_in)
91 :
92 : !--------------------------------------------------------------------- !
93 : ! !
94 : ! Purpose: Initialize constants for the liquid stratiform macrophysics !
95 : ! !
96 : ! Author: Sungsu Park, Dec.01.2009. !
97 : ! !
98 : !--------------------------------------------------------------------- !
99 :
100 : use cloud_fraction, only: cldfrc_getparams
101 : use cam_history, only: addfld
102 :
103 : integer, intent(in) :: rhminl_opt_in
104 : integer, intent(in) :: rhmini_opt_in
105 :
106 0 : i_rhminl = rhminl_opt_in
107 0 : i_rhmini = rhmini_opt_in
108 :
109 : call cldfrc_getparams(rhminl_out=rhminl_const, rhminl_adj_land_out=rhminl_adj_land_const, &
110 0 : rhminh_out=rhminh_const, premit_out=premit, premib_out=premib)
111 :
112 0 : if( masterproc ) then
113 0 : write(iulog,*) 'Park Macrophysics Parameters'
114 0 : write(iulog,*) ' rhminl = ', rhminl_const
115 0 : write(iulog,*) ' rhminl_adj_land = ', rhminl_adj_land_const
116 0 : write(iulog,*) ' rhminh = ', rhminh_const
117 0 : write(iulog,*) ' premit = ', premit
118 0 : write(iulog,*) ' premib = ', premib
119 0 : write(iulog,*) ' i_rhminl = ', i_rhminl
120 0 : write(iulog,*) ' i_rhmini = ', i_rhmini
121 : end if
122 :
123 :
124 0 : call addfld ('RHMIN_LIQ', (/ 'lev' /), 'A', 'fraction', 'Default critical RH for liquid-stratus')
125 0 : call addfld ('RHMIN_ICE', (/ 'lev' /), 'A', 'fraction', 'Default critical RH for ice-stratus')
126 0 : call addfld ('DRHMINPBL_LIQ', (/ 'lev' /), 'A', 'fraction', 'Drop of liquid-stratus critical RH by PBL turbulence')
127 0 : call addfld ('DRHMINPBL_ICE', (/ 'lev' /), 'A', 'fraction', 'Drop of ice-stratus critical RH by PBL turbulence')
128 0 : call addfld ('DRHMINDET_LIQ', (/ 'lev' /), 'A', 'fraction', 'Drop of liquid-stratus critical RH by convective detrainment')
129 0 : call addfld ('DRHMINDET_ICE', (/ 'lev' /), 'A', 'fraction', 'Drop of ice-stratus critical RH by convective detrainment')
130 :
131 0 : end subroutine ini_macro
132 :
133 : ! ------------------------------ !
134 : ! Stratiform Liquid Macrophysics !
135 : ! ------------------------------ !
136 :
137 : ! In the version, 'macro --> micro --> advective forcing --> macro...'
138 : ! A_...: only 'advective forcing' without 'microphysical tendency'
139 : ! C_...: only 'microphysical tendency'
140 : ! D_...: only 'detrainment of cumulus condensate'
141 : ! So, 'A' and 'C' are exclusive.
142 :
143 0 : subroutine mmacro_pcond( lchnk , ncol , dt , p , dp , &
144 : T0 , qv0 , ql0 , qi0 , nl0 , ni0 , &
145 : A_T , A_qv , A_ql , A_qi , A_nl , A_ni , &
146 : C_T , C_qv , C_ql , C_qi , C_nl , C_ni , C_qlst, &
147 : D_T , D_qv , D_ql , D_qi , D_nl , D_ni , &
148 : a_cud , a_cu0 , clrw_old , clri_old , landfrac , snowh , &
149 : tke , qtl_flx , qti_flx , cmfr_det , qlr_det , qir_det , &
150 : s_tendout , qv_tendout , ql_tendout , qi_tendout , nl_tendout , ni_tendout , &
151 : qme , qvadj , qladj , qiadj , qllim , qilim , &
152 : cld , al_st_star , ai_st_star , ql_st_star , qi_st_star , do_cldice )
153 :
154 0 : use constituents, only : qmin, cnst_get_ind
155 : use wv_saturation, only : findsp_vc
156 : use cam_history, only : outfld, hist_fld_active
157 :
158 : integer icol
159 : integer, intent(in) :: lchnk ! Chunk number
160 : integer, intent(in) :: ncol ! Number of active columns
161 :
162 : ! Input-Output variables
163 :
164 : real(r8), intent(inout) :: T0(pcols,pver) ! Temperature [K]
165 : real(r8), intent(inout) :: qv0(pcols,pver) ! Grid-mean water vapor specific humidity [kg/kg]
166 : real(r8), intent(inout) :: ql0(pcols,pver) ! Grid-mean liquid water content [kg/kg]
167 : real(r8), intent(inout) :: qi0(pcols,pver) ! Grid-mean ice water content [kg/kg]
168 : real(r8), intent(inout) :: nl0(pcols,pver) ! Grid-mean number concentration of cloud liquid droplet [#/kg]
169 : real(r8), intent(inout) :: ni0(pcols,pver) ! Grid-mean number concentration of cloud ice droplet [#/kg]
170 :
171 : ! Input variables
172 :
173 : real(r8), intent(in) :: dt ! Model integration time step [s]
174 : real(r8), intent(in) :: p(pcols,pver) ! Pressure at the layer mid-point [Pa]
175 : real(r8), intent(in) :: dp(pcols,pver) ! Pressure thickness [Pa] > 0
176 :
177 : real(r8), intent(in) :: A_T(pcols,pver) ! Non-microphysical advective external forcing of T [K/s]
178 : real(r8), intent(in) :: A_qv(pcols,pver) ! Non-microphysical advective external forcing of qv [kg/kg/s]
179 : real(r8), intent(in) :: A_ql(pcols,pver) ! Non-microphysical advective external forcing of ql [kg/kg/s]
180 : real(r8), intent(in) :: A_qi(pcols,pver) ! Non-microphysical advective external forcing of qi [kg/kg/s]
181 : real(r8), intent(in) :: A_nl(pcols,pver) ! Non-microphysical advective external forcing of nl [#/kg/s]
182 : real(r8), intent(in) :: A_ni(pcols,pver) ! Non-microphysical advective external forcing of ni [#/kg/s]
183 :
184 : real(r8), intent(in) :: C_T(pcols,pver) ! Microphysical advective external forcing of T [K/s]
185 : real(r8), intent(in) :: C_qv(pcols,pver) ! Microphysical advective external forcing of qv [kg/kg/s]
186 : real(r8), intent(in) :: C_ql(pcols,pver) ! Microphysical advective external forcing of ql [kg/kg/s]
187 : real(r8), intent(in) :: C_qi(pcols,pver) ! Microphysical advective external forcing of qi [kg/kg/s]
188 : real(r8), intent(in) :: C_nl(pcols,pver) ! Microphysical advective external forcing of nl [#/kg/s]
189 : real(r8), intent(in) :: C_ni(pcols,pver) ! Microphysical advective external forcing of ni [#/kg/s]
190 : real(r8), intent(in) :: C_qlst(pcols,pver) ! Microphysical advective external forcing of ql
191 : ! within liquid stratus [kg/kg/s]
192 :
193 : real(r8), intent(in) :: D_T(pcols,pver) ! Cumulus detrainment external forcing of T [K/s]
194 : real(r8), intent(in) :: D_qv(pcols,pver) ! Cumulus detrainment external forcing of qv [kg/kg/s]
195 : real(r8), intent(in) :: D_ql(pcols,pver) ! Cumulus detrainment external forcing of ql [kg/kg/s]
196 : real(r8), intent(in) :: D_qi(pcols,pver) ! Cumulus detrainment external forcing of qi [kg/kg/s]
197 : real(r8), intent(in) :: D_nl(pcols,pver) ! Cumulus detrainment external forcing of nl [#/kg/s]
198 : real(r8), intent(in) :: D_ni(pcols,pver) ! Cumulus detrainment external forcing of qi [#/kg/s]
199 :
200 : real(r8), intent(in) :: a_cud(pcols,pver) ! Old cumulus fraction before update
201 : real(r8), intent(in) :: a_cu0(pcols,pver) ! New cumulus fraction after update
202 :
203 : real(r8), intent(in) :: clrw_old(pcols,pver) ! Clear sky fraction at the previous time step for liquid stratus process
204 : real(r8), intent(in) :: clri_old(pcols,pver) ! Clear sky fraction at the previous time step for ice stratus process
205 : real(r8), pointer :: tke(:,:) ! (pcols,pverp) TKE from the PBL scheme
206 : real(r8), pointer :: qtl_flx(:,:) ! (pcols,pverp) overbar(w'qtl') from PBL scheme where qtl = qv + ql
207 : real(r8), pointer :: qti_flx(:,:) ! (pcols,pverp) overbar(w'qti') from PBL scheme where qti = qv + qi
208 : real(r8), pointer :: cmfr_det(:,:) ! (pcols,pver) Detrained mass flux from the convection scheme
209 : real(r8), pointer :: qlr_det(:,:) ! (pcols,pver) Detrained ql from the convection scheme
210 : real(r8), pointer :: qir_det(:,:) ! (pcols,pver) Detrained qi from the convection scheme
211 :
212 : real(r8), intent(in) :: landfrac(pcols) ! Land fraction
213 : real(r8), intent(in) :: snowh(pcols) ! Snow depth (liquid water equivalent)
214 : logical, intent(in) :: do_cldice ! Whether or not cldice should be prognosed
215 :
216 : ! Output variables
217 :
218 : real(r8), intent(out) :: s_tendout(pcols,pver) ! Net tendency of grid-mean s from 'Micro+Macro' processes [J/kg/s]
219 : real(r8), intent(out) :: qv_tendout(pcols,pver) ! Net tendency of grid-mean qv from 'Micro+Macro' processes [kg/kg/s]
220 : real(r8), intent(out) :: ql_tendout(pcols,pver) ! Net tendency of grid-mean ql from 'Micro+Macro' processes [kg/kg/s]
221 : real(r8), intent(out) :: qi_tendout(pcols,pver) ! Net tendency of grid-mean qi from 'Micro+Macro' processes [kg/kg/s]
222 : real(r8), intent(out) :: nl_tendout(pcols,pver) ! Net tendency of grid-mean nl from 'Micro+Macro' processes [#/kg/s]
223 : real(r8), intent(out) :: ni_tendout(pcols,pver) ! Net tendency of grid-mean ni from 'Micro+Macro' processes [#/kg/s]
224 :
225 : real(r8), intent(out) :: qme (pcols,pver) ! Net condensation rate [kg/kg/s]
226 : real(r8), intent(out) :: qvadj(pcols,pver) ! adjustment tendency from "positive_moisture" call (vapor)
227 : real(r8), intent(out) :: qladj(pcols,pver) ! adjustment tendency from "positive_moisture" call (liquid)
228 : real(r8), intent(out) :: qiadj(pcols,pver) ! adjustment tendency from "positive_moisture" call (ice)
229 : real(r8), intent(out) :: qllim(pcols,pver) ! tendency from "instratus_condensate" call (liquid)
230 : real(r8), intent(out) :: qilim(pcols,pver) ! tendency from "instratus_condensate" call (ice)
231 :
232 : real(r8), intent(out) :: cld(pcols,pver) ! Net cloud fraction ( 0 <= cld <= 1 )
233 : real(r8), intent(out) :: al_st_star(pcols,pver) ! Physical liquid stratus fraction
234 : real(r8), intent(out) :: ai_st_star(pcols,pver) ! Physical ice stratus fraction
235 : real(r8), intent(out) :: ql_st_star(pcols,pver) ! In-stratus LWC [kg/kg]
236 : real(r8), intent(out) :: qi_st_star(pcols,pver) ! In-stratus IWC [kg/kg]
237 :
238 : ! --------------- !
239 : ! Local variables !
240 : ! --------------- !
241 : integer :: ixcldliq, ixcldice
242 :
243 : integer :: i, j, k, iter, ii, jj ! Loop indexes
244 :
245 : ! Thermodynamic state variables
246 :
247 : real(r8) T(pcols,pver) ! Temperature of equilibrium reference state
248 : ! from which 'Micro & Macro' are computed [K]
249 : real(r8) T1(pcols,pver) ! Temperature after 'fice_force' on T01
250 : real(r8) T_0(pcols,pver) ! Temperature after 'instratus_condensate' on T1
251 : real(r8) T_05(pcols,pver) ! Temperature after 'advection' on T_0
252 : real(r8) T_prime0(pcols,pver) ! Temperature after 'Macrophysics (QQ)' on T_05star
253 : real(r8) T_dprime(pcols,pver) ! Temperature after 'fice_force' on T_prime
254 : real(r8) T_star(pcols,pver) ! Temperature after 'instratus_condensate' on T_dprime
255 :
256 : real(r8) qv(pcols,pver) ! Grid-mean qv of equilibrium reference state from which
257 : ! 'Micro & Macro' are computed [kg/kg]
258 : real(r8) qv1(pcols,pver) ! Grid-mean qv after 'fice_force' on qv01
259 : real(r8) qv_0(pcols,pver) ! Grid-mean qv after 'instratus_condensate' on qv1
260 : real(r8) qv_05(pcols,pver) ! Grid-mean qv after 'advection' on qv_0
261 : real(r8) qv_prime0(pcols,pver) ! Grid-mean qv after 'Macrophysics (QQ)' on qv_05star
262 : real(r8) qv_dprime(pcols,pver) ! Grid-mean qv after 'fice_force' on qv_prime
263 : real(r8) qv_star(pcols,pver) ! Grid-mean qv after 'instratus_condensate' on qv_dprime
264 :
265 : real(r8) ql(pcols,pver) ! Grid-mean ql of equilibrium reference state from which
266 : ! 'Micro & Macro' are computed [kg/kg]
267 : real(r8) ql1(pcols,pver) ! Grid-mean ql after 'fice_force' on ql01
268 : real(r8) ql_0(pcols,pver) ! Grid-mean ql after 'instratus_condensate' on ql1
269 : real(r8) ql_05(pcols,pver) ! Grid-mean ql after 'advection' on ql_0
270 : real(r8) ql_prime0(pcols,pver) ! Grid-mean ql after 'Macrophysics (QQ)' on ql_05star
271 : real(r8) ql_dprime(pcols,pver) ! Grid-mean ql after 'fice_force' on ql_prime
272 : real(r8) ql_star(pcols,pver) ! Grid-mean ql after 'instratus_condensate' on ql_dprime
273 :
274 : real(r8) qi(pcols,pver) ! Grid-mean qi of equilibrium reference state from which
275 : ! 'Micro & Macro' are computed [kg/kg]
276 : real(r8) qi1(pcols,pver) ! Grid-mean qi after 'fice_force' on qi01
277 : real(r8) qi_0(pcols,pver) ! Grid-mean qi after 'instratus_condensate' on qi1
278 : real(r8) qi_05(pcols,pver) ! Grid-mean qi after 'advection' on qi_0
279 : real(r8) qi_prime0(pcols,pver) ! Grid-mean qi after 'Macrophysics (QQ)' on qi_05star
280 : real(r8) qi_dprime(pcols,pver) ! Grid-mean qi after 'fice_force' on qi_prime
281 : real(r8) qi_star(pcols,pver) ! Grid-mean qi after 'instratus_condensate' on qi_dprime
282 :
283 : real(r8) nl(pcols,pver) ! Grid-mean nl of equilibrium reference state from which
284 : ! 'Micro & Macro' are computed [kg/kg]
285 : real(r8) nl1(pcols,pver) ! Grid-mean nl after 'fice_force' on nl01
286 : real(r8) nl_0(pcols,pver) ! Grid-mean nl after 'instratus_condensate' on nl1
287 : real(r8) nl_05(pcols,pver) ! Grid-mean nl after 'advection' on nl_0
288 : real(r8) nl_prime0(pcols,pver) ! Grid-mean nl after 'Macrophysics (QQ)' on nl_05star
289 : real(r8) nl_dprime(pcols,pver) ! Grid-mean nl after 'fice_force' on nl_prime
290 : real(r8) nl_star(pcols,pver) ! Grid-mean nl after 'instratus_condensate' on nl_dprime
291 :
292 : real(r8) ni(pcols,pver) ! Grid-mean ni of equilibrium reference state from which
293 : ! 'Micro & Macro' are computed [kg/kg]
294 : real(r8) ni1(pcols,pver) ! Grid-mean ni after 'fice_force' on ni01
295 : real(r8) ni_0(pcols,pver) ! Grid-mean ni after 'instratus_condensate' on ni1
296 : real(r8) ni_05(pcols,pver) ! Grid-mean ni after 'advection' on ni_0
297 : real(r8) ni_prime0(pcols,pver) ! Grid-mean ni after 'Macrophysics (QQ)' on ni_05star
298 : real(r8) ni_dprime(pcols,pver) ! Grid-mean ni after 'fice_force' on ni_prime
299 : real(r8) ni_star(pcols,pver) ! Grid-mean ni after 'instratus_condensate' on ni_dprime
300 :
301 : real(r8) a_st(pcols,pver) ! Stratus fraction of equilibrium reference state
302 : real(r8) a_st_0(pcols,pver) ! Stratus fraction at '_0' state
303 : real(r8) a_st_star(pcols,pver) ! Stratus fraction at '_star' state
304 :
305 : real(r8) al_st(pcols,pver) ! Liquid stratus fraction of equilibrium reference state
306 : real(r8) al_st_0(pcols,pver) ! Liquid stratus fraction at '_0' state
307 : real(r8) al_st_nc(pcols,pver) ! Non-physical liquid stratus fraction in the non-cumulus pixels
308 :
309 : real(r8) ai_st(pcols,pver) ! Ice stratus fraction of equilibrium reference state
310 : real(r8) ai_st_0(pcols,pver) ! Ice stratus fraction at '_0' state
311 : real(r8) ai_st_nc(pcols,pver) ! Non-physical ice stratus fraction in the non-cumulus pixels
312 :
313 : real(r8) ql_st(pcols,pver) ! In-stratus LWC of equilibrium reference state [kg/kg]
314 : real(r8) ql_st_0(pcols,pver) ! In-stratus LWC at '_0' state
315 :
316 : real(r8) qi_st(pcols,pver) ! In-stratus IWC of equilibrium reference state [kg/kg]
317 : real(r8) qi_st_0(pcols,pver) ! In-stratus IWC at '_0' state
318 :
319 : ! Cumulus properties
320 :
321 : real(r8) dacudt(pcols,pver)
322 : real(r8) a_cu(pcols,pver)
323 :
324 : ! Adjustment tendency in association with 'positive_moisture'
325 :
326 : real(r8) Tten_pwi1(pcols,pver) ! Pre-process T tendency of input equilibrium state [K/s]
327 : real(r8) qvten_pwi1(pcols,pver) ! Pre-process qv tendency of input equilibrium state [kg/kg/s]
328 : real(r8) qlten_pwi1(pcols,pver) ! Pre-process ql tendency of input equilibrium state [kg/kg/s]
329 : real(r8) qiten_pwi1(pcols,pver) ! Pre-process qi tendency of input equilibrium state [kg/kg/s]
330 : real(r8) nlten_pwi1(pcols,pver) ! Pre-process nl tendency of input equilibrium state [#/kg/s]
331 : real(r8) niten_pwi1(pcols,pver) ! Pre-process ni tendency of input equilibrium state [#/kg/s]
332 :
333 : real(r8) Tten_pwi2(pcols,pver) ! Post-process T tendency of provisional equilibrium state [K/s]
334 : real(r8) qvten_pwi2(pcols,pver) ! Post-process qv tendency of provisional equilibrium state [kg/kg/s]
335 : real(r8) qlten_pwi2(pcols,pver) ! Post-process ql tendency of provisional equilibrium state [kg/kg/s]
336 : real(r8) qiten_pwi2(pcols,pver) ! Post-process qi tendency of provisional equilibrium state [kg/kg/s]
337 : real(r8) nlten_pwi2(pcols,pver) ! Post-process nl tendency of provisoonal equilibrium state [#/kg/s]
338 : real(r8) niten_pwi2(pcols,pver) ! Post-process ni tendency of provisional equilibrium state [#/kg/s]
339 :
340 : real(r8) A_T_adj(pcols,pver) ! After applying external advective forcing [K/s]
341 : real(r8) A_qv_adj(pcols,pver) ! After applying external advective forcing [kg/kg/s]
342 : real(r8) A_ql_adj(pcols,pver) ! After applying external advective forcing [kg/kg/s]
343 : real(r8) A_qi_adj(pcols,pver) ! After applying external advective forcing [kg/kg/s]
344 : real(r8) A_nl_adj(pcols,pver) ! After applying external advective forcing [#/kg/s]
345 : real(r8) A_ni_adj(pcols,pver) ! After applying external advective forcing [#/kg/s]
346 :
347 : ! Adjustment tendency in association with 'instratus_condensate'
348 :
349 : real(r8) QQw1(pcols,pver) ! Effective adjustive condensation into water due to 'instratus_condensate' [kg/kg/s]
350 : real(r8) QQi1(pcols,pver) ! Effective adjustive condensation into ice due to 'instratus_condensate' [kg/kg/s]
351 : real(r8) QQw2(pcols,pver) ! Effective adjustive condensation into water due to 'instratus_condensate' [kg/kg/s]
352 : real(r8) QQi2(pcols,pver) ! Effective adjustive condensation into ice due to 'instratus_condensate' [kg/kg/s]
353 :
354 : real(r8) QQnl1(pcols,pver) ! Tendency of nl associated with QQw1 only when QQw1<0 (net evaporation) [#/kg/s]
355 : real(r8) QQni1(pcols,pver) ! Tendency of ni associated with QQi1 only when QQw1<0 (net evaporation) [#/kg/s]
356 : real(r8) QQnl2(pcols,pver) ! Tendency of nl associated with QQw2 only when QQw2<0 (net evaporation) [#/kg/s]
357 : real(r8) QQni2(pcols,pver) ! Tendency of ni associated with QQi2 only when QQw2<0 (net evaporation) [#/kg/s]
358 :
359 : ! Macrophysical process tendency variables
360 :
361 : real(r8) QQ(pcols,pver) ! Net condensation rate into water+ice [kg/kg/s]
362 : real(r8) QQw(pcols,pver) ! Net condensation rate into water [kg/kg/s]
363 : real(r8) QQi(pcols,pver) ! Net condensation rate into ice [kg/kg/s]
364 : real(r8) QQnl(pcols,pver) ! Tendency of nl associated with QQw both for condensation and evaporation [#/kg/s]
365 : real(r8) QQni(pcols,pver) ! Tendency of ni associated with QQi both for condensation and evaporation [#/kg/s]
366 : real(r8) ACnl(pcols,pver) ! Cloud liquid droplet (nl) activation tendency [#/kg/s]
367 : real(r8) ACni(pcols,pver) ! Cloud ice droplet (ni) activation tendency [#/kg/s]
368 :
369 : real(r8) QQw_prev(pcols,pver)
370 : real(r8) QQi_prev(pcols,pver)
371 : real(r8) QQnl_prev(pcols,pver)
372 : real(r8) QQni_prev(pcols,pver)
373 :
374 : real(r8) QQw_prog(pcols,pver)
375 : real(r8) QQi_prog(pcols,pver)
376 : real(r8) QQnl_prog(pcols,pver)
377 : real(r8) QQni_prog(pcols,pver)
378 :
379 : real(r8) QQ_final(pcols,pver)
380 : real(r8) QQw_final(pcols,pver)
381 : real(r8) QQi_final(pcols,pver)
382 : real(r8) QQn_final(pcols,pver)
383 : real(r8) QQnl_final(pcols,pver)
384 : real(r8) QQni_final(pcols,pver)
385 :
386 : real(r8) QQ_all(pcols,pver) ! QQw_all + QQi_all
387 : real(r8) QQw_all(pcols,pver) ! QQw_final + QQw1 + QQw2 + qlten_pwi1 + qlten_pwi2 + A_ql_adj [kg/kg/s]
388 : real(r8) QQi_all(pcols,pver) ! QQi_final + QQi1 + QQi2 + qiten_pwi1 + qiten_pwi2 + A_qi_adj [kg/kg/s]
389 : real(r8) QQn_all(pcols,pver) ! QQnl_all + QQni_all
390 : real(r8) QQnl_all(pcols,pver) ! QQnl_final + QQnl1 + QQnl2 + nlten_pwi1 + nlten_pwi2 + ACnl [#/kg/s]
391 : real(r8) QQni_all(pcols,pver) ! QQni_final + QQni1 + QQni2 + niten_pwi1 + niten_pwi2 + ACni [#/kg/s]
392 :
393 : ! Coefficient for computing QQ and related processes
394 :
395 : real(r8) U(pcols,pver) ! Grid-mean RH
396 : real(r8) U_nc(pcols,pver) ! Mean RH of non-cumulus pixels
397 : real(r8) G_nc(pcols,pver) ! d(U_nc)/d(a_st_nc)
398 : real(r8) F_nc(pcols,pver) ! A function of second parameter for a_st_nc
399 : real(r8) alpha ! = 1/qs
400 : real(r8) beta ! = (qv/qs**2)*dqsdT
401 : real(r8) betast ! = alpha*dqsdT
402 : real(r8) gammal ! = alpha + (latvap/cpair)*beta
403 : real(r8) gammai ! = alpha + ((latvap+latice)/cpair)*beta
404 : real(r8) gammaQ ! = alpha + (latvap/cpair)*beta
405 : real(r8) deltal ! = 1 + a_st*(latvap/cpair)*(betast/alpha)
406 : real(r8) deltai ! = 1 + a_st*((latvap+latice)/cpair)*(betast/alpha)
407 : real(r8) A_Tc ! Advective external forcing of Tc [K/s]
408 : real(r8) A_qt ! Advective external forcing of qt [kg/kg/s]
409 : real(r8) C_Tc ! Microphysical forcing of Tc [K/s]
410 : real(r8) C_qt ! Microphysical forcing of qt [kg/kg/s]
411 : real(r8) dTcdt ! d(Tc)/dt [K/s]
412 : real(r8) dqtdt ! d(qt)/dt [kg/kg/s]
413 : real(r8) dqtstldt ! d(qt_alst)/dt [kg/kg/s]
414 : real(r8) dqidt ! d(qi)/dt [kg/kg/s]
415 :
416 : real(r8) dqlstdt ! d(ql_st)/dt [kg/kg/s]
417 : real(r8) dalstdt ! d(al_st)/dt [1/s]
418 : real(r8) dastdt ! d(a_st)/dt [1/s]
419 :
420 : real(r8) anic ! Fractional area of non-cumulus and non-ice stratus fraction
421 : real(r8) GG ! G_nc(i,k)/anic
422 :
423 : real(r8) aa(2,2)
424 : real(r8) bb(2,1)
425 :
426 : real(r8) zeros(pcols,pver)
427 :
428 : real(r8) qmin1(pcols,pver)
429 : real(r8) qmin2(pcols,pver)
430 : real(r8) qmin3(pcols,pver)
431 :
432 : real(r8) esat_a(pcols) ! Saturation water vapor pressure [Pa]
433 : real(r8) qsat_a(pcols,pver) ! Saturation water vapor specific humidity [kg/kg]
434 : real(r8) Twb_aw(pcols) ! Wet-bulb temperature [K]
435 : real(r8) qvwb_aw(pcols,pver) ! Wet-bulb water vapor specific humidity [kg/kg]
436 :
437 : real(r8) esat_b(pcols)
438 : real(r8) qsat_b(pcols)
439 : real(r8) dqsdT_b(pcols)
440 :
441 : logical land
442 : real(r8) tmp
443 :
444 : real(r8) d_rhmin_liq_PBL(pcols,pver)
445 : real(r8) d_rhmin_ice_PBL(pcols,pver)
446 : real(r8) d_rhmin_liq_det(pcols,pver)
447 : real(r8) d_rhmin_ice_det(pcols,pver)
448 : real(r8) rhmaxi_arr(pcols,pver)
449 : real(r8) rhmini_arr(pcols,pver)
450 : real(r8) rhminl_arr(pcols,pver)
451 : real(r8) rhminl_adj_land_arr(pcols,pver)
452 : real(r8) rhminh_arr(pcols,pver)
453 : real(r8) rhmin_liq_diag(pcols,pver)
454 : real(r8) rhmin_ice_diag(pcols,pver)
455 :
456 : real(r8) QQmax,QQmin,QQwmin,QQimin ! For limiting QQ
457 : real(r8) cone ! Number close to but smaller than 1
458 :
459 0 : cone = 0.999_r8
460 0 : zeros(:ncol,:) = 0._r8
461 :
462 : ! ------------------------------------ !
463 : ! Global initialization of main output !
464 : ! ------------------------------------ !
465 :
466 0 : s_tendout(:ncol,:) = 0._r8
467 0 : qv_tendout(:ncol,:) = 0._r8
468 0 : ql_tendout(:ncol,:) = 0._r8
469 0 : qi_tendout(:ncol,:) = 0._r8
470 0 : nl_tendout(:ncol,:) = 0._r8
471 0 : ni_tendout(:ncol,:) = 0._r8
472 :
473 0 : qme(:ncol,:) = 0._r8
474 :
475 0 : cld(:ncol,:) = 0._r8
476 0 : al_st_star(:ncol,:) = 0._r8
477 0 : ai_st_star(:ncol,:) = 0._r8
478 0 : ql_st_star(:ncol,:) = 0._r8
479 0 : qi_st_star(:ncol,:) = 0._r8
480 :
481 : ! --------------------------------------- !
482 : ! Initialization of internal 2D variables !
483 : ! --------------------------------------- !
484 :
485 0 : T(:ncol,:) = 0._r8
486 0 : T1(:ncol,:) = 0._r8
487 0 : T_0(:ncol,:) = 0._r8
488 0 : T_05(:ncol,:) = 0._r8
489 0 : T_prime0(:ncol,:) = 0._r8
490 0 : T_dprime(:ncol,:) = 0._r8
491 0 : T_star(:ncol,:) = 0._r8
492 :
493 0 : qv(:ncol,:) = 0._r8
494 0 : qv1(:ncol,:) = 0._r8
495 0 : qv_0(:ncol,:) = 0._r8
496 0 : qv_05(:ncol,:) = 0._r8
497 0 : qv_prime0(:ncol,:) = 0._r8
498 0 : qv_dprime(:ncol,:) = 0._r8
499 0 : qv_star(:ncol,:) = 0._r8
500 :
501 0 : ql(:ncol,:) = 0._r8
502 0 : ql1(:ncol,:) = 0._r8
503 0 : ql_0(:ncol,:) = 0._r8
504 0 : ql_05(:ncol,:) = 0._r8
505 0 : ql_prime0(:ncol,:) = 0._r8
506 0 : ql_dprime(:ncol,:) = 0._r8
507 0 : ql_star(:ncol,:) = 0._r8
508 :
509 0 : qi(:ncol,:) = 0._r8
510 0 : qi1(:ncol,:) = 0._r8
511 0 : qi_0(:ncol,:) = 0._r8
512 0 : qi_05(:ncol,:) = 0._r8
513 0 : qi_prime0(:ncol,:) = 0._r8
514 0 : qi_dprime(:ncol,:) = 0._r8
515 0 : qi_star(:ncol,:) = 0._r8
516 :
517 0 : nl(:ncol,:) = 0._r8
518 0 : nl1(:ncol,:) = 0._r8
519 0 : nl_0(:ncol,:) = 0._r8
520 0 : nl_05(:ncol,:) = 0._r8
521 0 : nl_prime0(:ncol,:) = 0._r8
522 0 : nl_dprime(:ncol,:) = 0._r8
523 0 : nl_star(:ncol,:) = 0._r8
524 :
525 0 : ni(:ncol,:) = 0._r8
526 0 : ni1(:ncol,:) = 0._r8
527 0 : ni_0(:ncol,:) = 0._r8
528 0 : ni_05(:ncol,:) = 0._r8
529 0 : ni_prime0(:ncol,:) = 0._r8
530 0 : ni_dprime(:ncol,:) = 0._r8
531 0 : ni_star(:ncol,:) = 0._r8
532 :
533 0 : a_st(:ncol,:) = 0._r8
534 0 : a_st_0(:ncol,:) = 0._r8
535 0 : a_st_star(:ncol,:) = 0._r8
536 :
537 0 : al_st(:ncol,:) = 0._r8
538 0 : al_st_0(:ncol,:) = 0._r8
539 0 : al_st_nc(:ncol,:) = 0._r8
540 :
541 0 : ai_st(:ncol,:) = 0._r8
542 0 : ai_st_0(:ncol,:) = 0._r8
543 0 : ai_st_nc(:ncol,:) = 0._r8
544 :
545 0 : ql_st(:ncol,:) = 0._r8
546 0 : ql_st_0(:ncol,:) = 0._r8
547 :
548 0 : qi_st(:ncol,:) = 0._r8
549 0 : qi_st_0(:ncol,:) = 0._r8
550 :
551 : ! Cumulus properties
552 :
553 0 : dacudt(:ncol,:) = 0._r8
554 0 : a_cu(:ncol,:) = 0._r8
555 :
556 : ! Adjustment tendency in association with 'positive_moisture'
557 :
558 0 : Tten_pwi1(:ncol,:) = 0._r8
559 0 : qvten_pwi1(:ncol,:) = 0._r8
560 0 : qlten_pwi1(:ncol,:) = 0._r8
561 0 : qiten_pwi1(:ncol,:) = 0._r8
562 0 : nlten_pwi1(:ncol,:) = 0._r8
563 0 : niten_pwi1(:ncol,:) = 0._r8
564 :
565 0 : Tten_pwi2(:ncol,:) = 0._r8
566 0 : qvten_pwi2(:ncol,:) = 0._r8
567 0 : qlten_pwi2(:ncol,:) = 0._r8
568 0 : qiten_pwi2(:ncol,:) = 0._r8
569 0 : nlten_pwi2(:ncol,:) = 0._r8
570 0 : niten_pwi2(:ncol,:) = 0._r8
571 :
572 0 : A_T_adj(:ncol,:) = 0._r8
573 0 : A_qv_adj(:ncol,:) = 0._r8
574 0 : A_ql_adj(:ncol,:) = 0._r8
575 0 : A_qi_adj(:ncol,:) = 0._r8
576 0 : A_nl_adj(:ncol,:) = 0._r8
577 0 : A_ni_adj(:ncol,:) = 0._r8
578 :
579 0 : qvadj (:ncol,:) = 0._r8
580 0 : qladj (:ncol,:) = 0._r8
581 0 : qiadj (:ncol,:) = 0._r8
582 :
583 : ! Adjustment tendency in association with 'instratus_condensate'
584 :
585 0 : QQw1(:ncol,:) = 0._r8
586 0 : QQi1(:ncol,:) = 0._r8
587 0 : QQw2(:ncol,:) = 0._r8
588 0 : QQi2(:ncol,:) = 0._r8
589 :
590 0 : QQnl1(:ncol,:) = 0._r8
591 0 : QQni1(:ncol,:) = 0._r8
592 0 : QQnl2(:ncol,:) = 0._r8
593 0 : QQni2(:ncol,:) = 0._r8
594 :
595 0 : QQnl(:ncol,:) = 0._r8
596 0 : QQni(:ncol,:) = 0._r8
597 :
598 : ! Macrophysical process tendency variables
599 :
600 0 : QQ(:ncol,:) = 0._r8
601 0 : QQw(:ncol,:) = 0._r8
602 0 : QQi(:ncol,:) = 0._r8
603 0 : QQnl(:ncol,:) = 0._r8
604 0 : QQni(:ncol,:) = 0._r8
605 0 : ACnl(:ncol,:) = 0._r8
606 0 : ACni(:ncol,:) = 0._r8
607 :
608 0 : QQw_prev(:ncol,:) = 0._r8
609 0 : QQi_prev(:ncol,:) = 0._r8
610 0 : QQnl_prev(:ncol,:) = 0._r8
611 0 : QQni_prev(:ncol,:) = 0._r8
612 :
613 0 : QQw_prog(:ncol,:) = 0._r8
614 0 : QQi_prog(:ncol,:) = 0._r8
615 0 : QQnl_prog(:ncol,:) = 0._r8
616 0 : QQni_prog(:ncol,:) = 0._r8
617 :
618 0 : QQ_final(:ncol,:) = 0._r8
619 0 : QQw_final(:ncol,:) = 0._r8
620 0 : QQi_final(:ncol,:) = 0._r8
621 0 : QQn_final(:ncol,:) = 0._r8
622 0 : QQnl_final(:ncol,:) = 0._r8
623 0 : QQni_final(:ncol,:) = 0._r8
624 :
625 0 : QQ_all(:ncol,:) = 0._r8
626 0 : QQw_all(:ncol,:) = 0._r8
627 0 : QQi_all(:ncol,:) = 0._r8
628 0 : QQn_all(:ncol,:) = 0._r8
629 0 : QQnl_all(:ncol,:) = 0._r8
630 0 : QQni_all(:ncol,:) = 0._r8
631 :
632 : ! Coefficient for computing QQ and related processes
633 :
634 0 : U(:ncol,:) = 0._r8
635 0 : U_nc(:ncol,:) = 0._r8
636 0 : G_nc(:ncol,:) = 0._r8
637 0 : F_nc(:ncol,:) = 0._r8
638 :
639 : ! Other
640 :
641 0 : qmin1(:ncol,:) = 0._r8
642 0 : qmin2(:ncol,:) = 0._r8
643 0 : qmin3(:ncol,:) = 0._r8
644 :
645 : ! ---------------- !
646 : ! Main computation !
647 : ! ---------------- !
648 :
649 : ! Compute critical RH for stratus
650 0 : rhmaxi_arr(:ncol,:pver) = rhmaxi_const
651 : call rhcrit_calc( &
652 : ncol, dp, T0, p, &
653 : clrw_old, clri_old, tke, qtl_flx, &
654 : qti_flx, cmfr_det, qlr_det, qir_det, &
655 : rhmaxi_arr, rhmini_arr, rhminl_arr, rhminl_adj_land_arr, rhminh_arr, &
656 0 : d_rhmin_liq_PBL, d_rhmin_ice_PBL, d_rhmin_liq_det, d_rhmin_ice_det)
657 :
658 : ! ---------------------------------- !
659 : ! Compute cumulus-related properties !
660 : ! ---------------------------------- !
661 :
662 : dacudt(:ncol,top_lev:pver) = &
663 0 : (a_cu0(:ncol,top_lev:pver) - a_cud(:ncol,top_lev:pver))/dt
664 :
665 : ! ---------------------------------------------------------------------- !
666 : ! set to zero for levels above
667 : ! ---------------------------------------------------------------------- !
668 0 : ql0(:ncol,:top_lev-1) = 0._r8
669 0 : qi0(:ncol,:top_lev-1) = 0._r8
670 0 : nl0(:ncol,:top_lev-1) = 0._r8
671 0 : ni0(:ncol,:top_lev-1) = 0._r8
672 :
673 : ! ---------------------------------------------------------------------- !
674 : ! Check if input non-cumulus pixels satisfie a non-negative constraint. !
675 : ! If not, force all water vapor substances to be positive in all layers. !
676 : ! We should use 'old' cumulus properties for this routine. !
677 : ! ---------------------------------------------------------------------- !
678 :
679 0 : T1(:ncol,:) = T0(:ncol,:)
680 0 : qv1(:ncol,:) = qv0(:ncol,:)
681 0 : ql1(:ncol,:) = ql0(:ncol,:)
682 0 : qi1(:ncol,:) = qi0(:ncol,:)
683 0 : nl1(:ncol,:) = nl0(:ncol,:)
684 0 : ni1(:ncol,:) = ni0(:ncol,:)
685 :
686 :
687 0 : call cnst_get_ind( 'CLDLIQ', ixcldliq )
688 0 : call cnst_get_ind( 'CLDICE', ixcldice )
689 :
690 :
691 0 : qmin1(:ncol,:) = qmin(1)
692 0 : qmin2(:ncol,:) = qmin(ixcldliq)
693 0 : qmin3(:ncol,:) = qmin(ixcldice)
694 :
695 : call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, &
696 : qv1, ql1, qi1, T1, qvten_pwi1, qlten_pwi1, &
697 0 : qiten_pwi1, Tten_pwi1, do_cldice)
698 :
699 0 : do k = top_lev, pver
700 0 : do i = 1, ncol
701 0 : if( ql1(i,k) .lt. qsmall ) then
702 0 : nlten_pwi1(i,k) = -nl1(i,k)/dt
703 0 : nl1(i,k) = 0._r8
704 : endif
705 0 : if( qi1(i,k) .lt. qsmall ) then
706 0 : niten_pwi1(i,k) = -ni1(i,k)/dt
707 0 : ni1(i,k) = 0._r8
708 : endif
709 : enddo
710 : enddo
711 :
712 : ! ------------------------------------------------------------- !
713 : ! Impose 'in-stratus condensate amount constraint' !
714 : ! such that it is bounded by two limiting values. !
715 : ! This should also use 'old' cumulus properties since it is !
716 : ! before applying external forcings. !
717 : ! Below 'QQw1,QQi1' are effective adjustive condensation !
718 : ! Although this process also involves freezing of cloud !
719 : ! liquid into ice, they can be and only can be expressed !
720 : ! in terms of effective condensation. !
721 : ! ------------------------------------------------------------- !
722 :
723 0 : do k = top_lev, pver
724 : call instratus_condensate( lchnk, ncol, k, &
725 : p(:,k), T1(:,k), qv1(:,k), ql1(:,k), qi1(:,k), &
726 : ni1(:,k), &
727 : a_cud(:,k), zeros(:,k), zeros(:,k), &
728 : zeros(:,k), zeros(:,k), zeros(:,k), &
729 : landfrac, snowh, &
730 : rhmaxi_arr(:,k),rhmini_arr(:,k), rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k), &
731 : T_0(:,k), qv_0(:,k), ql_0(:,k), qi_0(:,k), &
732 0 : al_st_0(:,k), ai_st_0(:,k), ql_st_0(:,k), qi_st_0(:,k) )
733 0 : a_st_0(:ncol,k) = max(al_st_0(:ncol,k),ai_st_0(:ncol,k))
734 0 : QQw1(:ncol,k) = (ql_0(:ncol,k) - ql1(:ncol,k))/dt
735 0 : QQi1(:ncol,k) = (qi_0(:ncol,k) - qi1(:ncol,k))/dt
736 : ! -------------------------------------------------- !
737 : ! Reduce droplet concentration if evaporation occurs !
738 : ! Set a limit such that negative state not happens. !
739 : ! -------------------------------------------------- !
740 0 : do i = 1, ncol
741 0 : if( QQw1(i,k) .le. 0._r8 ) then
742 0 : if( ql1(i,k) .gt. qsmall ) then
743 0 : QQnl1(i,k) = QQw1(i,k)*nl1(i,k)/ql1(i,k)
744 0 : QQnl1(i,k) = min(0._r8,cone*max(QQnl1(i,k),-nl1(i,k)/dt))
745 : else
746 0 : QQnl1(i,k) = 0._r8
747 : endif
748 : endif
749 0 : if( QQi1(i,k) .le. 0._r8 ) then
750 0 : if( qi1(i,k) .gt. qsmall ) then
751 0 : QQni1(i,k) = QQi1(i,k)*ni1(i,k)/qi1(i,k)
752 0 : QQni1(i,k) = min(0._r8,cone*max(QQni1(i,k),-ni1(i,k)/dt))
753 : else
754 0 : QQni1(i,k) = 0._r8
755 : endif
756 : endif
757 : enddo
758 : enddo
759 0 : nl_0(:ncol,top_lev:) = max(0._r8,nl1(:ncol,top_lev:)+QQnl1(:ncol,top_lev:)*dt)
760 0 : ni_0(:ncol,top_lev:) = max(0._r8,ni1(:ncol,top_lev:)+QQni1(:ncol,top_lev:)*dt)
761 :
762 : ! ----------------------------------------------------------------------------- !
763 : ! Check if non-cumulus pixels of '_05' state satisfies non-negative constraint. !
764 : ! If not, force all water substances of '_05' state to be positive by imposing !
765 : ! adjustive advection. We should use 'new' cumulus properties for this routine. !
766 : ! ----------------------------------------------------------------------------- !
767 :
768 0 : T_05(:ncol,top_lev:) = T_0(:ncol,top_lev:) + ( A_T(:ncol,top_lev:) + C_T(:ncol,top_lev:) ) * dt
769 0 : qv_05(:ncol,top_lev:) = qv_0(:ncol,top_lev:) + ( A_qv(:ncol,top_lev:) + C_qv(:ncol,top_lev:) ) * dt
770 0 : ql_05(:ncol,top_lev:) = ql_0(:ncol,top_lev:) + ( A_ql(:ncol,top_lev:) + C_ql(:ncol,top_lev:) ) * dt
771 0 : qi_05(:ncol,top_lev:) = qi_0(:ncol,top_lev:) + ( A_qi(:ncol,top_lev:) + C_qi(:ncol,top_lev:) ) * dt
772 0 : nl_05(:ncol,top_lev:) = max(0._r8, nl_0(:ncol,top_lev:) + ( A_nl(:ncol,top_lev:) + C_nl(:ncol,top_lev:) ) * dt )
773 0 : ni_05(:ncol,top_lev:) = max(0._r8, ni_0(:ncol,top_lev:) + ( A_ni(:ncol,top_lev:) + C_ni(:ncol,top_lev:) ) * dt )
774 :
775 : call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, &
776 : qv_05, ql_05, qi_05, T_05, A_qv_adj, &
777 0 : A_ql_adj, A_qi_adj, A_T_adj, do_cldice)
778 :
779 : ! -------------------------------------------------------------- !
780 : ! Define reference state at the first iteration. This will be !
781 : ! continuously updated within the iteration loop below. !
782 : ! While equlibrium state properties are already output from the !
783 : ! 'instratus_condensate', they will be re-computed within the !
784 : ! each iteration process. At the first iteration, they will !
785 : ! produce exactly identical results. Note that except at the !
786 : ! very first iteration iter = 1, we must use updated cumulus !
787 : ! properties at all the other iteration processes. Even at the !
788 : ! first iteration, we should use updated cumulus properties !
789 : ! when computing limiters for (Q,P,E). !
790 : ! -------------------------------------------------------------- !
791 :
792 : ! -------------------------------------------------------------- !
793 : ! Define variables at the reference state of the first iteration !
794 : ! -------------------------------------------------------------- !
795 :
796 0 : T(:ncol,top_lev:) = T_0(:ncol,top_lev:)
797 0 : qv(:ncol,top_lev:) = qv_0(:ncol,top_lev:)
798 0 : ql(:ncol,top_lev:) = ql_0(:ncol,top_lev:)
799 0 : qi(:ncol,top_lev:) = qi_0(:ncol,top_lev:)
800 0 : al_st(:ncol,top_lev:) = al_st_0(:ncol,top_lev:)
801 0 : ai_st(:ncol,top_lev:) = ai_st_0(:ncol,top_lev:)
802 0 : a_st(:ncol,top_lev:) = a_st_0(:ncol,top_lev:)
803 0 : ql_st(:ncol,top_lev:) = ql_st_0(:ncol,top_lev:)
804 0 : qi_st(:ncol,top_lev:) = qi_st_0(:ncol,top_lev:)
805 0 : nl(:ncol,top_lev:) = nl_0(:ncol,top_lev:)
806 0 : ni(:ncol,top_lev:) = ni_0(:ncol,top_lev:)
807 :
808 : ! -------------------------- !
809 : ! Main iterative computation !
810 : ! -------------------------- !
811 :
812 0 : do k = top_lev, pver
813 : call findsp_vc(qv_05(:ncol,k), T_05(:ncol,k), p(:ncol,k), .false., &
814 0 : Twb_aw(:ncol), qvwb_aw(:ncol,k))
815 0 : call qsat_water(T_05(1:ncol,k), p(1:ncol,k), esat_a(1:ncol), qsat_a(1:ncol,k), ncol)
816 : enddo
817 :
818 0 : do iter = 1, niter
819 :
820 : ! ------------------------------------------ !
821 : ! Initialize array within the iteration loop !
822 : ! ------------------------------------------ !
823 :
824 0 : QQ(:,:) = 0._r8
825 0 : QQw(:,:) = 0._r8
826 0 : QQi(:,:) = 0._r8
827 0 : QQnl(:,:) = 0._r8
828 0 : QQni(:,:) = 0._r8
829 0 : QQw2(:,:) = 0._r8
830 0 : QQi2(:,:) = 0._r8
831 0 : QQnl2(:,:) = 0._r8
832 0 : QQni2(:,:) = 0._r8
833 0 : nlten_pwi2(:,:) = 0._r8
834 0 : niten_pwi2(:,:) = 0._r8
835 0 : ACnl(:,:) = 0._r8
836 0 : ACni(:,:) = 0._r8
837 0 : aa(:,:) = 0._r8
838 0 : bb(:,:) = 0._r8
839 :
840 0 : do k = top_lev, pver
841 0 : call qsat_water(T(1:ncol,k), p(1:ncol,k), esat_b(1:ncol), qsat_b(1:ncol), ncol, dqsdt=dqsdT_b(1:ncol))
842 0 : if( iter .eq. 1 ) then
843 0 : a_cu(:ncol,k) = a_cud(:ncol,k)
844 : else
845 0 : a_cu(:ncol,k) = a_cu0(:ncol,k)
846 : endif
847 0 : do i = 1, ncol
848 0 : U(i,k) = qv(i,k)/qsat_b(i)
849 0 : U_nc(i,k) = U(i,k)
850 : enddo
851 : if( CAMstfrac ) then
852 : call astG_RHU(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol,&
853 : rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k))
854 : else
855 : call astG_PDF(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol,&
856 0 : rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k))
857 : endif
858 : call aist_vector(qv(:,k),T(:,k),p(:,k),qi(:,k),ni(:,k),landfrac(:),snowh(:),ai_st_nc(:,k),ncol,&
859 0 : rhmaxi_arr(:,k), rhmini_arr(:,k), rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k))
860 :
861 0 : ai_st(:ncol,k) = (1._r8-a_cu(:ncol,k))*ai_st_nc(:ncol,k)
862 0 : al_st(:ncol,k) = (1._r8-a_cu(:ncol,k))*al_st_nc(:ncol,k)
863 0 : a_st(:ncol,k) = max(al_st(:ncol,k),ai_st(:ncol,k))
864 :
865 0 : do i = 1, ncol
866 :
867 : ! -------------------------------------------------------- !
868 : ! Compute basic thermodynamic coefficients for computing Q !
869 : ! -------------------------------------------------------- !
870 :
871 0 : alpha = 1._r8/qsat_b(i)
872 0 : beta = dqsdT_b(i)*(qv(i,k)/qsat_b(i)**2)
873 0 : betast = alpha*dqsdT_b(i)
874 0 : gammal = alpha + (latvap/cpair)*beta
875 0 : gammai = alpha + ((latvap+latice)/cpair)*beta
876 0 : gammaQ = alpha + (latvap/cpair)*beta
877 0 : deltal = 1._r8 + a_st(i,k)*(latvap/cpair)*(betast/alpha)
878 0 : deltai = 1._r8 + a_st(i,k)*((latvap+latice)/cpair)*(betast/alpha)
879 0 : A_Tc = A_T(i,k)+A_T_adj(i,k)-(latvap/cpair)*(A_ql(i,k)+A_ql_adj(i,k))-((latvap+latice)/cpair)*(A_qi(i,k)+A_qi_adj(i,k))
880 0 : 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)
881 0 : C_Tc = C_T(i,k) - (latvap/cpair)*C_ql(i,k) - ((latvap+latice)/cpair)*C_qi(i,k)
882 0 : C_qt = C_qv(i,k) + C_ql(i,k) + C_qi(i,k)
883 0 : dTcdt = A_Tc + C_Tc
884 0 : dqtdt = A_qt + C_qt
885 : ! dqtstldt = A_qt + C_ql(i,k)/max(1.e-2_r8,al_st(i,k)) ! Original
886 : ! 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.
887 0 : dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_qlst(i,k) ! New 2 on Dec.30.2009.
888 : ! dqtstldt = A_qt + C_qt ! Original Conservative treatment
889 : ! 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
890 0 : dqidt = A_qi(i,k) + A_qi_adj(i,k) + C_qi(i,k)
891 :
892 0 : anic = max(1.e-8_r8,(1._r8-a_cu(i,k)))
893 0 : GG = G_nc(i,k)/anic
894 0 : aa(1,1) = gammal*al_st(i,k)
895 0 : aa(1,2) = GG + gammal*cc*ql_st(i,k)
896 0 : aa(2,1) = alpha + (latvap/cpair)*betast*al_st(i,k)
897 0 : aa(2,2) = (latvap/cpair)*betast*cc*ql_st(i,k)
898 0 : bb(1,1) = alpha*dqtdt - beta*dTcdt - gammai*dqidt - GG*al_st_nc(i,k)*dacudt(i,k) + F_nc(i,k)
899 0 : bb(2,1) = alpha*dqtstldt - betast*(dTcdt + ((latvap+latice)/cpair)*dqidt)
900 0 : call gaussj(aa(1:2,1:2),2,2,bb(1:2,1),1,1)
901 0 : dqlstdt = bb(1,1)
902 0 : dalstdt = bb(2,1)
903 0 : 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) )
904 :
905 : ! ------------------------------------------------------------ !
906 : ! Limiter for QQ !
907 : ! Here, 'fice' should be from the reference equilibrium state !
908 : ! since QQ itself is computed from the reference state. !
909 : ! From the assumption used for derivation of QQ(i), it must be !
910 : ! that QQw(i) = QQ(i)*(1._r8-fice(i)), QQi(i) = QQ(i)*fice(i) !
911 : ! ------------------------------------------------------------ !
912 :
913 0 : if( QQ(i,k) .ge. 0._r8 ) then
914 0 : QQmax = (qv_05(i,k) - qmin(1))/dt ! For ghost cumulus & semi-ghost ice stratus
915 0 : QQmax = max(0._r8,QQmax)
916 0 : QQ(i,k) = min(QQ(i,k),QQmax)
917 0 : QQw(i,k) = QQ(i,k)
918 0 : QQi(i,k) = 0._r8
919 : else
920 0 : QQmin = 0._r8
921 0 : if( qv_05(i,k) .lt. qsat_a(i,k) ) QQmin = min(0._r8,cone*(qv_05(i,k)-qvwb_aw(i,k))/dt)
922 0 : QQ(i,k) = max(QQ(i,k),QQmin)
923 0 : QQw(i,k) = QQ(i,k)
924 0 : QQi(i,k) = 0._r8
925 0 : QQwmin = min(0._r8,-cone*ql_05(i,k)/dt)
926 0 : QQimin = min(0._r8,-cone*qi_05(i,k)/dt)
927 0 : QQw(i,k) = min(0._r8,max(QQw(i,k),QQwmin))
928 : QQi(i,k) = min(0._r8,max(QQi(i,k),QQimin))
929 : endif
930 :
931 : ! -------------------------------------------------- !
932 : ! Reduce droplet concentration if evaporation occurs !
933 : ! Note 'QQnl1,QQni1' are computed from the reference !
934 : ! equilibrium state but limiter is from 'nl_05'. !
935 : ! -------------------------------------------------- !
936 :
937 0 : if( QQw(i,k) .lt. 0._r8 ) then
938 0 : if( ql_05(i,k) .gt. qsmall ) then
939 0 : QQnl(i,k) = QQw(i,k)*nl_05(i,k)/ql_05(i,k)
940 0 : QQnl(i,k) = min(0._r8,cone*max(QQnl(i,k),-nl_05(i,k)/dt))
941 : else
942 0 : QQnl(i,k) = 0._r8
943 : endif
944 : endif
945 :
946 0 : if( QQi(i,k) .lt. 0._r8 ) then
947 0 : if( qi_05(i,k) .gt. qsmall ) then
948 0 : QQni(i,k) = QQi(i,k)*ni_05(i,k)/qi_05(i,k)
949 0 : QQni(i,k) = min(0._r8,cone*max(QQni(i,k),-ni_05(i,k)/dt))
950 : else
951 0 : QQni(i,k) = 0._r8
952 : endif
953 : endif
954 :
955 : enddo
956 : enddo
957 :
958 : ! -------------------------------------------------------------------- !
959 : ! Until now, we have finished computing all necessary tendencies !
960 : ! from the equilibrium input state (T_0). !
961 : ! If ramda = 0 : fully explicit scheme !
962 : ! ramda = 1 : fully implicit scheme !
963 : ! Note that 'ramda = 0.5 with niter = 2' can mimic !
964 : ! -------------------------------------------------------------------- !
965 :
966 0 : if( iter .eq. 1 ) then
967 0 : QQw_prev(:ncol,top_lev:) = QQw(:ncol,top_lev:)
968 0 : QQi_prev(:ncol,top_lev:) = QQi(:ncol,top_lev:)
969 0 : QQnl_prev(:ncol,top_lev:) = QQnl(:ncol,top_lev:)
970 0 : QQni_prev(:ncol,top_lev:) = QQni(:ncol,top_lev:)
971 : endif
972 :
973 0 : QQw_prog(:ncol,top_lev:) = ramda*QQw(:ncol,top_lev:) + (1._r8-ramda)*QQw_prev(:ncol,top_lev:)
974 0 : QQi_prog(:ncol,top_lev:) = ramda*QQi(:ncol,top_lev:) + (1._r8-ramda)*QQi_prev(:ncol,top_lev:)
975 0 : QQnl_prog(:ncol,top_lev:) = ramda*QQnl(:ncol,top_lev:) + (1._r8-ramda)*QQnl_prev(:ncol,top_lev:)
976 0 : QQni_prog(:ncol,top_lev:) = ramda*QQni(:ncol,top_lev:) + (1._r8-ramda)*QQni_prev(:ncol,top_lev:)
977 :
978 0 : QQw_prev(:ncol,top_lev:) = QQw_prog(:ncol,top_lev:)
979 0 : QQi_prev(:ncol,top_lev:) = QQi_prog(:ncol,top_lev:)
980 0 : QQnl_prev(:ncol,top_lev:) = QQnl_prog(:ncol,top_lev:)
981 0 : QQni_prev(:ncol,top_lev:) = QQni_prog(:ncol,top_lev:)
982 :
983 : ! -------------------------------------------------------- !
984 : ! Compute final prognostic state on which final diagnostic !
985 : ! in-stratus condensate adjustment is applied in the below.!
986 : ! Important : I must check whether there are any external !
987 : ! advective forcings of 'A_nl(i,k),A_ni(i,k)'. !
988 : ! Even they are (i.e., advection of aerosol), !
989 : ! actual droplet activation will be performd !
990 : ! in microphysics, so it will be completely !
991 : ! reasonable to 'A_nl(i,k)=A_ni(i,k)=0'. !
992 : ! -------------------------------------------------------- !
993 :
994 0 : do k = top_lev, pver
995 0 : do i = 1, ncol
996 0 : T_prime0(i,k) = T_0(i,k) + dt*( A_T(i,k) + A_T_adj(i,k) + C_T(i,k) + &
997 0 : (latvap*QQw_prog(i,k)+(latvap+latice)*QQi_prog(i,k))/cpair )
998 0 : 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) )
999 0 : 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) )
1000 0 : 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) )
1001 0 : nl_prime0(i,k) = max(0._r8,nl_0(i,k) + dt*( A_nl(i,k) + C_nl(i,k) + QQnl_prog(i,k) ))
1002 0 : ni_prime0(i,k) = max(0._r8,ni_0(i,k) + dt*( A_ni(i,k) + C_ni(i,k) + QQni_prog(i,k) ))
1003 0 : if( ql_prime0(i,k) .lt. qsmall ) nl_prime0(i,k) = 0._r8
1004 0 : if( qi_prime0(i,k) .lt. qsmall ) ni_prime0(i,k) = 0._r8
1005 : enddo
1006 : enddo
1007 :
1008 : ! -------------------------------------------------- !
1009 : ! Perform diagnostic 'positive_moisture' constraint. !
1010 : ! -------------------------------------------------- !
1011 :
1012 0 : T_dprime(:ncol,top_lev:) = T_prime0(:ncol,top_lev:)
1013 0 : qv_dprime(:ncol,top_lev:) = qv_prime0(:ncol,top_lev:)
1014 0 : ql_dprime(:ncol,top_lev:) = ql_prime0(:ncol,top_lev:)
1015 0 : qi_dprime(:ncol,top_lev:) = qi_prime0(:ncol,top_lev:)
1016 0 : nl_dprime(:ncol,top_lev:) = nl_prime0(:ncol,top_lev:)
1017 0 : ni_dprime(:ncol,top_lev:) = ni_prime0(:ncol,top_lev:)
1018 :
1019 : call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, &
1020 : qv_dprime, ql_dprime, qi_dprime, T_dprime, &
1021 0 : qvten_pwi2, qlten_pwi2, qiten_pwi2, Tten_pwi2, do_cldice)
1022 :
1023 0 : do k = top_lev, pver
1024 0 : do i = 1, ncol
1025 0 : if( ql_dprime(i,k) .lt. qsmall ) then
1026 0 : nlten_pwi2(i,k) = -nl_dprime(i,k)/dt
1027 0 : nl_dprime(i,k) = 0._r8
1028 : endif
1029 0 : if( qi_dprime(i,k) .lt. qsmall ) then
1030 0 : niten_pwi2(i,k) = -ni_dprime(i,k)/dt
1031 0 : ni_dprime(i,k) = 0._r8
1032 : endif
1033 : enddo
1034 : enddo
1035 :
1036 : ! -------------------------------------------------------------- !
1037 : ! Add tendency associated with detrainment of cumulus condensate !
1038 : ! This tendency is not used in computing Q !
1039 : ! Since D_ql,D_qi,D_nl,D_ni > 0, don't need to worry about !
1040 : ! negative scalar. !
1041 : ! This tendency is not reflected into Fzs2, which is OK. !
1042 : ! -------------------------------------------------------------- !
1043 :
1044 0 : T_dprime(:ncol,top_lev:) = T_dprime(:ncol,top_lev:) + D_T(:ncol,top_lev:) * dt
1045 0 : qv_dprime(:ncol,top_lev:) = qv_dprime(:ncol,top_lev:) + D_qv(:ncol,top_lev:) * dt
1046 0 : ql_dprime(:ncol,top_lev:) = ql_dprime(:ncol,top_lev:) + D_ql(:ncol,top_lev:) * dt
1047 0 : qi_dprime(:ncol,top_lev:) = qi_dprime(:ncol,top_lev:) + D_qi(:ncol,top_lev:) * dt
1048 0 : nl_dprime(:ncol,top_lev:) = nl_dprime(:ncol,top_lev:) + D_nl(:ncol,top_lev:) * dt
1049 0 : ni_dprime(:ncol,top_lev:) = ni_dprime(:ncol,top_lev:) + D_ni(:ncol,top_lev:) * dt
1050 :
1051 : ! ---------------------------------------------------------- !
1052 : ! Impose diagnostic upper and lower limits on the in-stratus !
1053 : ! condensate amount. This produces a final equilibrium state !
1054 : ! at the end of each iterative process. !
1055 : ! ---------------------------------------------------------- !
1056 :
1057 0 : do k = top_lev, pver
1058 : call instratus_condensate( lchnk , ncol , k , p(:,k) , &
1059 : T_dprime(:,k) , qv_dprime(:,k) , ql_dprime(:,k) , qi_dprime(:,k), &
1060 : ni_dprime(:,k) , &
1061 : a_cu0(:,k) , zeros(:,k) , zeros(:,k) , &
1062 : zeros(:,k) , zeros(:,k) , zeros(:,k) , &
1063 : landfrac , snowh , &
1064 : rhmaxi_arr(:,k),rhmini_arr(:,k), rhminl_arr(:,k), rhminl_adj_land_arr(:,k), rhminh_arr(:,k), &
1065 : T_star(:,k) , qv_star(:,k) , ql_star(:,k) , qi_star(:,k) , &
1066 0 : al_st_star(:,k), ai_st_star(:,k), ql_st_star(:,k), qi_st_star(:,k) )
1067 0 : a_st_star(:ncol,k) = max(al_st_star(:ncol,k),ai_st_star(:ncol,k))
1068 0 : QQw2(:ncol,k) = (ql_star(:ncol,k) - ql_dprime(:ncol,k))/dt
1069 0 : QQi2(:ncol,k) = (qi_star(:ncol,k) - qi_dprime(:ncol,k))/dt
1070 : ! -------------------------------------------------- !
1071 : ! Reduce droplet concentration if evaporation occurs !
1072 : ! -------------------------------------------------- !
1073 0 : do i = 1, ncol
1074 0 : if( QQw2(i,k) .le. 0._r8 ) then
1075 0 : if( ql_dprime(i,k) .ge. qsmall ) then
1076 0 : QQnl2(i,k) = QQw2(i,k)*nl_dprime(i,k)/ql_dprime(i,k)
1077 0 : QQnl2(i,k) = min(0._r8,cone*max(QQnl2(i,k),-nl_dprime(i,k)/dt))
1078 : else
1079 0 : QQnl2(i,k) = 0._r8
1080 : endif
1081 : endif
1082 0 : if( QQi2(i,k) .le. 0._r8 ) then
1083 0 : if( qi_dprime(i,k) .gt. qsmall ) then
1084 0 : QQni2(i,k) = QQi2(i,k)*ni_dprime(i,k)/qi_dprime(i,k)
1085 0 : QQni2(i,k) = min(0._r8,cone*max(QQni2(i,k),-ni_dprime(i,k)/dt))
1086 : else
1087 0 : QQni2(i,k) = 0._r8
1088 : endif
1089 : endif
1090 : enddo
1091 : enddo
1092 0 : nl_star(:ncol,top_lev:) = max(0._r8,nl_dprime(:ncol,top_lev:)+QQnl2(:ncol,top_lev:)*dt)
1093 0 : ni_star(:ncol,top_lev:) = max(0._r8,ni_dprime(:ncol,top_lev:)+QQni2(:ncol,top_lev:)*dt)
1094 :
1095 : ! ------------------------------------------ !
1096 : ! Final adjustment of droplet concentration. !
1097 : ! Set # to zero if there is no cloud. !
1098 : ! ------------------------------------------ !
1099 :
1100 0 : do k = top_lev, pver
1101 0 : do i = 1, ncol
1102 0 : if( ql_star(i,k) .lt. qsmall ) then
1103 0 : ACnl(i,k) = - nl_star(i,k)/dt
1104 0 : nl_star(i,k) = 0._r8
1105 : endif
1106 0 : if( qi_star(i,k) .lt. qsmall ) then
1107 0 : ACni(i,k) = - ni_star(i,k)/dt
1108 0 : ni_star(i,k) = 0._r8
1109 : endif
1110 : enddo
1111 : enddo
1112 :
1113 : ! ----------------------------------------------------- !
1114 : ! Define equilibrium reference state for next iteration !
1115 : ! ----------------------------------------------------- !
1116 :
1117 0 : T(:ncol,top_lev:) = T_star(:ncol,top_lev:)
1118 0 : qv(:ncol,top_lev:) = qv_star(:ncol,top_lev:)
1119 0 : ql(:ncol,top_lev:) = ql_star(:ncol,top_lev:)
1120 0 : qi(:ncol,top_lev:) = qi_star(:ncol,top_lev:)
1121 0 : al_st(:ncol,top_lev:) = al_st_star(:ncol,top_lev:)
1122 0 : ai_st(:ncol,top_lev:) = ai_st_star(:ncol,top_lev:)
1123 0 : a_st(:ncol,top_lev:) = a_st_star(:ncol,top_lev:)
1124 0 : ql_st(:ncol,top_lev:) = ql_st_star(:ncol,top_lev:)
1125 0 : qi_st(:ncol,top_lev:) = qi_st_star(:ncol,top_lev:)
1126 0 : nl(:ncol,top_lev:) = nl_star(:ncol,top_lev:)
1127 0 : ni(:ncol,top_lev:) = ni_star(:ncol,top_lev:)
1128 :
1129 : enddo ! End of 'iter' prognostic iterative computation
1130 :
1131 : ! ------------------------------------------------------------------------ !
1132 : ! Compute final tendencies of main output variables and diagnostic outputs !
1133 : ! Note that the very input state [T0,qv0,ql0,qi0] are !
1134 : ! marched to [T_star,qv_star,ql_star,qi_star] with equilibrium !
1135 : ! stratus informations of [a_st_star,ql_st_star,qi_st_star] by !
1136 : ! below final tendencies and [A_T,A_qv,A_ql,A_qi] !
1137 : ! ------------------------------------------------------------------------ !
1138 :
1139 : ! ------------------ !
1140 : ! Process tendencies !
1141 : ! ------------------ !
1142 :
1143 0 : QQw_final(:ncol,top_lev:) = QQw_prog(:ncol,top_lev:)
1144 0 : QQi_final(:ncol,top_lev:) = QQi_prog(:ncol,top_lev:)
1145 0 : QQ_final(:ncol,top_lev:) = QQw_final(:ncol,top_lev:) + QQi_final(:ncol,top_lev:)
1146 : QQw_all(:ncol,top_lev:) = QQw_prog(:ncol,top_lev:) + QQw1(:ncol,top_lev:) + QQw2(:ncol,top_lev:) + &
1147 0 : qlten_pwi1(:ncol,top_lev:) + qlten_pwi2(:ncol,top_lev:) + A_ql_adj(:ncol,top_lev:)
1148 : QQi_all(:ncol,top_lev:) = QQi_prog(:ncol,top_lev:) + QQi1(:ncol,top_lev:) + QQi2(:ncol,top_lev:) + &
1149 0 : qiten_pwi1(:ncol,top_lev:) + qiten_pwi2(:ncol,top_lev:) + A_qi_adj(:ncol,top_lev:)
1150 0 : QQ_all(:ncol,top_lev:) = QQw_all(:ncol,top_lev:) + QQi_all(:ncol,top_lev:)
1151 0 : QQnl_final(:ncol,top_lev:) = QQnl_prog(:ncol,top_lev:)
1152 0 : QQni_final(:ncol,top_lev:) = QQni_prog(:ncol,top_lev:)
1153 0 : QQn_final(:ncol,top_lev:) = QQnl_final(:ncol,top_lev:) + QQni_final(:ncol,top_lev:)
1154 : QQnl_all(:ncol,top_lev:) = QQnl_prog(:ncol,top_lev:) + QQnl1(:ncol,top_lev:) + QQnl2(:ncol,top_lev:) + &
1155 0 : nlten_pwi1(:ncol,top_lev:) + nlten_pwi2(:ncol,top_lev:) + ACnl(:ncol,top_lev:) + A_nl_adj(:ncol,top_lev:)
1156 : QQni_all(:ncol,top_lev:) = QQni_prog(:ncol,top_lev:) + QQni1(:ncol,top_lev:) + QQni2(:ncol,top_lev:) + &
1157 0 : niten_pwi1(:ncol,top_lev:) + niten_pwi2(:ncol,top_lev:) + ACni(:ncol,top_lev:) + A_ni_adj(:ncol,top_lev:)
1158 0 : QQn_all(:ncol,top_lev:) = QQnl_all(:ncol,top_lev:) + QQni_all(:ncol,top_lev:)
1159 0 : qme(:ncol,top_lev:) = QQ_final(:ncol,top_lev:)
1160 0 : qvadj(:ncol,top_lev:) = qvten_pwi1(:ncol,top_lev:) + qvten_pwi2(:ncol,top_lev:) + A_qv_adj(:ncol,top_lev:)
1161 0 : qladj(:ncol,top_lev:) = qlten_pwi1(:ncol,top_lev:) + qlten_pwi2(:ncol,top_lev:) + A_ql_adj(:ncol,top_lev:)
1162 0 : qiadj(:ncol,top_lev:) = qiten_pwi1(:ncol,top_lev:) + qiten_pwi2(:ncol,top_lev:) + A_qi_adj(:ncol,top_lev:)
1163 0 : qllim(:ncol,top_lev:) = QQw1 (:ncol,top_lev:) + QQw2 (:ncol,top_lev:)
1164 0 : qilim(:ncol,top_lev:) = QQi1 (:ncol,top_lev:) + QQi2 (:ncol,top_lev:)
1165 :
1166 : ! ----------------- !
1167 : ! Output tendencies !
1168 : ! ----------------- !
1169 :
1170 : s_tendout(:ncol,top_lev:) = cpair*( T_star(:ncol,top_lev:) - T0(:ncol,top_lev:) )/dt - &
1171 0 : cpair*(A_T(:ncol,top_lev:)+C_T(:ncol,top_lev:))
1172 : qv_tendout(:ncol,top_lev:) = ( qv_star(:ncol,top_lev:) - qv0(:ncol,top_lev:) )/dt - &
1173 0 : (A_qv(:ncol,top_lev:)+C_qv(:ncol,top_lev:))
1174 : ql_tendout(:ncol,top_lev:) = ( ql_star(:ncol,top_lev:) - ql0(:ncol,top_lev:) )/dt - &
1175 0 : (A_ql(:ncol,top_lev:)+C_ql(:ncol,top_lev:))
1176 : qi_tendout(:ncol,top_lev:) = ( qi_star(:ncol,top_lev:) - qi0(:ncol,top_lev:) )/dt - &
1177 0 : (A_qi(:ncol,top_lev:)+C_qi(:ncol,top_lev:))
1178 : nl_tendout(:ncol,top_lev:) = ( nl_star(:ncol,top_lev:) - nl0(:ncol,top_lev:) )/dt - &
1179 0 : (A_nl(:ncol,top_lev:)+C_nl(:ncol,top_lev:))
1180 : ni_tendout(:ncol,top_lev:) = ( ni_star(:ncol,top_lev:) - ni0(:ncol,top_lev:) )/dt - &
1181 0 : (A_ni(:ncol,top_lev:)+C_ni(:ncol,top_lev:))
1182 :
1183 0 : if (.not. do_cldice) then
1184 0 : do k = top_lev, pver
1185 0 : do i = 1, ncol
1186 :
1187 : ! Don't want either qi or ni tendencies, but the code above is somewhat convoluted and
1188 : ! is trying to adjust both (small numbers). Just force it to zero here.
1189 0 : qi_tendout(i,k) = 0._r8
1190 0 : ni_tendout(i,k) = 0._r8
1191 : end do
1192 : end do
1193 : end if
1194 :
1195 : ! ------------------ !
1196 : ! Net cloud fraction !
1197 : ! ------------------ !
1198 :
1199 0 : cld(:ncol,top_lev:) = a_st_star(:ncol,top_lev:) + a_cu0(:ncol,top_lev:)
1200 :
1201 : ! --------------------------------- !
1202 : ! Updated grid-mean state variables !
1203 : ! --------------------------------- !
1204 :
1205 0 : T0(:ncol,top_lev:) = T_star(:ncol,top_lev:)
1206 0 : qv0(:ncol,top_lev:) = qv_star(:ncol,top_lev:)
1207 0 : ql0(:ncol,top_lev:) = ql_star(:ncol,top_lev:)
1208 0 : qi0(:ncol,top_lev:) = qi_star(:ncol,top_lev:)
1209 0 : nl0(:ncol,top_lev:) = nl_star(:ncol,top_lev:)
1210 0 : ni0(:ncol,top_lev:) = ni_star(:ncol,top_lev:)
1211 :
1212 0 : if (hist_fld_active('RHMIN_LIQ')) then
1213 : ! Compute default critical RH as a function of height and surface type as in the current code.
1214 0 : rhmin_liq_diag(:,:) = 0._r8
1215 0 : do k = top_lev, pver
1216 0 : do i = 1, ncol
1217 0 : land = nint(landfrac(i)) == 1
1218 0 : if( p(i,k) .ge. premib ) then
1219 0 : if( land .and. (snowh(i).le.0.000001_r8) ) then
1220 0 : rhmin_liq_diag(i,k) = rhminl_const - rhminl_adj_land_const
1221 : else
1222 0 : rhmin_liq_diag(i,k) = rhminl_const
1223 : endif
1224 0 : elseif( p(i,k) .lt. premit ) then
1225 0 : rhmin_liq_diag(i,k) = rhminh_const
1226 : else
1227 0 : tmp = (premib-(max(p(i,k),premit)))/(premib-premit)
1228 0 : rhmin_liq_diag(i,k) = rhminh_const*tmp + rhminl_const*(1.0_r8-tmp)
1229 : endif
1230 : end do
1231 : end do
1232 0 : call outfld( 'RHMIN_LIQ', rhmin_liq_diag, pcols, lchnk )
1233 : end if
1234 :
1235 0 : rhmin_ice_diag(:,:) = rhminh_const
1236 0 : call outfld( 'RHMIN_ICE', rhmin_ice_diag, pcols, lchnk )
1237 :
1238 0 : call outfld( 'DRHMINPBL_LIQ', d_rhmin_liq_PBL, pcols, lchnk )
1239 0 : call outfld( 'DRHMINPBL_ICE', d_rhmin_ice_PBL, pcols, lchnk )
1240 0 : call outfld( 'DRHMINDET_LIQ', d_rhmin_liq_det, pcols, lchnk )
1241 0 : call outfld( 'DRHMINDET_ICE', d_rhmin_ice_det, pcols, lchnk )
1242 :
1243 0 : end subroutine mmacro_pcond
1244 :
1245 :
1246 : !=======================================================================================================
1247 :
1248 0 : subroutine rhcrit_calc( &
1249 : ncol, dp, T0, p, &
1250 : clrw_old, clri_old, tke, qtl_flx, &
1251 : qti_flx, cmfr_det, qlr_det, qir_det, &
1252 : rhmaxi_arr, rhmini_arr, rhminl_arr, rhminl_adj_land_arr, rhminh_arr, &
1253 : d_rhmin_liq_PBL, d_rhmin_ice_PBL, d_rhmin_liq_det, d_rhmin_ice_det)
1254 :
1255 : ! ------------------------------------------------- !
1256 : ! Compute a drop of critical RH for stratus by !
1257 : ! (1) PBL turbulence, and !
1258 : ! (2) convective detrainment. !
1259 : ! Note that all of 'd_rhmin...' terms are positive. !
1260 : ! ------------------------------------------------- !
1261 :
1262 : integer, intent(in) :: ncol ! Number of active columns
1263 : real(r8), intent(in) :: dp(pcols,pver) ! Pressure thickness [Pa] > 0
1264 : real(r8), intent(in) :: T0(pcols,pver) ! Temperature [K]
1265 : real(r8), intent(in) :: p(pcols,pver) ! Pressure at the layer mid-point [Pa]
1266 : real(r8), intent(in) :: clrw_old(pcols,pver) ! Clear sky fraction at the previous time step for liquid stratus process
1267 : real(r8), intent(in) :: clri_old(pcols,pver) ! Clear sky fraction at the previous time step for ice stratus process
1268 : real(r8), pointer :: tke(:,:) ! (pcols,pverp) TKE from the PBL scheme
1269 : real(r8), pointer :: qtl_flx(:,:) ! (pcols,pverp) overbar(w'qtl') from PBL scheme where qtl = qv + ql
1270 : real(r8), pointer :: qti_flx(:,:) ! (pcols,pverp) overbar(w'qti') from PBL scheme where qti = qv + qi
1271 : real(r8), pointer :: cmfr_det(:,:) ! (pcols,pver) Detrained mass flux from the convection scheme
1272 : real(r8), pointer :: qlr_det(:,:) ! (pcols,pver) Detrained ql from the convection scheme
1273 : real(r8), pointer :: qir_det(:,:) ! (pcols,pver) Detrained qi from the convection scheme
1274 :
1275 : real(r8), intent(in) :: rhmaxi_arr(pcols,pver)
1276 :
1277 : real(r8), intent(out) :: rhmini_arr(pcols,pver)
1278 : real(r8), intent(out) :: rhminl_arr(pcols,pver)
1279 : real(r8), intent(out) :: rhminl_adj_land_arr(pcols,pver)
1280 : real(r8), intent(out) :: rhminh_arr(pcols,pver)
1281 : real(r8), intent(out) :: d_rhmin_liq_PBL(pcols,pver)
1282 : real(r8), intent(out) :: d_rhmin_ice_PBL(pcols,pver)
1283 : real(r8), intent(out) :: d_rhmin_liq_det(pcols,pver)
1284 : real(r8), intent(out) :: d_rhmin_ice_det(pcols,pver)
1285 :
1286 : ! local variables
1287 :
1288 : integer :: i, k
1289 :
1290 : real(r8) :: esat_tmp(pcols) ! Dummy for saturation vapor pressure calc.
1291 : real(r8) :: qsat_tmp(pcols) ! Saturation water vapor specific humidity [kg/kg]
1292 : real(r8) :: sig_tmp
1293 : !---------------------------------------------------------------------------------------------------
1294 :
1295 :
1296 :
1297 : ! ---------------------------------- !
1298 : ! Calc critical RH for ice stratus !
1299 : ! ---------------------------------- !
1300 :
1301 0 : rhmini_arr(:,:) = rhmini_const
1302 :
1303 0 : if (i_rhmini > 0) then
1304 :
1305 : ! Compute the drop of critical RH by convective detrainment of cloud condensate
1306 :
1307 0 : do k = top_lev, pver
1308 0 : do i = 1, ncol
1309 0 : d_rhmin_ice_det(i,k) = tau_deti*(gravit/dp(i,k))*cmfr_det(i,k)*clri_old(i,k)*qir_det(i,k)*3.6e6_r8
1310 0 : d_rhmin_ice_det(i,k) = max(0._r8,min(0.5_r8,d_rhmin_ice_det(i,k)))
1311 : end do
1312 : end do
1313 :
1314 0 : if (i_rhmini == 1) then
1315 0 : rhmini_arr(:ncol,:) = rhmini_const - d_rhmin_ice_det(:ncol,:)
1316 : end if
1317 :
1318 : end if
1319 :
1320 0 : if (i_rhmini == 2) then
1321 :
1322 : ! Compute the drop of critical RH by the variability induced by PBL turbulence
1323 0 : do k = top_lev, pver
1324 0 : call qsat_ice(T0(1:ncol,k), p(1:ncol,k), esat_tmp(1:ncol), qsat_tmp(1:ncol), ncol)
1325 0 : do i = 1, ncol
1326 0 : sig_tmp = 0.5_r8 * ( qti_flx(i,k) / sqrt(max(qsmall,tke(i,k))) + &
1327 0 : qti_flx(i,k+1) / sqrt(max(qsmall,tke(i,k+1))) )
1328 0 : d_rhmin_ice_PBL(i,k) = c_aniso*sig_tmp/max(qsmall,qsat_tmp(i))
1329 0 : d_rhmin_ice_PBL(i,k) = max(0._r8,min(0.5_r8,d_rhmin_ice_PBL(i,k)))
1330 :
1331 0 : rhmini_arr(i,k) = 1._r8 - d_rhmin_ice_PBL(i,k) - d_rhmin_ice_det(i,k)
1332 : end do
1333 : end do
1334 : end if
1335 :
1336 0 : if (i_rhmini > 0) then
1337 0 : do k = top_lev, pver
1338 0 : do i = 1, ncol
1339 0 : rhmini_arr(i,k) = max(0._r8,min(rhmaxi_arr(i,k),rhmini_arr(i,k)))
1340 : end do
1341 : end do
1342 : end if
1343 :
1344 : ! ------------------------------------- !
1345 : ! Choose critical RH for liquid stratus !
1346 : ! ------------------------------------- !
1347 :
1348 0 : rhminl_arr(:,:) = rhminl_const
1349 0 : rhminl_adj_land_arr(:,:) = rhminl_adj_land_const
1350 0 : rhminh_arr(:,:) = rhminh_const
1351 :
1352 0 : if (i_rhminl > 0) then
1353 :
1354 : ! Compute the drop of critical RH by convective detrainment of cloud condensate
1355 :
1356 0 : do k = top_lev, pver
1357 0 : do i = 1, ncol
1358 0 : d_rhmin_liq_det(i,k) = tau_detw*(gravit/dp(i,k))*cmfr_det(i,k)*clrw_old(i,k)*qlr_det(i,k)*3.6e6_r8
1359 0 : d_rhmin_liq_det(i,k) = max(0._r8,min(0.5_r8,d_rhmin_liq_det(i,k)))
1360 : end do
1361 : end do
1362 :
1363 0 : if (i_rhminl == 1) then
1364 0 : rhminl_arr(:ncol,top_lev:) = rhminl_const - d_rhmin_liq_det(:ncol,top_lev:)
1365 0 : rhminh_arr(:ncol,top_lev:) = rhminh_const - d_rhmin_liq_det(:ncol,top_lev:)
1366 : end if
1367 :
1368 : end if
1369 :
1370 0 : if (i_rhminl == 2) then
1371 :
1372 : ! Compute the drop of critical RH by the variability induced by PBL turbulence
1373 0 : do k = top_lev, pver
1374 0 : call qsat_water(T0(1:ncol,k), p(1:ncol,k), esat_tmp(1:ncol), qsat_tmp(1:ncol), ncol)
1375 0 : do i = 1, ncol
1376 0 : sig_tmp = 0.5_r8 * ( qtl_flx(i,k) / sqrt(max(qsmall,tke(i,k))) + &
1377 0 : qtl_flx(i,k+1) / sqrt(max(qsmall,tke(i,k+1))) )
1378 0 : d_rhmin_liq_PBL(i,k) = c_aniso*sig_tmp/max(qsmall,qsat_tmp(i))
1379 0 : d_rhmin_liq_PBL(i,k) = max(0._r8,min(0.5_r8,d_rhmin_liq_PBL(i,k)))
1380 :
1381 0 : rhminl_arr(i,k) = 1._r8 - d_rhmin_liq_PBL(i,k) - d_rhmin_liq_det(i,k)
1382 0 : rhminl_adj_land_arr(i,k) = 0._r8
1383 0 : rhminh_arr(i,k) = rhminl_arr(i,k)
1384 : end do
1385 : end do
1386 : end if
1387 :
1388 0 : if (i_rhminl > 0) then
1389 0 : do k = top_lev, pver
1390 0 : do i = 1, ncol
1391 0 : rhminl_arr(i,k) = max(rhminl_adj_land_arr(i,k),min(1._r8,rhminl_arr(i,k)))
1392 0 : rhminh_arr(i,k) = max(0._r8,min(1._r8,rhminh_arr(i,k)))
1393 : end do
1394 : end do
1395 : end if
1396 :
1397 0 : end subroutine rhcrit_calc
1398 :
1399 : !=======================================================================================================
1400 :
1401 0 : subroutine instratus_condensate( lchnk, ncol, k, &
1402 : p_in, T0_in, qv0_in, ql0_in, qi0_in, &
1403 : ni0_in, &
1404 : a_dc_in, ql_dc_in, qi_dc_in, &
1405 : a_sc_in, ql_sc_in, qi_sc_in, &
1406 : landfrac, snowh, &
1407 : rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, &
1408 : T_out, qv_out, ql_out, qi_out, &
1409 : al_st_out, ai_st_out, ql_st_out, qi_st_out )
1410 :
1411 : ! ------------------------------------------------------- !
1412 : ! Diagnostically force in-stratus condensate to be !
1413 : ! in the range of 'qlst_min < qc_st < qlst_max' !
1414 : ! whenever stratus exists in the equilibrium state !
1415 : ! ------------------------------------------------------- !
1416 :
1417 : integer, intent(in) :: lchnk ! Chunk identifier
1418 : integer, intent(in) :: ncol ! Number of atmospheric columns
1419 : integer, intent(in) :: k ! Layer index
1420 :
1421 : real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa]
1422 : real(r8), intent(in) :: T0_in(pcols) ! Temperature [K]
1423 : real(r8), intent(in) :: qv0_in(pcols) ! Grid-mean water vapor [kg/kg]
1424 : real(r8), intent(in) :: ql0_in(pcols) ! Grid-mean LWC [kg/kg]
1425 : real(r8), intent(in) :: qi0_in(pcols) ! Grid-mean IWC [kg/kg]
1426 : real(r8), intent(in) :: ni0_in(pcols)
1427 :
1428 : real(r8), intent(in) :: a_dc_in(pcols) ! Deep cumulus cloud fraction
1429 : real(r8), intent(in) :: ql_dc_in(pcols) ! In-deep cumulus LWC [kg/kg]
1430 : real(r8), intent(in) :: qi_dc_in(pcols) ! In-deep cumulus IWC [kg/kg]
1431 : real(r8), intent(in) :: a_sc_in(pcols) ! Shallow cumulus cloud fraction
1432 : real(r8), intent(in) :: ql_sc_in(pcols) ! In-shallow cumulus LWC [kg/kg]
1433 : real(r8), intent(in) :: qi_sc_in(pcols) ! In-shallow cumulus IWC [kg/kg]
1434 :
1435 : real(r8), intent(in) :: landfrac(pcols) ! Land fraction
1436 : real(r8), intent(in) :: snowh(pcols) ! Snow depth (liquid water equivalent)
1437 :
1438 : real(r8), intent(in) :: rhmaxi_in(pcols)
1439 : real(r8), intent(in) :: rhmini_in(pcols)
1440 : real(r8), intent(in) :: rhminl_in(pcols)
1441 : real(r8), intent(in) :: rhminl_adj_land_in(pcols)
1442 : real(r8), intent(in) :: rhminh_in(pcols)
1443 :
1444 : real(r8), intent(out) :: T_out(pcols) ! Temperature [K]
1445 : real(r8), intent(out) :: qv_out(pcols) ! Grid-mean water vapor [kg/kg]
1446 : real(r8), intent(out) :: ql_out(pcols) ! Grid-mean LWC [kg/kg]
1447 : real(r8), intent(out) :: qi_out(pcols) ! Grid-mean IWC [kg/kg]
1448 :
1449 : real(r8), intent(out) :: al_st_out(pcols) ! Liquid stratus fraction
1450 : real(r8), intent(out) :: ai_st_out(pcols) ! Ice stratus fraction
1451 : real(r8), intent(out) :: ql_st_out(pcols) ! In-stratus LWC [kg/kg]
1452 : real(r8), intent(out) :: qi_st_out(pcols) ! In-stratus IWC [kg/kg]
1453 :
1454 : ! Local variables
1455 :
1456 : integer i ! Column index
1457 :
1458 : real(r8) p
1459 : real(r8) T0
1460 : real(r8) qv0
1461 : real(r8) ql0
1462 : real(r8) qi0
1463 : real(r8) a_dc
1464 : real(r8) ql_dc
1465 : real(r8) qi_dc
1466 : real(r8) a_sc
1467 : real(r8) ql_sc
1468 : real(r8) qi_sc
1469 : real(r8) esat0
1470 : real(r8) qsat0
1471 : real(r8) U0
1472 : real(r8) U0_nc
1473 : real(r8) G0_nc
1474 : real(r8) al0_st_nc
1475 : real(r8) al0_st
1476 : real(r8) ai0_st_nc
1477 : real(r8) ai0_st
1478 : real(r8) a0_st
1479 : real(r8) ql0_nc
1480 : real(r8) qi0_nc
1481 : real(r8) qc0_nc
1482 : real(r8) ql0_st
1483 : real(r8) qi0_st
1484 : real(r8) qc0_st
1485 : real(r8) T
1486 : real(r8) qv
1487 : real(r8) ql
1488 : real(r8) qi
1489 : real(r8) ql_st
1490 : real(r8) qi_st
1491 : real(r8) es
1492 : real(r8) qs
1493 : real(r8) esat_in(pcols)
1494 : real(r8) qsat_in(pcols)
1495 : real(r8) U0_in(pcols)
1496 : real(r8) al0_st_nc_in(pcols)
1497 : real(r8) ai0_st_nc_in(pcols)
1498 : real(r8) G0_nc_in(pcols)
1499 : integer idxmod
1500 : real(r8) U
1501 : real(r8) U_nc
1502 : real(r8) al_st_nc
1503 : real(r8) ai_st_nc
1504 : real(r8) G_nc
1505 : real(r8) a_st
1506 : real(r8) al_st
1507 : real(r8) ai_st
1508 : real(r8) Tmin0
1509 : real(r8) Tmax0
1510 : real(r8) Tmin
1511 : real(r8) Tmax
1512 : integer caseid
1513 :
1514 : real(r8) rhmaxi
1515 : real(r8) rhmini
1516 : real(r8) rhminl
1517 : real(r8) rhminl_adj_land
1518 : real(r8) rhminh
1519 :
1520 : ! ---------------- !
1521 : ! Main Computation !
1522 : ! ---------------- !
1523 :
1524 0 : call qsat_water(T0_in(1:ncol), p_in(1:ncol), esat_in(1:ncol), qsat_in(1:ncol), ncol)
1525 0 : U0_in(:ncol) = qv0_in(:ncol)/qsat_in(:ncol)
1526 0 : if( CAMstfrac ) then
1527 : call astG_RHU(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol,&
1528 : rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:))
1529 : else
1530 : call astG_PDF(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol,&
1531 : rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:))
1532 : endif
1533 : call aist_vector(qv0_in(:),T0_in(:),p_in(:),qi0_in(:),ni0_in(:),landfrac(:),snowh(:),ai0_st_nc_in(:),ncol,&
1534 0 : rhmaxi_in(:), rhmini_in(:), rhminl_in(:), rhminl_adj_land_in(:), rhminh_in(:))
1535 :
1536 0 : do i = 1, ncol
1537 :
1538 : ! ---------------------- !
1539 : ! Define local variables !
1540 : ! ---------------------- !
1541 :
1542 0 : p = p_in(i)
1543 :
1544 0 : T0 = T0_in(i)
1545 0 : qv0 = qv0_in(i)
1546 0 : ql0 = ql0_in(i)
1547 0 : qi0 = qi0_in(i)
1548 :
1549 0 : a_dc = a_dc_in(i)
1550 0 : ql_dc = ql_dc_in(i)
1551 0 : qi_dc = qi_dc_in(i)
1552 :
1553 0 : a_sc = a_sc_in(i)
1554 0 : ql_sc = ql_sc_in(i)
1555 0 : qi_sc = qi_sc_in(i)
1556 :
1557 0 : ql_dc = 0._r8
1558 0 : qi_dc = 0._r8
1559 0 : ql_sc = 0._r8
1560 0 : qi_sc = 0._r8
1561 :
1562 0 : es = esat_in(i)
1563 0 : qs = qsat_in(i)
1564 :
1565 0 : rhmaxi = rhmaxi_in(i)
1566 0 : rhmini = rhmini_in(i)
1567 0 : rhminl = rhminl_in(i)
1568 0 : rhminl_adj_land = rhminl_adj_land_in(i)
1569 0 : rhminh = rhminh_in(i)
1570 :
1571 0 : idxmod = 0
1572 0 : caseid = -1
1573 :
1574 : ! ------------------------------------------------------------ !
1575 : ! Force the grid-mean RH to be smaller than 1 if oversaturated !
1576 : ! In order to be compatible with reduced 3x3 QQ, condensation !
1577 : ! should occur only into the liquid in gridmean_RH. !
1578 : ! ------------------------------------------------------------ !
1579 :
1580 0 : if( qv0 .gt. qs ) then
1581 : call gridmean_RH( lchnk, i, k, p, T0, qv0, ql0, qi0, &
1582 : a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, &
1583 : landfrac(i), snowh(i) )
1584 0 : call qsat_water(T0, p, esat0, qsat0)
1585 0 : U0 = (qv0/qsat0)
1586 0 : U0_nc = U0
1587 : if( CAMstfrac ) then
1588 : call astG_RHU_single(U0_nc, p, qv0, landfrac(i), snowh(i), al0_st_nc, G0_nc, &
1589 : rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
1590 : else
1591 0 : call astG_PDF_single(U0_nc, p, qv0, landfrac(i), snowh(i), al0_st_nc, G0_nc, &
1592 0 : rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
1593 : endif
1594 0 : call aist_single(qv0,T0,p,qi0,landfrac(i),snowh(i),ai0_st_nc,&
1595 0 : rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh)
1596 0 : ai0_st = (1._r8-a_dc-a_sc)*ai0_st_nc
1597 0 : al0_st = (1._r8-a_dc-a_sc)*al0_st_nc
1598 0 : a0_st = max(ai0_st,al0_st)
1599 0 : idxmod = 1
1600 : else
1601 0 : ai0_st = (1._r8-a_dc-a_sc)*ai0_st_nc_in(i)
1602 0 : al0_st = (1._r8-a_dc-a_sc)*al0_st_nc_in(i)
1603 : endif
1604 0 : a0_st = max(ai0_st,al0_st)
1605 :
1606 : ! ----------------------- !
1607 : ! Handling of input state !
1608 : ! ----------------------- !
1609 :
1610 0 : ql0_nc = max(0._r8,ql0-a_dc*ql_dc-a_sc*ql_sc)
1611 0 : qi0_nc = max(0._r8,qi0-a_dc*qi_dc-a_sc*qi_sc)
1612 0 : qc0_nc = ql0_nc + qi0_nc
1613 :
1614 0 : Tmin0 = T0 - (latvap/cpair)*ql0
1615 0 : Tmax0 = T0 + ((latvap+latice)/cpair)*qv0
1616 :
1617 : ! ------------------------------------------------------------- !
1618 : ! Do nothing and just exit if generalized in-stratus condensate !
1619 : ! condition is satisfied. This includes the case I. !
1620 : ! For 4x4 liquid stratus, a0_st --> al0_st. !
1621 : ! ------------------------------------------------------------- !
1622 0 : if( ( ql0_nc .ge. qlst_min*al0_st ) .and. ( ql0_nc .le. qlst_max*al0_st ) ) then
1623 :
1624 : ! ------------------ !
1625 : ! This is the case I !
1626 : ! ------------------ !
1627 0 : T = T0
1628 0 : qv = qv0
1629 0 : ql = ql0
1630 0 : qi = qi0
1631 0 : caseid = 0
1632 0 : goto 10
1633 : else
1634 : ! ----------------------------- !
1635 : ! This is case II : Dense Cloud !
1636 : ! ----------------------------- !
1637 0 : if( al0_st .eq. 0._r8 .and. ql0_nc .gt. 0._r8 ) then
1638 : ! ------------------------------------- !
1639 : ! Compute hypothetical full evaporation !
1640 : ! ------------------------------------- !
1641 0 : T = Tmin0
1642 0 : qv = qv0 + ql0
1643 0 : call qsat_water(T, p, es, qs)
1644 0 : U = qv/qs
1645 0 : U_nc = U
1646 : if( CAMstfrac ) then
1647 : call astG_RHU_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, &
1648 : rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
1649 : else
1650 0 : call astG_PDF_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, &
1651 0 : rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
1652 : endif
1653 0 : al_st = (1._r8-a_dc-a_sc)*al_st_nc
1654 0 : caseid = 0
1655 :
1656 0 : if( al_st .eq. 0._r8 ) then
1657 0 : ql = 0._r8
1658 0 : qi = qi0
1659 0 : idxmod = 1
1660 0 : caseid = 1
1661 0 : goto 10
1662 : else
1663 : ! ------------------------------------------- !
1664 : ! Evaporate until qc_st decreases to qlst_max !
1665 : ! ------------------------------------------- !
1666 0 : Tmin = Tmin0
1667 0 : Tmax = T0
1668 : call instratus_core( lchnk, i, k, p, &
1669 : T0, qv0, ql0, 0._r8, &
1670 : a_dc, ql_dc, qi_dc, &
1671 : a_sc, ql_sc, qi_sc, ai0_st, &
1672 0 : qlst_max, Tmin, Tmax, landfrac(i), snowh(i), &
1673 : rhminl, rhminl_adj_land, rhminh, &
1674 0 : T, qv, ql, qi )
1675 0 : idxmod = 1
1676 0 : caseid = 2
1677 0 : goto 10
1678 : endif
1679 : ! ------------------------------ !
1680 : ! This is case III : Empty Cloud !
1681 : ! ------------------------------ !
1682 0 : elseif( al0_st .gt. 0._r8 .and. ql0_nc .eq. 0._r8 ) then
1683 : ! ------------------------------------------ !
1684 : ! Condense until qc_st increases to qlst_min !
1685 : ! ------------------------------------------ !
1686 0 : Tmin = Tmin0
1687 0 : Tmax = Tmax0
1688 : call instratus_core( lchnk, i, k, p, &
1689 : T0, qv0, ql0, 0._r8, &
1690 : a_dc, ql_dc, qi_dc, &
1691 : a_sc, ql_sc, qi_sc, ai0_st, &
1692 0 : qlst_min, Tmin, Tmax, landfrac(i), snowh(i), &
1693 : rhminl, rhminl_adj_land, rhminh, &
1694 0 : T, qv, ql, qi )
1695 0 : idxmod = 1
1696 0 : caseid = 3
1697 0 : goto 10
1698 : ! --------------- !
1699 : ! This is case IV !
1700 : ! --------------- !
1701 0 : elseif( al0_st .gt. 0._r8 .and. ql0_nc .gt. 0._r8 ) then
1702 :
1703 0 : if( ql0_nc .gt. qlst_max*al0_st ) then
1704 : ! --------------------------------------- !
1705 : ! Evaporate until qc_st drops to qlst_max !
1706 : ! --------------------------------------- !
1707 0 : Tmin = Tmin0
1708 0 : Tmax = Tmax0
1709 : call instratus_core( lchnk, i, k, p, &
1710 : T0, qv0, ql0, 0._r8, &
1711 : a_dc, ql_dc, qi_dc, &
1712 : a_sc, ql_sc, qi_sc, ai0_st, &
1713 0 : qlst_max, Tmin, Tmax, landfrac(i), snowh(i), &
1714 : rhminl, rhminl_adj_land, rhminh, &
1715 0 : T, qv, ql, qi )
1716 0 : idxmod = 1
1717 0 : caseid = 4
1718 0 : goto 10
1719 0 : elseif( ql0_nc .lt. qlst_min*al0_st ) then
1720 : ! -------------------------------------------- !
1721 : ! Condensate until qc_st increases to qlst_min !
1722 : ! -------------------------------------------- !
1723 0 : Tmin = Tmin0
1724 0 : Tmax = Tmax0
1725 : call instratus_core( lchnk, i, k, p, &
1726 : T0, qv0, ql0, 0._r8, &
1727 : a_dc, ql_dc, qi_dc, &
1728 : a_sc, ql_sc, qi_sc, ai0_st, &
1729 0 : qlst_min, Tmin, Tmax, landfrac(i), snowh(i), &
1730 : rhminl, rhminl_adj_land, rhminh, &
1731 0 : T, qv, ql, qi )
1732 0 : idxmod = 1
1733 0 : caseid = 5
1734 0 : goto 10
1735 : else
1736 : ! ------------------------------------------------ !
1737 : ! This case should not happen. Issue error message !
1738 : ! ------------------------------------------------ !
1739 0 : write(iulog,*) 'Impossible case1 in instratus_condensate'
1740 0 : call endrun
1741 : endif
1742 : ! ------------------------------------------------ !
1743 : ! This case should not happen. Issue error message !
1744 : ! ------------------------------------------------ !
1745 : else
1746 0 : write(iulog,*) 'Impossible case2 in instratus_condensate'
1747 0 : write(iulog,*) al0_st, a_sc, a_dc
1748 0 : write(iulog,*) 1000*ql0_nc, 1000*(ql0+qi0)
1749 0 : call endrun
1750 : endif
1751 : endif
1752 :
1753 : 10 continue
1754 :
1755 : ! -------------------------------------------------- !
1756 : ! Force final energy-moisture conserving consistency !
1757 : ! -------------------------------------------------- !
1758 :
1759 0 : qi = qi0
1760 :
1761 0 : if( idxmod .eq. 1 ) then
1762 0 : call aist_single(qv,T,p,qi,landfrac(i),snowh(i),ai_st_nc,&
1763 0 : rhmaxi, rhmini, rhminl, rhminl_adj_land, rhminh)
1764 0 : ai_st = (1._r8-a_dc-a_sc)*ai_st_nc
1765 0 : call qsat_water(T, p, es, qs)
1766 0 : U = (qv/qs)
1767 0 : U_nc = U
1768 : if( CAMstfrac ) then
1769 : call astG_RHU_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, &
1770 : rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
1771 : else
1772 0 : call astG_PDF_single(U_nc, p, qv, landfrac(i), snowh(i), al_st_nc, G_nc, &
1773 0 : rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
1774 : endif
1775 0 : al_st = (1._r8-a_dc-a_sc)*al_st_nc
1776 : else
1777 0 : ai_st = (1._r8-a_dc-a_sc)*ai0_st_nc_in(i)
1778 0 : al_st = (1._r8-a_dc-a_sc)*al0_st_nc_in(i)
1779 : endif
1780 :
1781 0 : a_st = max(ai_st,al_st)
1782 :
1783 0 : if( al_st .eq. 0._r8 ) then
1784 : ql_st = 0._r8
1785 : else
1786 0 : ql_st = ql/al_st
1787 0 : ql_st = min(qlst_max,max(qlst_min,ql_st)) ! PJR
1788 : endif
1789 0 : if( ai_st .eq. 0._r8 ) then
1790 : qi_st = 0._r8
1791 : else
1792 0 : qi_st = qi/ai_st
1793 : endif
1794 :
1795 0 : qi = ai_st*qi_st
1796 0 : ql = al_st*ql_st
1797 :
1798 0 : T = T0 - (latvap/cpair)*(ql0-ql) - ((latvap+latice)/cpair)*(qi0-qi)
1799 0 : qv = qv0 + ql0 - ql + qi0 - qi
1800 :
1801 : ! -------------- !
1802 : ! Send to output !
1803 : ! -------------- !
1804 :
1805 0 : T_out(i) = T
1806 0 : qv_out(i) = qv
1807 0 : ql_out(i) = ql
1808 0 : qi_out(i) = qi
1809 0 : al_st_out(i) = al_st
1810 0 : ai_st_out(i) = ai_st
1811 0 : ql_st_out(i) = ql_st
1812 0 : qi_st_out(i) = qi_st
1813 :
1814 : enddo
1815 :
1816 0 : return
1817 : end subroutine instratus_condensate
1818 :
1819 : ! ----------------- !
1820 : ! End of subroutine !
1821 : ! ----------------- !
1822 :
1823 0 : subroutine instratus_core( lchnk, icol, k, p, &
1824 : T0, qv0, ql0, qi0, &
1825 : a_dc, ql_dc, qi_dc, &
1826 : a_sc, ql_sc, qi_sc, ai_st, &
1827 : qcst_crit, Tmin, Tmax, landfrac, snowh, &
1828 : rhminl, rhminl_adj_land, rhminh, &
1829 : T, qv, ql, qi )
1830 :
1831 : ! ------------------------------------------------------ !
1832 : ! Subroutine to find saturation equilibrium state using !
1833 : ! a Newton iteration method, so that 'qc_st = qcst_crit' !
1834 : ! is satisfied. !
1835 : ! ------------------------------------------------------ !
1836 :
1837 : integer, intent(in) :: lchnk ! Chunk identifier
1838 : integer, intent(in) :: icol ! Number of atmospheric columns
1839 : integer, intent(in) :: k ! Layer index
1840 :
1841 : real(r8), intent(in) :: p ! Pressure [Pa]
1842 : real(r8), intent(in) :: T0 ! Temperature [K]
1843 : real(r8), intent(in) :: qv0 ! Grid-mean water vapor [kg/kg]
1844 : real(r8), intent(in) :: ql0 ! Grid-mean LWC [kg/kg]
1845 : real(r8), intent(in) :: qi0 ! Grid-mean IWC [kg/kg]
1846 :
1847 : real(r8), intent(in) :: a_dc ! Deep cumulus cloud fraction
1848 : real(r8), intent(in) :: ql_dc ! In-deep cumulus LWC [kg/kg]
1849 : real(r8), intent(in) :: qi_dc ! In-deep cumulus IWC [kg/kg]
1850 : real(r8), intent(in) :: a_sc ! Shallow cumulus cloud fraction
1851 : real(r8), intent(in) :: ql_sc ! In-shallow cumulus LWC [kg/kg]
1852 : real(r8), intent(in) :: qi_sc ! In-shallow cumulus IWC [kg/kg]
1853 :
1854 : real(r8), intent(in) :: ai_st ! Ice stratus fraction (fixed)
1855 :
1856 : real(r8), intent(in) :: Tmin ! Minimum temperature system can have [K]
1857 : real(r8), intent(in) :: Tmax ! Maximum temperature system can have [K]
1858 : real(r8), intent(in) :: qcst_crit ! Critical in-stratus condensate [kg/kg]
1859 : real(r8), intent(in) :: landfrac ! Land fraction
1860 : real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
1861 :
1862 : real(r8), intent(in) :: rhminl
1863 : real(r8), intent(in) :: rhminl_adj_land
1864 : real(r8), intent(in) :: rhminh
1865 :
1866 : real(r8), intent(out) :: T ! Temperature [K]
1867 : real(r8), intent(out) :: qv ! Grid-mean water vapor [kg/kg]
1868 : real(r8), intent(out) :: ql ! Grid-mean LWC [kg/kg]
1869 : real(r8), intent(out) :: qi ! Grid-mean IWC [kg/kg]
1870 :
1871 : ! Local variables
1872 :
1873 : integer i ! Iteration index
1874 :
1875 : real(r8) muQ0, muQ
1876 : real(r8) ql_nc0, qi_nc0, qc_nc0, qc_nc
1877 : real(r8) fice0, fice
1878 : real(r8) ficeg0, ficeg
1879 : real(r8) esat0
1880 : real(r8) qsat0
1881 : real(r8) dqcncdt, dastdt, dUdt
1882 : real(r8) alpha, beta
1883 : real(r8) U, U_nc
1884 : real(r8) al_st_nc, G_nc
1885 : real(r8) al_st
1886 :
1887 : ! Variables for root-finding algorithm
1888 :
1889 : integer j
1890 : real(r8) x1, x2
1891 : real(r8) rtsafe
1892 : real(r8) df, dx, dxold, f, fh, fl, temp, xh, xl
1893 : real(r8), parameter :: xacc = 1.e-3_r8
1894 :
1895 : ! ---------------- !
1896 : ! Main computation !
1897 : ! ---------------- !
1898 :
1899 0 : ql_nc0 = max(0._r8,ql0-a_dc*ql_dc-a_sc*ql_sc)
1900 0 : qi_nc0 = max(0._r8,qi0-a_dc*qi_dc-a_sc*qi_sc)
1901 0 : qc_nc0 = max(0._r8,ql0+qi0-a_dc*(ql_dc+qi_dc)-a_sc*(ql_sc+qi_sc))
1902 0 : fice0 = 0._r8
1903 0 : ficeg0 = 0._r8
1904 0 : muQ0 = 1._r8
1905 :
1906 : ! ------------ !
1907 : ! Root finding !
1908 : ! ------------ !
1909 :
1910 0 : x1 = Tmin
1911 0 : x2 = Tmax
1912 : call funcd_instratus( x1, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
1913 : a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
1914 : qcst_crit, landfrac, snowh, &
1915 : rhminl, rhminl_adj_land, rhminh, &
1916 0 : fl, df, qc_nc, fice, al_st )
1917 : call funcd_instratus( x2, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
1918 : a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
1919 : qcst_crit, landfrac, snowh, &
1920 : rhminl, rhminl_adj_land, rhminh, &
1921 0 : fh, df, qc_nc, fice, al_st )
1922 0 : if((fl > 0._r8 .and. fh > 0._r8) .or. (fl < 0._r8 .and. fh < 0._r8)) then
1923 : call funcd_instratus( T0, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
1924 : a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
1925 : qcst_crit, landfrac, snowh, &
1926 : rhminl, rhminl_adj_land, rhminh, &
1927 0 : fl, df, qc_nc, fice, al_st )
1928 0 : rtsafe = T0
1929 0 : goto 10
1930 : endif
1931 0 : if( fl == 0._r8) then
1932 0 : rtsafe = x1
1933 0 : goto 10
1934 0 : elseif( fh == 0._r8) then
1935 0 : rtsafe = x2
1936 0 : goto 10
1937 0 : elseif( fl < 0._r8) then
1938 0 : xl = x1
1939 0 : xh = x2
1940 : else
1941 0 : xh = x1
1942 0 : xl = x2
1943 : end if
1944 0 : rtsafe = 0.5_r8*(x1+x2)
1945 0 : dxold = abs(x2-x1)
1946 0 : dx = dxold
1947 : call funcd_instratus( rtsafe, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
1948 : a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
1949 : qcst_crit, landfrac, snowh, &
1950 : rhminl, rhminl_adj_land, rhminh, &
1951 0 : f, df, qc_nc, fice, al_st )
1952 0 : do j = 1, 20
1953 0 : if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) > 0._r8 .or. abs(2.0_r8*f) > abs(dxold*df) ) then
1954 0 : dxold = dx
1955 0 : dx = 0.5_r8*(xh-xl)
1956 0 : rtsafe = xl + dx
1957 0 : if(xl == rtsafe) goto 10
1958 : else
1959 0 : dxold = dx
1960 0 : dx = f/df
1961 0 : temp = rtsafe
1962 0 : rtsafe = rtsafe - dx
1963 0 : if (temp == rtsafe) goto 10
1964 : end if
1965 : ! if(abs(dx) < xacc) goto 10
1966 : call funcd_instratus( rtsafe, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
1967 : a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
1968 : qcst_crit, landfrac, snowh, &
1969 : rhminl, rhminl_adj_land, rhminh, &
1970 0 : f, df, qc_nc, fice, al_st )
1971 : ! Sep.21.2010. Sungsu modified to enhance convergence and guarantee 'qlst_min < qlst < qlst_max'.
1972 0 : if( qcst_crit < 0.5_r8 * ( qlst_min + qlst_max ) ) then
1973 0 : if( ( qc_nc*(1._r8-fice) .gt. qlst_min*al_st .and. &
1974 : qc_nc*(1._r8-fice) .lt. 1.1_r8 * qlst_min*al_st ) ) goto 10
1975 : else
1976 0 : if( ( qc_nc*(1._r8-fice) .gt. 0.9_r8 * qlst_max*al_st .and. &
1977 : qc_nc*(1._r8-fice) .lt. qlst_max*al_st ) ) goto 10
1978 : endif
1979 0 : if(f < 0._r8) then
1980 0 : xl = rtsafe
1981 : else
1982 0 : xh = rtsafe
1983 : endif
1984 :
1985 : enddo
1986 :
1987 : 10 continue
1988 :
1989 : ! ------------------------------------------- !
1990 : ! Final safety check before sending to output !
1991 : ! ------------------------------------------- !
1992 :
1993 0 : qc_nc = max(0._r8,qc_nc)
1994 :
1995 0 : T = rtsafe
1996 0 : ql = qc_nc*(1._r8-fice) + a_dc*ql_dc + a_sc*ql_sc
1997 0 : qi = qc_nc*fice + a_dc*qi_dc + a_sc*qi_sc
1998 0 : qv = qv0 + ql0 + qi0 - (qc_nc + a_dc*(ql_dc+qi_dc) + a_sc*(ql_sc+qi_sc))
1999 0 : qv = max(qv,1.e-12_r8)
2000 :
2001 0 : return
2002 : end subroutine instratus_core
2003 :
2004 : ! ----------------- !
2005 : ! End of subroutine !
2006 : ! ----------------- !
2007 :
2008 0 : subroutine funcd_instratus( T, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
2009 : a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
2010 : qcst_crit, landfrac, snowh, &
2011 : rhminl, rhminl_adj_land, rhminh, &
2012 : f, fg, qc_nc, fice, al_st )
2013 :
2014 : ! --------------------------------------------------- !
2015 : ! Subroutine to find function value and gradient at T !
2016 : ! --------------------------------------------------- !
2017 :
2018 : implicit none
2019 :
2020 : real(r8), intent(in) :: T ! Iteration temperature [K]
2021 :
2022 : real(r8), intent(in) :: p ! Pressure [Pa]
2023 : real(r8), intent(in) :: T0 ! Initial temperature [K]
2024 : real(r8), intent(in) :: qv0 ! Grid-mean water vapor [kg/kg]
2025 : real(r8), intent(in) :: ql0 ! Grid-mean LWC [kg/kg]
2026 : real(r8), intent(in) :: qi0 ! Grid-mean IWC [kg/kg]
2027 : real(r8), intent(in) :: fice0 !
2028 : real(r8), intent(in) :: muQ0 !
2029 : real(r8), intent(in) :: qc_nc0 !
2030 :
2031 : real(r8), intent(in) :: a_dc ! Deep cumulus cloud fraction
2032 : real(r8), intent(in) :: ql_dc ! In-deep cumulus LWC [kg/kg]
2033 : real(r8), intent(in) :: qi_dc ! In-deep cumulus IWC [kg/kg]
2034 : real(r8), intent(in) :: a_sc ! Shallow cumulus cloud fraction
2035 : real(r8), intent(in) :: ql_sc ! In-shallow cumulus LWC [kg/kg]
2036 : real(r8), intent(in) :: qi_sc ! In-shallow cumulus IWC [kg/kg]
2037 :
2038 : real(r8), intent(in) :: ai_st ! Ice stratus fraction (fixed)
2039 :
2040 : real(r8), intent(in) :: qcst_crit ! Critical in-stratus condensate [kg/kg]
2041 : real(r8), intent(in) :: landfrac ! Land fraction
2042 : real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
2043 :
2044 : real(r8), intent(in) :: rhminl
2045 : real(r8), intent(in) :: rhminl_adj_land
2046 : real(r8), intent(in) :: rhminh
2047 :
2048 : real(r8), intent(out) :: f ! Value of minimization function at T
2049 : real(r8), intent(out) :: fg ! Gradient of minimization function
2050 : real(r8), intent(out) :: qc_nc !
2051 : real(r8), intent(out) :: al_st !
2052 : real(r8), intent(out) :: fice !
2053 :
2054 : ! Local variables
2055 :
2056 : real(r8) es
2057 : real(r8) qs
2058 : real(r8) dqsdT
2059 : real(r8) dqcncdt
2060 : real(r8) alpha
2061 : real(r8) beta
2062 : real(r8) U
2063 : real(r8) U_nc
2064 : real(r8) al_st_nc
2065 : real(r8) G_nc
2066 : real(r8) dUdt
2067 : real(r8) dalstdt
2068 : real(r8) qv
2069 :
2070 : ! ---------------- !
2071 : ! Main computation !
2072 : ! ---------------- !
2073 :
2074 0 : call qsat_water(T, p, es, qs, dqsdt=dqsdT)
2075 :
2076 0 : fice = fice0
2077 0 : qc_nc = (cpair/latvap)*(T-T0)+muQ0*qc_nc0
2078 0 : dqcncdt = (cpair/latvap)
2079 0 : qv = (qv0 + ql0 + qi0 - (qc_nc + a_dc*(ql_dc+qi_dc) + a_sc*(ql_sc+qi_sc)))
2080 0 : alpha = (1._r8/qs)
2081 0 : beta = (qv/qs**2._r8)*dqsdT
2082 :
2083 0 : U = (qv/qs)
2084 0 : U_nc = U
2085 : if( CAMstfrac ) then
2086 : call astG_RHU_single(U_nc, p, qv, landfrac, snowh, al_st_nc, G_nc, &
2087 : rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
2088 : else
2089 : call astG_PDF_single(U_nc, p, qv, landfrac, snowh, al_st_nc, G_nc, &
2090 0 : rhminl_in=rhminl, rhminl_adj_land_in=rhminl_adj_land, rhminh_in=rhminh)
2091 : endif
2092 0 : al_st = (1._r8-a_dc-a_sc)*al_st_nc
2093 0 : dUdt = -(alpha*dqcncdt+beta)
2094 0 : dalstdt = (1._r8/G_nc)*dUdt
2095 0 : if( U_nc .eq. 1._r8 ) dalstdt = 0._r8
2096 :
2097 0 : f = qc_nc - qcst_crit*al_st
2098 0 : fg = dqcncdt - qcst_crit*dalstdt
2099 :
2100 0 : return
2101 : end subroutine funcd_instratus
2102 :
2103 : ! ----------------- !
2104 : ! End of subroutine !
2105 : ! ----------------- !
2106 :
2107 0 : subroutine gridmean_RH( lchnk, icol, k, p, T, qv, ql, qi, &
2108 : a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, &
2109 : landfrac, snowh )
2110 :
2111 : ! ------------------------------------------------------------- !
2112 : ! Subroutine to force grid-mean RH = 1 when RH > 1 !
2113 : ! This is condensation process similar to instratus_condensate. !
2114 : ! During condensation, we assume 'fice' is maintained in this !
2115 : ! verison for MG not for RK. !
2116 : ! ------------------------------------------------------------- !
2117 :
2118 : integer, intent(in) :: lchnk ! Chunk identifier
2119 : integer, intent(in) :: icol ! Number of atmospheric columns
2120 : integer, intent(in) :: k ! Layer index
2121 :
2122 : real(r8), intent(in) :: p ! Pressure [Pa]
2123 : real(r8), intent(inout) :: T ! Temperature [K]
2124 : real(r8), intent(inout) :: qv ! Grid-mean water vapor [kg/kg]
2125 : real(r8), intent(inout) :: ql ! Grid-mean LWC [kg/kg]
2126 : real(r8), intent(inout) :: qi ! Grid-mean IWC [kg/kg]
2127 :
2128 : real(r8), intent(in) :: a_dc ! Deep cumulus cloud fraction
2129 : real(r8), intent(in) :: ql_dc ! In-deep cumulus LWC [kg/kg]
2130 : real(r8), intent(in) :: qi_dc ! In-deep cumulus IWC [kg/kg]
2131 : real(r8), intent(in) :: a_sc ! Shallow cumulus cloud fraction
2132 : real(r8), intent(in) :: ql_sc ! In-shallow cumulus LWC [kg/kg]
2133 : real(r8), intent(in) :: qi_sc ! In-shallow cumulus IWC [kg/kg]
2134 :
2135 : real(r8), intent(in) :: landfrac ! Land fraction
2136 : real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
2137 :
2138 : ! Local variables
2139 :
2140 : integer m ! Iteration index
2141 :
2142 : real(r8) ql_nc0, qi_nc0, qc_nc0
2143 : real(r8) Tscale
2144 : real(r8) Tc, qt, qc, dqcdt, qc_nc
2145 : real(r8) es, qs, dqsdT
2146 : real(r8) al_st_nc, G_nc
2147 : real(r8) f, fg
2148 : real(r8), parameter :: xacc = 1.e-3_r8
2149 :
2150 : ! ---------------- !
2151 : ! Main computation !
2152 : ! ---------------- !
2153 :
2154 0 : ql_nc0 = max(0._r8,ql-a_dc*ql_dc-a_sc*ql_sc)
2155 0 : qi_nc0 = max(0._r8,qi-a_dc*qi_dc-a_sc*qi_sc)
2156 0 : qc_nc0 = max(0._r8,ql+qi-a_dc*(ql_dc+qi_dc)-a_sc*(ql_sc+qi_sc))
2157 0 : Tc = T - (latvap/cpair)*ql
2158 0 : qt = qv + ql
2159 :
2160 0 : do m = 1, 20
2161 0 : call qsat_water(T, p, es, qs, dqsdt=dqsdT)
2162 0 : Tscale = latvap/cpair
2163 0 : qc = (T-Tc)/Tscale
2164 0 : dqcdt = 1._r8/Tscale
2165 0 : f = qs + qc - qt
2166 0 : fg = dqsdT + dqcdt
2167 0 : fg = sign(1._r8,fg)*max(1.e-10_r8,abs(fg))
2168 : ! Sungsu modified convergence criteria to speed up convergence and guarantee RH <= 1.
2169 0 : if( qc .ge. 0._r8 .and. ( qt - qc ) .ge. 0.999_r8*qs .and. ( qt - qc ) .le. 1._r8*qs ) then
2170 : goto 10
2171 : endif
2172 0 : T = T - f/fg
2173 : enddo
2174 : ! write(iulog,*) 'Convergence in gridmean_RH is not reached. RH = ', ( qt - qc ) / qs
2175 : 10 continue
2176 :
2177 0 : call qsat_water(T, p, es, qs)
2178 : ! Sungsu modified 'qv = qs' in consistent with the modified convergence criteria above.
2179 0 : qv = min(qt,qs) ! Modified
2180 0 : ql = qt - qv
2181 0 : T = Tc + (latvap/cpair)*ql
2182 :
2183 0 : return
2184 : end subroutine gridmean_RH
2185 :
2186 : ! ----------------- !
2187 : ! End of subroutine !
2188 : ! ----------------- !
2189 :
2190 0 : subroutine positive_moisture( ncol, dt, qvmin, qlmin, qimin, dp, &
2191 : qv, ql, qi, t, qvten, &
2192 : qlten, qiten, tten, do_cldice)
2193 :
2194 : ! ------------------------------------------------------------------------------- !
2195 : ! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, !
2196 : ! force them to be larger than minimum value by (1) condensating water vapor !
2197 : ! into liquid or ice, and (2) by transporting water vapor from the very lower !
2198 : ! layer. '2._r8' is multiplied to the minimum values for safety. !
2199 : ! Update final state variables and tendencies associated with this correction. !
2200 : ! If any condensation happens, update (s,t) too. !
2201 : ! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
2202 : ! input tendencies. !
2203 : ! Be careful the order of k : '1': top layer, 'pver' : near-surface layer !
2204 : ! ------------------------------------------------------------------------------- !
2205 :
2206 : implicit none
2207 : integer, intent(in) :: ncol
2208 : real(r8), intent(in) :: dt
2209 : real(r8), intent(in) :: dp(pcols,pver), qvmin(pcols,pver), qlmin(pcols,pver), qimin(pcols,pver)
2210 : real(r8), intent(inout) :: qv(pcols,pver), ql(pcols,pver), qi(pcols,pver), t(pcols,pver)
2211 : real(r8), intent(out) :: qvten(pcols,pver), qlten(pcols,pver), qiten(pcols,pver), tten(pcols,pver)
2212 : logical, intent(in) :: do_cldice
2213 : integer i, k
2214 : real(r8) dql, dqi, dqv, sum, aa, dum
2215 :
2216 0 : tten(:ncol,:pver) = 0._r8
2217 0 : qvten(:ncol,:pver) = 0._r8
2218 0 : qlten(:ncol,:pver) = 0._r8
2219 0 : qiten(:ncol,:pver) = 0._r8
2220 :
2221 0 : do i = 1, ncol
2222 0 : do k = top_lev, pver
2223 0 : 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
2224 : goto 10
2225 : endif
2226 : enddo
2227 0 : goto 11
2228 : 10 continue
2229 0 : do k = top_lev, pver ! From the top to the 1st (lowest) layer from the surface
2230 0 : dql = max(0._r8,1._r8*qlmin(i,k)-ql(i,k))
2231 :
2232 0 : if (do_cldice) then
2233 0 : dqi = max(0._r8,1._r8*qimin(i,k)-qi(i,k))
2234 : else
2235 : dqi = 0._r8
2236 : end if
2237 :
2238 0 : qlten(i,k) = qlten(i,k) + dql/dt
2239 0 : qiten(i,k) = qiten(i,k) + dqi/dt
2240 0 : qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
2241 0 : tten(i,k) = tten(i,k) + (latvap/cpair)*(dql/dt) + ((latvap+latice)/cpair)*(dqi/dt)
2242 0 : ql(i,k) = ql(i,k) + dql
2243 0 : qi(i,k) = qi(i,k) + dqi
2244 0 : qv(i,k) = qv(i,k) - dql - dqi
2245 0 : t(i,k) = t(i,k) + (latvap * dql + (latvap+latice) * dqi)/cpair
2246 0 : dqv = max(0._r8,1._r8*qvmin(i,k)-qv(i,k))
2247 0 : qvten(i,k) = qvten(i,k) + dqv/dt
2248 0 : qv(i,k) = qv(i,k) + dqv
2249 0 : if( k .ne. pver ) then
2250 0 : qv(i,k+1) = qv(i,k+1) - dqv*dp(i,k)/dp(i,k+1)
2251 0 : qvten(i,k+1) = qvten(i,k+1) - dqv*dp(i,k)/dp(i,k+1)/dt
2252 : endif
2253 0 : qv(i,k) = max(qv(i,k),qvmin(i,k))
2254 0 : ql(i,k) = max(ql(i,k),qlmin(i,k))
2255 0 : qi(i,k) = max(qi(i,k),qimin(i,k))
2256 : end do
2257 : ! Extra moisture used to satisfy 'qv(i,pver)=qvmin' is proportionally
2258 : ! extracted from all the layers that has 'qv > 2*qvmin'. This fully
2259 : ! preserves column moisture.
2260 0 : if( dqv .gt. 1.e-20_r8 ) then
2261 : sum = 0._r8
2262 0 : do k = top_lev, pver
2263 0 : if( qv(i,k) .gt. 2._r8*qvmin(i,k) ) sum = sum + qv(i,k)*dp(i,k)
2264 : enddo
2265 0 : aa = dqv*dp(i,pver)/max(1.e-20_r8,sum)
2266 0 : if( aa .lt. 0.5_r8 ) then
2267 0 : do k = top_lev, pver
2268 0 : if( qv(i,k) .gt. 2._r8*qvmin(i,k) ) then
2269 0 : dum = aa*qv(i,k)
2270 0 : qv(i,k) = qv(i,k) - dum
2271 0 : qvten(i,k) = qvten(i,k) - dum/dt
2272 : endif
2273 : enddo
2274 : else
2275 0 : write(iulog,*) 'Full positive_moisture is impossible in Park Macro'
2276 : endif
2277 : endif
2278 0 : 11 continue
2279 : enddo
2280 0 : return
2281 :
2282 : end subroutine positive_moisture
2283 :
2284 : ! ----------------- !
2285 : ! End of subroutine !
2286 : ! ----------------- !
2287 :
2288 0 : SUBROUTINE gaussj(a,n,np,b,m,mp)
2289 : INTEGER m,mp,n,np,NMAX
2290 : real(r8) a(np,np),b(np,mp)
2291 0 : real(r8) aa(np,np),bb(np,mp)
2292 : PARAMETER (NMAX=50)
2293 : INTEGER i,icol,irow,j,k,l,ll,ii,jj,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
2294 : real(r8) big,dum,pivinv
2295 :
2296 0 : aa(:,:) = a(:,:)
2297 0 : bb(:,:) = b(:,:)
2298 :
2299 0 : do 11 j=1,n
2300 0 : ipiv(j)=0
2301 0 : 11 continue
2302 0 : do 22 i=1,n
2303 0 : big=0._r8
2304 0 : do 13 j=1,n
2305 0 : if(ipiv(j).ne.1)then
2306 0 : do 12 k=1,n
2307 0 : if (ipiv(k).eq.0) then
2308 0 : if (abs(a(j,k)).ge.big)then
2309 0 : big=abs(a(j,k))
2310 0 : irow=j
2311 0 : icol=k
2312 : endif
2313 0 : else if (ipiv(k).gt.1) then
2314 0 : write(iulog,*) 'singular matrix in gaussj 1'
2315 0 : do ii = 1, np
2316 0 : do jj = 1, np
2317 0 : write(iulog,*) ii, jj, aa(ii,jj), bb(ii,1)
2318 : end do
2319 : end do
2320 0 : call endrun
2321 : endif
2322 0 : 12 continue
2323 : endif
2324 0 : 13 continue
2325 0 : ipiv(icol)=ipiv(icol)+1
2326 0 : if (irow.ne.icol) then
2327 0 : do 14 l=1,n
2328 0 : dum=a(irow,l)
2329 0 : a(irow,l)=a(icol,l)
2330 0 : a(icol,l)=dum
2331 0 : 14 continue
2332 0 : do 15 l=1,m
2333 0 : dum=b(irow,l)
2334 0 : b(irow,l)=b(icol,l)
2335 0 : b(icol,l)=dum
2336 0 : 15 continue
2337 : endif
2338 0 : indxr(i)=irow
2339 0 : indxc(i)=icol
2340 0 : if (a(icol,icol).eq.0._r8) then
2341 0 : write(iulog,*) 'singular matrix in gaussj 2'
2342 0 : do ii = 1, np
2343 0 : do jj = 1, np
2344 0 : write(iulog,*) ii, jj, aa(ii,jj), bb(ii,1)
2345 : end do
2346 : end do
2347 0 : call endrun
2348 : endif
2349 0 : pivinv=1._r8/a(icol,icol)
2350 0 : a(icol,icol)=1._r8
2351 0 : do 16 l=1,n
2352 0 : a(icol,l)=a(icol,l)*pivinv
2353 0 : 16 continue
2354 0 : do 17 l=1,m
2355 0 : b(icol,l)=b(icol,l)*pivinv
2356 0 : 17 continue
2357 0 : do 21 ll=1,n
2358 0 : if(ll.ne.icol)then
2359 0 : dum=a(ll,icol)
2360 0 : a(ll,icol)=0._r8
2361 0 : do 18 l=1,n
2362 0 : a(ll,l)=a(ll,l)-a(icol,l)*dum
2363 0 : 18 continue
2364 0 : do 19 l=1,m
2365 0 : b(ll,l)=b(ll,l)-b(icol,l)*dum
2366 0 : 19 continue
2367 : endif
2368 0 : 21 continue
2369 0 : 22 continue
2370 0 : do 24 l=n,1,-1
2371 0 : if(indxr(l).ne.indxc(l))then
2372 0 : do 23 k=1,n
2373 0 : dum=a(k,indxr(l))
2374 0 : a(k,indxr(l))=a(k,indxc(l))
2375 0 : a(k,indxc(l))=dum
2376 0 : 23 continue
2377 : endif
2378 0 : 24 continue
2379 :
2380 0 : return
2381 : end subroutine gaussj
2382 :
2383 : ! ----------------- !
2384 : ! End of subroutine !
2385 : ! ----------------- !
2386 :
2387 : end module cldwat2m_macro
|
