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