1 |
|
|
|
2 |
|
|
! $Header$ |
3 |
|
|
|
4 |
|
|
SUBROUTINE conccm(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, rain, snow, & |
5 |
|
|
kbascm, ktopcm) |
6 |
|
|
|
7 |
|
|
USE dimphy |
8 |
|
|
IMPLICIT NONE |
9 |
|
|
! ====================================================================== |
10 |
|
|
! Auteur(s): Z.X. Li (LMD/CNRS) date: le 14 mars 1996 |
11 |
|
|
! Objet: Schema simple (avec flux de masse) pour la convection |
12 |
|
|
! (schema standard du modele NCAR CCM2) |
13 |
|
|
! ====================================================================== |
14 |
|
|
include "YOMCST.h" |
15 |
|
|
include "YOETHF.h" |
16 |
|
|
|
17 |
|
|
! Entree: |
18 |
|
|
REAL dtime ! pas d'integration |
19 |
|
|
REAL paprs(klon, klev+1) ! pression inter-couche (Pa) |
20 |
|
|
REAL pplay(klon, klev) ! pression au milieu de couche (Pa) |
21 |
|
|
REAL t(klon, klev) ! temperature (K) |
22 |
|
|
REAL q(klon, klev) ! humidite specifique (g/g) |
23 |
|
|
REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s) |
24 |
|
|
! Sortie: |
25 |
|
|
REAL d_t(klon, klev) ! incrementation temperature |
26 |
|
|
REAL d_q(klon, klev) ! incrementation vapeur |
27 |
|
|
REAL rain(klon) ! pluie (mm/s) |
28 |
|
|
REAL snow(klon) ! neige (mm/s) |
29 |
|
|
INTEGER kbascm(klon) ! niveau du bas de convection |
30 |
|
|
INTEGER ktopcm(klon) ! niveau du haut de convection |
31 |
|
|
|
32 |
|
|
REAL pt(klon, klev) |
33 |
|
|
REAL pq(klon, klev) |
34 |
|
|
REAL pres(klon, klev) |
35 |
|
|
REAL dp(klon, klev) |
36 |
|
|
REAL zgeom(klon, klev) |
37 |
|
|
REAL cmfprs(klon) |
38 |
|
|
REAL cmfprt(klon) |
39 |
|
|
INTEGER ntop(klon) |
40 |
|
|
INTEGER nbas(klon) |
41 |
|
|
INTEGER i, k |
42 |
|
|
REAL zlvdcp, zlsdcp, zdelta, zz, za, zb |
43 |
|
|
|
44 |
|
|
LOGICAL usekuo ! utiliser convection profonde (schema Kuo) |
45 |
|
|
PARAMETER (usekuo=.TRUE.) |
46 |
|
|
|
47 |
|
|
REAL d_t_bis(klon, klev) |
48 |
|
|
REAL d_q_bis(klon, klev) |
49 |
|
|
REAL rain_bis(klon) |
50 |
|
|
REAL snow_bis(klon) |
51 |
|
|
INTEGER ibas_bis(klon) |
52 |
|
|
INTEGER itop_bis(klon) |
53 |
|
|
REAL d_ql_bis(klon, klev) |
54 |
|
|
REAL rneb_bis(klon, klev) |
55 |
|
|
|
56 |
|
|
! initialiser les variables de sortie (pour securite) |
57 |
|
|
DO i = 1, klon |
58 |
|
|
rain(i) = 0.0 |
59 |
|
|
snow(i) = 0.0 |
60 |
|
|
kbascm(i) = 0 |
61 |
|
|
ktopcm(i) = 0 |
62 |
|
|
END DO |
63 |
|
|
DO k = 1, klev |
64 |
|
|
DO i = 1, klon |
65 |
|
|
d_t(i, k) = 0.0 |
66 |
|
|
d_q(i, k) = 0.0 |
67 |
|
|
END DO |
68 |
|
|
END DO |
69 |
|
|
|
70 |
|
|
! preparer les variables d'entree (attention: l'ordre des niveaux |
71 |
|
|
! verticaux augmente du haut vers le bas) |
72 |
|
|
DO k = 1, klev |
73 |
|
|
DO i = 1, klon |
74 |
|
|
pt(i, k) = t(i, klev-k+1) |
75 |
|
|
pq(i, k) = q(i, klev-k+1) |
76 |
|
|
pres(i, k) = pplay(i, klev-k+1) |
77 |
|
|
dp(i, k) = paprs(i, klev+1-k) - paprs(i, klev+1-k+1) |
78 |
|
|
END DO |
79 |
|
|
END DO |
80 |
|
|
DO i = 1, klon |
81 |
|
|
zgeom(i, klev) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i, & |
82 |
|
|
1)))*(paprs(i,1)-pplay(i,1)) |
83 |
|
|
END DO |
84 |
|
|
DO k = 2, klev |
85 |
|
|
DO i = 1, klon |
86 |
|
|
zgeom(i, klev+1-k) = zgeom(i, klev+1-k+1) + rd*0.5*(t(i,k-1)+t(i,k))/ & |
87 |
|
|
paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) |
88 |
|
|
END DO |
89 |
|
|
END DO |
90 |
|
|
|
91 |
|
|
CALL cmfmca(dtime, pres, dp, zgeom, pt, pq, cmfprt, cmfprs, ntop, nbas) |
92 |
|
|
|
93 |
|
|
DO k = 1, klev |
94 |
|
|
DO i = 1, klon |
95 |
|
|
d_q(i, klev+1-k) = pq(i, k) - q(i, klev+1-k) |
96 |
|
|
d_t(i, klev+1-k) = pt(i, k) - t(i, klev+1-k) |
97 |
|
|
END DO |
98 |
|
|
END DO |
99 |
|
|
|
100 |
|
|
DO i = 1, klon |
101 |
|
|
rain(i) = cmfprt(i)*rhoh2o |
102 |
|
|
snow(i) = cmfprs(i)*rhoh2o |
103 |
|
|
kbascm(i) = klev + 1 - nbas(i) |
104 |
|
|
ktopcm(i) = klev + 1 - ntop(i) |
105 |
|
|
END DO |
106 |
|
|
|
107 |
|
|
IF (usekuo) THEN |
108 |
|
|
CALL conkuo(dtime, paprs, pplay, t, q, conv_q, d_t_bis, d_q_bis, & |
109 |
|
|
d_ql_bis, rneb_bis, rain_bis, snow_bis, ibas_bis, itop_bis) |
110 |
|
|
DO k = 1, klev |
111 |
|
|
DO i = 1, klon |
112 |
|
|
d_t(i, k) = d_t(i, k) + d_t_bis(i, k) |
113 |
|
|
d_q(i, k) = d_q(i, k) + d_q_bis(i, k) |
114 |
|
|
END DO |
115 |
|
|
END DO |
116 |
|
|
DO i = 1, klon |
117 |
|
|
rain(i) = rain(i) + rain_bis(i) |
118 |
|
|
snow(i) = snow(i) + snow_bis(i) |
119 |
|
|
kbascm(i) = min(kbascm(i), ibas_bis(i)) |
120 |
|
|
ktopcm(i) = max(ktopcm(i), itop_bis(i)) |
121 |
|
|
END DO |
122 |
|
|
DO k = 1, klev ! eau liquide convective est |
123 |
|
|
DO i = 1, klon ! dispersee dans l'air |
124 |
|
|
zlvdcp = rlvtt/rcpd/(1.0+rvtmp2*q(i,k)) |
125 |
|
|
zlsdcp = rlstt/rcpd/(1.0+rvtmp2*q(i,k)) |
126 |
|
|
zdelta = max(0., sign(1.,rtt-t(i,k))) |
127 |
|
|
zz = d_ql_bis(i, k) ! re-evap. de l'eau liquide |
128 |
|
|
zb = max(0.0, zz) |
129 |
|
|
za = -max(0.0, zz)*(zlvdcp*(1.-zdelta)+zlsdcp*zdelta) |
130 |
|
|
d_t(i, k) = d_t(i, k) + za |
131 |
|
|
d_q(i, k) = d_q(i, k) + zb |
132 |
|
|
END DO |
133 |
|
|
END DO |
134 |
|
|
END IF |
135 |
|
|
|
136 |
|
|
RETURN |
137 |
|
|
END SUBROUTINE conccm |
138 |
|
|
SUBROUTINE cmfmca(deltat, p, dp, gz, tb, shb, cmfprt, cmfprs, cnt, cnb) |
139 |
|
|
USE dimphy |
140 |
|
|
IMPLICIT NONE |
141 |
|
|
! ----------------------------------------------------------------------- |
142 |
|
|
! Moist convective mass flux procedure: |
143 |
|
|
! If stratification is unstable to nonentraining parcel ascent, |
144 |
|
|
! complete an adjustment making use of a simple cloud model |
145 |
|
|
|
146 |
|
|
! Code generalized to allow specification of parcel ("updraft") |
147 |
|
|
! properties, as well as convective transport of an arbitrary |
148 |
|
|
! number of passive constituents (see cmrb array). |
149 |
|
|
! ----------------------------Code History------------------------------- |
150 |
|
|
! Original version: J. J. Hack, March 22, 1990 |
151 |
|
|
! Standardized: J. Rosinski, June 1992 |
152 |
|
|
! Reviewed: J. Hack, G. Taylor, August 1992 |
153 |
|
|
! Adaptation au LMD: Z.X. Li, mars 1996 (reference: Hack 1994, |
154 |
|
|
! J. Geophys. Res. vol 99, D3, 5551-5568). J'ai |
155 |
|
|
! introduit les constantes et les fonctions thermo- |
156 |
|
|
! dynamiques du Centre Europeen. J'ai elimine le |
157 |
|
|
! re-indicage du code en esperant que cela pourra |
158 |
|
|
! simplifier la lecture et la comprehension. |
159 |
|
|
! ----------------------------------------------------------------------- |
160 |
|
|
INTEGER pcnst ! nombre de traceurs passifs |
161 |
|
|
PARAMETER (pcnst=1) |
162 |
|
|
! ------------------------------Arguments-------------------------------- |
163 |
|
|
! Input arguments |
164 |
|
|
|
165 |
|
|
REAL deltat ! time step (seconds) |
166 |
|
|
REAL p(klon, klev) ! pressure |
167 |
|
|
REAL dp(klon, klev) ! delta-p |
168 |
|
|
REAL gz(klon, klev) ! geopotential (a partir du sol) |
169 |
|
|
|
170 |
|
|
REAL thtap(klon) ! PBL perturbation theta |
171 |
|
|
REAL shp(klon) ! PBL perturbation specific humidity |
172 |
|
|
REAL pblh(klon) ! PBL height (provided by PBL routine) |
173 |
|
|
REAL cmrp(klon, pcnst) ! constituent perturbations in PBL |
174 |
|
|
|
175 |
|
|
! Updated arguments: |
176 |
|
|
|
177 |
|
|
REAL tb(klon, klev) ! temperature (t bar) |
178 |
|
|
REAL shb(klon, klev) ! specific humidity (sh bar) |
179 |
|
|
REAL cmrb(klon, klev, pcnst) ! constituent mixing ratios (cmr bar) |
180 |
|
|
|
181 |
|
|
! Output arguments |
182 |
|
|
|
183 |
|
|
REAL cmfdt(klon, klev) ! dT/dt due to moist convection |
184 |
|
|
REAL cmfdq(klon, klev) ! dq/dt due to moist convection |
185 |
|
|
REAL cmfmc(klon, klev) ! moist convection cloud mass flux |
186 |
|
|
REAL cmfdqr(klon, klev) ! dq/dt due to convective rainout |
187 |
|
|
REAL cmfsl(klon, klev) ! convective lw static energy flux |
188 |
|
|
REAL cmflq(klon, klev) ! convective total water flux |
189 |
|
|
REAL cmfprt(klon) ! convective precipitation rate |
190 |
|
|
REAL cmfprs(klon) ! convective snowfall rate |
191 |
|
|
REAL qc(klon, klev) ! dq/dt due to rainout terms |
192 |
|
|
INTEGER cnt(klon) ! top level of convective activity |
193 |
|
|
INTEGER cnb(klon) ! bottom level of convective activity |
194 |
|
|
! ------------------------------Parameters------------------------------- |
195 |
|
|
REAL c0 ! rain water autoconversion coefficient |
196 |
|
|
PARAMETER (c0=1.0E-4) |
197 |
|
|
REAL dzmin ! minimum convective depth for precipitation |
198 |
|
|
PARAMETER (dzmin=0.0) |
199 |
|
|
REAL betamn ! minimum overshoot parameter |
200 |
|
|
PARAMETER (betamn=0.10) |
201 |
|
|
REAL cmftau ! characteristic adjustment time scale |
202 |
|
|
PARAMETER (cmftau=3600.) |
203 |
|
|
INTEGER limcnv ! top interface level limit for convection |
204 |
|
|
PARAMETER (limcnv=1) |
205 |
|
|
REAL tpmax ! maximum acceptable t perturbation (degrees C) |
206 |
|
|
PARAMETER (tpmax=1.50) |
207 |
|
|
REAL shpmax ! maximum acceptable q perturbation (g/g) |
208 |
|
|
PARAMETER (shpmax=1.50E-3) |
209 |
|
|
REAL tiny ! arbitrary small num used in transport estimates |
210 |
|
|
PARAMETER (tiny=1.0E-36) |
211 |
|
|
REAL eps ! convergence criteria (machine dependent) |
212 |
|
|
PARAMETER (eps=1.0E-13) |
213 |
|
|
REAL tmelt ! freezing point of water(req'd for rain vs snow) |
214 |
|
|
PARAMETER (tmelt=273.15) |
215 |
|
|
REAL ssfac ! supersaturation bound (detrained air) |
216 |
|
|
PARAMETER (ssfac=1.001) |
217 |
|
|
|
218 |
|
|
! ---------------------------Local workspace----------------------------- |
219 |
|
|
REAL gam(klon, klev) ! L/cp (d(qsat)/dT) |
220 |
|
|
REAL sb(klon, klev) ! dry static energy (s bar) |
221 |
|
|
REAL hb(klon, klev) ! moist static energy (h bar) |
222 |
|
|
REAL shbs(klon, klev) ! sat. specific humidity (sh bar star) |
223 |
|
|
REAL hbs(klon, klev) ! sat. moist static energy (h bar star) |
224 |
|
|
REAL shbh(klon, klev+1) ! specific humidity on interfaces |
225 |
|
|
REAL sbh(klon, klev+1) ! s bar on interfaces |
226 |
|
|
REAL hbh(klon, klev+1) ! h bar on interfaces |
227 |
|
|
REAL cmrh(klon, klev+1) ! interface constituent mixing ratio |
228 |
|
|
REAL prec(klon) ! instantaneous total precipitation |
229 |
|
|
REAL dzcld(klon) ! depth of convective layer (m) |
230 |
|
|
REAL beta(klon) ! overshoot parameter (fraction) |
231 |
|
|
REAL betamx ! local maximum on overshoot |
232 |
|
|
REAL eta(klon) ! convective mass flux (kg/m^2 s) |
233 |
|
|
REAL etagdt ! eta*grav*deltat |
234 |
|
|
REAL cldwtr(klon) ! cloud water (mass) |
235 |
|
|
REAL rnwtr(klon) ! rain water (mass) |
236 |
|
|
REAL sc(klon) ! dry static energy ("in-cloud") |
237 |
|
|
REAL shc(klon) ! specific humidity ("in-cloud") |
238 |
|
|
REAL hc(klon) ! moist static energy ("in-cloud") |
239 |
|
|
REAL cmrc(klon) ! constituent mix rat ("in-cloud") |
240 |
|
|
REAL dq1(klon) ! shb convective change (lower lvl) |
241 |
|
|
REAL dq2(klon) ! shb convective change (mid level) |
242 |
|
|
REAL dq3(klon) ! shb convective change (upper lvl) |
243 |
|
|
REAL ds1(klon) ! sb convective change (lower lvl) |
244 |
|
|
REAL ds2(klon) ! sb convective change (mid level) |
245 |
|
|
REAL ds3(klon) ! sb convective change (upper lvl) |
246 |
|
|
REAL dcmr1(klon) ! cmrb convective change (lower lvl) |
247 |
|
|
REAL dcmr2(klon) ! cmrb convective change (mid level) |
248 |
|
|
REAL dcmr3(klon) ! cmrb convective change (upper lvl) |
249 |
|
|
REAL flotab(klon) ! hc - hbs (mesure d'instabilite) |
250 |
|
|
LOGICAL ldcum(klon) ! .true. si la convection existe |
251 |
|
|
LOGICAL etagt0 ! true if eta > 0.0 |
252 |
|
|
REAL dt ! current 2 delta-t (model time step) |
253 |
|
|
REAL cats ! modified characteristic adj. time |
254 |
|
|
REAL rdt ! 1./dt |
255 |
|
|
REAL qprime ! modified specific humidity pert. |
256 |
|
|
REAL tprime ! modified thermal perturbation |
257 |
|
|
REAL pblhgt ! bounded pbl height (max[pblh,1m]) |
258 |
|
|
REAL fac1 ! intermediate scratch variable |
259 |
|
|
REAL shprme ! intermediate specific humidity pert. |
260 |
|
|
REAL qsattp ! saturation mixing ratio for |
261 |
|
|
! ! thermally perturbed PBL parcels |
262 |
|
|
REAL dz ! local layer depth |
263 |
|
|
REAL b1 ! bouyancy measure in detrainment lvl |
264 |
|
|
REAL b2 ! bouyancy measure in condensation lvl |
265 |
|
|
REAL g ! bounded vertical gradient of hb |
266 |
|
|
REAL tmass ! total mass available for convective exchange |
267 |
|
|
REAL denom ! intermediate scratch variable |
268 |
|
|
REAL qtest1 ! used in negative q test (middle lvl) |
269 |
|
|
REAL qtest2 ! used in negative q test (lower lvl) |
270 |
|
|
REAL fslkp ! flux lw static energy (bot interface) |
271 |
|
|
REAL fslkm ! flux lw static energy (top interface) |
272 |
|
|
REAL fqlkp ! flux total water (bottom interface) |
273 |
|
|
REAL fqlkm ! flux total water (top interface) |
274 |
|
|
REAL botflx ! bottom constituent mixing ratio flux |
275 |
|
|
REAL topflx ! top constituent mixing ratio flux |
276 |
|
|
REAL efac1 ! ratio cmrb to convectively induced change (bot lvl) |
277 |
|
|
REAL efac2 ! ratio cmrb to convectively induced change (mid lvl) |
278 |
|
|
REAL efac3 ! ratio cmrb to convectively induced change (top lvl) |
279 |
|
|
|
280 |
|
|
INTEGER i, k ! indices horizontal et vertical |
281 |
|
|
INTEGER km1 ! k-1 (index offset) |
282 |
|
|
INTEGER kp1 ! k+1 (index offset) |
283 |
|
|
INTEGER m ! constituent index |
284 |
|
|
INTEGER ktp ! temporary index used to track top |
285 |
|
|
INTEGER is ! nombre de points a ajuster |
286 |
|
|
|
287 |
|
|
REAL tmp1, tmp2, tmp3, tmp4 |
288 |
|
|
REAL zx_t, zx_p, zx_q, zx_qs, zx_gam |
289 |
|
|
REAL zcor, zdelta, zcvm5 |
290 |
|
|
|
291 |
|
|
REAL qhalf, sh1, sh2, shbs1, shbs2 |
292 |
|
|
include "YOMCST.h" |
293 |
|
|
include "YOETHF.h" |
294 |
|
|
include "FCTTRE.h" |
295 |
|
|
qhalf(sh1, sh2, shbs1, shbs2) = min(max(sh1,sh2), & |
296 |
|
|
(shbs2*sh1+shbs1*sh2)/(shbs1+shbs2)) |
297 |
|
|
|
298 |
|
|
! ----------------------------------------------------------------------- |
299 |
|
|
! pas de traceur pour l'instant |
300 |
|
|
DO m = 1, pcnst |
301 |
|
|
DO k = 1, klev |
302 |
|
|
DO i = 1, klon |
303 |
|
|
cmrb(i, k, m) = 0.0 |
304 |
|
|
END DO |
305 |
|
|
END DO |
306 |
|
|
END DO |
307 |
|
|
|
308 |
|
|
! Les perturbations de la couche limite sont zero pour l'instant |
309 |
|
|
|
310 |
|
|
DO m = 1, pcnst |
311 |
|
|
DO i = 1, klon |
312 |
|
|
cmrp(i, m) = 0.0 |
313 |
|
|
END DO |
314 |
|
|
END DO |
315 |
|
|
DO i = 1, klon |
316 |
|
|
thtap(i) = 0.0 |
317 |
|
|
shp(i) = 0.0 |
318 |
|
|
pblh(i) = 1.0 |
319 |
|
|
END DO |
320 |
|
|
|
321 |
|
|
! Ensure that characteristic adjustment time scale (cmftau) assumed |
322 |
|
|
! in estimate of eta isn't smaller than model time scale (deltat) |
323 |
|
|
|
324 |
|
|
dt = deltat |
325 |
|
|
cats = max(dt, cmftau) |
326 |
|
|
rdt = 1.0/dt |
327 |
|
|
|
328 |
|
|
! Compute sb,hb,shbs,hbs |
329 |
|
|
|
330 |
|
|
DO k = 1, klev |
331 |
|
|
DO i = 1, klon |
332 |
|
|
zx_t = tb(i, k) |
333 |
|
|
zx_p = p(i, k) |
334 |
|
|
zx_q = shb(i, k) |
335 |
|
|
zdelta = max(0., sign(1.,rtt-zx_t)) |
336 |
|
|
zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta |
337 |
|
|
zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zx_q) |
338 |
|
|
zx_qs = r2es*foeew(zx_t, zdelta)/zx_p |
339 |
|
|
zx_qs = min(0.5, zx_qs) |
340 |
|
|
zcor = 1./(1.-retv*zx_qs) |
341 |
|
|
zx_qs = zx_qs*zcor |
342 |
|
|
zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) |
343 |
|
|
shbs(i, k) = zx_qs |
344 |
|
|
gam(i, k) = zx_gam |
345 |
|
|
END DO |
346 |
|
|
END DO |
347 |
|
|
|
348 |
|
|
DO k = limcnv, klev |
349 |
|
|
DO i = 1, klon |
350 |
|
|
sb(i, k) = rcpd*tb(i, k) + gz(i, k) |
351 |
|
|
hb(i, k) = sb(i, k) + rlvtt*shb(i, k) |
352 |
|
|
hbs(i, k) = sb(i, k) + rlvtt*shbs(i, k) |
353 |
|
|
END DO |
354 |
|
|
END DO |
355 |
|
|
|
356 |
|
|
! Compute sbh, shbh |
357 |
|
|
|
358 |
|
|
DO k = limcnv + 1, klev |
359 |
|
|
km1 = k - 1 |
360 |
|
|
DO i = 1, klon |
361 |
|
|
sbh(i, k) = 0.5*(sb(i,km1)+sb(i,k)) |
362 |
|
|
shbh(i, k) = qhalf(shb(i,km1), shb(i,k), shbs(i,km1), shbs(i,k)) |
363 |
|
|
hbh(i, k) = sbh(i, k) + rlvtt*shbh(i, k) |
364 |
|
|
END DO |
365 |
|
|
END DO |
366 |
|
|
|
367 |
|
|
! Specify properties at top of model (not used, but filling anyway) |
368 |
|
|
|
369 |
|
|
DO i = 1, klon |
370 |
|
|
sbh(i, limcnv) = sb(i, limcnv) |
371 |
|
|
shbh(i, limcnv) = shb(i, limcnv) |
372 |
|
|
hbh(i, limcnv) = hb(i, limcnv) |
373 |
|
|
END DO |
374 |
|
|
|
375 |
|
|
! Zero vertically independent control, tendency & diagnostic arrays |
376 |
|
|
|
377 |
|
|
DO i = 1, klon |
378 |
|
|
prec(i) = 0.0 |
379 |
|
|
dzcld(i) = 0.0 |
380 |
|
|
cnb(i) = 0 |
381 |
|
|
cnt(i) = klev + 1 |
382 |
|
|
END DO |
383 |
|
|
|
384 |
|
|
DO k = 1, klev |
385 |
|
|
DO i = 1, klon |
386 |
|
|
cmfdt(i, k) = 0. |
387 |
|
|
cmfdq(i, k) = 0. |
388 |
|
|
cmfdqr(i, k) = 0. |
389 |
|
|
cmfmc(i, k) = 0. |
390 |
|
|
cmfsl(i, k) = 0. |
391 |
|
|
cmflq(i, k) = 0. |
392 |
|
|
END DO |
393 |
|
|
END DO |
394 |
|
|
|
395 |
|
|
! Begin moist convective mass flux adjustment procedure. |
396 |
|
|
! Formalism ensures that negative cloud liquid water can never occur |
397 |
|
|
|
398 |
|
|
DO k = klev - 1, limcnv + 1, -1 |
399 |
|
|
km1 = k - 1 |
400 |
|
|
kp1 = k + 1 |
401 |
|
|
DO i = 1, klon |
402 |
|
|
eta(i) = 0.0 |
403 |
|
|
beta(i) = 0.0 |
404 |
|
|
ds1(i) = 0.0 |
405 |
|
|
ds2(i) = 0.0 |
406 |
|
|
ds3(i) = 0.0 |
407 |
|
|
dq1(i) = 0.0 |
408 |
|
|
dq2(i) = 0.0 |
409 |
|
|
dq3(i) = 0.0 |
410 |
|
|
|
411 |
|
|
! Specification of "cloud base" conditions |
412 |
|
|
|
413 |
|
|
qprime = 0.0 |
414 |
|
|
tprime = 0.0 |
415 |
|
|
|
416 |
|
|
! Assign tprime within the PBL to be proportional to the quantity |
417 |
|
|
! thtap (which will be bounded by tpmax), passed to this routine by |
418 |
|
|
! the PBL routine. Don't allow perturbation to produce a dry |
419 |
|
|
! adiabatically unstable parcel. Assign qprime within the PBL to be |
420 |
|
|
! an appropriately modified value of the quantity shp (which will be |
421 |
|
|
! bounded by shpmax) passed to this routine by the PBL routine. The |
422 |
|
|
! quantity qprime should be less than the local saturation value |
423 |
|
|
! (qsattp=qsat[t+tprime,p]). In both cases, thtap and shp are |
424 |
|
|
! linearly reduced toward zero as the PBL top is approached. |
425 |
|
|
|
426 |
|
|
pblhgt = max(pblh(i), 1.0) |
427 |
|
|
IF (gz(i,kp1)/rg<=pblhgt .AND. dzcld(i)==0.0) THEN |
428 |
|
|
fac1 = max(0.0, 1.0-gz(i,kp1)/rg/pblhgt) |
429 |
|
|
tprime = min(thtap(i), tpmax)*fac1 |
430 |
|
|
qsattp = shbs(i, kp1) + rcpd/rlvtt*gam(i, kp1)*tprime |
431 |
|
|
shprme = min(min(shp(i),shpmax)*fac1, max(qsattp-shb(i,kp1),0.0)) |
432 |
|
|
qprime = max(qprime, shprme) |
433 |
|
|
ELSE |
434 |
|
|
tprime = 0.0 |
435 |
|
|
qprime = 0.0 |
436 |
|
|
END IF |
437 |
|
|
|
438 |
|
|
! Specify "updraft" (in-cloud) thermodynamic properties |
439 |
|
|
|
440 |
|
|
sc(i) = sb(i, kp1) + rcpd*tprime |
441 |
|
|
shc(i) = shb(i, kp1) + qprime |
442 |
|
|
hc(i) = sc(i) + rlvtt*shc(i) |
443 |
|
|
flotab(i) = hc(i) - hbs(i, k) |
444 |
|
|
dz = dp(i, k)*rd*tb(i, k)/rg/p(i, k) |
445 |
|
|
IF (flotab(i)>0.0) THEN |
446 |
|
|
dzcld(i) = dzcld(i) + dz |
447 |
|
|
ELSE |
448 |
|
|
dzcld(i) = 0.0 |
449 |
|
|
END IF |
450 |
|
|
END DO |
451 |
|
|
|
452 |
|
|
! Check on moist convective instability |
453 |
|
|
|
454 |
|
|
is = 0 |
455 |
|
|
DO i = 1, klon |
456 |
|
|
IF (flotab(i)>0.0) THEN |
457 |
|
|
ldcum(i) = .TRUE. |
458 |
|
|
is = is + 1 |
459 |
|
|
ELSE |
460 |
|
|
ldcum(i) = .FALSE. |
461 |
|
|
END IF |
462 |
|
|
END DO |
463 |
|
|
|
464 |
|
|
IF (is==0) THEN |
465 |
|
|
DO i = 1, klon |
466 |
|
|
dzcld(i) = 0.0 |
467 |
|
|
END DO |
468 |
|
|
GO TO 70 |
469 |
|
|
END IF |
470 |
|
|
|
471 |
|
|
! Current level just below top level => no overshoot |
472 |
|
|
|
473 |
|
|
IF (k<=limcnv+1) THEN |
474 |
|
|
DO i = 1, klon |
475 |
|
|
IF (ldcum(i)) THEN |
476 |
|
|
cldwtr(i) = sb(i, k) - sc(i) + flotab(i)/(1.0+gam(i,k)) |
477 |
|
|
cldwtr(i) = max(0.0, cldwtr(i)) |
478 |
|
|
beta(i) = 0.0 |
479 |
|
|
END IF |
480 |
|
|
END DO |
481 |
|
|
GO TO 20 |
482 |
|
|
END IF |
483 |
|
|
|
484 |
|
|
! First guess at overshoot parameter using crude buoyancy closure |
485 |
|
|
! 10% overshoot assumed as a minimum and 1-c0*dz maximum to start |
486 |
|
|
! If pre-existing supersaturation in detrainment layer, beta=0 |
487 |
|
|
! cldwtr is temporarily equal to RLVTT*l (l=> liquid water) |
488 |
|
|
|
489 |
|
|
DO i = 1, klon |
490 |
|
|
IF (ldcum(i)) THEN |
491 |
|
|
cldwtr(i) = sb(i, k) - sc(i) + flotab(i)/(1.0+gam(i,k)) |
492 |
|
|
cldwtr(i) = max(0.0, cldwtr(i)) |
493 |
|
|
betamx = 1.0 - c0*max(0.0, (dzcld(i)-dzmin)) |
494 |
|
|
b1 = (hc(i)-hbs(i,km1))*dp(i, km1) |
495 |
|
|
b2 = (hc(i)-hbs(i,k))*dp(i, k) |
496 |
|
|
beta(i) = max(betamn, min(betamx,1.0+b1/b2)) |
497 |
|
|
IF (hbs(i,km1)<=hb(i,km1)) beta(i) = 0.0 |
498 |
|
|
END IF |
499 |
|
|
END DO |
500 |
|
|
|
501 |
|
|
! Bound maximum beta to ensure physically realistic solutions |
502 |
|
|
|
503 |
|
|
! First check constrains beta so that eta remains positive |
504 |
|
|
! (assuming that eta is already positive for beta equal zero) |
505 |
|
|
! La premiere contrainte de beta est que le flux eta doit etre positif. |
506 |
|
|
|
507 |
|
|
DO i = 1, klon |
508 |
|
|
IF (ldcum(i)) THEN |
509 |
|
|
tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) - & |
510 |
|
|
(hbh(i,kp1)-hc(i))*dp(i, k)/dp(i, kp1) |
511 |
|
|
tmp2 = (1.0+gam(i,k))*(sc(i)-sbh(i,k)) |
512 |
|
|
IF ((beta(i)*tmp2-tmp1)>0.0) THEN |
513 |
|
|
betamx = 0.99*(tmp1/tmp2) |
514 |
|
|
beta(i) = max(0.0, min(betamx,beta(i))) |
515 |
|
|
END IF |
516 |
|
|
|
517 |
|
|
! Second check involves supersaturation of "detrainment layer" |
518 |
|
|
! small amount of supersaturation acceptable (by ssfac factor) |
519 |
|
|
! La 2e contrainte est que la convection ne doit pas sursaturer |
520 |
|
|
! la "detrainment layer", Neanmoins, une petite sursaturation |
521 |
|
|
! est acceptee (facteur ssfac). |
522 |
|
|
|
523 |
|
|
IF (hb(i,km1)<hbs(i,km1)) THEN |
524 |
|
|
tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) - & |
525 |
|
|
(hbh(i,kp1)-hc(i))*dp(i, k)/dp(i, kp1) |
526 |
|
|
tmp1 = tmp1/dp(i, k) |
527 |
|
|
tmp2 = gam(i, km1)*(sbh(i,k)-sc(i)+cldwtr(i)) - hbh(i, k) + hc(i) - & |
528 |
|
|
sc(i) + sbh(i, k) |
529 |
|
|
tmp3 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i, k) |
530 |
|
|
tmp4 = (dt/cats)*(hc(i)-hbs(i,k))*tmp2/(dp(i,km1)*(hbs(i,km1)-hb(i, & |
531 |
|
|
km1))) + tmp3 |
532 |
|
|
IF ((beta(i)*tmp4-tmp1)>0.0) THEN |
533 |
|
|
betamx = ssfac*(tmp1/tmp4) |
534 |
|
|
beta(i) = max(0.0, min(betamx,beta(i))) |
535 |
|
|
END IF |
536 |
|
|
ELSE |
537 |
|
|
beta(i) = 0.0 |
538 |
|
|
END IF |
539 |
|
|
|
540 |
|
|
! Third check to avoid introducing 2 delta x thermodynamic |
541 |
|
|
! noise in the vertical ... constrain adjusted h (or theta e) |
542 |
|
|
! so that the adjustment doesn't contribute to "kinks" in h |
543 |
|
|
|
544 |
|
|
g = min(0.0, hb(i,k)-hb(i,km1)) |
545 |
|
|
tmp3 = (hb(i,k)-hb(i,km1)-g)*(cats/dt)/(hc(i)-hbs(i,k)) |
546 |
|
|
tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) - & |
547 |
|
|
(hbh(i,kp1)-hc(i))*dp(i, k)/dp(i, kp1) |
548 |
|
|
tmp1 = tmp1/dp(i, k) |
549 |
|
|
tmp1 = tmp3*tmp1 + (hc(i)-hbh(i,kp1))/dp(i, k) |
550 |
|
|
tmp2 = tmp3*(1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i, k) + & |
551 |
|
|
(hc(i)-hbh(i,k)-cldwtr(i))*(1.0/dp(i,k)+1.0/dp(i,kp1)) |
552 |
|
|
IF ((beta(i)*tmp2-tmp1)>0.0) THEN |
553 |
|
|
betamx = 0.0 |
554 |
|
|
IF (tmp2/=0.0) betamx = tmp1/tmp2 |
555 |
|
|
beta(i) = max(0.0, min(betamx,beta(i))) |
556 |
|
|
END IF |
557 |
|
|
END IF |
558 |
|
|
END DO |
559 |
|
|
|
560 |
|
|
! Calculate mass flux required for stabilization. |
561 |
|
|
|
562 |
|
|
! Ensure that the convective mass flux, eta, is positive by |
563 |
|
|
! setting negative values of eta to zero.. |
564 |
|
|
! Ensure that estimated mass flux cannot move more than the |
565 |
|
|
! minimum of total mass contained in either layer k or layer k+1. |
566 |
|
|
! Also test for other pathological cases that result in non- |
567 |
|
|
! physical states and adjust eta accordingly. |
568 |
|
|
|
569 |
|
|
20 CONTINUE |
570 |
|
|
DO i = 1, klon |
571 |
|
|
IF (ldcum(i)) THEN |
572 |
|
|
beta(i) = max(0.0, beta(i)) |
573 |
|
|
tmp1 = hc(i) - hbs(i, k) |
574 |
|
|
tmp2 = ((1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i))-beta(i)*(1.0+gam( & |
575 |
|
|
i,k))*(sc(i)-sbh(i,k)))/dp(i, k) - (hbh(i,kp1)-hc(i))/dp(i, kp1) |
576 |
|
|
eta(i) = tmp1/(tmp2*rg*cats) |
577 |
|
|
tmass = min(dp(i,k), dp(i,kp1))/rg |
578 |
|
|
IF (eta(i)>tmass*rdt .OR. eta(i)<=0.0) eta(i) = 0.0 |
579 |
|
|
|
580 |
|
|
! Check on negative q in top layer (bound beta) |
581 |
|
|
|
582 |
|
|
IF (shc(i)-shbh(i,k)<0.0 .AND. beta(i)*eta(i)/=0.0) THEN |
583 |
|
|
denom = eta(i)*rg*dt*(shc(i)-shbh(i,k))/dp(i, km1) |
584 |
|
|
beta(i) = max(0.0, min(-0.999*shb(i,km1)/denom,beta(i))) |
585 |
|
|
END IF |
586 |
|
|
|
587 |
|
|
! Check on negative q in middle layer (zero eta) |
588 |
|
|
|
589 |
|
|
qtest1 = shb(i, k) + eta(i)*rg*dt*((shc(i)-shbh(i, & |
590 |
|
|
kp1))-(1.0-beta(i))*cldwtr(i)/rlvtt-beta(i)*(shc(i)-shbh(i, & |
591 |
|
|
k)))/dp(i, k) |
592 |
|
|
IF (qtest1<=0.0) eta(i) = 0.0 |
593 |
|
|
|
594 |
|
|
! Check on negative q in lower layer (bound eta) |
595 |
|
|
|
596 |
|
|
fac1 = -(shbh(i,kp1)-shc(i))/dp(i, kp1) |
597 |
|
|
qtest2 = shb(i, kp1) - eta(i)*rg*dt*fac1 |
598 |
|
|
IF (qtest2<0.0) THEN |
599 |
|
|
eta(i) = 0.99*shb(i, kp1)/(rg*dt*fac1) |
600 |
|
|
END IF |
601 |
|
|
END IF |
602 |
|
|
END DO |
603 |
|
|
|
604 |
|
|
|
605 |
|
|
! Calculate cloud water, rain water, and thermodynamic changes |
606 |
|
|
|
607 |
|
|
DO i = 1, klon |
608 |
|
|
IF (ldcum(i)) THEN |
609 |
|
|
etagdt = eta(i)*rg*dt |
610 |
|
|
cldwtr(i) = etagdt*cldwtr(i)/rlvtt/rg |
611 |
|
|
rnwtr(i) = (1.0-beta(i))*cldwtr(i) |
612 |
|
|
ds1(i) = etagdt*(sbh(i,kp1)-sc(i))/dp(i, kp1) |
613 |
|
|
dq1(i) = etagdt*(shbh(i,kp1)-shc(i))/dp(i, kp1) |
614 |
|
|
ds2(i) = (etagdt*(sc(i)-sbh(i,kp1))+rlvtt*rg*cldwtr(i)-beta(i)*etagdt & |
615 |
|
|
*(sc(i)-sbh(i,k)))/dp(i, k) |
616 |
|
|
dq2(i) = (etagdt*(shc(i)-shbh(i,kp1))-rg*rnwtr(i)-beta(i)*etagdt*(shc & |
617 |
|
|
(i)-shbh(i,k)))/dp(i, k) |
618 |
|
|
ds3(i) = beta(i)*(etagdt*(sc(i)-sbh(i,k))-rlvtt*rg*cldwtr(i))/dp(i, & |
619 |
|
|
km1) |
620 |
|
|
dq3(i) = beta(i)*etagdt*(shc(i)-shbh(i,k))/dp(i, km1) |
621 |
|
|
|
622 |
|
|
! Isolate convective fluxes for later diagnostics |
623 |
|
|
|
624 |
|
|
fslkp = eta(i)*(sc(i)-sbh(i,kp1)) |
625 |
|
|
fslkm = beta(i)*(eta(i)*(sc(i)-sbh(i,k))-rlvtt*cldwtr(i)*rdt) |
626 |
|
|
fqlkp = eta(i)*(shc(i)-shbh(i,kp1)) |
627 |
|
|
fqlkm = beta(i)*eta(i)*(shc(i)-shbh(i,k)) |
628 |
|
|
|
629 |
|
|
|
630 |
|
|
! Update thermodynamic profile (update sb, hb, & hbs later) |
631 |
|
|
|
632 |
|
|
tb(i, kp1) = tb(i, kp1) + ds1(i)/rcpd |
633 |
|
|
tb(i, k) = tb(i, k) + ds2(i)/rcpd |
634 |
|
|
tb(i, km1) = tb(i, km1) + ds3(i)/rcpd |
635 |
|
|
shb(i, kp1) = shb(i, kp1) + dq1(i) |
636 |
|
|
shb(i, k) = shb(i, k) + dq2(i) |
637 |
|
|
shb(i, km1) = shb(i, km1) + dq3(i) |
638 |
|
|
prec(i) = prec(i) + rnwtr(i)/rhoh2o |
639 |
|
|
|
640 |
|
|
! Update diagnostic information for final budget |
641 |
|
|
! Tracking temperature & specific humidity tendencies, |
642 |
|
|
! rainout term, convective mass flux, convective liquid |
643 |
|
|
! water static energy flux, and convective total water flux |
644 |
|
|
|
645 |
|
|
cmfdt(i, kp1) = cmfdt(i, kp1) + ds1(i)/rcpd*rdt |
646 |
|
|
cmfdt(i, k) = cmfdt(i, k) + ds2(i)/rcpd*rdt |
647 |
|
|
cmfdt(i, km1) = cmfdt(i, km1) + ds3(i)/rcpd*rdt |
648 |
|
|
cmfdq(i, kp1) = cmfdq(i, kp1) + dq1(i)*rdt |
649 |
|
|
cmfdq(i, k) = cmfdq(i, k) + dq2(i)*rdt |
650 |
|
|
cmfdq(i, km1) = cmfdq(i, km1) + dq3(i)*rdt |
651 |
|
|
cmfdqr(i, k) = cmfdqr(i, k) + (rg*rnwtr(i)/dp(i,k))*rdt |
652 |
|
|
cmfmc(i, kp1) = cmfmc(i, kp1) + eta(i) |
653 |
|
|
cmfmc(i, k) = cmfmc(i, k) + beta(i)*eta(i) |
654 |
|
|
cmfsl(i, kp1) = cmfsl(i, kp1) + fslkp |
655 |
|
|
cmfsl(i, k) = cmfsl(i, k) + fslkm |
656 |
|
|
cmflq(i, kp1) = cmflq(i, kp1) + rlvtt*fqlkp |
657 |
|
|
cmflq(i, k) = cmflq(i, k) + rlvtt*fqlkm |
658 |
|
|
qc(i, k) = (rg*rnwtr(i)/dp(i,k))*rdt |
659 |
|
|
END IF |
660 |
|
|
END DO |
661 |
|
|
|
662 |
|
|
! Next, convectively modify passive constituents |
663 |
|
|
|
664 |
|
|
DO m = 1, pcnst |
665 |
|
|
DO i = 1, klon |
666 |
|
|
IF (ldcum(i)) THEN |
667 |
|
|
|
668 |
|
|
! If any of the reported values of the constituent is negative in |
669 |
|
|
! the three adjacent levels, nothing will be done to the profile |
670 |
|
|
|
671 |
|
|
IF ((cmrb(i,kp1,m)<0.0) .OR. (cmrb(i,k,m)<0.0) .OR. (cmrb(i,km1, & |
672 |
|
|
m)<0.0)) GO TO 40 |
673 |
|
|
|
674 |
|
|
! Specify constituent interface values (linear interpolation) |
675 |
|
|
|
676 |
|
|
cmrh(i, k) = 0.5*(cmrb(i,km1,m)+cmrb(i,k,m)) |
677 |
|
|
cmrh(i, kp1) = 0.5*(cmrb(i,k,m)+cmrb(i,kp1,m)) |
678 |
|
|
|
679 |
|
|
! Specify perturbation properties of constituents in PBL |
680 |
|
|
|
681 |
|
|
pblhgt = max(pblh(i), 1.0) |
682 |
|
|
IF (gz(i,kp1)/rg<=pblhgt .AND. dzcld(i)==0.) THEN |
683 |
|
|
fac1 = max(0.0, 1.0-gz(i,kp1)/rg/pblhgt) |
684 |
|
|
cmrc(i) = cmrb(i, kp1, m) + cmrp(i, m)*fac1 |
685 |
|
|
ELSE |
686 |
|
|
cmrc(i) = cmrb(i, kp1, m) |
687 |
|
|
END IF |
688 |
|
|
|
689 |
|
|
! Determine fluxes, flux divergence => changes due to convection |
690 |
|
|
! Logic must be included to avoid producing negative values. A bit |
691 |
|
|
! messy since there are no a priori assumptions about profiles. |
692 |
|
|
! Tendency is modified (reduced) when pending disaster detected. |
693 |
|
|
|
694 |
|
|
etagdt = eta(i)*rg*dt |
695 |
|
|
botflx = etagdt*(cmrc(i)-cmrh(i,kp1)) |
696 |
|
|
topflx = beta(i)*etagdt*(cmrc(i)-cmrh(i,k)) |
697 |
|
|
dcmr1(i) = -botflx/dp(i, kp1) |
698 |
|
|
efac1 = 1.0 |
699 |
|
|
efac2 = 1.0 |
700 |
|
|
efac3 = 1.0 |
701 |
|
|
|
702 |
|
|
IF (cmrb(i,kp1,m)+dcmr1(i)<0.0) THEN |
703 |
|
|
efac1 = max(tiny, abs(cmrb(i,kp1,m)/dcmr1(i))-eps) |
704 |
|
|
END IF |
705 |
|
|
|
706 |
|
|
IF (efac1==tiny .OR. efac1>1.0) efac1 = 0.0 |
707 |
|
|
dcmr1(i) = -efac1*botflx/dp(i, kp1) |
708 |
|
|
dcmr2(i) = (efac1*botflx-topflx)/dp(i, k) |
709 |
|
|
|
710 |
|
|
IF (cmrb(i,k,m)+dcmr2(i)<0.0) THEN |
711 |
|
|
efac2 = max(tiny, abs(cmrb(i,k,m)/dcmr2(i))-eps) |
712 |
|
|
END IF |
713 |
|
|
|
714 |
|
|
IF (efac2==tiny .OR. efac2>1.0) efac2 = 0.0 |
715 |
|
|
dcmr2(i) = (efac1*botflx-efac2*topflx)/dp(i, k) |
716 |
|
|
dcmr3(i) = efac2*topflx/dp(i, km1) |
717 |
|
|
|
718 |
|
|
IF (cmrb(i,km1,m)+dcmr3(i)<0.0) THEN |
719 |
|
|
efac3 = max(tiny, abs(cmrb(i,km1,m)/dcmr3(i))-eps) |
720 |
|
|
END IF |
721 |
|
|
|
722 |
|
|
IF (efac3==tiny .OR. efac3>1.0) efac3 = 0.0 |
723 |
|
|
efac3 = min(efac2, efac3) |
724 |
|
|
dcmr2(i) = (efac1*botflx-efac3*topflx)/dp(i, k) |
725 |
|
|
dcmr3(i) = efac3*topflx/dp(i, km1) |
726 |
|
|
|
727 |
|
|
cmrb(i, kp1, m) = cmrb(i, kp1, m) + dcmr1(i) |
728 |
|
|
cmrb(i, k, m) = cmrb(i, k, m) + dcmr2(i) |
729 |
|
|
cmrb(i, km1, m) = cmrb(i, km1, m) + dcmr3(i) |
730 |
|
|
END IF |
731 |
|
|
40 END DO |
732 |
|
|
END DO ! end of m=1,pcnst loop |
733 |
|
|
|
734 |
|
|
IF (k==limcnv+1) GO TO 60 ! on ne pourra plus glisser |
735 |
|
|
|
736 |
|
|
! Dans la procedure de glissage ascendant, les variables thermo- |
737 |
|
|
! dynamiques des couches k et km1 servent au calcul des couches |
738 |
|
|
! superieures. Elles ont donc besoin d'une mise-a-jour. |
739 |
|
|
|
740 |
|
|
DO i = 1, klon |
741 |
|
|
IF (ldcum(i)) THEN |
742 |
|
|
zx_t = tb(i, k) |
743 |
|
|
zx_p = p(i, k) |
744 |
|
|
zx_q = shb(i, k) |
745 |
|
|
zdelta = max(0., sign(1.,rtt-zx_t)) |
746 |
|
|
zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta |
747 |
|
|
zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zx_q) |
748 |
|
|
zx_qs = r2es*foeew(zx_t, zdelta)/zx_p |
749 |
|
|
zx_qs = min(0.5, zx_qs) |
750 |
|
|
zcor = 1./(1.-retv*zx_qs) |
751 |
|
|
zx_qs = zx_qs*zcor |
752 |
|
|
zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) |
753 |
|
|
shbs(i, k) = zx_qs |
754 |
|
|
gam(i, k) = zx_gam |
755 |
|
|
|
756 |
|
|
zx_t = tb(i, km1) |
757 |
|
|
zx_p = p(i, km1) |
758 |
|
|
zx_q = shb(i, km1) |
759 |
|
|
zdelta = max(0., sign(1.,rtt-zx_t)) |
760 |
|
|
zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta |
761 |
|
|
zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zx_q) |
762 |
|
|
zx_qs = r2es*foeew(zx_t, zdelta)/zx_p |
763 |
|
|
zx_qs = min(0.5, zx_qs) |
764 |
|
|
zcor = 1./(1.-retv*zx_qs) |
765 |
|
|
zx_qs = zx_qs*zcor |
766 |
|
|
zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) |
767 |
|
|
shbs(i, km1) = zx_qs |
768 |
|
|
gam(i, km1) = zx_gam |
769 |
|
|
|
770 |
|
|
sb(i, k) = sb(i, k) + ds2(i) |
771 |
|
|
sb(i, km1) = sb(i, km1) + ds3(i) |
772 |
|
|
hb(i, k) = sb(i, k) + rlvtt*shb(i, k) |
773 |
|
|
hb(i, km1) = sb(i, km1) + rlvtt*shb(i, km1) |
774 |
|
|
hbs(i, k) = sb(i, k) + rlvtt*shbs(i, k) |
775 |
|
|
hbs(i, km1) = sb(i, km1) + rlvtt*shbs(i, km1) |
776 |
|
|
|
777 |
|
|
sbh(i, k) = 0.5*(sb(i,k)+sb(i,km1)) |
778 |
|
|
shbh(i, k) = qhalf(shb(i,km1), shb(i,k), shbs(i,km1), shbs(i,k)) |
779 |
|
|
hbh(i, k) = sbh(i, k) + rlvtt*shbh(i, k) |
780 |
|
|
sbh(i, km1) = 0.5*(sb(i,km1)+sb(i,k-2)) |
781 |
|
|
shbh(i, km1) = qhalf(shb(i,k-2), shb(i,km1), shbs(i,k-2), & |
782 |
|
|
shbs(i,km1)) |
783 |
|
|
hbh(i, km1) = sbh(i, km1) + rlvtt*shbh(i, km1) |
784 |
|
|
END IF |
785 |
|
|
END DO |
786 |
|
|
|
787 |
|
|
! Ensure that dzcld is reset if convective mass flux zero |
788 |
|
|
! specify the current vertical extent of the convective activity |
789 |
|
|
! top of convective layer determined by size of overshoot param. |
790 |
|
|
|
791 |
|
|
60 CONTINUE |
792 |
|
|
DO i = 1, klon |
793 |
|
|
etagt0 = eta(i) > 0.0 |
794 |
|
|
IF (.NOT. etagt0) dzcld(i) = 0.0 |
795 |
|
|
IF (etagt0 .AND. beta(i)>betamn) THEN |
796 |
|
|
ktp = km1 |
797 |
|
|
ELSE |
798 |
|
|
ktp = k |
799 |
|
|
END IF |
800 |
|
|
IF (etagt0) THEN |
801 |
|
|
cnt(i) = min(cnt(i), ktp) |
802 |
|
|
cnb(i) = max(cnb(i), k) |
803 |
|
|
END IF |
804 |
|
|
END DO |
805 |
|
|
70 END DO ! end of k loop |
806 |
|
|
|
807 |
|
|
! determine whether precipitation, prec, is frozen (snow) or not |
808 |
|
|
|
809 |
|
|
DO i = 1, klon |
810 |
|
|
IF (tb(i,klev)<tmelt .AND. tb(i,klev-1)<tmelt) THEN |
811 |
|
|
cmfprs(i) = prec(i)*rdt |
812 |
|
|
ELSE |
813 |
|
|
cmfprt(i) = prec(i)*rdt |
814 |
|
|
END IF |
815 |
|
|
END DO |
816 |
|
|
|
817 |
|
|
RETURN ! we're all done ... return to calling procedure |
818 |
|
|
END SUBROUTINE cmfmca |