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

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

Correction du bug sur les precip negatives. O.Boucher et J.Quaas
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.2 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
835            zzzmb = MAX(CMFCMIN, -pvervel(i,k)/RG)
836            zmfmax = (paph(i,k)-paph(i,k-1))/(RG*pdtime)
837            pmfub(i) = MIN(zzzmb,zmfmax)
838            pmfu(i,k+1) = pmfub(i)
839            pmfus(i,k+1) = pmfub(i)*(RCPD*ptu(i,k+1)+pgeoh(i,k+1))
840            pmfuq(i,k+1) = pmfub(i)*pqu(i,k+1)
841            pmful(i,k+1) = 0.0
842            pdmfup(i,k+1) = 0.0
843            kcbot(i) = k
844            klab(i,k+1) = 1
845            ktype(i) = 3
846            pentr(i) = ENTRMID
847         ENDIF
848         ENDDO
849      ENDIF
850c
851      is = 0
852      DO i = 1, klon
853         is = is + klab(i,k+1)
854         IF (klab(i,k+1) .EQ. 0) klab(i,k) = 0
855         llflag(i) = .FALSE.
856         IF (klab(i,k+1) .GT. 0) llflag(i) = .TRUE.
857      ENDDO
858      IF (is .EQ. 0) GOTO 480
859c
860c calculer le taux d'entrainement et de detrainement
861c
862      DO i = 1, klon
863         pen_u(i,k) = 0.0
864         pde_u(i,k) = 0.0
865         zrho(i)=paph(i,k+1)/(RD*ptenh(i,k+1))
866         zpbot(i)=paph(i,kcbot(i))
867         zptop(i)=paph(i,kctop0(i))
868      ENDDO
869c
870      DO 125 i = 1, klon
871      IF(ldcum(i)) THEN
872         zdprho=(paph(i,k+1)-paph(i,k))/(RG*zrho(i))
873         zentr=pentr(i)*pmfu(i,k+1)*zdprho
874         llo1=k.LT.kcbot(i)
875         IF(llo1) pde_u(i,k)=zentr
876         zpmid=0.5*(zpbot(i)+zptop(i))
877         llo2=llo1.AND.ktype(i).EQ.2.AND.
878     .        (zpbot(i)-paph(i,k).LT.0.2E5.OR.
879     .         paph(i,k).GT.zpmid)
880         IF(llo2) pen_u(i,k)=zentr
881         llo2=llo1.AND.(ktype(i).EQ.1.OR.ktype(i).EQ.3).AND.
882     .        (k.GE.MAX(klwmin(i),kctop0(i)+2).OR.pap(i,k).GT.zpmid)
883         IF(llo2) pen_u(i,k)=zentr
884         llo1=pen_u(i,k).GT.0..AND.(ktype(i).EQ.1.OR.ktype(i).EQ.2)
885         IF(llo1) THEN
886            zentr=zentr*(1.+3.*(1.-MIN(1.,(zpbot(i)-pap(i,k))/1.5E4)))
887            pen_u(i,k)=pen_u(i,k)*(1.+3.*(1.-MIN(1.,
888     .                 (zpbot(i)-pap(i,k))/1.5E4)))
889            pde_u(i,k)=pde_u(i,k)*(1.+3.*(1.-MIN(1.,
890     .                 (zpbot(i)-pap(i,k))/1.5E4)))
891         ENDIF
892         IF(llo2.AND.pqenh(i,k+1).GT.1.E-5)
893     .   pen_u(i,k)=zentr+MAX(pqte(i,k),0.)/pqenh(i,k+1)*
894     .              zrho(i)*zdprho
895      ENDIF
896  125 CONTINUE
897c
898C----------------------------------------------------------------------
899c DO ADIABATIC ASCENT FOR ENTRAINING/DETRAINING PLUME
900C----------------------------------------------------------------------
901c
902      DO 420 i = 1, klon
903      IF (llflag(i)) THEN
904         IF (k.LT.kcbot(i)) THEN
905            zmftest = pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)
906            zmfmax = MIN(zmftest,(paph(i,k)-paph(i,k-1))/(RG*pdtime))
907            pen_u(i,k)=MAX(pen_u(i,k)-MAX(0.0,zmftest-zmfmax),0.0)
908         ENDIF
909         pde_u(i,k)=MIN(pde_u(i,k),0.75*pmfu(i,k+1))
910c calculer le flux de masse du niveau k a partir de celui du k+1
911         pmfu(i,k)=pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)
912c calculer les valeurs Su, Qu et l du niveau k dans le panache montant
913         zqeen=pqenh(i,k+1)*pen_u(i,k)
914         zseen=(RCPD*ptenh(i,k+1)+pgeoh(i,k+1))*pen_u(i,k)
915         zscde=(RCPD*ptu(i,k+1)+pgeoh(i,k+1))*pde_u(i,k)
916         zqude=pqu(i,k+1)*pde_u(i,k)
917         plude(i,k)=plu(i,k+1)*pde_u(i,k)
918         zmfusk=pmfus(i,k+1)+zseen-zscde
919         zmfuqk=pmfuq(i,k+1)+zqeen-zqude
920         zmfulk=pmful(i,k+1)    -plude(i,k)
921         plu(i,k)=zmfulk*(1./MAX(CMFCMIN,pmfu(i,k)))
922         pqu(i,k)=zmfuqk*(1./MAX(CMFCMIN,pmfu(i,k)))
923         ptu(i,k)=(zmfusk*(1./MAX(CMFCMIN,pmfu(i,k)))-
924     1               pgeoh(i,k))/RCPD
925         ptu(i,k)=MAX(100.,ptu(i,k))
926         ptu(i,k)=MIN(400.,ptu(i,k))
927         zqold(i)=pqu(i,k)
928      ELSE
929         zqold(i)=0.0
930      ENDIF
931  420 CONTINUE
932c
933C----------------------------------------------------------------------
934c DO CORRECTIONS FOR MOIST ASCENT BY ADJUSTING T,Q AND L
935C----------------------------------------------------------------------
936c
937      icall = 1
938      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
939C
940      DO 440 i = 1, klon
941      IF(llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN
942         klab(i,k) = 2
943         plu(i,k) = plu(i,k)+zqold(i)-pqu(i,k)
944         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
945     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))
946         IF (klab(i,k+1).EQ.1) zbuo=zbuo+0.5
947         IF (zbuo.GT.0..AND.pmfu(i,k).GE.0.1*pmfub(i)) THEN
948            kctop(i) = k
949            ldcum(i) = .TRUE.
950            zdnoprc = 1.5E4
951            IF (ldland(i)) zdnoprc = zdland(i)
952            zprcon = CPRCON
953            IF ((zpbot(i)-paph(i,k)).LT.zdnoprc) zprcon = 0.0
954            zlnew=plu(i,k)/(1.+zprcon*(pgeoh(i,k)-pgeoh(i,k+1)))
955            pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))
956            plu(i,k)=zlnew
957         ELSE
958            klab(i,k)=0
959            pmfu(i,k)=0.
960         ENDIF
961      ENDIF
962  440 CONTINUE
963      DO 455 i = 1, klon
964      IF (llflag(i)) THEN
965         pmful(i,k)=plu(i,k)*pmfu(i,k)
966         pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)
967         pmfuq(i,k)=pqu(i,k)*pmfu(i,k)
968      ENDIF
969  455 CONTINUE
970C
971  480 CONTINUE
972C----------------------------------------------------------------------
973C DETERMINE CONVECTIVE FLUXES ABOVE NON-BUOYANCY LEVEL
974C    (NOTE: CLOUD VARIABLES LIKE T,Q AND L ARE NOT
975C           AFFECTED BY DETRAINMENT AND ARE ALREADY KNOWN
976C           FROM PREVIOUS CALCULATIONS ABOVE)
977C----------------------------------------------------------------------
978      DO i = 1, klon
979         IF (kctop(i).EQ.klev-1) ldcum(i) = .FALSE.
980         kcbot(i) = MAX(kcbot(i),kctop(i))
981      ENDDO
982c
983      ldcum(1)=ldcum(1)
984c
985      is = 0
986      DO i = 1, klon
987         if (ldcum(i)) is = is + 1
988      ENDDO
989      kcum = is
990      IF (is.EQ.0) GOTO 800
991c
992      DO 530 i = 1, klon
993      IF (ldcum(i)) THEN
994         k=kctop(i)-1
995         pde_u(i,k)=(1.-CMFCTOP)*pmfu(i,k+1)
996         plude(i,k)=pde_u(i,k)*plu(i,k+1)
997         pmfu(i,k)=pmfu(i,k+1)-pde_u(i,k)
998         zlnew=plu(i,k)
999         pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))
1000         plu(i,k)=zlnew
1001         pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)
1002         pmfuq(i,k)=pqu(i,k)*pmfu(i,k)
1003         pmful(i,k)=plu(i,k)*pmfu(i,k)
1004         plude(i,k-1)=pmful(i,k)
1005      ENDIF
1006  530 CONTINUE
1007C
1008  800 CONTINUE
1009      RETURN
1010      END
1011      SUBROUTINE flxflux(pdtime, pqen, pqsen, ptenh, pqenh, pap
1012     .  ,  paph, ldland, pgeoh, kcbot, kctop, lddraf, kdtop
1013     .  ,  ktype, ldcum, pmfu, pmfd, pmfus, pmfds
1014     .  ,  pmfuq, pmfdq, pmful, plude, pdmfup, pdmfdp
1015     .  ,  pten, prfl, psfl, pdpmel, ktopm2
1016     .  ,  pmflxr, pmflxs)
1017      IMPLICIT none
1018C----------------------------------------------------------------------
1019C THIS ROUTINE DOES THE FINAL CALCULATION OF CONVECTIVE
1020C FLUXES IN THE CLOUD LAYER AND IN THE SUBCLOUD LAYER
1021C----------------------------------------------------------------------
1022#include "dimensions.h"
1023#include "dimphy.h"
1024#include "YOMCST.h"
1025#include "YOETHF.h"
1026#include "YOECUMF.h"
1027C
1028      REAL cevapcu(klev)
1029C     -----------------------------------------------------------------
1030      REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)
1031      REAL pten(klon,klev), ptenh(klon,klev)
1032      REAL paph(klon,klev+1), pgeoh(klon,klev)
1033c
1034      REAL pap(klon,klev)
1035      REAL ztmsmlt, zdelta, zqsat
1036C
1037      REAL pmfu(klon,klev), pmfus(klon,klev)
1038      REAL pmfd(klon,klev), pmfds(klon,klev)
1039      REAL pmfuq(klon,klev), pmful(klon,klev)
1040      REAL pmfdq(klon,klev)
1041      REAL plude(klon,klev)
1042      REAL pdmfup(klon,klev), pdpmel(klon,klev)
1043cjq The variable maxpdmfdp(klon) has been introduced by Olivier Boucher
1044cjq 14/11/00 to fix the problem with the negative precipitation.     
1045      REAL pdmfdp(klon,klev), maxpdmfdp(klon)
1046      REAL prfl(klon), psfl(klon)
1047      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
1048      INTEGER  kcbot(klon), kctop(klon), ktype(klon)
1049      LOGICAL  ldland(klon), ldcum(klon)
1050      INTEGER k, i
1051      REAL zcons1, zcons2, zcucov, ztmelp2
1052      REAL pdtime, zdp, zzp, zfac, zsnmlt, zrfl, zrnew
1053      REAL zrmin, zrfln, zdrfl
1054      REAL zpds, zpdr, zdenom
1055      INTEGER ktopm2, itop, ikb
1056c
1057      LOGICAL lddraf(klon)
1058      INTEGER kdtop(klon)
1059c
1060#include "FCTTRE.h"
1061c
1062      DO 101 k=1,klev
1063      CEVAPCU(k)=1.93E-6*261.*SQRT(1.E3/(38.3*0.293)
1064     1 *SQRT(0.5*(paph(1,k)+paph(1,k+1))/paph(1,klev+1)) ) * 0.5/RG
1065 101  CONTINUE
1066c
1067c SPECIFY CONSTANTS
1068c
1069      zcons1 = RCPD/(RLMLT*RG*pdtime)
1070      zcons2 = 1./(RG*pdtime)
1071      zcucov = 0.05
1072      ztmelp2 = RTT + 2.
1073c
1074c DETERMINE FINAL CONVECTIVE FLUXES
1075c
1076      itop=klev
1077      DO 110 i = 1, klon
1078         itop=MIN(itop,kctop(i))
1079         IF (.NOT.ldcum(i) .OR. kdtop(i).LT.kctop(i)) lddraf(i)=.FALSE.
1080         IF(.NOT.ldcum(i)) ktype(i)=0
1081  110 CONTINUE
1082c
1083      ktopm2=itop-2
1084      DO 120 k=ktopm2,klev
1085      DO 115 i = 1, klon
1086      IF(ldcum(i).AND.k.GE.kctop(i)-1) THEN
1087         pmfus(i,k)=pmfus(i,k)-pmfu(i,k)*
1088     .                (RCPD*ptenh(i,k)+pgeoh(i,k))
1089         pmfuq(i,k)=pmfuq(i,k)-pmfu(i,k)*pqenh(i,k)
1090         zdp = 1.5E4
1091         IF ( ldland(i) ) zdp = 3.E4
1092c
1093c        l'eau liquide detrainee est precipitee quand certaines
1094c        conditions sont reunies (sinon, elle est consideree
1095c        evaporee dans l'environnement)
1096c
1097         IF(paph(i,kcbot(i))-paph(i,kctop(i)).GE.zdp.AND.
1098     .      pqen(i,k-1).GT.0.8*pqsen(i,k-1))
1099     .      pdmfup(i,k-1)=pdmfup(i,k-1)+plude(i,k-1)
1100c
1101         IF(lddraf(i).AND.k.GE.kdtop(i)) THEN
1102            pmfds(i,k)=pmfds(i,k)-pmfd(i,k)*
1103     .                   (RCPD*ptenh(i,k)+pgeoh(i,k))
1104            pmfdq(i,k)=pmfdq(i,k)-pmfd(i,k)*pqenh(i,k)
1105         ELSE
1106            pmfd(i,k)=0.
1107            pmfds(i,k)=0.
1108            pmfdq(i,k)=0.
1109            pdmfdp(i,k-1)=0.
1110         END IF
1111      ELSE
1112         pmfu(i,k)=0.
1113         pmfus(i,k)=0.
1114         pmfuq(i,k)=0.
1115         pmful(i,k)=0.
1116         pdmfup(i,k-1)=0.
1117         plude(i,k-1)=0.
1118         pmfd(i,k)=0.
1119         pmfds(i,k)=0.
1120         pmfdq(i,k)=0.
1121         pdmfdp(i,k-1)=0.
1122      ENDIF
1123  115 CONTINUE
1124  120 CONTINUE
1125c
1126      DO 130 k=ktopm2,klev
1127      DO 125 i = 1, klon
1128      IF(ldcum(i).AND.k.GT.kcbot(i)) THEN
1129         ikb=kcbot(i)
1130         zzp=((paph(i,klev+1)-paph(i,k))/
1131     .        (paph(i,klev+1)-paph(i,ikb)))
1132         IF (ktype(i).EQ.3) zzp = zzp**2
1133         pmfu(i,k)=pmfu(i,ikb)*zzp
1134         pmfus(i,k)=pmfus(i,ikb)*zzp
1135         pmfuq(i,k)=pmfuq(i,ikb)*zzp
1136         pmful(i,k)=pmful(i,ikb)*zzp
1137      ENDIF
1138  125 CONTINUE
1139  130 CONTINUE
1140c
1141c CALCULATE RAIN/SNOW FALL RATES
1142c CALCULATE MELTING OF SNOW
1143c CALCULATE EVAPORATION OF PRECIP
1144c
1145      DO k = 1, klev+1
1146      DO i = 1, klon
1147         pmflxr(i,k) = 0.0
1148         pmflxs(i,k) = 0.0
1149      ENDDO
1150      ENDDO
1151      DO k = ktopm2, klev
1152      DO i = 1, klon
1153      IF (ldcum(i)) THEN
1154         IF (pmflxs(i,k).GT.0.0 .AND. pten(i,k).GT.ztmelp2) THEN
1155            zfac=zcons1*(paph(i,k+1)-paph(i,k))
1156            zsnmlt=MIN(pmflxs(i,k),zfac*(pten(i,k)-ztmelp2))
1157            pdpmel(i,k)=zsnmlt
1158            ztmsmlt=pten(i,k)-zsnmlt/zfac
1159            zdelta=MAX(0.,SIGN(1.,RTT-ztmsmlt))
1160            zqsat=R2ES*FOEEW(ztmsmlt, zdelta) / pap(i,k)
1161            zqsat=MIN(0.5,zqsat)
1162            zqsat=zqsat/(1.-RETV  *zqsat)
1163            pqsen(i,k) = zqsat
1164         ENDIF
1165         IF (pten(i,k).GT.RTT) THEN
1166         pmflxr(i,k+1)=pmflxr(i,k)+pdmfup(i,k)+pdmfdp(i,k)+pdpmel(i,k)
1167         pmflxs(i,k+1)=pmflxs(i,k)-pdpmel(i,k)
1168         ELSE
1169           pmflxs(i,k+1)=pmflxs(i,k)+pdmfup(i,k)+pdmfdp(i,k)
1170           pmflxr(i,k+1)=pmflxr(i,k)
1171         ENDIF
1172c        si la precipitation est negative, on ajuste le plux du
1173c        panache descendant pour eliminer la negativite
1174         IF ((pmflxr(i,k+1)+pmflxs(i,k+1)).LT.0.0) THEN
1175cjq            write (*,*) 'total precip negative. rain ',pmflxr(i,k+1),
1176cjq     .                  'snow ',pmflxs(i,k+1)
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 i = 1, klon
1195         maxpdmfdp(i)=0.0
1196      ENDDO
1197      DO k = ktopm2, klev
1198         DO i = 1, klon
1199         if (k.GE.kcbot(i)) then
1200         maxpdmfdp(i)=maxpdmfdp(i)+pdmfdp(i,k)
1201         endif
1202         ENDDO
1203      ENDDO
1204cjq End of initialization
1205c     
1206      DO k = ktopm2, klev
1207      DO i = 1, klon
1208      IF (ldcum(i) .AND. k.GE.kcbot(i)) THEN
1209         zrfl = pmflxr(i,k) + pmflxs(i,k)
1210         IF (zrfl.GT.1.0E-20) THEN
1211            zrnew=(MAX(0.,SQRT(zrfl/zcucov)-
1212     .            CEVAPCU(k)*(paph(i,k+1)-paph(i,k))*
1213     .            MAX(0.,pqsen(i,k)-pqen(i,k))))**2*zcucov
1214            zrmin=zrfl-zcucov*MAX(0.,0.8*pqsen(i,k)-pqen(i,k))
1215     .            *zcons2*(paph(i,k+1)-paph(i,k))
1216            zrnew=MAX(zrnew,zrmin)
1217            zrfln=MAX(zrnew,0.)
1218            zdrfl=MIN(0.,zrfln-zrfl)
1219cjq At least the amount of precipiation needed to feed the downdraft
1220cjq with humidity below the base of convection has to be left and can't
1221cjq be evaporated (surely the evaporation can't be positive):           
1222            zdrfl=MAX(zdrfl,
1223     .            MIN(-pmflxr(i,k)-pmflxs(i,k)-maxpdmfdp(i),0.0))
1224cjq End of insertion
1225c           
1226            zdenom=1.0/MAX(1.0E-20,pmflxr(i,k)+pmflxs(i,k))
1227            IF (pten(i,k).GT.RTT) THEN
1228               zpdr = pdmfdp(i,k)
1229               zpds = 0.0
1230            ELSE
1231               zpdr = 0.0
1232               zpds = pdmfdp(i,k)
1233            ENDIF
1234            pmflxr(i,k+1) = pmflxr(i,k) + zpdr + pdpmel(i,k)
1235     .                    + zdrfl*pmflxr(i,k)*zdenom
1236            pmflxs(i,k+1) = pmflxs(i,k) + zpds - pdpmel(i,k)
1237     .                    + zdrfl*pmflxs(i,k)*zdenom
1238            pdmfup(i,k) = pdmfup(i,k) + zdrfl
1239         ELSE
1240            pmflxr(i,k+1) = 0.0
1241            pmflxs(i,k+1) = 0.0
1242            pdmfdp(i,k) = 0.0
1243            pdpmel(i,k) = 0.0
1244         ENDIF
1245      ENDIF
1246      ENDDO
1247      ENDDO
1248c
1249      DO 210 i = 1, klon
1250         prfl(i) = pmflxr(i,klev+1)
1251         psfl(i) = pmflxs(i,klev+1)
1252  210 CONTINUE
1253c
1254      RETURN
1255      END
1256      SUBROUTINE flxdtdq(pdtime, ktopm2, paph, ldcum, pten
1257     .  ,  pmfus, pmfds, pmfuq, pmfdq, pmful, pdmfup, pdmfdp
1258     .  ,  pdpmel, dt_con, dq_con)
1259      IMPLICIT none
1260c----------------------------------------------------------------------
1261c calculer les tendances T et Q
1262c----------------------------------------------------------------------
1263#include "dimensions.h"
1264#include "dimphy.h"
1265#include "YOMCST.h"
1266#include "YOETHF.h"
1267#include "YOECUMF.h"
1268C     -----------------------------------------------------------------
1269      LOGICAL  llo1
1270C
1271      REAL pten(klon,klev), paph(klon,klev+1)
1272      REAL pmfus(klon,klev), pmfuq(klon,klev), pmful(klon,klev)
1273      REAL pmfds(klon,klev), pmfdq(klon,klev)
1274      REAL pdmfup(klon,klev)
1275      REAL pdmfdp(klon,klev)
1276      REAL pdpmel(klon,klev)
1277      LOGICAL ldcum(klon)
1278      REAL dt_con(klon,klev), dq_con(klon,klev)
1279c
1280      INTEGER ktopm2
1281      REAL pdtime
1282c
1283      INTEGER i, k
1284      REAL zalv, zdtdt, zdqdt
1285c
1286      DO 210 k=ktopm2,klev-1
1287      DO 220 i = 1, klon
1288      IF (ldcum(i)) THEN
1289         llo1 = (pten(i,k)-RTT).GT.0.
1290         zalv = RLSTT
1291         IF (llo1) zalv = RLVTT
1292         zdtdt=RG/(paph(i,k+1)-paph(i,k))/RCPD
1293     .        *(pmfus(i,k+1)-pmfus(i,k)
1294     .         +pmfds(i,k+1)-pmfds(i,k)
1295     .          -RLMLT*pdpmel(i,k)
1296     .          -zalv*(pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1297     .         )
1298         dt_con(i,k)=zdtdt
1299         zdqdt=RG/(paph(i,k+1)-paph(i,k))
1300     .        *(pmfuq(i,k+1)-pmfuq(i,k)
1301     .         +pmfdq(i,k+1)-pmfdq(i,k)
1302     .          +pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1303         dq_con(i,k)=zdqdt
1304      ENDIF
1305  220 CONTINUE
1306  210 CONTINUE
1307C
1308      k = klev
1309      DO 230 i = 1, klon
1310      IF (ldcum(i)) THEN
1311         llo1 = (pten(i,k)-RTT).GT.0.
1312         zalv = RLSTT
1313         IF (llo1) zalv = RLVTT
1314         zdtdt=-RG/(paph(i,k+1)-paph(i,k))/RCPD
1315     .         *(pmfus(i,k)+pmfds(i,k)+RLMLT*pdpmel(i,k)
1316     .           -zalv*(pmful(i,k)+pdmfup(i,k)+pdmfdp(i,k)))
1317         dt_con(i,k)=zdtdt
1318         zdqdt=-RG/(paph(i,k+1)-paph(i,k))
1319     .            *(pmfuq(i,k)+pmfdq(i,k)+pmful(i,k)
1320     .             +pdmfup(i,k)+pdmfdp(i,k))
1321         dq_con(i,k)=zdqdt
1322      ENDIF
1323  230 CONTINUE
1324C
1325      RETURN
1326      END
1327      SUBROUTINE flxdlfs(ptenh, pqenh, pgeoh, paph, ptu, pqu,
1328     .     ldcum, kcbot, kctop, pmfub, prfl, ptd, pqd,
1329     .     pmfd, pmfds, pmfdq, pdmfdp, kdtop, lddraf)
1330      IMPLICIT none
1331C
1332C----------------------------------------------------------------------
1333C THIS ROUTINE CALCULATES LEVEL OF FREE SINKING FOR
1334C CUMULUS DOWNDRAFTS AND SPECIFIES T,Q,U AND V VALUES
1335C
1336C TO PRODUCE LFS-VALUES FOR CUMULUS DOWNDRAFTS
1337C FOR MASSFLUX CUMULUS PARAMETERIZATION
1338C
1339C INPUT ARE ENVIRONMENTAL VALUES OF T,Q,U,V,P,PHI
1340C AND UPDRAFT VALUES T,Q,U AND V AND ALSO
1341C CLOUD BASE MASSFLUX AND CU-PRECIPITATION RATE.
1342C IT RETURNS T,Q,U AND V VALUES AND MASSFLUX AT LFS.
1343C
1344C CHECK FOR NEGATIVE BUOYANCY OF AIR OF EQUAL PARTS OF
1345C MOIST ENVIRONMENTAL AIR AND CLOUD AIR.
1346C----------------------------------------------------------------------
1347#include "dimensions.h"
1348#include "dimphy.h"
1349#include "YOMCST.h"
1350#include "YOETHF.h"
1351#include "YOECUMF.h"
1352C
1353      REAL ptenh(klon,klev)
1354      REAL pqenh(klon,klev)
1355      REAL pgeoh(klon,klev), paph(klon,klev+1)
1356      REAL ptu(klon,klev), pqu(klon,klev)
1357      REAL pmfub(klon)
1358      REAL prfl(klon)
1359C
1360      REAL ptd(klon,klev), pqd(klon,klev)
1361      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1362      REAL pdmfdp(klon,klev)
1363      INTEGER  kcbot(klon), kctop(klon), kdtop(klon)
1364      LOGICAL  ldcum(klon), lddraf(klon)
1365C
1366      REAL ztenwb(klon,klev), zqenwb(klon,klev), zcond(klon)
1367      REAL zttest, zqtest, zbuo, zmftop
1368      LOGICAL  llo2(klon)
1369      INTEGER i, k, is, icall
1370C----------------------------------------------------------------------
1371      DO i= 1, klon
1372         lddraf(i)=.FALSE.
1373         kdtop(i)=klev+1
1374      ENDDO
1375C
1376C----------------------------------------------------------------------
1377C DETERMINE LEVEL OF FREE SINKING BY
1378C DOING A SCAN FROM TOP TO BASE OF CUMULUS CLOUDS
1379C
1380C FOR EVERY POINT AND PROCEED AS FOLLOWS:
1381C     (1) DETEMINE WET BULB ENVIRONMENTAL T AND Q
1382C     (2) DO MIXING WITH CUMULUS CLOUD AIR
1383C     (3) CHECK FOR NEGATIVE BUOYANCY
1384C
1385C THE ASSUMPTION IS THAT AIR OF DOWNDRAFTS IS MIXTURE
1386C OF 50% CLOUD AIR + 50% ENVIRONMENTAL AIR AT WET BULB
1387C TEMPERATURE (I.E. WHICH BECAME SATURATED DUE TO
1388C EVAPORATION OF RAIN AND CLOUD WATER)
1389C----------------------------------------------------------------------
1390C
1391      DO 290 k = 3, klev-3
1392C
1393      is=0
1394      DO 212 i= 1, klon
1395         ztenwb(i,k)=ptenh(i,k)
1396         zqenwb(i,k)=pqenh(i,k)
1397         llo2(i) = ldcum(i).AND.prfl(i).GT.0.
1398     .             .AND..NOT.lddraf(i)
1399     .             .AND.(k.LT.kcbot(i).AND.k.GT.kctop(i))
1400         IF ( llo2(i) ) is = is + 1
1401  212 CONTINUE
1402      IF(is.EQ.0) GO TO 290
1403C
1404      icall=2
1405      CALL flxadjtq(paph(1,k), ztenwb(1,k), zqenwb(1,k), llo2, icall)
1406C
1407C----------------------------------------------------------------------
1408C DO MIXING OF CUMULUS AND ENVIRONMENTAL AIR
1409C AND CHECK FOR NEGATIVE BUOYANCY.
1410C THEN SET VALUES FOR DOWNDRAFT AT LFS.
1411C----------------------------------------------------------------------
1412      DO 222 i= 1, klon
1413      IF (llo2(i)) THEN
1414         zttest=0.5*(ptu(i,k)+ztenwb(i,k))
1415         zqtest=0.5*(pqu(i,k)+zqenwb(i,k))
1416         zbuo=zttest*(1.+RETV*zqtest)-
1417     .        ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1418         zcond(i)=pqenh(i,k)-zqenwb(i,k)
1419         zmftop=-CMFDEPS*pmfub(i)
1420         IF (zbuo.LT.0..AND.prfl(i).GT.10.*zmftop*zcond(i)) THEN
1421            kdtop(i)=k
1422            lddraf(i)=.TRUE.
1423            ptd(i,k)=zttest
1424            pqd(i,k)=zqtest
1425            pmfd(i,k)=zmftop
1426            pmfds(i,k)=pmfd(i,k)*(RCPD*ptd(i,k)+pgeoh(i,k))
1427            pmfdq(i,k)=pmfd(i,k)*pqd(i,k)
1428            pdmfdp(i,k-1)=-0.5*pmfd(i,k)*zcond(i)
1429            prfl(i)=prfl(i)+pdmfdp(i,k-1)
1430         ENDIF
1431      ENDIF
1432  222 CONTINUE
1433c
1434  290 CONTINUE
1435C
1436      RETURN
1437      END
1438      SUBROUTINE flxddraf(ptenh, pqenh, pgeoh, paph, prfl,
1439     .           ptd, pqd, pmfd, pmfds, pmfdq, pdmfdp,
1440     .           lddraf, pen_d, pde_d)
1441      IMPLICIT none
1442C
1443C----------------------------------------------------------------------
1444C          THIS ROUTINE CALCULATES CUMULUS DOWNDRAFT DESCENT
1445C
1446C          TO PRODUCE THE VERTICAL PROFILES FOR CUMULUS DOWNDRAFTS
1447C          (I.E. T,Q,U AND V AND FLUXES)
1448C
1449C          INPUT IS T,Q,P,PHI,U,V AT HALF LEVELS.
1450C          IT RETURNS FLUXES OF S,Q AND EVAPORATION RATE
1451C          AND U,V AT LEVELS WHERE DOWNDRAFT OCCURS
1452C
1453C          CALCULATE MOIST DESCENT FOR ENTRAINING/DETRAINING PLUME BY
1454C          A) MOVING AIR DRY-ADIABATICALLY TO NEXT LEVEL BELOW AND
1455C          B) CORRECTING FOR EVAPORATION TO OBTAIN SATURATED STATE.
1456C
1457C----------------------------------------------------------------------
1458#include "dimensions.h"
1459#include "dimphy.h"
1460#include "YOMCST.h"
1461#include "YOETHF.h"
1462#include "YOECUMF.h"
1463C
1464      REAL ptenh(klon,klev), pqenh(klon,klev)
1465      REAL pgeoh(klon,klev), paph(klon,klev+1)
1466C
1467      REAL ptd(klon,klev), pqd(klon,klev)
1468      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1469      REAL pdmfdp(klon,klev)
1470      REAL prfl(klon)
1471      LOGICAL lddraf(klon)
1472C
1473      REAL pen_d(klon,klev), pde_d(klon,klev), zcond(klon)
1474      LOGICAL llo2(klon), llo1
1475      INTEGER i, k, is, icall, itopde
1476      REAL zentr, zseen, zqeen, zsdde, zqdde, zmfdsk, zmfdqk, zdmfdp
1477      REAL zbuo
1478C----------------------------------------------------------------------
1479C CALCULATE MOIST DESCENT FOR CUMULUS DOWNDRAFT BY
1480C       (A) CALCULATING ENTRAINMENT RATES, ASSUMING
1481C           LINEAR DECREASE OF MASSFLUX IN PBL
1482C       (B) DOING MOIST DESCENT - EVAPORATIVE COOLING
1483C           AND MOISTENING IS CALCULATED IN *flxadjtq*
1484C       (C) CHECKING FOR NEGATIVE BUOYANCY AND
1485C           SPECIFYING FINAL T,Q,U,V AND DOWNWARD FLUXES
1486C
1487      DO 180 k = 3, klev
1488c
1489      is = 0
1490      DO i = 1, klon
1491         llo2(i)=lddraf(i).AND.pmfd(i,k-1).LT.0.
1492         IF (llo2(i)) is = is + 1
1493      ENDDO
1494      IF (is.EQ.0) GOTO 180
1495c
1496      DO i = 1, klon
1497      IF (llo2(i)) THEN
1498         zentr = ENTRDD*pmfd(i,k-1)*RD*ptenh(i,k-1)/
1499     .           (RG*paph(i,k-1))*(paph(i,k)-paph(i,k-1))
1500         pen_d(i,k) = zentr
1501         pde_d(i,k) = zentr
1502      ENDIF
1503      ENDDO
1504c
1505      itopde = klev-2
1506      IF (k.GT.itopde) THEN
1507         DO i = 1, klon
1508         IF (llo2(i)) THEN
1509            pen_d(i,k)=0.
1510            pde_d(i,k)=pmfd(i,itopde)*
1511     .      (paph(i,k)-paph(i,k-1))/(paph(i,klev+1)-paph(i,itopde))
1512         ENDIF
1513         ENDDO
1514      ENDIF
1515C
1516      DO i = 1, klon
1517      IF (llo2(i)) THEN
1518         pmfd(i,k) = pmfd(i,k-1)+pen_d(i,k)-pde_d(i,k)
1519         zseen = (RCPD*ptenh(i,k-1)+pgeoh(i,k-1))*pen_d(i,k)
1520         zqeen = pqenh(i,k-1)*pen_d(i,k)
1521         zsdde = (RCPD*ptd(i,k-1)+pgeoh(i,k-1))*pde_d(i,k)
1522         zqdde = pqd(i,k-1)*pde_d(i,k)
1523         zmfdsk = pmfds(i,k-1)+zseen-zsdde
1524         zmfdqk = pmfdq(i,k-1)+zqeen-zqdde
1525         pqd(i,k) = zmfdqk*(1./MIN(-CMFCMIN,pmfd(i,k)))
1526         ptd(i,k) = (zmfdsk*(1./MIN(-CMFCMIN,pmfd(i,k)))-
1527     .               pgeoh(i,k))/RCPD
1528         ptd(i,k) = MIN(400.,ptd(i,k))
1529         ptd(i,k) = MAX(100.,ptd(i,k))
1530         zcond(i) = pqd(i,k)
1531      ENDIF
1532      ENDDO
1533C
1534      icall = 2
1535      CALL flxadjtq(paph(1,k), ptd(1,k), pqd(1,k), llo2, icall)
1536C
1537      DO i = 1, klon
1538      IF (llo2(i)) THEN
1539         zcond(i) = zcond(i)-pqd(i,k)
1540         zbuo = ptd(i,k)*(1.+RETV  *pqd(i,k))-
1541     .          ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1542         llo1 = zbuo.LT.0..AND.(prfl(i)-pmfd(i,k)*zcond(i).GT.0.)
1543         IF (.not.llo1) pmfd(i,k) = 0.0
1544         pmfds(i,k) = (RCPD*ptd(i,k)+pgeoh(i,k))*pmfd(i,k)
1545         pmfdq(i,k) = pqd(i,k)*pmfd(i,k)
1546         zdmfdp = -pmfd(i,k)*zcond(i)
1547         pdmfdp(i,k-1) = zdmfdp
1548         prfl(i) = prfl(i)+zdmfdp
1549      ENDIF
1550      ENDDO
1551c
1552  180 CONTINUE
1553      RETURN
1554      END
1555      SUBROUTINE flxadjtq(pp, pt, pq, ldflag, kcall)
1556      IMPLICIT none
1557c======================================================================
1558c Objet: ajustement entre T et Q
1559c======================================================================
1560C NOTE: INPUT PARAMETER kcall DEFINES CALCULATION AS
1561C        kcall=0    ENV. T AND QS IN*CUINI*
1562C        kcall=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)
1563C        kcall=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1564C
1565#include "dimensions.h"
1566#include "dimphy.h"
1567#include "YOMCST.h"
1568C
1569      REAL pt(klon), pq(klon), pp(klon)
1570      LOGICAL ldflag(klon)
1571      INTEGER kcall
1572c
1573      REAL zcond(klon), zcond1
1574      REAL Z5alvcp, z5alscp, zalvdcp, zalsdcp
1575      REAL zdelta, zcvm5, zldcp, zqsat, zcor
1576      INTEGER is, i
1577#include "YOETHF.h"
1578#include "FCTTRE.h"
1579C
1580      z5alvcp = r5les*RLVTT/RCPD
1581      z5alscp = r5ies*RLSTT/RCPD
1582      zalvdcp = rlvtt/RCPD
1583      zalsdcp = rlstt/RCPD
1584C
1585
1586      DO i = 1, klon
1587         zcond(i) = 0.0
1588      ENDDO
1589
1590      DO 210 i =1, klon
1591      IF (ldflag(i)) THEN
1592         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1593         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1594         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1595         zqsat = R2ES*FOEEW(pt(i),zdelta) / pp(i)
1596         zqsat = MIN(0.5,zqsat)
1597         zcor = 1./(1.-RETV*zqsat)
1598         zqsat = zqsat*zcor
1599         zcond(i) = (pq(i)-zqsat)
1600     .     / (1. + FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor))
1601         IF (kcall.EQ.1) zcond(i) = MAX(zcond(i),0.)
1602         IF (kcall.EQ.2) zcond(i) = MIN(zcond(i),0.)
1603         pt(i) = pt(i) + zldcp*zcond(i)
1604         pq(i) = pq(i) - zcond(i)
1605      ENDIF
1606  210 CONTINUE
1607C
1608      is = 0
1609      DO i =1, klon
1610         IF (zcond(i).NE.0.) is = is + 1
1611      ENDDO
1612      IF (is.EQ.0) GOTO 230
1613C
1614      DO 220 i = 1, klon
1615      IF(ldflag(i).AND.zcond(i).NE.0.) THEN
1616         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1617         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1618         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1619         zqsat = R2ES* FOEEW(pt(i),zdelta) / pp(i)
1620         zqsat = MIN(0.5,zqsat)
1621         zcor = 1./(1.-RETV*zqsat)
1622         zqsat = zqsat*zcor
1623         zcond1 = (pq(i)-zqsat)
1624     .     / (1. + FOEDE(pt(i),zdelta,zcvm5,zqsat,zcor))
1625         pt(i) = pt(i) + zldcp*zcond1
1626         pq(i) = pq(i) - zcond1
1627      ENDIF
1628  220 CONTINUE
1629C
1630  230 CONTINUE
1631      RETURN
1632      END
1633      SUBROUTINE flxsetup
1634      IMPLICIT none
1635C
1636C     THIS ROUTINE DEFINES DISPOSABLE PARAMETERS FOR MASSFLUX SCHEME
1637C
1638#include "YOECUMF.h"
1639C
1640      ENTRPEN=1.0E-4  ! ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
1641      ENTRSCV=3.0E-4  ! ENTRAINMENT RATE FOR SHALLOW CONVECTION
1642      ENTRMID=1.0E-4  ! ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
1643      ENTRDD =2.0E-4  ! ENTRAINMENT RATE FOR DOWNDRAFTS
1644      CMFCTOP=0.33  ! RELATIVE CLOUD MASSFLUX AT LEVEL ABOVE NONBUO LEVEL
1645      CMFCMAX=1.0  ! MAXIMUM MASSFLUX VALUE ALLOWED FOR UPDRAFTS ETC
1646      CMFCMIN=1.E-10  ! MINIMUM MASSFLUX VALUE (FOR SAFETY)
1647      CMFDEPS=0.3  ! FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
1648      CPRCON =2.0E-4  ! CONVERSION FROM CLOUD WATER TO RAIN
1649      RHCDD=1.  ! RELATIVE SATURATION IN DOWNDRAFRS (NO LONGER USED)
1650c                 (FORMULATION IMPLIES SATURATION)
1651      LMFPEN = .TRUE.
1652      LMFSCV = .TRUE.
1653      LMFMID = .TRUE.
1654      LMFDD = .TRUE.
1655      LMFDUDV = .TRUE.
1656c
1657      RETURN
1658      END
Note: See TracBrowser for help on using the repository browser.