source: LMDZ.3.3/trunk/libf/phylmd/conflx.F @ 286

Last change on this file since 286 was 286, checked in by lmdz, 23 years ago

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