source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/phylmd/conccm.F @ 5095

Last change on this file since 5095 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

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