GCC Code Coverage Report


Directory: ./
File: phys/conccm.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 386 0.0%
Branches: 0 184 0.0%

Line Branch Exec Source
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
819