source: LMDZ6/branches/LMDZ-tracers/libf/phylmd/conccm.F90 @ 3871

Last change on this file since 3871 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM
  • 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.6 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  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
137END SUBROUTINE conccm
138SUBROUTINE 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
56920  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
73140    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
79160  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
80570 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
818END SUBROUTINE cmfmca
Note: See TracBrowser for help on using the repository browser.