source: LMDZ5/trunk/libf/phylmd/conflx.F @ 1923

Last change on this file since 1923 was 1923, checked in by lguez, 10 years ago

Replace test on level index (which was probably tailored for a vertical grid with 19 levels) by test on sigma = pressure / surface_pressure.
  • 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: 55.9 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE conflx (dtime,pres_h,pres_f,
5     e                   t, q, con_t, con_q, pqhfl, w,
6     s                   d_t, d_q, rain, snow,
7     s                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
8     s                   kcbot, kctop, kdtop, pmflxr, pmflxs)
9c
10      USE dimphy
11      IMPLICIT none
12c======================================================================
13c Auteur(s): Z.X. Li (LMD/CNRS) date: 19941014
14c Objet: Schema flux de masse pour la convection
15c        (schema de Tiedtke avec qqs modifications mineures)
16c Dec.97: Prise en compte des modifications introduites par
17c         Olivier Boucher et Alexandre Armengaud pour melange
18c         et lessivage des traceurs passifs.
19c======================================================================
20c     yminclude "dimensions.h"
21c     yminclude "dimphy.h"
22      include "YOMCST.h"
23      include "YOETHF.h"
24c Entree:
25      REAL dtime            ! pas d'integration (s)
26      REAL pres_h(klon,klev+1) ! pression half-level (Pa)
27      REAL pres_f(klon,klev)! pression full-level (Pa)
28      REAL t(klon,klev)     ! temperature (K)
29      REAL q(klon,klev)     ! humidite specifique (g/g)
30      REAL w(klon,klev)     ! vitesse verticale (Pa/s)
31      REAL con_t(klon,klev) ! convergence de temperature (K/s)
32      REAL con_q(klon,klev) ! convergence de l'eau vapeur (g/g/s)
33      REAL pqhfl(klon)      ! evaporation (negative vers haut) mm/s
34c Sortie:
35      REAL d_t(klon,klev)   ! incrementation de temperature
36      REAL d_q(klon,klev)   ! incrementation d'humidite
37      REAL pmfu(klon,klev)  ! flux masse (kg/m2/s) panache ascendant
38      REAL pmfd(klon,klev)  ! flux masse (kg/m2/s) panache descendant
39      REAL pen_u(klon,klev)
40      REAL pen_d(klon,klev)
41      REAL pde_u(klon,klev)
42      REAL pde_d(klon,klev)
43      REAL rain(klon)       ! pluie (mm/s)
44      REAL snow(klon)       ! neige (mm/s)
45      REAL pmflxr(klon,klev+1)
46      REAL pmflxs(klon,klev+1)
47      INTEGER kcbot(klon)  ! niveau du bas de la convection
48      INTEGER kctop(klon)  ! niveau du haut de la convection
49      INTEGER kdtop(klon)  ! niveau du haut des downdrafts
50c Local:
51      REAL pt(klon,klev)
52      REAL pq(klon,klev)
53      REAL pqs(klon,klev)
54      REAL pvervel(klon,klev)
55      LOGICAL land(klon)
56c
57      REAL d_t_bis(klon,klev)
58      REAL d_q_bis(klon,klev)
59      REAL paprs(klon,klev+1)
60      REAL paprsf(klon,klev)
61      REAL zgeom(klon,klev)
62      REAL zcvgq(klon,klev)
63      REAL zcvgt(klon,klev)
64cAA
65      REAL zmfu(klon,klev)
66      REAL zmfd(klon,klev)
67      REAL zen_u(klon,klev)
68      REAL zen_d(klon,klev)
69      REAL zde_u(klon,klev)
70      REAL zde_d(klon,klev)
71      REAL zmflxr(klon,klev+1)
72      REAL zmflxs(klon,klev+1)
73cAA
74
75c
76      INTEGER i, k
77      REAL zdelta, zqsat
78c
79      include "FCTTRE.h"
80c
81c initialiser les variables de sortie (pour securite)
82      DO i = 1, klon
83         rain(i) = 0.0
84         snow(i) = 0.0
85         kcbot(i) = 0
86         kctop(i) = 0
87         kdtop(i) = 0
88      ENDDO
89      DO k = 1, klev
90      DO i = 1, klon
91         d_t(i,k) = 0.0
92         d_q(i,k) = 0.0
93         pmfu(i,k) = 0.0
94         pmfd(i,k) = 0.0
95         pen_u(i,k) = 0.0
96         pde_u(i,k) = 0.0
97         pen_d(i,k) = 0.0
98         pde_d(i,k) = 0.0
99         zmfu(i,k) = 0.0
100         zmfd(i,k) = 0.0
101         zen_u(i,k) = 0.0
102         zde_u(i,k) = 0.0
103         zen_d(i,k) = 0.0
104         zde_d(i,k) = 0.0
105      ENDDO
106      ENDDO
107      DO k = 1, klev+1
108      DO i = 1, klon
109         zmflxr(i,k) = 0.0
110         zmflxs(i,k) = 0.0
111      ENDDO
112      ENDDO
113c
114c calculer la nature du sol (pour l'instant, ocean partout)
115      DO i = 1, klon
116         land(i) = .FALSE.
117      ENDDO
118c
119c preparer les variables d'entree (attention: l'ordre des niveaux
120c verticaux augmente du haut vers le bas)
121      DO k = 1, klev
122      DO i = 1, klon
123         pt(i,k) = t(i,klev-k+1)
124         pq(i,k) = q(i,klev-k+1)
125         paprsf(i,k) = pres_f(i,klev-k+1)
126         paprs(i,k) = pres_h(i,klev+1-k+1)
127         pvervel(i,k) = w(i,klev+1-k)
128         zcvgt(i,k) = con_t(i,klev-k+1)
129         zcvgq(i,k) = con_q(i,klev-k+1)
130c
131         zdelta=MAX(0.,SIGN(1.,RTT-pt(i,k)))
132         zqsat=R2ES*FOEEW ( pt(i,k), zdelta ) / paprsf(i,k)
133         zqsat=MIN(0.5,zqsat)
134         zqsat=zqsat/(1.-RETV  *zqsat)
135         pqs(i,k) = zqsat
136      ENDDO
137      ENDDO
138      DO i = 1, klon
139         paprs(i,klev+1) = pres_h(i,1)
140         zgeom(i,klev) = RD * pt(i,klev)
141     .                   / (0.5*(paprs(i,klev+1)+paprsf(i,klev)))
142     .                   * (paprs(i,klev+1)-paprsf(i,klev))
143      ENDDO
144      DO k = klev-1, 1, -1
145      DO i = 1, klon
146         zgeom(i,k) = zgeom(i,k+1)
147     .              + RD * 0.5*(pt(i,k+1)+pt(i,k)) / paprs(i,k+1)
148     .                   * (paprsf(i,k+1)-paprsf(i,k))
149      ENDDO
150      ENDDO
151c
152c appeler la routine principale
153c
154      CALL flxmain(dtime, pt, pq, pqs, pqhfl,
155     .             paprsf, paprs, zgeom, land, zcvgt, zcvgq, pvervel,
156     .             rain, snow, kcbot, kctop, kdtop,
157     .             zmfu, zmfd, zen_u, zde_u, zen_d, zde_d,
158     .             d_t_bis, d_q_bis, zmflxr, zmflxs)
159C
160cAA--------------------------------------------------------
161cAA rem : De la meme facon que l'on effectue le reindicage
162cAA       pour la temperature t et le champ q
163cAA       on reindice les flux necessaires a la convection
164cAA       des traceurs
165cAA--------------------------------------------------------
166      DO k = 1, klev
167      DO i = 1, klon
168         d_q(i,klev+1-k) = dtime*d_q_bis(i,k)
169         d_t(i,klev+1-k) = dtime*d_t_bis(i,k)
170      ENDDO
171      ENDDO
172c
173      DO i = 1, klon
174         pmfu(i,1)= 0.
175         pmfd(i,1)= 0.
176         pen_d(i,1)= 0.
177         pde_d(i,1)= 0.
178      ENDDO
179     
180      DO k = 2, klev
181      DO i = 1, klon
182         pmfu(i,klev+2-k)= zmfu(i,k)
183         pmfd(i,klev+2-k)= zmfd(i,k)
184      ENDDO
185      ENDDO
186c
187      DO k = 1, klev
188      DO i = 1, klon
189         pen_u(i,klev+1-k)=  zen_u(i,k)
190         pde_u(i,klev+1-k)=  zde_u(i,k)
191      ENDDO
192      ENDDO
193c
194      DO k = 1, klev-1
195      DO i = 1, klon
196         pen_d(i,klev+1-k)= -zen_d(i,k+1)
197         pde_d(i,klev+1-k)= -zde_d(i,k+1)
198      ENDDO
199      ENDDO
200
201      DO k = 1, klev+1
202      DO i = 1, klon
203         pmflxr(i,klev+2-k)= zmflxr(i,k)
204         pmflxs(i,klev+2-k)= zmflxs(i,k)
205      ENDDO
206      ENDDO
207
208      RETURN
209      END subroutine conflx
210c--------------------------------------------------------------------
211      SUBROUTINE flxmain(pdtime, pten, pqen, pqsen, pqhfl, pap, paph,
212     .                   pgeo, ldland, ptte, pqte, pvervel,
213     .                   prsfc, pssfc, kcbot, kctop, kdtop,
214c     *                   ldcum, ktype,
215     .                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
216     .                   dt_con, dq_con, pmflxr, pmflxs)
217      USE dimphy
218      IMPLICIT none
219C     ------------------------------------------------------------------
220c     yminclude "dimensions.h"
221c     yminclude "dimphy.h"
222      include "YOMCST.h"
223      include "YOETHF.h"
224      include "YOECUMF.h"
225C     ----------------------------------------------------------------
226      REAL pten(klon,klev), pqen(klon,klev), pqsen(klon,klev)
227      REAL ptte(klon,klev)
228      REAL pqte(klon,klev)
229      REAL pvervel(klon,klev)
230      REAL pgeo(klon,klev), pap(klon,klev), paph(klon,klev+1)
231      REAL pqhfl(klon)
232c
233      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
234      REAL plude(klon,klev)
235      REAL pmfu(klon,klev)
236      REAL prsfc(klon), pssfc(klon)
237      INTEGER  kcbot(klon), kctop(klon), ktype(klon)
238      LOGICAL  ldland(klon), ldcum(klon)
239c
240      REAL ztenh(klon,klev), zqenh(klon,klev), zqsenh(klon,klev)
241      REAL zgeoh(klon,klev)
242      REAL zmfub(klon), zmfub1(klon)
243      REAL zmfus(klon,klev), zmfuq(klon,klev), zmful(klon,klev)
244      REAL zdmfup(klon,klev), zdpmel(klon,klev)
245      REAL zentr(klon), zhcbase(klon)
246      REAL zdqpbl(klon), zdqcv(klon), zdhpbl(klon)
247      REAL zrfl(klon)
248      REAL pmflxr(klon,klev+1)
249      REAL pmflxs(klon,klev+1)
250      INTEGER  ilab(klon,klev), ictop0(klon)
251      LOGICAL  llo1
252      REAL dt_con(klon,klev), dq_con(klon,klev)
253      REAL zmfmax, zdh
254      REAL pdtime, zqumqe, zdqmin, zalvdcp, zhsat, zzz
255      REAL zhhat, zpbmpt, zgam, zeps, zfac
256      INTEGER i, k, ikb, itopm2, kcum
257c
258      REAL pen_u(klon,klev), pde_u(klon,klev)
259      REAL pen_d(klon,klev), pde_d(klon,klev)
260c
261      REAL ptd(klon,klev), pqd(klon,klev), pmfd(klon,klev)
262      REAL zmfds(klon,klev), zmfdq(klon,klev), zdmfdp(klon,klev)
263      INTEGER kdtop(klon)
264      LOGICAL lddraf(klon)
265C---------------------------------------------------------------------
266      LOGICAL firstcal
267      SAVE firstcal
268      DATA firstcal / .TRUE. /
269c$OMP THREADPRIVATE(firstcal)
270C---------------------------------------------------------------------
271      IF (firstcal) THEN
272         CALL flxsetup
273         firstcal = .FALSE.
274      ENDIF
275C---------------------------------------------------------------------
276      DO i = 1, klon
277         ldcum(i) = .FALSE.
278      ENDDO
279      DO k = 1, klev
280      DO i = 1, klon
281         dt_con(i,k) = 0.0
282         dq_con(i,k) = 0.0
283      ENDDO
284      ENDDO
285c----------------------------------------------------------------------
286c initialiser les variables et faire l'interpolation verticale
287c----------------------------------------------------------------------
288      CALL flxini(pten, pqen, pqsen, pgeo,
289     .     paph, zgeoh, ztenh, zqenh, zqsenh,
290     .     ptu, pqu, ptd, pqd, pmfd, zmfds, zmfdq, zdmfdp,
291     .     pmfu, zmfus, zmfuq, zdmfup,
292     .     zdpmel, plu, plude, ilab, pen_u, pde_u, pen_d, pde_d)
293c---------------------------------------------------------------------
294c determiner les valeurs au niveau de base de la tour convective
295c---------------------------------------------------------------------
296      CALL flxbase(ztenh, zqenh, zgeoh, paph,
297     *            ptu, pqu, plu, ldcum, kcbot, ilab)
298c---------------------------------------------------------------------
299c calculer la convergence totale de l'humidite et celle en provenance
300c de la couche limite, plus precisement, la convergence integree entre
301c le sol et la base de la convection. Cette derniere convergence est
302c comparee avec l'evaporation obtenue dans la couche limite pour
303c determiner le type de la convection
304c---------------------------------------------------------------------
305      k=1
306      DO i = 1, klon
307         zdqcv(i) = pqte(i,k)*(paph(i,k+1)-paph(i,k))
308         zdhpbl(i) = 0.0
309         zdqpbl(i) = 0.0
310      ENDDO
311c
312      DO k=2,klev
313      DO i = 1, klon
314          zdqcv(i)=zdqcv(i)+pqte(i,k)*(paph(i,k+1)-paph(i,k))
315          IF (k.GE.kcbot(i)) THEN
316             zdqpbl(i)=zdqpbl(i)+pqte(i,k)*(paph(i,k+1)-paph(i,k))
317             zdhpbl(i)=zdhpbl(i)+(RCPD*ptte(i,k)+RLVTT*pqte(i,k))
318     .                          *(paph(i,k+1)-paph(i,k))
319          ENDIF
320      ENDDO
321      ENDDO
322c
323      DO i = 1, klon
324         ktype(i) = 2
325         if (zdqcv(i).GT.MAX(0.,-1.5*pqhfl(i)*RG)) ktype(i) = 1
326ccc         if (zdqcv(i).GT.MAX(0.,-1.1*pqhfl(i)*RG)) ktype(i) = 1
327      ENDDO
328c
329c---------------------------------------------------------------------
330c determiner le flux de masse entrant a travers la base.
331c on ignore, pour l'instant, l'effet du panache descendant
332c---------------------------------------------------------------------
333      DO i = 1, klon
334         ikb=kcbot(i)
335         zqumqe=pqu(i,ikb)+plu(i,ikb)-zqenh(i,ikb)
336         zdqmin=MAX(0.01*zqenh(i,ikb),1.E-10)
337         IF (zdqpbl(i).GT.0..AND.zqumqe.GT.zdqmin.AND.ldcum(i)) THEN
338            zmfub(i) = zdqpbl(i)/(RG*MAX(zqumqe,zdqmin))
339         ELSE
340            zmfub(i) = 0.01
341            ldcum(i)=.FALSE.
342         ENDIF
343         IF (ktype(i).EQ.2) THEN
344            zdh = RCPD*(ptu(i,ikb)-ztenh(i,ikb)) + RLVTT*zqumqe
345            zdh = RG * MAX(zdh,1.0E5*zdqmin)
346            IF (zdhpbl(i).GT.0..AND.ldcum(i))zmfub(i)=zdhpbl(i)/zdh
347         ENDIF
348         zmfmax = (paph(i,ikb)-paph(i,ikb-1)) / (RG*pdtime)
349         zmfub(i) = MIN(zmfub(i),zmfmax)
350         zentr(i) = ENTRSCV
351         IF (ktype(i).EQ.1) zentr(i) = ENTRPEN
352      ENDDO
353C-----------------------------------------------------------------------
354C DETERMINE CLOUD ASCENT FOR ENTRAINING PLUME
355C-----------------------------------------------------------------------
356c (A) calculer d'abord la hauteur "theorique" de la tour convective sans
357c     considerer l'entrainement ni le detrainement du panache, sachant
358c     ces derniers peuvent abaisser la hauteur theorique.
359c
360      DO i = 1, klon
361         ikb=kcbot(i)
362         zhcbase(i)=RCPD*ptu(i,ikb)+zgeoh(i,ikb)+RLVTT*pqu(i,ikb)
363         ictop0(i)=kcbot(i)-1
364      ENDDO
365c
366      zalvdcp=RLVTT/RCPD
367      DO k=klev-1,3,-1
368      DO i = 1, klon
369         zhsat=RCPD*ztenh(i,k)+zgeoh(i,k)+RLVTT*zqsenh(i,k)
370         zgam=R5LES*zalvdcp*zqsenh(i,k)/
371     .        ((1.-RETV  *zqsenh(i,k))*(ztenh(i,k)-R4LES)**2)
372         zzz=RCPD*ztenh(i,k)*0.608
373         zhhat=zhsat-(zzz+zgam*zzz)/(1.+zgam*zzz/RLVTT)*
374     .               MAX(zqsenh(i,k)-zqenh(i,k),0.)
375         IF(k.LT.ictop0(i).AND.zhcbase(i).GT.zhhat) ictop0(i)=k
376      ENDDO
377      ENDDO
378c
379c (B) calculer le panache ascendant
380c
381      CALL flxasc(pdtime,ztenh, zqenh, pten, pqen, pqsen,
382     .     pgeo, zgeoh, pap, paph, pqte, pvervel,
383     .     ldland, ldcum, ktype, ilab,
384     .     ptu, pqu, plu, pmfu, zmfub, zentr,
385     .     zmfus, zmfuq, zmful, plude, zdmfup,
386     .     kcbot, kctop, ictop0, kcum, pen_u, pde_u)
387      IF (kcum.EQ.0) GO TO 1000
388C
389C verifier l'epaisseur de la convection et changer eventuellement
390c le taux d'entrainement/detrainement
391C
392      DO i = 1, klon
393         zpbmpt=paph(i,kcbot(i))-paph(i,kctop(i))
394         IF(ldcum(i).AND.ktype(i).EQ.1.AND.zpbmpt.LT.2.E4)ktype(i)=2
395         IF(ldcum(i)) ictop0(i)=kctop(i)
396         IF(ktype(i).EQ.2) zentr(i)=ENTRSCV
397      ENDDO
398c
399      IF (lmfdd) THEN  ! si l'on considere le panache descendant
400c
401c calculer la precipitation issue du panache ascendant pour
402c determiner l'existence du panache descendant dans la convection
403      DO i = 1, klon
404         zrfl(i)=zdmfup(i,1)
405      ENDDO
406      DO k=2,klev
407      DO i = 1, klon
408         zrfl(i)=zrfl(i)+zdmfup(i,k)
409      ENDDO
410      ENDDO
411c
412c determiner le LFS (level of free sinking: niveau de plonge libre)
413      CALL flxdlfs(ztenh, zqenh, zgeoh, paph, ptu, pqu,
414     *     ldcum,    kcbot,    kctop,    zmfub,    zrfl,
415     *     ptd,      pqd,
416     *     pmfd,     zmfds,    zmfdq,    zdmfdp,
417     *     kdtop,    lddraf)
418c
419c calculer le panache descendant
420      CALL flxddraf(ztenh,    zqenh,
421     *     zgeoh,    paph,     zrfl,
422     *     ptd,      pqd,
423     *     pmfd,     zmfds,    zmfdq,    zdmfdp,
424     *     lddraf, pen_d, pde_d)
425c
426c calculer de nouveau le flux de masse entrant a travers la base
427c de la convection, sachant qu'il a ete modifie par le panache
428c descendant
429      DO i = 1, klon
430      IF (lddraf(i)) THEN
431         ikb = kcbot(i)
432         llo1 = PMFD(i,ikb).LT.0.
433         zeps = 0.
434         IF ( llo1 ) zeps = CMFDEPS
435         zqumqe = pqu(i,ikb)+plu(i,ikb)-
436     .            zeps*pqd(i,ikb)-(1.-zeps)*zqenh(i,ikb)
437         zdqmin = MAX(0.01*zqenh(i,ikb),1.E-10)
438         zmfmax = (paph(i,ikb)-paph(i,ikb-1)) / (RG*pdtime)
439         IF (zdqpbl(i).GT.0..AND.zqumqe.GT.zdqmin.AND.ldcum(i)
440     .       .AND.zmfub(i).LT.zmfmax) THEN
441            zmfub1(i) = zdqpbl(i) / (RG*MAX(zqumqe,zdqmin))
442         ELSE
443            zmfub1(i) = zmfub(i)
444         ENDIF
445         IF (ktype(i).EQ.2) THEN
446            zdh = RCPD*(ptu(i,ikb)-zeps*ptd(i,ikb)-
447     .            (1.-zeps)*ztenh(i,ikb))+RLVTT*zqumqe
448            zdh = RG * MAX(zdh,1.0E5*zdqmin)
449            IF (zdhpbl(i).GT.0..AND.ldcum(i))zmfub1(i)=zdhpbl(i)/zdh
450         ENDIF
451         IF ( .NOT.((ktype(i).EQ.1.OR.ktype(i).EQ.2).AND.
452     .              ABS(zmfub1(i)-zmfub(i)).LT.0.2*zmfub(i)) )
453     .      zmfub1(i) = zmfub(i)
454      ENDIF
455      ENDDO
456      DO k = 1, klev
457      DO i = 1, klon
458      IF (lddraf(i)) THEN
459         zfac = zmfub1(i)/MAX(zmfub(i),1.E-10)
460         pmfd(i,k) = pmfd(i,k)*zfac
461         zmfds(i,k) = zmfds(i,k)*zfac
462         zmfdq(i,k) = zmfdq(i,k)*zfac
463         zdmfdp(i,k) = zdmfdp(i,k)*zfac
464         pen_d(i,k) = pen_d(i,k)*zfac
465         pde_d(i,k) = pde_d(i,k)*zfac
466      ENDIF
467      ENDDO
468      ENDDO
469      DO i = 1, klon
470         IF (lddraf(i)) zmfub(i)=zmfub1(i)
471      ENDDO
472c
473      ENDIF   ! fin de test sur lmfdd
474c
475c-----------------------------------------------------------------------
476c calculer de nouveau le panache ascendant
477c-----------------------------------------------------------------------
478      CALL flxasc(pdtime,ztenh, zqenh, pten, pqen, pqsen,
479     .     pgeo, zgeoh, pap, paph, pqte, pvervel,
480     .     ldland, ldcum, ktype, ilab,
481     .     ptu, pqu, plu, pmfu, zmfub, zentr,
482     .     zmfus, zmfuq, zmful, plude, zdmfup,
483     .     kcbot, kctop, ictop0, kcum, pen_u, pde_u)
484c
485c-----------------------------------------------------------------------
486c determiner les flux convectifs en forme finale, ainsi que
487c la quantite des precipitations
488c-----------------------------------------------------------------------
489      CALL flxflux(pdtime, pqen, pqsen, ztenh, zqenh, pap, paph,
490     .     ldland, zgeoh, kcbot, kctop, lddraf, kdtop, ktype, ldcum,
491     .     pmfu, pmfd, zmfus, zmfds, zmfuq, zmfdq, zmful, plude,
492     .     zdmfup, zdmfdp, pten, prsfc, pssfc, zdpmel, itopm2,
493     .     pmflxr, pmflxs)
494c
495c----------------------------------------------------------------------
496c calculer les tendances pour T et Q
497c----------------------------------------------------------------------
498      CALL flxdtdq(pdtime, itopm2, paph, ldcum, pten,
499     e     zmfus, zmfds, zmfuq, zmfdq, zmful, zdmfup, zdmfdp, zdpmel,
500     s     dt_con,dq_con)
501c
502 1000 CONTINUE
503      RETURN
504      END SUBROUTINE flxmain
505      SUBROUTINE flxini(pten, pqen, pqsen, pgeo, paph, pgeoh, ptenh,
506     .           pqenh, pqsenh, ptu, pqu, ptd, pqd, pmfd, pmfds, pmfdq,
507     .           pdmfdp, pmfu, pmfus, pmfuq, pdmfup, pdpmel, plu, plude,
508     .           klab,pen_u, pde_u, pen_d, pde_d)
509      USE dimphy
510      IMPLICIT none
511C----------------------------------------------------------------------
512C THIS ROUTINE INTERPOLATES LARGE-SCALE FIELDS OF T,Q ETC.
513C TO HALF LEVELS (I.E. GRID FOR MASSFLUX SCHEME),
514C AND INITIALIZES VALUES FOR UPDRAFTS
515C----------------------------------------------------------------------
516c     yminclude "dimensions.h"
517c     yminclude "dimphy.h"
518      include "YOMCST.h"
519      include "YOETHF.h"
520C
521      REAL pten(klon,klev)   ! temperature (environnement)
522      REAL pqen(klon,klev)   ! humidite (environnement)
523      REAL pqsen(klon,klev)  ! humidite saturante (environnement)
524      REAL pgeo(klon,klev)   ! geopotentiel (g * metre)
525      REAL pgeoh(klon,klev)  ! geopotentiel aux demi-niveaux
526      REAL paph(klon,klev+1) ! pression aux demi-niveaux
527      REAL ptenh(klon,klev)  ! temperature aux demi-niveaux
528      REAL pqenh(klon,klev)  ! humidite aux demi-niveaux
529      REAL pqsenh(klon,klev) ! humidite saturante aux demi-niveaux
530C
531      REAL ptu(klon,klev)    ! temperature du panache ascendant (p-a)
532      REAL pqu(klon,klev)    ! humidite du p-a
533      REAL plu(klon,klev)    ! eau liquide du p-a
534      REAL pmfu(klon,klev)   ! flux de masse du p-a
535      REAL pmfus(klon,klev)  ! flux de l'energie seche dans le p-a
536      REAL pmfuq(klon,klev)  ! flux de l'humidite dans le p-a
537      REAL pdmfup(klon,klev) ! quantite de l'eau precipitee dans p-a
538      REAL plude(klon,klev)  ! quantite de l'eau liquide jetee du
539c                              p-a a l'environnement
540      REAL pdpmel(klon,klev) ! quantite de neige fondue
541c
542      REAL ptd(klon,klev)    ! temperature du panache descendant (p-d)
543      REAL pqd(klon,klev)    ! humidite du p-d
544      REAL pmfd(klon,klev)   ! flux de masse du p-d
545      REAL pmfds(klon,klev)  ! flux de l'energie seche dans le p-d
546      REAL pmfdq(klon,klev)  ! flux de l'humidite dans le p-d
547      REAL pdmfdp(klon,klev) ! quantite de precipitation dans p-d
548c
549      REAL pen_u(klon,klev) ! quantite de masse entrainee pour p-a
550      REAL pde_u(klon,klev) ! quantite de masse detrainee pour p-a
551      REAL pen_d(klon,klev) ! quantite de masse entrainee pour p-d
552      REAL pde_d(klon,klev) ! quantite de masse detrainee pour p-d
553C
554      INTEGER  klab(klon,klev)
555      LOGICAL  llflag(klon)
556      INTEGER k, i, icall
557      REAL zzs
558C----------------------------------------------------------------------
559C SPECIFY LARGE SCALE PARAMETERS AT HALF LEVELS
560C ADJUST TEMPERATURE FIELDS IF STATICLY UNSTABLE
561C----------------------------------------------------------------------
562      DO 130 k = 2, klev
563c
564      DO i = 1, klon
565         pgeoh(i,k)=pgeo(i,k)+(pgeo(i,k-1)-pgeo(i,k))*0.5
566         ptenh(i,k)=(MAX(RCPD*pten(i,k-1)+pgeo(i,k-1),
567     .             RCPD*pten(i,k)+pgeo(i,k))-pgeoh(i,k))/RCPD
568         pqsenh(i,k)=pqsen(i,k-1)
569         llflag(i)=.TRUE.
570      ENDDO
571c
572      icall=0
573      CALL flxadjtq(paph(1,k),ptenh(1,k),pqsenh(1,k),llflag,icall)
574c
575      DO i = 1, klon
576         pqenh(i,k)=MIN(pqen(i,k-1),pqsen(i,k-1))
577     .               +(pqsenh(i,k)-pqsen(i,k-1))
578         pqenh(i,k)=MAX(pqenh(i,k),0.)
579      ENDDO
580c
581  130 CONTINUE
582C
583      DO 140 i = 1, klon
584         ptenh(i,klev)=(RCPD*pten(i,klev)+pgeo(i,klev)-
585     1                   pgeoh(i,klev))/RCPD
586         pqenh(i,klev)=pqen(i,klev)
587         ptenh(i,1)=pten(i,1)
588         pqenh(i,1)=pqen(i,1)
589         pgeoh(i,1)=pgeo(i,1)
590  140 CONTINUE
591c
592      DO 160 k = klev-1, 2, -1
593      DO 150 i = 1, klon
594         zzs = MAX(RCPD*ptenh(i,k)+pgeoh(i,k),
595     .             RCPD*ptenh(i,k+1)+pgeoh(i,k+1))
596         ptenh(i,k) = (zzs-pgeoh(i,k))/RCPD
597  150 CONTINUE
598  160 CONTINUE
599C
600C-----------------------------------------------------------------------
601C INITIALIZE VALUES FOR UPDRAFTS AND DOWNDRAFTS
602C-----------------------------------------------------------------------
603      DO k = 1, klev
604      DO i = 1, klon
605         ptu(i,k) = ptenh(i,k)
606         pqu(i,k) = pqenh(i,k)
607         plu(i,k) = 0.
608         pmfu(i,k) = 0.
609         pmfus(i,k) = 0.
610         pmfuq(i,k) = 0.
611         pdmfup(i,k) = 0.
612         pdpmel(i,k) = 0.
613         plude(i,k) = 0.
614c
615         klab(i,k) = 0
616c
617         ptd(i,k) = ptenh(i,k)
618         pqd(i,k) = pqenh(i,k)
619         pmfd(i,k) = 0.0
620         pmfds(i,k) = 0.0
621         pmfdq(i,k) = 0.0
622         pdmfdp(i,k) = 0.0
623c
624         pen_u(i,k) = 0.0
625         pde_u(i,k) = 0.0
626         pen_d(i,k) = 0.0
627         pde_d(i,k) = 0.0
628      ENDDO
629      ENDDO
630C
631      RETURN
632      END SUBROUTINE flxini
633      SUBROUTINE flxbase(ptenh, pqenh, pgeoh, paph,
634     *     ptu, pqu, plu, ldcum, kcbot, klab)
635      USE dimphy
636      IMPLICIT none
637C----------------------------------------------------------------------
638C THIS ROUTINE CALCULATES CLOUD BASE VALUES (T AND Q)
639C
640C INPUT ARE ENVIRONM. VALUES OF T,Q,P,PHI AT HALF LEVELS.
641C IT RETURNS CLOUD BASE VALUES AND FLAGS AS FOLLOWS;
642C   klab=1 FOR SUBCLOUD LEVELS
643C   klab=2 FOR CONDENSATION LEVEL
644C
645C LIFT SURFACE AIR DRY-ADIABATICALLY TO CLOUD BASE
646C (NON ENTRAINING PLUME,I.E.CONSTANT MASSFLUX)
647C----------------------------------------------------------------------
648c     yminclude "dimensions.h"
649c     yminclude "dimphy.h"
650      include "YOMCST.h"
651      include "YOETHF.h"
652C       ----------------------------------------------------------------
653      REAL ptenh(klon,klev), pqenh(klon,klev)
654      REAL pgeoh(klon,klev), paph(klon,klev+1)
655C
656      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
657      INTEGER  klab(klon,klev), kcbot(klon)
658C
659      LOGICAL llflag(klon), ldcum(klon)
660      INTEGER i, k, icall, is
661      REAL zbuo, zqold(klon)
662C----------------------------------------------------------------------
663C INITIALIZE VALUES AT LIFTING LEVEL
664C----------------------------------------------------------------------
665      DO i = 1, klon
666         klab(i,klev)=1
667         kcbot(i)=klev-1
668         ldcum(i)=.FALSE.
669      ENDDO
670C----------------------------------------------------------------------
671C DO ASCENT IN SUBCLOUD LAYER,
672C CHECK FOR EXISTENCE OF CONDENSATION LEVEL,
673C ADJUST T,Q AND L ACCORDINGLY
674C CHECK FOR BUOYANCY AND SET FLAGS
675C----------------------------------------------------------------------
676      DO 290 k = klev-1, 2, -1
677c
678      is = 0
679      DO i = 1, klon
680         IF (klab(i,k+1).EQ.1) is = is + 1
681         llflag(i) = .FALSE.
682         IF (klab(i,k+1).EQ.1) llflag(i) = .TRUE.
683      ENDDO
684      IF (is.EQ.0) GOTO 290
685c
686      DO i = 1, klon
687      IF(llflag(i)) THEN
688         pqu(i,k) = pqu(i,k+1)
689         ptu(i,k) = ptu(i,k+1)+(pgeoh(i,k+1)-pgeoh(i,k))/RCPD
690         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
691     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))+0.5
692         IF (zbuo.GT.0.) klab(i,k) = 1
693         zqold(i) = pqu(i,k)
694      ENDIF
695      ENDDO
696c
697      icall=1
698      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
699c
700      DO i = 1, klon
701      IF (llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN
702         klab(i,k) = 2
703         plu(i,k) = plu(i,k) + zqold(i)-pqu(i,k)
704         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
705     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))+0.5
706         IF (zbuo.GT.0.) kcbot(i) = k
707         IF (zbuo.GT.0.) ldcum(i) = .TRUE.
708      ENDIF
709      ENDDO
710c
711  290 CONTINUE
712c
713      RETURN
714      END SUBROUTINE flxbase
715      SUBROUTINE flxasc(pdtime, ptenh, pqenh, pten, pqen, pqsen,
716     .     pgeo, pgeoh, pap, paph, pqte, pvervel,
717     .     ldland, ldcum, ktype, klab, ptu, pqu, plu,
718     .     pmfu, pmfub, pentr, pmfus, pmfuq,
719     .     pmful, plude, pdmfup, kcbot, kctop, kctop0, kcum,
720     .     pen_u, pde_u)
721      USE dimphy
722      IMPLICIT none
723C----------------------------------------------------------------------
724C THIS ROUTINE DOES THE CALCULATIONS FOR CLOUD ASCENTS
725C FOR CUMULUS PARAMETERIZATION
726C----------------------------------------------------------------------
727c     yminclude "dimensions.h"
728c     yminclude "dimphy.h"
729      include "YOMCST.h"
730      include "YOETHF.h"
731      include "YOECUMF.h"
732C
733      REAL pdtime
734      REAL pten(klon,klev), ptenh(klon,klev)
735      REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)
736      REAL pgeo(klon,klev), pgeoh(klon,klev)
737      REAL pap(klon,klev), paph(klon,klev+1)
738      REAL pqte(klon,klev)
739      REAL pvervel(klon,klev) ! vitesse verticale en Pa/s
740C
741      REAL pmfub(klon), pentr(klon)
742      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
743      REAL plude(klon,klev)
744      REAL pmfu(klon,klev), pmfus(klon,klev)
745      REAL pmfuq(klon,klev), pmful(klon,klev)
746      REAL pdmfup(klon,klev)
747      INTEGER ktype(klon), klab(klon,klev), kcbot(klon), kctop(klon)
748      INTEGER kctop0(klon)
749      LOGICAL ldland(klon), ldcum(klon)
750C
751      REAL pen_u(klon,klev), pde_u(klon,klev)
752      REAL zqold(klon)
753      REAL zdland(klon)
754      LOGICAL llflag(klon)
755      INTEGER k, i, is, icall, kcum
756      REAL ztglace, zdphi, zqeen, zseen, zscde, zqude
757      REAL zmfusk, zmfuqk, zmfulk, zbuo, zdnoprc, zprcon, zlnew
758c
759      REAL zpbot(klon), zptop(klon), zrho(klon)
760      REAL zdprho, zentr, zpmid, zmftest, zmfmax
761      LOGICAL llo1, llo2
762c
763      REAL zwmax(klon), zzzmb
764      INTEGER klwmin(klon) ! level of maximum vertical velocity
765      real fact
766C----------------------------------------------------------------------
767      ztglace = RTT - 13.
768c
769c Chercher le niveau ou la vitesse verticale est maximale:
770      DO i = 1, klon
771         klwmin(i) = klev
772         zwmax(i) = 0.0
773      ENDDO
774      DO k = klev, 3, -1
775      DO i = 1, klon
776      IF (pvervel(i,k).LT.zwmax(i)) THEN
777         zwmax(i) = pvervel(i,k)
778         klwmin(i) = k
779      ENDIF
780      ENDDO
781      ENDDO
782C----------------------------------------------------------------------
783C SET DEFAULT VALUES
784C----------------------------------------------------------------------
785      DO i = 1, klon
786         IF (.NOT.ldcum(i)) ktype(i)=0
787      ENDDO
788c
789      DO k=1,klev
790      DO i = 1, klon
791         plu(i,k)=0.
792         pmfu(i,k)=0.
793         pmfus(i,k)=0.
794         pmfuq(i,k)=0.
795         pmful(i,k)=0.
796         plude(i,k)=0.
797         pdmfup(i,k)=0.
798         IF(.NOT.ldcum(i).OR.ktype(i).EQ.3) klab(i,k)=0
799         IF(.NOT.ldcum(i).AND.paph(i,k).LT.4.E4) kctop0(i)=k
800      ENDDO
801      ENDDO
802c
803      DO i = 1, klon
804      IF (ldland(i)) THEN
805         zdland(i)=3.0E4
806         zdphi=pgeoh(i,kctop0(i))-pgeoh(i,kcbot(i))
807         IF (ptu(i,kctop0(i)).GE.ztglace) zdland(i)=zdphi
808         zdland(i)=MAX(3.0E4,zdland(i))
809         zdland(i)=MIN(5.0E4,zdland(i))
810      ENDIF
811      ENDDO
812C
813C Initialiser les valeurs au niveau d'ascendance
814C
815      DO i = 1, klon
816         kctop(i) = klev-1
817         IF (.NOT.ldcum(i)) THEN
818            kcbot(i) = klev-1
819            pmfub(i) = 0.
820            pqu(i,klev) = 0.
821         ENDIF
822         pmfu(i,klev) = pmfub(i)
823         pmfus(i,klev) = pmfub(i)*(RCPD*ptu(i,klev)+pgeoh(i,klev))
824         pmfuq(i,klev) = pmfub(i)*pqu(i,klev)
825      ENDDO
826c
827      DO i = 1, klon
828         ldcum(i) = .FALSE.
829      ENDDO
830C----------------------------------------------------------------------
831C  DO ASCENT: SUBCLOUD LAYER (klab=1) ,CLOUDS (klab=2)
832C  BY DOING FIRST DRY-ADIABATIC ASCENT AND THEN
833C  BY ADJUSTING T,Q AND L ACCORDINGLY IN *flxadjtq*,
834C  THEN CHECK FOR BUOYANCY AND SET FLAGS ACCORDINGLY
835C----------------------------------------------------------------------
836      DO 480 k = klev-1,3,-1
837c
838      IF (LMFMID .AND. k.LT.klev-1) THEN
839         DO i = 1, klon
840         IF (.NOT.ldcum(i) .AND. klab(i,k+1).EQ.0 .AND.
841     .       pqen(i,k).GT.0.9*pqsen(i,k)
842     $           .and. pap(i, k) / paph(i, klev + 1)> 0.4) THEN
843            ptu(i,k+1) = pten(i,k) +(pgeo(i,k)-pgeoh(i,k+1))/RCPD
844            pqu(i,k+1) = pqen(i,k)
845            plu(i,k+1) = 0.0
846            zzzmb = MAX(CMFCMIN, -pvervel(i,k)/RG)
847            zmfmax = (paph(i,k)-paph(i,k-1))/(RG*pdtime)
848            pmfub(i) = MIN(zzzmb,zmfmax)
849            pmfu(i,k+1) = pmfub(i)
850            pmfus(i,k+1) = pmfub(i)*(RCPD*ptu(i,k+1)+pgeoh(i,k+1))
851            pmfuq(i,k+1) = pmfub(i)*pqu(i,k+1)
852            pmful(i,k+1) = 0.0
853            pdmfup(i,k+1) = 0.0
854            kcbot(i) = k
855            klab(i,k+1) = 1
856            ktype(i) = 3
857            pentr(i) = ENTRMID
858         ENDIF
859         ENDDO
860      ENDIF
861c
862      is = 0
863      DO i = 1, klon
864         is = is + klab(i,k+1)
865         IF (klab(i,k+1) .EQ. 0) klab(i,k) = 0
866         llflag(i) = .FALSE.
867         IF (klab(i,k+1) .GT. 0) llflag(i) = .TRUE.
868      ENDDO
869      IF (is .EQ. 0) GOTO 480
870c
871c calculer le taux d'entrainement et de detrainement
872c
873      DO i = 1, klon
874         pen_u(i,k) = 0.0
875         pde_u(i,k) = 0.0
876         zrho(i)=paph(i,k+1)/(RD*ptenh(i,k+1))
877         zpbot(i)=paph(i,kcbot(i))
878         zptop(i)=paph(i,kctop0(i))
879      ENDDO
880c
881      DO 125 i = 1, klon
882      IF(ldcum(i)) THEN
883         zdprho=(paph(i,k+1)-paph(i,k))/(RG*zrho(i))
884         zentr=pentr(i)*pmfu(i,k+1)*zdprho
885         llo1=k.LT.kcbot(i)
886         IF(llo1) pde_u(i,k)=zentr
887         zpmid=0.5*(zpbot(i)+zptop(i))
888         llo2=llo1.AND.ktype(i).EQ.2.AND.
889     .        (zpbot(i)-paph(i,k).LT.0.2E5.OR.
890     .         paph(i,k).GT.zpmid)
891         IF(llo2) pen_u(i,k)=zentr
892         llo2=llo1.AND.(ktype(i).EQ.1.OR.ktype(i).EQ.3).AND.
893     .        (k.GE.MAX(klwmin(i),kctop0(i)+2).OR.pap(i,k).GT.zpmid)
894         IF(llo2) pen_u(i,k)=zentr
895         llo1=pen_u(i,k).GT.0..AND.(ktype(i).EQ.1.OR.ktype(i).EQ.2)
896         IF(llo1) THEN
897            fact=1.+3.*(1.-MIN(1.,(zpbot(i)-pap(i,k))/1.5E4))
898            zentr=zentr*fact
899            pen_u(i,k)=pen_u(i,k)*fact
900            pde_u(i,k)=pde_u(i,k)*fact
901         ENDIF
902         IF(llo2.AND.pqenh(i,k+1).GT.1.E-5)
903     .   pen_u(i,k)=zentr+MAX(pqte(i,k),0.)/pqenh(i,k+1)*
904     .              zrho(i)*zdprho
905      ENDIF
906  125 CONTINUE
907c
908C----------------------------------------------------------------------
909c DO ADIABATIC ASCENT FOR ENTRAINING/DETRAINING PLUME
910C----------------------------------------------------------------------
911c
912      DO 420 i = 1, klon
913      IF (llflag(i)) THEN
914         IF (k.LT.kcbot(i)) THEN
915            zmftest = pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)
916            zmfmax = MIN(zmftest,(paph(i,k)-paph(i,k-1))/(RG*pdtime))
917            pen_u(i,k)=MAX(pen_u(i,k)-MAX(0.0,zmftest-zmfmax),0.0)
918         ENDIF
919         pde_u(i,k)=MIN(pde_u(i,k),0.75*pmfu(i,k+1))
920c calculer le flux de masse du niveau k a partir de celui du k+1
921         pmfu(i,k)=pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)
922c calculer les valeurs Su, Qu et l du niveau k dans le panache montant
923         zqeen=pqenh(i,k+1)*pen_u(i,k)
924         zseen=(RCPD*ptenh(i,k+1)+pgeoh(i,k+1))*pen_u(i,k)
925         zscde=(RCPD*ptu(i,k+1)+pgeoh(i,k+1))*pde_u(i,k)
926         zqude=pqu(i,k+1)*pde_u(i,k)
927         plude(i,k)=plu(i,k+1)*pde_u(i,k)
928         zmfusk=pmfus(i,k+1)+zseen-zscde
929         zmfuqk=pmfuq(i,k+1)+zqeen-zqude
930         zmfulk=pmful(i,k+1)    -plude(i,k)
931         plu(i,k)=zmfulk*(1./MAX(CMFCMIN,pmfu(i,k)))
932         pqu(i,k)=zmfuqk*(1./MAX(CMFCMIN,pmfu(i,k)))
933         ptu(i,k)=(zmfusk*(1./MAX(CMFCMIN,pmfu(i,k)))-
934     1               pgeoh(i,k))/RCPD
935         ptu(i,k)=MAX(100.,ptu(i,k))
936         ptu(i,k)=MIN(400.,ptu(i,k))
937         zqold(i)=pqu(i,k)
938      ELSE
939         zqold(i)=0.0
940      ENDIF
941  420 CONTINUE
942c
943C----------------------------------------------------------------------
944c DO CORRECTIONS FOR MOIST ASCENT BY ADJUSTING T,Q AND L
945C----------------------------------------------------------------------
946c
947      icall = 1
948      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
949C
950      DO 440 i = 1, klon
951      IF(llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN
952         klab(i,k) = 2
953         plu(i,k) = plu(i,k)+zqold(i)-pqu(i,k)
954         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
955     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))
956         IF (klab(i,k+1).EQ.1) zbuo=zbuo+0.5
957         IF (zbuo.GT.0..AND.pmfu(i,k).GE.0.1*pmfub(i)) THEN
958            kctop(i) = k
959            ldcum(i) = .TRUE.
960            zdnoprc = 1.5E4
961            IF (ldland(i)) zdnoprc = zdland(i)
962            zprcon = CPRCON
963            IF ((zpbot(i)-paph(i,k)).LT.zdnoprc) zprcon = 0.0
964            zlnew=plu(i,k)/(1.+zprcon*(pgeoh(i,k)-pgeoh(i,k+1)))
965            pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))
966            plu(i,k)=zlnew
967         ELSE
968            klab(i,k)=0
969            pmfu(i,k)=0.
970         ENDIF
971      ENDIF
972  440 CONTINUE
973      DO 455 i = 1, klon
974      IF (llflag(i)) THEN
975         pmful(i,k)=plu(i,k)*pmfu(i,k)
976         pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)
977         pmfuq(i,k)=pqu(i,k)*pmfu(i,k)
978      ENDIF
979  455 CONTINUE
980C
981  480 CONTINUE
982C----------------------------------------------------------------------
983C DETERMINE CONVECTIVE FLUXES ABOVE NON-BUOYANCY LEVEL
984C    (NOTE: CLOUD VARIABLES LIKE T,Q AND L ARE NOT
985C           AFFECTED BY DETRAINMENT AND ARE ALREADY KNOWN
986C           FROM PREVIOUS CALCULATIONS ABOVE)
987C----------------------------------------------------------------------
988      DO i = 1, klon
989         IF (kctop(i).EQ.klev-1) ldcum(i) = .FALSE.
990         kcbot(i) = MAX(kcbot(i),kctop(i))
991      ENDDO
992c
993      ldcum(1)=ldcum(1)
994c
995      is = 0
996      DO i = 1, klon
997         if (ldcum(i)) is = is + 1
998      ENDDO
999      kcum = is
1000      IF (is.EQ.0) GOTO 800
1001c
1002      DO 530 i = 1, klon
1003      IF (ldcum(i)) THEN
1004         k=kctop(i)-1
1005         pde_u(i,k)=(1.-CMFCTOP)*pmfu(i,k+1)
1006         plude(i,k)=pde_u(i,k)*plu(i,k+1)
1007         pmfu(i,k)=pmfu(i,k+1)-pde_u(i,k)
1008         zlnew=plu(i,k)
1009         pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))
1010         plu(i,k)=zlnew
1011         pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)
1012         pmfuq(i,k)=pqu(i,k)*pmfu(i,k)
1013         pmful(i,k)=plu(i,k)*pmfu(i,k)
1014         plude(i,k-1)=pmful(i,k)
1015      ENDIF
1016  530 CONTINUE
1017C
1018  800 CONTINUE
1019      RETURN
1020      END SUBROUTINE flxasc
1021      SUBROUTINE flxflux(pdtime, pqen, pqsen, ptenh, pqenh, pap
1022     .  ,  paph, ldland, pgeoh, kcbot, kctop, lddraf, kdtop
1023     .  ,  ktype, ldcum, pmfu, pmfd, pmfus, pmfds
1024     .  ,  pmfuq, pmfdq, pmful, plude, pdmfup, pdmfdp
1025     .  ,  pten, prfl, psfl, pdpmel, ktopm2
1026     .  ,  pmflxr, pmflxs)
1027      USE dimphy
1028      IMPLICIT none
1029C----------------------------------------------------------------------
1030C THIS ROUTINE DOES THE FINAL CALCULATION OF CONVECTIVE
1031C FLUXES IN THE CLOUD LAYER AND IN THE SUBCLOUD LAYER
1032C----------------------------------------------------------------------
1033c     yminclude "dimensions.h"
1034c     yminclude "dimphy.h"
1035      include "YOMCST.h"
1036      include "YOETHF.h"
1037      include "YOECUMF.h"
1038      include "iniprint.h"
1039C
1040      REAL cevapcu(klon,klev)
1041C     -----------------------------------------------------------------
1042      REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)
1043      REAL pten(klon,klev), ptenh(klon,klev)
1044      REAL paph(klon,klev+1), pgeoh(klon,klev)
1045c
1046      REAL pap(klon,klev)
1047      REAL ztmsmlt, zdelta, zqsat
1048C
1049      REAL pmfu(klon,klev), pmfus(klon,klev)
1050      REAL pmfd(klon,klev), pmfds(klon,klev)
1051      REAL pmfuq(klon,klev), pmful(klon,klev)
1052      REAL pmfdq(klon,klev)
1053      REAL plude(klon,klev)
1054      REAL pdmfup(klon,klev), pdpmel(klon,klev)
1055cjq The variable maxpdmfdp(klon) has been introduced by Olivier Boucher
1056cjq 14/11/00 to fix the problem with the negative precipitation.     
1057      REAL pdmfdp(klon,klev), maxpdmfdp(klon,klev)
1058      REAL prfl(klon), psfl(klon)
1059      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
1060      INTEGER  kcbot(klon), kctop(klon), ktype(klon)
1061      LOGICAL  ldland(klon), ldcum(klon)
1062      INTEGER k, kp, i
1063      REAL zcons1, zcons2, zcucov, ztmelp2
1064      REAL pdtime, zdp, zzp, zfac, zsnmlt, zrfl, zrnew
1065      REAL zrmin, zrfln, zdrfl
1066      REAL zpds, zpdr, zdenom
1067      INTEGER ktopm2, itop, ikb
1068c
1069      LOGICAL lddraf(klon)
1070      INTEGER kdtop(klon)
1071c
1072      include "FCTTRE.h"
1073c
1074      DO 101 k=1,klev
1075        DO i=1,klon
1076      CEVAPCU(i,k)=1.93E-6*261.*SQRT(1.E3/(38.3*0.293)
1077     1 *SQRT(0.5*(paph(i,k)+paph(i,k+1))/paph(i,klev+1)) ) * 0.5/RG
1078        ENDDO
1079 101  CONTINUE
1080c
1081c SPECIFY CONSTANTS
1082c
1083      zcons1 = RCPD/(RLMLT*RG*pdtime)
1084      zcons2 = 1./(RG*pdtime)
1085      zcucov = 0.05
1086      ztmelp2 = RTT + 2.
1087c
1088c DETERMINE FINAL CONVECTIVE FLUXES
1089c
1090      itop=klev
1091      DO 110 i = 1, klon
1092         itop=MIN(itop,kctop(i))
1093         IF (.NOT.ldcum(i) .OR. kdtop(i).LT.kctop(i)) lddraf(i)=.FALSE.
1094         IF(.NOT.ldcum(i)) ktype(i)=0
1095  110 CONTINUE
1096c
1097      ktopm2=itop-2
1098      DO 120 k=ktopm2,klev
1099      DO 115 i = 1, klon
1100      IF(ldcum(i).AND.k.GE.kctop(i)-1) THEN
1101         pmfus(i,k)=pmfus(i,k)-pmfu(i,k)*
1102     .                (RCPD*ptenh(i,k)+pgeoh(i,k))
1103         pmfuq(i,k)=pmfuq(i,k)-pmfu(i,k)*pqenh(i,k)
1104         zdp = 1.5E4
1105         IF ( ldland(i) ) zdp = 3.E4
1106c
1107c        l'eau liquide detrainee est precipitee quand certaines
1108c        conditions sont reunies (sinon, elle est consideree
1109c        evaporee dans l'environnement)
1110c
1111         IF(paph(i,kcbot(i))-paph(i,kctop(i)).GE.zdp.AND.
1112     .      pqen(i,k-1).GT.0.8*pqsen(i,k-1))
1113     .      pdmfup(i,k-1)=pdmfup(i,k-1)+plude(i,k-1)
1114c
1115         IF(lddraf(i).AND.k.GE.kdtop(i)) THEN
1116            pmfds(i,k)=pmfds(i,k)-pmfd(i,k)*
1117     .                   (RCPD*ptenh(i,k)+pgeoh(i,k))
1118            pmfdq(i,k)=pmfdq(i,k)-pmfd(i,k)*pqenh(i,k)
1119         ELSE
1120            pmfd(i,k)=0.
1121            pmfds(i,k)=0.
1122            pmfdq(i,k)=0.
1123            pdmfdp(i,k-1)=0.
1124         END IF
1125      ELSE
1126         pmfu(i,k)=0.
1127         pmfus(i,k)=0.
1128         pmfuq(i,k)=0.
1129         pmful(i,k)=0.
1130         pdmfup(i,k-1)=0.
1131         plude(i,k-1)=0.
1132         pmfd(i,k)=0.
1133         pmfds(i,k)=0.
1134         pmfdq(i,k)=0.
1135         pdmfdp(i,k-1)=0.
1136      ENDIF
1137  115 CONTINUE
1138  120 CONTINUE
1139c
1140      DO 130 k=ktopm2,klev
1141      DO 125 i = 1, klon
1142      IF(ldcum(i).AND.k.GT.kcbot(i)) THEN
1143         ikb=kcbot(i)
1144         zzp=((paph(i,klev+1)-paph(i,k))/
1145     .        (paph(i,klev+1)-paph(i,ikb)))
1146         IF (ktype(i).EQ.3) zzp = zzp**2
1147         pmfu(i,k)=pmfu(i,ikb)*zzp
1148         pmfus(i,k)=pmfus(i,ikb)*zzp
1149         pmfuq(i,k)=pmfuq(i,ikb)*zzp
1150         pmful(i,k)=pmful(i,ikb)*zzp
1151      ENDIF
1152  125 CONTINUE
1153  130 CONTINUE
1154c
1155c CALCULATE RAIN/SNOW FALL RATES
1156c CALCULATE MELTING OF SNOW
1157c CALCULATE EVAPORATION OF PRECIP
1158c
1159      DO k = 1, klev+1
1160      DO i = 1, klon
1161         pmflxr(i,k) = 0.0
1162         pmflxs(i,k) = 0.0
1163      ENDDO
1164      ENDDO
1165      DO k = ktopm2, klev
1166      DO i = 1, klon
1167      IF (ldcum(i)) THEN
1168         IF (pmflxs(i,k).GT.0.0 .AND. pten(i,k).GT.ztmelp2) THEN
1169            zfac=zcons1*(paph(i,k+1)-paph(i,k))
1170            zsnmlt=MIN(pmflxs(i,k),zfac*(pten(i,k)-ztmelp2))
1171            pdpmel(i,k)=zsnmlt
1172            ztmsmlt=pten(i,k)-zsnmlt/zfac
1173            zdelta=MAX(0.,SIGN(1.,RTT-ztmsmlt))
1174            zqsat=R2ES*FOEEW(ztmsmlt, zdelta) / pap(i,k)
1175            zqsat=MIN(0.5,zqsat)
1176            zqsat=zqsat/(1.-RETV  *zqsat)
1177            pqsen(i,k) = zqsat
1178         ENDIF
1179         IF (pten(i,k).GT.RTT) THEN
1180         pmflxr(i,k+1)=pmflxr(i,k)+pdmfup(i,k)+pdmfdp(i,k)+pdpmel(i,k)
1181         pmflxs(i,k+1)=pmflxs(i,k)-pdpmel(i,k)
1182         ELSE
1183           pmflxs(i,k+1)=pmflxs(i,k)+pdmfup(i,k)+pdmfdp(i,k)
1184           pmflxr(i,k+1)=pmflxr(i,k)
1185         ENDIF
1186c        si la precipitation est negative, on ajuste le plux du
1187c        panache descendant pour eliminer la negativite
1188         IF ((pmflxr(i,k+1)+pmflxs(i,k+1)).LT.0.0) THEN
1189            pdmfdp(i,k) = -pmflxr(i,k)-pmflxs(i,k)-pdmfup(i,k)
1190            pmflxr(i,k+1) = 0.0
1191            pmflxs(i,k+1) = 0.0
1192            pdpmel(i,k) = 0.0
1193         ENDIF
1194      ENDIF
1195      ENDDO
1196      ENDDO
1197c
1198cjq The new variable is initialized here.
1199cjq It contains the humidity which is fed to the downdraft
1200cjq by evaporation of precipitation in the column below the base
1201cjq of convection.
1202cjq
1203cjq In the former version, this term has been subtracted from precip
1204cjq as well as the evaporation.
1205cjq     
1206      DO k = 1, klev
1207      DO i = 1, klon
1208         maxpdmfdp(i,k)=0.0
1209      ENDDO
1210      ENDDO
1211      DO k = 1, klev
1212       DO kp = k, klev
1213        DO i = 1, klon
1214         maxpdmfdp(i,k)=maxpdmfdp(i,k)+pdmfdp(i,kp)
1215        ENDDO
1216       ENDDO
1217      ENDDO
1218cjq End of initialization
1219c     
1220      DO k = ktopm2, klev
1221      DO i = 1, klon
1222      IF (ldcum(i) .AND. k.GE.kcbot(i)) THEN
1223         zrfl = pmflxr(i,k) + pmflxs(i,k)
1224         IF (zrfl.GT.1.0E-20) THEN
1225            zrnew=(MAX(0.,SQRT(zrfl/zcucov)-
1226     .            CEVAPCU(i,k)*(paph(i,k+1)-paph(i,k))*
1227     .            MAX(0.,pqsen(i,k)-pqen(i,k))))**2*zcucov
1228            zrmin=zrfl-zcucov*MAX(0.,0.8*pqsen(i,k)-pqen(i,k))
1229     .            *zcons2*(paph(i,k+1)-paph(i,k))
1230            zrnew=MAX(zrnew,zrmin)
1231            zrfln=MAX(zrnew,0.)
1232            zdrfl=MIN(0.,zrfln-zrfl)
1233cjq At least the amount of precipiation needed to feed the downdraft
1234cjq with humidity below the base of convection has to be left and can't
1235cjq be evaporated (surely the evaporation can't be positive):           
1236            zdrfl=MAX(zdrfl,
1237     .            MIN(-pmflxr(i,k)-pmflxs(i,k)-maxpdmfdp(i,k),0.0))
1238cjq End of insertion
1239c           
1240            zdenom=1.0/MAX(1.0E-20,pmflxr(i,k)+pmflxs(i,k))
1241            IF (pten(i,k).GT.RTT) THEN
1242               zpdr = pdmfdp(i,k)
1243               zpds = 0.0
1244            ELSE
1245               zpdr = 0.0
1246               zpds = pdmfdp(i,k)
1247            ENDIF
1248            pmflxr(i,k+1) = pmflxr(i,k) + zpdr + pdpmel(i,k)
1249     .                    + zdrfl*pmflxr(i,k)*zdenom
1250            pmflxs(i,k+1) = pmflxs(i,k) + zpds - pdpmel(i,k)
1251     .                    + zdrfl*pmflxs(i,k)*zdenom
1252            pdmfup(i,k) = pdmfup(i,k) + zdrfl
1253         ELSE
1254            pmflxr(i,k+1) = 0.0
1255            pmflxs(i,k+1) = 0.0
1256            pdmfdp(i,k) = 0.0
1257            pdpmel(i,k) = 0.0
1258         ENDIF         
1259         if (pmflxr(i,k) + pmflxs(i,k) < -1.e-26 .and. prt_level >= 1)
1260     .    write(*,*) 'precip. < 1e-16 ',pmflxr(i,k) + pmflxs(i,k)
1261      ENDIF
1262      ENDDO
1263      ENDDO
1264c
1265      DO 210 i = 1, klon
1266         prfl(i) = pmflxr(i,klev+1)
1267         psfl(i) = pmflxs(i,klev+1)
1268  210 CONTINUE
1269c
1270      RETURN
1271      END SUBROUTINE flxflux
1272      SUBROUTINE flxdtdq(pdtime, ktopm2, paph, ldcum, pten
1273     .  ,  pmfus, pmfds, pmfuq, pmfdq, pmful, pdmfup, pdmfdp
1274     .  ,  pdpmel, dt_con, dq_con)
1275      USE dimphy
1276      IMPLICIT none
1277c----------------------------------------------------------------------
1278c calculer les tendances T et Q
1279c----------------------------------------------------------------------
1280c     yminclude "dimensions.h"
1281c     yminclude "dimphy.h"
1282      include "YOMCST.h"
1283      include "YOETHF.h"
1284      include "YOECUMF.h"
1285C     -----------------------------------------------------------------
1286      LOGICAL  llo1
1287C
1288      REAL pten(klon,klev), paph(klon,klev+1)
1289      REAL pmfus(klon,klev), pmfuq(klon,klev), pmful(klon,klev)
1290      REAL pmfds(klon,klev), pmfdq(klon,klev)
1291      REAL pdmfup(klon,klev)
1292      REAL pdmfdp(klon,klev)
1293      REAL pdpmel(klon,klev)
1294      LOGICAL ldcum(klon)
1295      REAL dt_con(klon,klev), dq_con(klon,klev)
1296c
1297      INTEGER ktopm2
1298      REAL pdtime
1299c
1300      INTEGER i, k
1301      REAL zalv, zdtdt, zdqdt
1302c
1303      DO 210 k=ktopm2,klev-1
1304      DO 220 i = 1, klon
1305      IF (ldcum(i)) THEN
1306         llo1 = (pten(i,k)-RTT).GT.0.
1307         zalv = RLSTT
1308         IF (llo1) zalv = RLVTT
1309         zdtdt=RG/(paph(i,k+1)-paph(i,k))/RCPD
1310     .        *(pmfus(i,k+1)-pmfus(i,k)
1311     .         +pmfds(i,k+1)-pmfds(i,k)
1312     .          -RLMLT*pdpmel(i,k)
1313     .          -zalv*(pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1314     .         )
1315         dt_con(i,k)=zdtdt
1316         zdqdt=RG/(paph(i,k+1)-paph(i,k))
1317     .        *(pmfuq(i,k+1)-pmfuq(i,k)
1318     .         +pmfdq(i,k+1)-pmfdq(i,k)
1319     .          +pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1320         dq_con(i,k)=zdqdt
1321      ENDIF
1322  220 CONTINUE
1323  210 CONTINUE
1324C
1325      k = klev
1326      DO 230 i = 1, klon
1327      IF (ldcum(i)) THEN
1328         llo1 = (pten(i,k)-RTT).GT.0.
1329         zalv = RLSTT
1330         IF (llo1) zalv = RLVTT
1331         zdtdt=-RG/(paph(i,k+1)-paph(i,k))/RCPD
1332     .         *(pmfus(i,k)+pmfds(i,k)+RLMLT*pdpmel(i,k)
1333     .           -zalv*(pmful(i,k)+pdmfup(i,k)+pdmfdp(i,k)))
1334         dt_con(i,k)=zdtdt
1335         zdqdt=-RG/(paph(i,k+1)-paph(i,k))
1336     .            *(pmfuq(i,k)+pmfdq(i,k)+pmful(i,k)
1337     .             +pdmfup(i,k)+pdmfdp(i,k))
1338         dq_con(i,k)=zdqdt
1339      ENDIF
1340  230 CONTINUE
1341C
1342      RETURN
1343      END SUBROUTINE flxdtdq
1344      SUBROUTINE flxdlfs(ptenh, pqenh, pgeoh, paph, ptu, pqu,
1345     .     ldcum, kcbot, kctop, pmfub, prfl, ptd, pqd,
1346     .     pmfd, pmfds, pmfdq, pdmfdp, kdtop, lddraf)
1347      USE dimphy
1348      IMPLICIT none
1349C
1350C----------------------------------------------------------------------
1351C THIS ROUTINE CALCULATES LEVEL OF FREE SINKING FOR
1352C CUMULUS DOWNDRAFTS AND SPECIFIES T,Q,U AND V VALUES
1353C
1354C TO PRODUCE LFS-VALUES FOR CUMULUS DOWNDRAFTS
1355C FOR MASSFLUX CUMULUS PARAMETERIZATION
1356C
1357C INPUT ARE ENVIRONMENTAL VALUES OF T,Q,U,V,P,PHI
1358C AND UPDRAFT VALUES T,Q,U AND V AND ALSO
1359C CLOUD BASE MASSFLUX AND CU-PRECIPITATION RATE.
1360C IT RETURNS T,Q,U AND V VALUES AND MASSFLUX AT LFS.
1361C
1362C CHECK FOR NEGATIVE BUOYANCY OF AIR OF EQUAL PARTS OF
1363C MOIST ENVIRONMENTAL AIR AND CLOUD AIR.
1364C----------------------------------------------------------------------
1365c     yminclude "dimensions.h"
1366c     yminclude "dimphy.h"
1367      include "YOMCST.h"
1368      include "YOETHF.h"
1369      include "YOECUMF.h"
1370C
1371      REAL ptenh(klon,klev)
1372      REAL pqenh(klon,klev)
1373      REAL pgeoh(klon,klev), paph(klon,klev+1)
1374      REAL ptu(klon,klev), pqu(klon,klev)
1375      REAL pmfub(klon)
1376      REAL prfl(klon)
1377C
1378      REAL ptd(klon,klev), pqd(klon,klev)
1379      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1380      REAL pdmfdp(klon,klev)
1381      INTEGER  kcbot(klon), kctop(klon), kdtop(klon)
1382      LOGICAL  ldcum(klon), lddraf(klon)
1383C
1384      REAL ztenwb(klon,klev), zqenwb(klon,klev), zcond(klon)
1385      REAL zttest, zqtest, zbuo, zmftop
1386      LOGICAL  llo2(klon)
1387      INTEGER i, k, is, icall
1388C----------------------------------------------------------------------
1389      DO i= 1, klon
1390         lddraf(i)=.FALSE.
1391         kdtop(i)=klev+1
1392      ENDDO
1393C
1394C----------------------------------------------------------------------
1395C DETERMINE LEVEL OF FREE SINKING BY
1396C DOING A SCAN FROM TOP TO BASE OF CUMULUS CLOUDS
1397C
1398C FOR EVERY POINT AND PROCEED AS FOLLOWS:
1399C     (1) DETEMINE WET BULB ENVIRONMENTAL T AND Q
1400C     (2) DO MIXING WITH CUMULUS CLOUD AIR
1401C     (3) CHECK FOR NEGATIVE BUOYANCY
1402C
1403C THE ASSUMPTION IS THAT AIR OF DOWNDRAFTS IS MIXTURE
1404C OF 50% CLOUD AIR + 50% ENVIRONMENTAL AIR AT WET BULB
1405C TEMPERATURE (I.E. WHICH BECAME SATURATED DUE TO
1406C EVAPORATION OF RAIN AND CLOUD WATER)
1407C----------------------------------------------------------------------
1408C
1409      DO 290 k = 3, klev-3
1410C
1411      is=0
1412      DO 212 i= 1, klon
1413         ztenwb(i,k)=ptenh(i,k)
1414         zqenwb(i,k)=pqenh(i,k)
1415         llo2(i) = ldcum(i).AND.prfl(i).GT.0.
1416     .             .AND..NOT.lddraf(i)
1417     .             .AND.(k.LT.kcbot(i).AND.k.GT.kctop(i))
1418         IF ( llo2(i) ) is = is + 1
1419  212 CONTINUE
1420      IF(is.EQ.0) GO TO 290
1421C
1422      icall=2
1423      CALL flxadjtq(paph(1,k), ztenwb(1,k), zqenwb(1,k), llo2, icall)
1424C
1425C----------------------------------------------------------------------
1426C DO MIXING OF CUMULUS AND ENVIRONMENTAL AIR
1427C AND CHECK FOR NEGATIVE BUOYANCY.
1428C THEN SET VALUES FOR DOWNDRAFT AT LFS.
1429C----------------------------------------------------------------------
1430      DO 222 i= 1, klon
1431      IF (llo2(i)) THEN
1432         zttest=0.5*(ptu(i,k)+ztenwb(i,k))
1433         zqtest=0.5*(pqu(i,k)+zqenwb(i,k))
1434         zbuo=zttest*(1.+RETV*zqtest)-
1435     .        ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1436         zcond(i)=pqenh(i,k)-zqenwb(i,k)
1437         zmftop=-CMFDEPS*pmfub(i)
1438         IF (zbuo.LT.0..AND.prfl(i).GT.10.*zmftop*zcond(i)) THEN
1439            kdtop(i)=k
1440            lddraf(i)=.TRUE.
1441            ptd(i,k)=zttest
1442            pqd(i,k)=zqtest
1443            pmfd(i,k)=zmftop
1444            pmfds(i,k)=pmfd(i,k)*(RCPD*ptd(i,k)+pgeoh(i,k))
1445            pmfdq(i,k)=pmfd(i,k)*pqd(i,k)
1446            pdmfdp(i,k-1)=-0.5*pmfd(i,k)*zcond(i)
1447            prfl(i)=prfl(i)+pdmfdp(i,k-1)
1448         ENDIF
1449      ENDIF
1450  222 CONTINUE
1451c
1452  290 CONTINUE
1453C
1454      RETURN
1455      END SUBROUTINE flxdlfs
1456      SUBROUTINE flxddraf(ptenh, pqenh, pgeoh, paph, prfl,
1457     .           ptd, pqd, pmfd, pmfds, pmfdq, pdmfdp,
1458     .           lddraf, pen_d, pde_d)
1459      USE dimphy
1460      IMPLICIT none
1461C
1462C----------------------------------------------------------------------
1463C          THIS ROUTINE CALCULATES CUMULUS DOWNDRAFT DESCENT
1464C
1465C          TO PRODUCE THE VERTICAL PROFILES FOR CUMULUS DOWNDRAFTS
1466C          (I.E. T,Q,U AND V AND FLUXES)
1467C
1468C          INPUT IS T,Q,P,PHI,U,V AT HALF LEVELS.
1469C          IT RETURNS FLUXES OF S,Q AND EVAPORATION RATE
1470C          AND U,V AT LEVELS WHERE DOWNDRAFT OCCURS
1471C
1472C          CALCULATE MOIST DESCENT FOR ENTRAINING/DETRAINING PLUME BY
1473C          A) MOVING AIR DRY-ADIABATICALLY TO NEXT LEVEL BELOW AND
1474C          B) CORRECTING FOR EVAPORATION TO OBTAIN SATURATED STATE.
1475C
1476C----------------------------------------------------------------------
1477c     yminclude "dimensions.h"
1478c     yminclude "dimphy.h"
1479      include "YOMCST.h"
1480      include "YOETHF.h"
1481      include "YOECUMF.h"
1482C
1483      REAL ptenh(klon,klev), pqenh(klon,klev)
1484      REAL pgeoh(klon,klev), paph(klon,klev+1)
1485C
1486      REAL ptd(klon,klev), pqd(klon,klev)
1487      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1488      REAL pdmfdp(klon,klev)
1489      REAL prfl(klon)
1490      LOGICAL lddraf(klon)
1491C
1492      REAL pen_d(klon,klev), pde_d(klon,klev), zcond(klon)
1493      LOGICAL llo2(klon), llo1
1494      INTEGER i, k, is, icall, itopde
1495      REAL zentr, zseen, zqeen, zsdde, zqdde, zmfdsk, zmfdqk, zdmfdp
1496      REAL zbuo
1497C----------------------------------------------------------------------
1498C CALCULATE MOIST DESCENT FOR CUMULUS DOWNDRAFT BY
1499C       (A) CALCULATING ENTRAINMENT RATES, ASSUMING
1500C           LINEAR DECREASE OF MASSFLUX IN PBL
1501C       (B) DOING MOIST DESCENT - EVAPORATIVE COOLING
1502C           AND MOISTENING IS CALCULATED IN *flxadjtq*
1503C       (C) CHECKING FOR NEGATIVE BUOYANCY AND
1504C           SPECIFYING FINAL T,Q,U,V AND DOWNWARD FLUXES
1505C
1506      DO 180 k = 3, klev
1507c
1508      is = 0
1509      DO i = 1, klon
1510         llo2(i)=lddraf(i).AND.pmfd(i,k-1).LT.0.
1511         IF (llo2(i)) is = is + 1
1512      ENDDO
1513      IF (is.EQ.0) GOTO 180
1514c
1515      DO i = 1, klon
1516      IF (llo2(i)) THEN
1517         zentr = ENTRDD*pmfd(i,k-1)*RD*ptenh(i,k-1)/
1518     .           (RG*paph(i,k-1))*(paph(i,k)-paph(i,k-1))
1519         pen_d(i,k) = zentr
1520         pde_d(i,k) = zentr
1521      ENDIF
1522      ENDDO
1523c
1524      itopde = klev-2
1525      IF (k.GT.itopde) THEN
1526         DO i = 1, klon
1527         IF (llo2(i)) THEN
1528            pen_d(i,k)=0.
1529            pde_d(i,k)=pmfd(i,itopde)*
1530     .      (paph(i,k)-paph(i,k-1))/(paph(i,klev+1)-paph(i,itopde))
1531         ENDIF
1532         ENDDO
1533      ENDIF
1534C
1535      DO i = 1, klon
1536      IF (llo2(i)) THEN
1537         pmfd(i,k) = pmfd(i,k-1)+pen_d(i,k)-pde_d(i,k)
1538         zseen = (RCPD*ptenh(i,k-1)+pgeoh(i,k-1))*pen_d(i,k)
1539         zqeen = pqenh(i,k-1)*pen_d(i,k)
1540         zsdde = (RCPD*ptd(i,k-1)+pgeoh(i,k-1))*pde_d(i,k)
1541         zqdde = pqd(i,k-1)*pde_d(i,k)
1542         zmfdsk = pmfds(i,k-1)+zseen-zsdde
1543         zmfdqk = pmfdq(i,k-1)+zqeen-zqdde
1544         pqd(i,k) = zmfdqk*(1./MIN(-CMFCMIN,pmfd(i,k)))
1545         ptd(i,k) = (zmfdsk*(1./MIN(-CMFCMIN,pmfd(i,k)))-
1546     .               pgeoh(i,k))/RCPD
1547         ptd(i,k) = MIN(400.,ptd(i,k))
1548         ptd(i,k) = MAX(100.,ptd(i,k))
1549         zcond(i) = pqd(i,k)
1550      ENDIF
1551      ENDDO
1552C
1553      icall = 2
1554      CALL flxadjtq(paph(1,k), ptd(1,k), pqd(1,k), llo2, icall)
1555C
1556      DO i = 1, klon
1557      IF (llo2(i)) THEN
1558         zcond(i) = zcond(i)-pqd(i,k)
1559         zbuo = ptd(i,k)*(1.+RETV  *pqd(i,k))-
1560     .          ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1561         llo1 = zbuo.LT.0..AND.(prfl(i)-pmfd(i,k)*zcond(i).GT.0.)
1562         IF (.not.llo1) pmfd(i,k) = 0.0
1563         pmfds(i,k) = (RCPD*ptd(i,k)+pgeoh(i,k))*pmfd(i,k)
1564         pmfdq(i,k) = pqd(i,k)*pmfd(i,k)
1565         zdmfdp = -pmfd(i,k)*zcond(i)
1566         pdmfdp(i,k-1) = zdmfdp
1567         prfl(i) = prfl(i)+zdmfdp
1568      ENDIF
1569      ENDDO
1570c
1571  180 CONTINUE
1572      RETURN
1573      END SUBROUTINE flxddraf
1574      SUBROUTINE flxadjtq(pp, pt, pq, ldflag, kcall)
1575      USE dimphy
1576      IMPLICIT none
1577c======================================================================
1578c Objet: ajustement entre T et Q
1579c======================================================================
1580C NOTE: INPUT PARAMETER kcall DEFINES CALCULATION AS
1581C        kcall=0    ENV. T AND QS IN*CUINI*
1582C        kcall=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)
1583C        kcall=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1584C
1585c     yminclude "dimensions.h"
1586c     yminclude "dimphy.h"
1587      include "YOMCST.h"
1588C
1589      REAL pt(klon), pq(klon), pp(klon)
1590      LOGICAL ldflag(klon)
1591      INTEGER kcall
1592c
1593      REAL zcond(klon), zcond1
1594      REAL Z5alvcp, z5alscp, zalvdcp, zalsdcp
1595      REAL zdelta, zcvm5, zldcp, zqsat, zcor
1596      INTEGER is, i
1597      include "YOETHF.h"
1598      include "FCTTRE.h"
1599C
1600      z5alvcp = r5les*RLVTT/RCPD
1601      z5alscp = r5ies*RLSTT/RCPD
1602      zalvdcp = rlvtt/RCPD
1603      zalsdcp = rlstt/RCPD
1604C
1605
1606      DO i = 1, klon
1607         zcond(i) = 0.0
1608      ENDDO
1609
1610      DO 210 i =1, klon
1611      IF (ldflag(i)) THEN
1612         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1613         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1614         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1615         zqsat = R2ES*FOEEW(pt(i),zdelta) / pp(i)
1616         zqsat = MIN(0.5,zqsat)
1617         zcor = 1./(1.-RETV*zqsat)
1618         zqsat = zqsat*zcor
1619         zcond(i) = (pq(i)-zqsat)
1620     .     / (1. + FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor))
1621         IF (kcall.EQ.1) zcond(i) = MAX(zcond(i),0.)
1622         IF (kcall.EQ.2) zcond(i) = MIN(zcond(i),0.)
1623         pt(i) = pt(i) + zldcp*zcond(i)
1624         pq(i) = pq(i) - zcond(i)
1625      ENDIF
1626  210 CONTINUE
1627C
1628      is = 0
1629      DO i =1, klon
1630         IF (zcond(i).NE.0.) is = is + 1
1631      ENDDO
1632      IF (is.EQ.0) GOTO 230
1633C
1634      DO 220 i = 1, klon
1635      IF(ldflag(i).AND.zcond(i).NE.0.) THEN
1636         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1637         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1638         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1639         zqsat = R2ES* FOEEW(pt(i),zdelta) / pp(i)
1640         zqsat = MIN(0.5,zqsat)
1641         zcor = 1./(1.-RETV*zqsat)
1642         zqsat = zqsat*zcor
1643         zcond1 = (pq(i)-zqsat)
1644     .     / (1. + FOEDE(pt(i),zdelta,zcvm5,zqsat,zcor))
1645         pt(i) = pt(i) + zldcp*zcond1
1646         pq(i) = pq(i) - zcond1
1647      ENDIF
1648  220 CONTINUE
1649C
1650  230 CONTINUE
1651      RETURN
1652      END SUBROUTINE flxadjtq
1653      SUBROUTINE flxsetup
1654      IMPLICIT none
1655C
1656C     THIS ROUTINE DEFINES DISPOSABLE PARAMETERS FOR MASSFLUX SCHEME
1657C
1658      include "YOECUMF.h"
1659C
1660      ENTRPEN=1.0E-4  ! ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
1661      ENTRSCV=3.0E-4  ! ENTRAINMENT RATE FOR SHALLOW CONVECTION
1662      ENTRMID=1.0E-4  ! ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
1663      ENTRDD =2.0E-4  ! ENTRAINMENT RATE FOR DOWNDRAFTS
1664      CMFCTOP=0.33  ! RELATIVE CLOUD MASSFLUX AT LEVEL ABOVE NONBUO LEVEL
1665      CMFCMAX=1.0  ! MAXIMUM MASSFLUX VALUE ALLOWED FOR UPDRAFTS ETC
1666      CMFCMIN=1.E-10  ! MINIMUM MASSFLUX VALUE (FOR SAFETY)
1667      CMFDEPS=0.3  ! FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
1668      CPRCON =2.0E-4  ! CONVERSION FROM CLOUD WATER TO RAIN
1669      RHCDD=1.  ! RELATIVE SATURATION IN DOWNDRAFRS (NO LONGER USED)
1670c                 (FORMULATION IMPLIES SATURATION)
1671      LMFPEN = .TRUE.
1672      LMFSCV = .TRUE.
1673      LMFMID = .TRUE.
1674      LMFDD = .TRUE.
1675      LMFDUDV = .TRUE.
1676c
1677      RETURN
1678      END SUBROUTINE flxsetup
Note: See TracBrowser for help on using the repository browser.