source: LMDZ.3.3/branches/LF/libf/phylmd/conccm.F @ 3

Last change on this file since 3 was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 30.8 KB
Line 
1      SUBROUTINE conccm (dtime,paprs,pplay,t,q,conv_q,
2     s                   d_t, d_q, rain, snow, kbascm, ktopcm)
3c
4      IMPLICIT none
5c======================================================================
6c Auteur(s): Z.X. Li (LMD/CNRS) date: le 14 mars 1996
7c Objet: Schema simple (avec flux de masse) pour la convection
8c        (schema standard du modele NCAR CCM2)
9c======================================================================
10#include "dimensions.h"
11#include "dimphy.h"
12#include "YOMCST.h"
13#include "YOETHF.h"
14c
15c Entree:
16      REAL dtime              ! pas d'integration
17      REAL paprs(klon,klev+1) ! pression inter-couche (Pa)
18      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
19      REAL t(klon,klev)       ! temperature (K)
20      REAL q(klon,klev)       ! humidite specifique (g/g)
21      REAL conv_q(klon,klev)  ! taux de convergence humidite (g/g/s)
22c Sortie:
23      REAL d_t(klon,klev)     ! incrementation temperature
24      REAL d_q(klon,klev)     ! incrementation vapeur
25      REAL rain(klon)         ! pluie (mm/s)
26      REAL snow(klon)         ! neige (mm/s)
27      INTEGER kbascm(klon)    ! niveau du bas de convection
28      INTEGER ktopcm(klon)    ! niveau du haut de convection
29c
30      REAL pt(klon,klev)
31      REAL pq(klon,klev)
32      REAL pres(klon,klev)
33      REAL dp(klon,klev)
34      REAL zgeom(klon,klev)
35      REAL cmfprs(klon)
36      REAL cmfprt(klon)
37      INTEGER ntop(klon)
38      INTEGER nbas(klon)
39      INTEGER i, k
40      REAL zlvdcp, zlsdcp, zdelta, zz, za, zb
41c
42      LOGICAL usekuo ! utiliser convection profonde (schema Kuo)
43      PARAMETER (usekuo=.TRUE.)
44c
45      REAL d_t_bis(klon,klev)
46      REAL d_q_bis(klon,klev)
47      REAL rain_bis(klon)
48      REAL snow_bis(klon)
49      INTEGER ibas_bis(klon)
50      INTEGER itop_bis(klon)
51      REAL d_ql_bis(klon,klev)
52      REAL rneb_bis(klon,klev)
53c
54c initialiser les variables de sortie (pour securite)
55      DO i = 1, klon
56         rain(i) = 0.0
57         snow(i) = 0.0
58         kbascm(i) = 0
59         ktopcm(i) = 0
60      ENDDO
61      DO k = 1, klev
62      DO i = 1, klon
63         d_t(i,k) = 0.0
64         d_q(i,k) = 0.0
65      ENDDO
66      ENDDO
67c
68c preparer les variables d'entree (attention: l'ordre des niveaux
69c verticaux augmente du haut vers le bas)
70      DO k = 1, klev
71      DO i = 1, klon
72         pt(i,k) = t(i,klev-k+1)
73         pq(i,k) = q(i,klev-k+1)
74         pres(i,k) = pplay(i,klev-k+1)
75         dp(i,k) = paprs(i,klev+1-k)-paprs(i,klev+1-k+1)
76      ENDDO
77      ENDDO
78      DO i = 1, klon
79         zgeom(i,klev) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
80     .                      * (paprs(i,1)-pplay(i,1))
81      ENDDO
82      DO k = 2, klev
83      DO i = 1, klon
84         zgeom(i,klev+1-k) = zgeom(i,klev+1-k+1)
85     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
86     .                   * (pplay(i,k-1)-pplay(i,k))
87      ENDDO
88      ENDDO
89c
90      CALL cmfmca(dtime, pres, dp, zgeom, pt, pq,
91     $                  cmfprt, cmfprs, ntop, nbas)
92C
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      ENDDO
98      ENDDO
99c
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      ENDDO
106c
107      IF (usekuo) THEN
108      CALL conkuo(dtime, paprs, pplay, t, q, conv_q,
109     s            d_t_bis, d_q_bis, d_ql_bis, rneb_bis,
110     s            rain_bis, snow_bis, ibas_bis, itop_bis)
111      DO k = 1, klev
112      DO i = 1, klon
113         d_t(i,k) = d_t(i,k) + d_t_bis(i,k)
114         d_q(i,k) = d_q(i,k) + d_q_bis(i,k)
115      ENDDO
116      ENDDO
117      DO i = 1, klon
118         rain(i) = rain(i) + rain_bis(i)
119         snow(i) = snow(i) + snow_bis(i)
120         kbascm(i) = MIN(kbascm(i),ibas_bis(i))
121         ktopcm(i) = MAX(ktopcm(i),itop_bis(i))
122      ENDDO
123      DO k = 1, klev ! eau liquide convective est
124      DO i = 1, klon ! dispersee dans l'air
125         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q(i,k))
126         zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q(i,k))
127         zdelta = MAX(0.,SIGN(1.,RTT-t(i,k)))
128         zz = d_ql_bis(i,k) ! re-evap. de l'eau liquide
129         zb = MAX(0.0,zz)
130         za = - MAX(0.0,zz) * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
131         d_t(i,k) = d_t(i,k) + za
132         d_q(i,k) = d_q(i,k) + zb
133      ENDDO
134      ENDDO
135      ENDIF
136c
137      RETURN
138      END
139      SUBROUTINE cmfmca(deltat, p, dp, gz,
140     $                  tb, shb,
141     $                  cmfprt, cmfprs, cnt, cnb)
142      IMPLICIT none
143C-----------------------------------------------------------------------
144C Moist convective mass flux procedure:
145C If stratification is unstable to nonentraining parcel ascent,
146C complete an adjustment making use of a simple cloud model
147C
148C Code generalized to allow specification of parcel ("updraft")
149C properties, as well as convective transport of an arbitrary
150C number of passive constituents (see cmrb array).
151C----------------------------Code History-------------------------------
152C Original version:  J. J. Hack, March 22, 1990
153C Standardized:      J. Rosinski, June 1992
154C Reviewed:          J. Hack, G. Taylor, August 1992
155c Adaptation au LMD: Z.X. Li, mars 1996 (reference: Hack 1994,
156c                    J. Geophys. Res. vol 99, D3, 5551-5568). J'ai
157c                    introduit les constantes et les fonctions thermo-
158c                    dynamiques du Centre Europeen. J'ai elimine le
159c                    re-indicage du code en esperant que cela pourra
160c                    simplifier la lecture et la comprehension.
161C-----------------------------------------------------------------------
162#include "dimensions.h"
163#include "dimphy.h"
164      INTEGER pcnst ! nombre de traceurs passifs
165      PARAMETER (pcnst=1)
166C------------------------------Arguments--------------------------------
167C Input arguments
168C
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)
173c
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
178c
179c Updated arguments:
180c
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)
184C
185C Output arguments
186C
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
198C------------------------------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)
221C
222C---------------------------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
265C                   !  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)
283c
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
290C
291      REAL tmp1, tmp2, tmp3, tmp4
292      REAL zx_t, zx_p, zx_q, zx_qs, zx_gam
293      REAL zcor, zdelta, zcvm5
294C
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))
301C
302C-----------------------------------------------------------------------
303c 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      ENDDO
309      ENDDO
310      ENDDO
311c
312c Les perturbations de la couche limite sont zero pour l'instant
313c
314      DO m = 1, pcnst
315      DO i = 1, klon
316         cmrp(i,m) = 0.0
317      ENDDO
318      ENDDO
319      DO i = 1, klon
320         thtap(i) = 0.0
321         shp(i) = 0.0
322         pblh(i) = 1.0
323      ENDDO
324C
325C Ensure that characteristic adjustment time scale (cmftau) assumed
326C in estimate of eta isn't smaller than model time scale (deltat)
327C
328      dt   = deltat
329      cats = MAX(dt,cmftau)
330      rdt  = 1.0/dt
331C
332C Compute sb,hb,shbs,hbs
333C
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      ENDDO
350      ENDDO
351C
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         ENDDO
358      ENDDO
359C
360C Compute sbh, shbh
361C
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         ENDDO
369      ENDDO
370C
371C Specify properties at top of model (not used, but filling anyway)
372C
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      ENDDO
378C
379C Zero vertically independent control, tendency & diagnostic arrays
380C
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      ENDDO
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        ENDDO
397      ENDDO
398C
399C Begin moist convective mass flux adjustment procedure.
400C Formalism ensures that negative cloud liquid water can never occur
401C
402      DO 70 k=klev-1,limcnv+1,-1
403         km1 = k - 1
404         kp1 = k + 1
405         DO 10 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
414C
415C Specification of "cloud base" conditions
416C
417            qprime    = 0.0
418            tprime    = 0.0
419C
420C Assign tprime within the PBL to be proportional to the quantity
421C thtap (which will be bounded by tpmax), passed to this routine by
422C the PBL routine.  Don't allow perturbation to produce a dry
423C adiabatically unstable parcel.  Assign qprime within the PBL to be
424C an appropriately modified value of the quantity shp (which will be
425C bounded by shpmax) passed to this routine by the PBL routine.  The
426C quantity qprime should be less than the local saturation value
427C (qsattp=qsat[t+tprime,p]).  In both cases, thtap and shp are
428C linearly reduced toward zero as the PBL top is approached.
429C
430            pblhgt = MAX(pblh(i),1.0)
431            IF (gz(i,kp1)/RG.LE.pblhgt .AND. dzcld(i).EQ.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,
436     $                        MAX(qsattp-shb(i,kp1),0.0))
437               qprime = MAX(qprime,shprme)
438            ELSE
439               tprime = 0.0
440               qprime = 0.0
441            ENDIF
442C
443C Specify "updraft" (in-cloud) thermodynamic properties
444C
445            sc (i)    = sb (i,kp1) + RCPD*tprime
446            shc(i)    = shb(i,kp1) + qprime
447            hc (i)    = sc (i    ) + RLVTT*shc(i)
448            flotab(i) = hc(i) - hbs(i,k)
449            dz        = dp(i,k)*RD*tb(i,k)/RG/p(i,k)
450            IF (flotab(i).gt.0.0) THEN
451               dzcld(i) = dzcld(i) + dz
452            ELSE
453               dzcld(i) = 0.0
454            ENDIF
455   10    CONTINUE
456C
457C Check on moist convective instability
458C
459         is = 0
460         DO i = 1, klon
461            IF (flotab(i).GT.0.0) THEN
462               ldcum(i) = .TRUE.
463               is = is + 1
464            ELSE
465               ldcum(i) = .FALSE.
466            ENDIF
467         ENDDO
468C
469         IF (is.EQ.0) THEN
470            DO i=1,klon
471               dzcld(i) = 0.0
472            ENDDO
473            GOTO 70
474         ENDIF
475C
476C Current level just below top level => no overshoot
477C
478         IF (k.le.limcnv+1) THEN
479            DO i=1,klon
480            IF (ldcum(i)) THEN
481               cldwtr(i) = sb(i,k)-sc(i)+flotab(i)/(1.0+gam(i,k))
482               cldwtr(i) = MAX(0.0,cldwtr(i))
483               beta(i)   = 0.0
484            ENDIF
485            ENDDO
486            GOTO 20
487         ENDIF
488C
489C First guess at overshoot parameter using crude buoyancy closure
490C 10% overshoot assumed as a minimum and 1-c0*dz maximum to start
491C If pre-existing supersaturation in detrainment layer, beta=0
492C cldwtr is temporarily equal to RLVTT*l (l=> liquid water)
493C
494         DO i=1,klon
495         IF (ldcum(i)) THEN
496            cldwtr(i) = sb(i,k)-sc(i)+flotab(i)/(1.0+gam(i,k))
497            cldwtr(i) = MAX(0.0,cldwtr(i))
498            betamx = 1.0 - c0*MAX(0.0,(dzcld(i)-dzmin))
499            b1        = (hc(i) - hbs(i,km1))*dp(i,km1)
500            b2        = (hc(i) - hbs(i,k  ))*dp(i,k  )
501            beta(i)   = MAX(betamn,MIN(betamx, 1.0+b1/b2))
502            IF (hbs(i,km1).le.hb(i,km1)) beta(i) = 0.0
503         ENDIF
504         ENDDO
505C
506C Bound maximum beta to ensure physically realistic solutions
507C
508C First check constrains beta so that eta remains positive
509C (assuming that eta is already positive for beta equal zero)
510c La premiere contrainte de beta est que le flux eta doit etre positif.
511C
512         DO i=1,klon
513         IF (ldcum(i)) THEN
514            tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
515     $            - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
516            tmp2 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))
517            IF ((beta(i)*tmp2-tmp1).GT.0.0) THEN
518               betamx = 0.99*(tmp1/tmp2)
519               beta(i) = MAX(0.0,MIN(betamx,beta(i)))
520            ENDIF
521C
522C Second check involves supersaturation of "detrainment layer"
523C small amount of supersaturation acceptable (by ssfac factor)
524c La 2e contrainte est que la convection ne doit pas sursaturer
525c la "detrainment layer", Neanmoins, une petite sursaturation
526c est acceptee (facteur ssfac).
527C
528            IF (hb(i,km1).lt.hbs(i,km1)) THEN
529               tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
530     $               - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
531               tmp1 = tmp1/dp(i,k)
532               tmp2 = gam(i,km1)*(sbh(i,k)-sc(i) + cldwtr(i)) -
533     $                 hbh(i,k) + hc(i) - sc(i) + sbh(i,k)
534               tmp3 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i,k)
535               tmp4 = (dt/cats)*(hc(i)-hbs(i,k))*tmp2
536     $               / (dp(i,km1)*(hbs(i,km1)-hb(i,km1))) + tmp3
537               IF ((beta(i)*tmp4-tmp1).GT.0.0) THEN
538                  betamx = ssfac*(tmp1/tmp4)
539                  beta(i)   = MAX(0.0,MIN(betamx,beta(i)))
540               ENDIF
541            ELSE
542               beta(i) = 0.0
543            ENDIF
544C
545C Third check to avoid introducing 2 delta x thermodynamic
546C noise in the vertical ... constrain adjusted h (or theta e)
547C so that the adjustment doesn't contribute to "kinks" in h
548C
549            g = MIN(0.0,hb(i,k)-hb(i,km1))
550            tmp3 = (hb(i,k)-hb(i,km1)-g)*(cats/dt) / (hc(i)-hbs(i,k))
551            tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
552     $            - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
553            tmp1 = tmp1/dp(i,k)
554            tmp1 = tmp3*tmp1 + (hc(i) - hbh(i,kp1))/dp(i,k)
555            tmp2 = tmp3*(1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i,k)
556     $            + (hc(i)-hbh(i,k)-cldwtr(i))
557     $             *(1.0/dp(i,k)+1.0/dp(i,kp1))
558            IF ((beta(i)*tmp2-tmp1).GT.0.0) THEN
559               betamx = 0.0
560               IF (tmp2.NE.0.0) betamx = tmp1/tmp2
561               beta(i) = MAX(0.0,MIN(betamx,beta(i)))
562            ENDIF
563         ENDIF
564         ENDDO
565C
566C Calculate mass flux required for stabilization.
567C
568C Ensure that the convective mass flux, eta, is positive by
569C setting negative values of eta to zero..
570C Ensure that estimated mass flux cannot move more than the
571C minimum of total mass contained in either layer k or layer k+1.
572C Also test for other pathological cases that result in non-
573C physical states and adjust eta accordingly.
574C
575   20    CONTINUE
576         DO i=1,klon
577         IF (ldcum(i)) THEN
578            beta(i) = MAX(0.0,beta(i))
579            tmp1 = hc(i) - hbs(i,k)
580            tmp2 = ((1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) -
581     $               beta(i)*(1.0+gam(i,k))*(sc(i)-sbh(i,k)))/dp(i,k) -
582     $              (hbh(i,kp1)-hc(i))/dp(i,kp1)
583            eta(i) = tmp1/(tmp2*RG*cats)
584            tmass = MIN(dp(i,k),dp(i,kp1))/RG
585            IF (eta(i).GT.tmass*rdt .OR. eta(i).LE.0.0) eta(i) = 0.0
586C
587C Check on negative q in top layer (bound beta)
588C
589            IF(shc(i)-shbh(i,k).LT.0.0 .AND. beta(i)*eta(i).NE.0.0)THEN
590               denom = eta(i)*RG*dt*(shc(i) - shbh(i,k))/dp(i,km1)
591               beta(i) = MAX(0.0,MIN(-0.999*shb(i,km1)/denom,beta(i)))
592            ENDIF
593C
594C Check on negative q in middle layer (zero eta)
595C
596            qtest1 = shb(i,k) + eta(i)*RG*dt*((shc(i) - shbh(i,kp1)) -
597     $               (1.0 - beta(i))*cldwtr(i)/RLVTT -
598     $               beta(i)*(shc(i) - shbh(i,k)))/dp(i,k)
599            IF (qtest1.le.0.0) eta(i) = 0.0
600C
601C Check on negative q in lower layer (bound eta)
602C
603            fac1 = -(shbh(i,kp1) - shc(i))/dp(i,kp1)
604            qtest2 = shb(i,kp1) - eta(i)*RG*dt*fac1
605            IF (qtest2 .lt. 0.0) THEN
606               eta(i) = 0.99*shb(i,kp1)/(RG*dt*fac1)
607            ENDIF
608         ENDIF
609         ENDDO
610C
611C
612C Calculate cloud water, rain water, and thermodynamic changes
613C
614         DO 30 i=1,klon
615         IF (ldcum(i)) THEN
616            etagdt = eta(i)*RG*dt
617            cldwtr(i) = etagdt*cldwtr(i)/RLVTT/RG
618            rnwtr(i) = (1.0 - beta(i))*cldwtr(i)
619            ds1(i) = etagdt*(sbh(i,kp1) - sc(i))/dp(i,kp1)
620            dq1(i) = etagdt*(shbh(i,kp1) - shc(i))/dp(i,kp1)
621            ds2(i) = (etagdt*(sc(i) - sbh(i,kp1)) +
622     $                RLVTT*RG*cldwtr(i) - beta(i)*etagdt*
623     $                (sc(i) - sbh(i,k)))/dp(i,k)
624            dq2(i) = (etagdt*(shc(i) - shbh(i,kp1)) -
625     $                RG*rnwtr(i) - beta(i)*etagdt*
626     $                (shc(i) - shbh(i,k)))/dp(i,k)
627            ds3(i) = beta(i)*(etagdt*(sc(i) - sbh(i,k)) -
628     $               RLVTT*RG*cldwtr(i))/dp(i,km1)
629            dq3(i) = beta(i)*etagdt*(shc(i) - shbh(i,k))/dp(i,km1)
630C
631C Isolate convective fluxes for later diagnostics
632C
633            fslkp = eta(i)*(sc(i) - sbh(i,kp1))
634            fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) -
635     $                       RLVTT*cldwtr(i)*rdt)
636            fqlkp = eta(i)*(shc(i) - shbh(i,kp1))
637            fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k))
638C
639C
640C Update thermodynamic profile (update sb, hb, & hbs later)
641C
642            tb (i,kp1) = tb(i,kp1)  + ds1(i) / RCPD
643            tb (i,k  ) = tb(i,k  )  + ds2(i) / RCPD
644            tb (i,km1) = tb(i,km1)  + ds3(i) / RCPD
645            shb(i,kp1) = shb(i,kp1) + dq1(i)
646            shb(i,k  ) = shb(i,k  ) + dq2(i)
647            shb(i,km1) = shb(i,km1) + dq3(i)
648            prec(i)    = prec(i)    + rnwtr(i)/rhoh2o
649C
650C Update diagnostic information for final budget
651C Tracking temperature & specific humidity tendencies,
652C rainout term, convective mass flux, convective liquid
653C water static energy flux, and convective total water flux
654C
655            cmfdt (i,kp1) = cmfdt (i,kp1) + ds1(i)/RCPD*rdt
656            cmfdt (i,k  ) = cmfdt (i,k  ) + ds2(i)/RCPD*rdt
657            cmfdt (i,km1) = cmfdt (i,km1) + ds3(i)/RCPD*rdt
658            cmfdq (i,kp1) = cmfdq (i,kp1) + dq1(i)*rdt
659            cmfdq (i,k  ) = cmfdq (i,k  ) + dq2(i)*rdt
660            cmfdq (i,km1) = cmfdq (i,km1) + dq3(i)*rdt
661            cmfdqr(i,k  ) = cmfdqr(i,k  ) + (RG*rnwtr(i)/dp(i,k))*rdt
662            cmfmc (i,kp1) = cmfmc (i,kp1) + eta(i)
663            cmfmc (i,k  ) = cmfmc (i,k  ) + beta(i)*eta(i)
664            cmfsl (i,kp1) = cmfsl (i,kp1) + fslkp
665            cmfsl (i,k  ) = cmfsl (i,k  ) + fslkm
666            cmflq (i,kp1) = cmflq (i,kp1) + RLVTT*fqlkp
667            cmflq (i,k  ) = cmflq (i,k  ) + RLVTT*fqlkm
668            qc    (i,k  ) =                (RG*rnwtr(i)/dp(i,k))*rdt
669         ENDIF
670   30    CONTINUE
671C
672C Next, convectively modify passive constituents
673C
674         DO 50 m=1,pcnst
675         DO 40 i=1,klon
676         IF (ldcum(i)) THEN
677C
678C If any of the reported values of the constituent is negative in
679C the three adjacent levels, nothing will be done to the profile
680C
681            IF ((cmrb(i,kp1,m).LT.0.0) .OR.
682     $          (cmrb(i,k,m).LT.0.0) .OR.
683     $          (cmrb(i,km1,m).LT.0.0)) GOTO 40
684C
685C Specify constituent interface values (linear interpolation)
686C
687            cmrh(i,k  ) = 0.5*(cmrb(i,km1,m) + cmrb(i,k  ,m))
688            cmrh(i,kp1) = 0.5*(cmrb(i,k  ,m) + cmrb(i,kp1,m))
689C
690C Specify perturbation properties of constituents in PBL
691C
692            pblhgt = MAX(pblh(i),1.0)
693            IF (gz(i,kp1)/RG.LE.pblhgt .AND. dzcld(i).EQ.0.) THEN
694               fac1 = MAX(0.0,1.0-gz(i,kp1)/RG/pblhgt)
695               cmrc(i) = cmrb(i,kp1,m) + cmrp(i,m)*fac1
696            ELSE
697               cmrc(i) = cmrb(i,kp1,m)
698            ENDIF
699C
700C Determine fluxes, flux divergence => changes due to convection
701C Logic must be included to avoid producing negative values. A bit
702C messy since there are no a priori assumptions about profiles.
703C Tendency is modified (reduced) when pending disaster detected.
704C
705            etagdt = eta(i)*RG*dt
706            botflx   = etagdt*(cmrc(i) - cmrh(i,kp1))
707            topflx   = beta(i)*etagdt*(cmrc(i)-cmrh(i,k))
708            dcmr1(i) = -botflx/dp(i,kp1)
709            efac1    = 1.0
710            efac2    = 1.0
711            efac3    = 1.0
712C
713            IF (cmrb(i,kp1,m)+dcmr1(i) .LT. 0.0) THEN
714               efac1 = MAX(tiny,ABS(cmrb(i,kp1,m)/dcmr1(i)) - eps)
715            ENDIF
716C
717            IF (efac1.EQ.tiny .OR. efac1.GT.1.0) efac1 = 0.0
718            dcmr1(i) = -efac1*botflx/dp(i,kp1)
719            dcmr2(i) = (efac1*botflx - topflx)/dp(i,k)
720
721            IF (cmrb(i,k,m)+dcmr2(i) .LT. 0.0) THEN
722               efac2 = MAX(tiny,ABS(cmrb(i,k  ,m)/dcmr2(i)) - eps)
723            ENDIF
724C
725            IF (efac2.EQ.tiny .OR. efac2.GT.1.0) efac2 = 0.0
726            dcmr2(i) = (efac1*botflx - efac2*topflx)/dp(i,k)
727            dcmr3(i) = efac2*topflx/dp(i,km1)
728C
729            IF (cmrb(i,km1,m)+dcmr3(i) .LT. 0.0) THEN
730               efac3 = MAX(tiny,ABS(cmrb(i,km1,m)/dcmr3(i)) - eps)
731            ENDIF
732C
733            IF (efac3.EQ.tiny .OR. efac3.GT.1.0) efac3 = 0.0
734            efac3    = MIN(efac2,efac3)
735            dcmr2(i) = (efac1*botflx - efac3*topflx)/dp(i,k)
736            dcmr3(i) = efac3*topflx/dp(i,km1)
737C
738            cmrb(i,kp1,m) = cmrb(i,kp1,m) + dcmr1(i)
739            cmrb(i,k  ,m) = cmrb(i,k  ,m) + dcmr2(i)
740            cmrb(i,km1,m) = cmrb(i,km1,m) + dcmr3(i)
741         ENDIF
742   40    CONTINUE
743   50    CONTINUE              ! end of m=1,pcnst loop
744C
745         IF (k.EQ.limcnv+1) GOTO 60 ! on ne pourra plus glisser
746c
747c Dans la procedure de glissage ascendant, les variables thermo-
748c dynamiques des couches k et km1 servent au calcul des couches
749c superieures. Elles ont donc besoin d'une mise-a-jour.
750C
751         DO i = 1, klon
752         IF (ldcum(i)) THEN
753            zx_t = tb(i,k)
754            zx_p = p(i,k)
755            zx_q = shb(i,k)
756              zdelta=MAX(0.,SIGN(1.,RTT-zx_t))
757              zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
758              zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q)
759              zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p
760              zx_qs=MIN(0.5,zx_qs)
761              zcor=1./(1.-retv*zx_qs)
762              zx_qs=zx_qs*zcor
763              zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor)
764            shbs(i,k) = zx_qs
765            gam(i,k) = zx_gam
766c
767            zx_t = tb(i,km1)
768            zx_p = p(i,km1)
769            zx_q = shb(i,km1)
770              zdelta=MAX(0.,SIGN(1.,RTT-zx_t))
771              zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
772              zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q)
773              zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p
774              zx_qs=MIN(0.5,zx_qs)
775              zcor=1./(1.-retv*zx_qs)
776              zx_qs=zx_qs*zcor
777              zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor)
778            shbs(i,km1) = zx_qs
779            gam(i,km1) = zx_gam
780C
781            sb (i,k  ) = sb(i,k  ) + ds2(i)
782            sb (i,km1) = sb(i,km1) + ds3(i)
783            hb (i,k  ) = sb(i,k  ) + RLVTT*shb(i,k)
784            hb (i,km1) = sb(i,km1) + RLVTT*shb(i,km1)
785            hbs(i,k  ) = sb(i,k  ) + RLVTT*shbs(i,k  )
786            hbs(i,km1) = sb(i,km1) + RLVTT*shbs(i,km1)
787C
788            sbh (i,k) = 0.5*(sb(i,k) + sb(i,km1))
789            shbh(i,k) = qhalf(shb(i,km1),shb(i,k)
790     $                       ,shbs(i,km1),shbs(i,k))
791            hbh (i,k) = sbh(i,k) + RLVTT*shbh(i,k)
792            sbh (i,km1) = 0.5*(sb(i,km1) + sb(i,k-2))
793            shbh(i,km1) = qhalf(shb(i,k-2),shb(i,km1),
794     $                    shbs(i,k-2),shbs(i,km1))
795            hbh (i,km1) = sbh(i,km1) + RLVTT*shbh(i,km1)
796         ENDIF
797         ENDDO
798C
799C Ensure that dzcld is reset if convective mass flux zero
800C specify the current vertical extent of the convective activity
801C top of convective layer determined by size of overshoot param.
802C
803   60    CONTINUE
804         DO i=1,klon
805            etagt0 = eta(i).gt.0.0
806            IF (.not.etagt0) dzcld(i) = 0.0
807            IF (etagt0 .and. beta(i).gt.betamn) THEN
808               ktp = km1
809            ELSE
810               ktp = k
811            ENDIF
812            IF (etagt0) THEN
813               cnt(i) = MIN(cnt(i),ktp)
814               cnb(i) = MAX(cnb(i),k)
815            ENDIF
816         ENDDO
817   70 CONTINUE        ! end of k loop
818C
819C determine whether precipitation, prec, is frozen (snow) or not
820C
821      DO i=1,klon
822         IF (tb(i,klev).LT.tmelt .AND. tb(i,klev-1).lt.tmelt) THEN
823             cmfprs(i) = prec(i)*rdt
824         ELSE
825             cmfprt(i) = prec(i)*rdt
826         ENDIF
827      ENDDO
828C
829      RETURN  ! we're all done ... return to calling procedure
830      END
Note: See TracBrowser for help on using the repository browser.