source: LMDZ5/trunk/libf/phylmd/conccm.F90 @ 2237

Last change on this file since 2237 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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