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

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

la tendance de l'eau doit etre inversee dans conflx D. Hauglustaine
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 53.8 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         pmflxr(i,k) = 0.0
106         pmflxs(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.1*pqhfl(i)*RG)) ktype(i) = 1
320      ENDDO
321c
322c---------------------------------------------------------------------
323c determiner le flux de masse entrant a travers la base.
324c on ignore, pour l'instant, l'effet du panache descendant
325c---------------------------------------------------------------------
326      DO i = 1, klon
327         ikb=kcbot(i)
328         zqumqe=pqu(i,ikb)+plu(i,ikb)-zqenh(i,ikb)
329         zdqmin=MAX(0.01*zqenh(i,ikb),1.E-10)
330         IF (zdqpbl(i).GT.0..AND.zqumqe.GT.zdqmin.AND.ldcum(i)) THEN
331            zmfub(i) = zdqpbl(i)/(RG*MAX(zqumqe,zdqmin))
332         ELSE
333            zmfub(i) = 0.01
334            ldcum(i)=.FALSE.
335         ENDIF
336         IF (ktype(i).EQ.2) THEN
337            zdh = RCPD*(ptu(i,ikb)-ztenh(i,ikb)) + RLVTT*zqumqe
338            zdh = RG * MAX(zdh,1.0E5*zdqmin)
339            IF (zdhpbl(i).GT.0..AND.ldcum(i))zmfub(i)=zdhpbl(i)/zdh
340         ENDIF
341         zmfmax = (paph(i,ikb)-paph(i,ikb-1)) / (RG*pdtime)
342         zmfub(i) = MIN(zmfub(i),zmfmax)
343         zentr(i) = ENTRSCV
344         IF (ktype(i).EQ.1) zentr(i) = ENTRPEN
345      ENDDO
346C-----------------------------------------------------------------------
347C DETERMINE CLOUD ASCENT FOR ENTRAINING PLUME
348C-----------------------------------------------------------------------
349c (A) calculer d'abord la hauteur "theorique" de la tour convective sans
350c     considerer l'entrainement ni le detrainement du panache, sachant
351c     ces derniers peuvent abaisser la hauteur theorique.
352c
353      DO i = 1, klon
354         ikb=kcbot(i)
355         zhcbase(i)=RCPD*ptu(i,ikb)+zgeoh(i,ikb)+RLVTT*pqu(i,ikb)
356         ictop0(i)=kcbot(i)-1
357      ENDDO
358c
359      zalvdcp=RLVTT/RCPD
360      DO k=klev-1,3,-1
361      DO i = 1, klon
362         zhsat=RCPD*ztenh(i,k)+zgeoh(i,k)+RLVTT*zqsenh(i,k)
363         zgam=R5LES*zalvdcp*zqsenh(i,k)/
364     .        ((1.-RETV  *zqsenh(i,k))*(ztenh(i,k)-R4LES)**2)
365         zzz=RCPD*ztenh(i,k)*0.608
366         zhhat=zhsat-(zzz+zgam*zzz)/(1.+zgam*zzz/RLVTT)*
367     .               MAX(zqsenh(i,k)-zqenh(i,k),0.)
368         IF(k.LT.ictop0(i).AND.zhcbase(i).GT.zhhat) ictop0(i)=k
369      ENDDO
370      ENDDO
371c
372c (B) calculer le panache ascendant
373c
374      CALL flxasc(pdtime,ztenh, zqenh, pten, pqen, pqsen,
375     .     pgeo, zgeoh, pap, paph, pqte, pvervel,
376     .     ldland, ldcum, ktype, ilab,
377     .     ptu, pqu, plu, pmfu, zmfub, zentr,
378     .     zmfus, zmfuq, zmful, plude, zdmfup,
379     .     kcbot, kctop, ictop0, kcum, pen_u, pde_u)
380      IF (kcum.EQ.0) GO TO 1000
381C
382C verifier l'epaisseur de la convection et changer eventuellement
383c le taux d'entrainement/detrainement
384C
385      DO i = 1, klon
386         zpbmpt=paph(i,kcbot(i))-paph(i,kctop(i))
387         IF(ldcum(i).AND.ktype(i).EQ.1.AND.zpbmpt.LT.2.E4)ktype(i)=2
388         IF(ldcum(i)) ictop0(i)=kctop(i)
389         IF(ktype(i).EQ.2) zentr(i)=ENTRSCV
390      ENDDO
391c
392      IF (lmfdd) THEN  ! si l'on considere le panache descendant
393c
394c calculer la precipitation issue du panache ascendant pour
395c determiner l'existence du panache descendant dans la convection
396      DO i = 1, klon
397         zrfl(i)=zdmfup(i,1)
398      ENDDO
399      DO k=2,klev
400      DO i = 1, klon
401         zrfl(i)=zrfl(i)+zdmfup(i,k)
402      ENDDO
403      ENDDO
404c
405c determiner le LFS (level of free sinking: niveau de plonge libre)
406      CALL flxdlfs(ztenh, zqenh, zgeoh, paph, ptu, pqu,
407     *     ldcum,    kcbot,    kctop,    zmfub,    zrfl,
408     *     ptd,      pqd,
409     *     pmfd,     zmfds,    zmfdq,    zdmfdp,
410     *     kdtop,    lddraf)
411c
412c calculer le panache descendant
413      CALL flxddraf(ztenh,    zqenh,
414     *     zgeoh,    paph,     zrfl,
415     *     ptd,      pqd,
416     *     pmfd,     zmfds,    zmfdq,    zdmfdp,
417     *     lddraf, pen_d, pde_d)
418c
419c calculer de nouveau le flux de masse entrant a travers la base
420c de la convection, sachant qu'il a ete modifie par le panache
421c descendant
422      DO i = 1, klon
423      IF (lddraf(i)) THEN
424         ikb = kcbot(i)
425         llo1 = PMFD(i,ikb).LT.0.
426         zeps = 0.
427         IF ( llo1 ) zeps = CMFDEPS
428         zqumqe = pqu(i,ikb)+plu(i,ikb)-
429     .            zeps*pqd(i,ikb)-(1.-zeps)*zqenh(i,ikb)
430         zdqmin = MAX(0.01*zqenh(i,ikb),1.E-10)
431         zmfmax = (paph(i,ikb)-paph(i,ikb-1)) / (RG*pdtime)
432         IF (zdqpbl(i).GT.0..AND.zqumqe.GT.zdqmin.AND.ldcum(i)
433     .       .AND.zmfub(i).LT.zmfmax) THEN
434            zmfub1(i) = zdqpbl(i) / (RG*MAX(zqumqe,zdqmin))
435         ELSE
436            zmfub1(i) = zmfub(i)
437         ENDIF
438         IF (ktype(i).EQ.2) THEN
439            zdh = RCPD*(ptu(i,ikb)-zeps*ptd(i,ikb)-
440     .            (1.-zeps)*ztenh(i,ikb))+RLVTT*zqumqe
441            zdh = RG * MAX(zdh,1.0E5*zdqmin)
442            IF (zdhpbl(i).GT.0..AND.ldcum(i))zmfub1(i)=zdhpbl(i)/zdh
443         ENDIF
444         IF ( .NOT.((ktype(i).EQ.1.OR.ktype(i).EQ.2).AND.
445     .              ABS(zmfub1(i)-zmfub(i)).LT.0.2*zmfub(i)) )
446     .      zmfub1(i) = zmfub(i)
447      ENDIF
448      ENDDO
449      DO k = 1, klev
450      DO i = 1, klon
451      IF (lddraf(i)) THEN
452         zfac = zmfub1(i)/MAX(zmfub(i),1.E-10)
453         pmfd(i,k) = pmfd(i,k)*zfac
454         zmfds(i,k) = zmfds(i,k)*zfac
455         zmfdq(i,k) = zmfdq(i,k)*zfac
456         zdmfdp(i,k) = zdmfdp(i,k)*zfac
457         pen_d(i,k) = pen_d(i,k)*zfac
458         pde_d(i,k) = pde_d(i,k)*zfac
459      ENDIF
460      ENDDO
461      ENDDO
462      DO i = 1, klon
463         IF (lddraf(i)) zmfub(i)=zmfub1(i)
464      ENDDO
465c
466      ENDIF   ! fin de test sur lmfdd
467c
468c-----------------------------------------------------------------------
469c calculer de nouveau le panache ascendant
470c-----------------------------------------------------------------------
471      CALL flxasc(pdtime,ztenh, zqenh, pten, pqen, pqsen,
472     .     pgeo, zgeoh, pap, paph, pqte, pvervel,
473     .     ldland, ldcum, ktype, ilab,
474     .     ptu, pqu, plu, pmfu, zmfub, zentr,
475     .     zmfus, zmfuq, zmful, plude, zdmfup,
476     .     kcbot, kctop, ictop0, kcum, pen_u, pde_u)
477c
478c-----------------------------------------------------------------------
479c determiner les flux convectifs en forme finale, ainsi que
480c la quantite des precipitations
481c-----------------------------------------------------------------------
482      CALL flxflux(pdtime, pqen, pqsen, ztenh, zqenh, pap, paph,
483     .     ldland, zgeoh, kcbot, kctop, lddraf, kdtop, ktype, ldcum,
484     .     pmfu, pmfd, zmfus, zmfds, zmfuq, zmfdq, zmful, plude,
485     .     zdmfup, zdmfdp, pten, prsfc, pssfc, zdpmel, itopm2,
486     .     pmflxr, pmflxs)
487c
488c----------------------------------------------------------------------
489c calculer les tendances pour T et Q
490c----------------------------------------------------------------------
491      CALL flxdtdq(pdtime, itopm2, paph, ldcum, pten,
492     e     zmfus, zmfds, zmfuq, zmfdq, zmful, zdmfup, zdmfdp, zdpmel,
493     s     dt_con,dq_con)
494c
495 1000 CONTINUE
496      RETURN
497      END
498      SUBROUTINE flxini(pten, pqen, pqsen, pgeo, paph, pgeoh, ptenh,
499     .           pqenh, pqsenh, ptu, pqu, ptd, pqd, pmfd, pmfds, pmfdq,
500     .           pdmfdp, pmfu, pmfus, pmfuq, pdmfup, pdpmel, plu, plude,
501     .           klab,pen_u, pde_u, pen_d, pde_d)
502      IMPLICIT none
503C----------------------------------------------------------------------
504C THIS ROUTINE INTERPOLATES LARGE-SCALE FIELDS OF T,Q ETC.
505C TO HALF LEVELS (I.E. GRID FOR MASSFLUX SCHEME),
506C AND INITIALIZES VALUES FOR UPDRAFTS
507C----------------------------------------------------------------------
508#include "dimensions.h"
509#include "dimphy.h"
510#include "YOMCST.h"
511#include "YOETHF.h"
512C
513      REAL pten(klon,klev)   ! temperature (environnement)
514      REAL pqen(klon,klev)   ! humidite (environnement)
515      REAL pqsen(klon,klev)  ! humidite saturante (environnement)
516      REAL pgeo(klon,klev)   ! geopotentiel (g * metre)
517      REAL pgeoh(klon,klev)  ! geopotentiel aux demi-niveaux
518      REAL paph(klon,klev+1) ! pression aux demi-niveaux
519      REAL ptenh(klon,klev)  ! temperature aux demi-niveaux
520      REAL pqenh(klon,klev)  ! humidite aux demi-niveaux
521      REAL pqsenh(klon,klev) ! humidite saturante aux demi-niveaux
522C
523      REAL ptu(klon,klev)    ! temperature du panache ascendant (p-a)
524      REAL pqu(klon,klev)    ! humidite du p-a
525      REAL plu(klon,klev)    ! eau liquide du p-a
526      REAL pmfu(klon,klev)   ! flux de masse du p-a
527      REAL pmfus(klon,klev)  ! flux de l'energie seche dans le p-a
528      REAL pmfuq(klon,klev)  ! flux de l'humidite dans le p-a
529      REAL pdmfup(klon,klev) ! quantite de l'eau precipitee dans p-a
530      REAL plude(klon,klev)  ! quantite de l'eau liquide jetee du
531c                              p-a a l'environnement
532      REAL pdpmel(klon,klev) ! quantite de neige fondue
533c
534      REAL ptd(klon,klev)    ! temperature du panache descendant (p-d)
535      REAL pqd(klon,klev)    ! humidite du p-d
536      REAL pmfd(klon,klev)   ! flux de masse du p-d
537      REAL pmfds(klon,klev)  ! flux de l'energie seche dans le p-d
538      REAL pmfdq(klon,klev)  ! flux de l'humidite dans le p-d
539      REAL pdmfdp(klon,klev) ! quantite de precipitation dans p-d
540c
541      REAL pen_u(klon,klev) ! quantite de masse entrainee pour p-a
542      REAL pde_u(klon,klev) ! quantite de masse detrainee pour p-a
543      REAL pen_d(klon,klev) ! quantite de masse entrainee pour p-d
544      REAL pde_d(klon,klev) ! quantite de masse detrainee pour p-d
545C
546      INTEGER  klab(klon,klev)
547      LOGICAL  llflag(klon)
548      INTEGER k, i, icall
549      REAL zzs
550C----------------------------------------------------------------------
551C SPECIFY LARGE SCALE PARAMETERS AT HALF LEVELS
552C ADJUST TEMPERATURE FIELDS IF STATICLY UNSTABLE
553C----------------------------------------------------------------------
554      DO 130 k = 2, klev
555c
556      DO i = 1, klon
557         pgeoh(i,k)=pgeo(i,k)+(pgeo(i,k-1)-pgeo(i,k))*0.5
558         ptenh(i,k)=(MAX(RCPD*pten(i,k-1)+pgeo(i,k-1),
559     .             RCPD*pten(i,k)+pgeo(i,k))-pgeoh(i,k))/RCPD
560         pqsenh(i,k)=pqsen(i,k-1)
561         llflag(i)=.TRUE.
562      ENDDO
563c
564      icall=0
565      CALL flxadjtq(paph(1,k),ptenh(1,k),pqsenh(1,k),llflag,icall)
566c
567      DO i = 1, klon
568         pqenh(i,k)=MIN(pqen(i,k-1),pqsen(i,k-1))
569     .               +(pqsenh(i,k)-pqsen(i,k-1))
570         pqenh(i,k)=MAX(pqenh(i,k),0.)
571      ENDDO
572c
573  130 CONTINUE
574C
575      DO 140 i = 1, klon
576         ptenh(i,klev)=(RCPD*pten(i,klev)+pgeo(i,klev)-
577     1                   pgeoh(i,klev))/RCPD
578         pqenh(i,klev)=pqen(i,klev)
579         ptenh(i,1)=pten(i,1)
580         pqenh(i,1)=pqen(i,1)
581         pgeoh(i,1)=pgeo(i,1)
582  140 CONTINUE
583c
584      DO 160 k = klev-1, 2, -1
585      DO 150 i = 1, klon
586         zzs = MAX(RCPD*ptenh(i,k)+pgeoh(i,k),
587     .             RCPD*ptenh(i,k+1)+pgeoh(i,k+1))
588         ptenh(i,k) = (zzs-pgeoh(i,k))/RCPD
589  150 CONTINUE
590  160 CONTINUE
591C
592C-----------------------------------------------------------------------
593C INITIALIZE VALUES FOR UPDRAFTS AND DOWNDRAFTS
594C-----------------------------------------------------------------------
595      DO k = 1, klev
596      DO i = 1, klon
597         ptu(i,k) = ptenh(i,k)
598         pqu(i,k) = pqenh(i,k)
599         plu(i,k) = 0.
600         pmfu(i,k) = 0.
601         pmfus(i,k) = 0.
602         pmfuq(i,k) = 0.
603         pdmfup(i,k) = 0.
604         pdpmel(i,k) = 0.
605         plude(i,k) = 0.
606c
607         klab(i,k) = 0
608c
609         ptd(i,k) = ptenh(i,k)
610         pqd(i,k) = pqenh(i,k)
611         pmfd(i,k) = 0.0
612         pmfds(i,k) = 0.0
613         pmfdq(i,k) = 0.0
614         pdmfdp(i,k) = 0.0
615c
616         pen_u(i,k) = 0.0
617         pde_u(i,k) = 0.0
618         pen_d(i,k) = 0.0
619         pde_d(i,k) = 0.0
620      ENDDO
621      ENDDO
622C
623      RETURN
624      END
625      SUBROUTINE flxbase(ptenh, pqenh, pgeoh, paph,
626     *     ptu, pqu, plu, ldcum, kcbot, klab)
627      IMPLICIT none
628C----------------------------------------------------------------------
629C THIS ROUTINE CALCULATES CLOUD BASE VALUES (T AND Q)
630C
631C INPUT ARE ENVIRONM. VALUES OF T,Q,P,PHI AT HALF LEVELS.
632C IT RETURNS CLOUD BASE VALUES AND FLAGS AS FOLLOWS;
633C   klab=1 FOR SUBCLOUD LEVELS
634C   klab=2 FOR CONDENSATION LEVEL
635C
636C LIFT SURFACE AIR DRY-ADIABATICALLY TO CLOUD BASE
637C (NON ENTRAINING PLUME,I.E.CONSTANT MASSFLUX)
638C----------------------------------------------------------------------
639#include "dimensions.h"
640#include "dimphy.h"
641#include "YOMCST.h"
642#include "YOETHF.h"
643C       ----------------------------------------------------------------
644      REAL ptenh(klon,klev), pqenh(klon,klev)
645      REAL pgeoh(klon,klev), paph(klon,klev+1)
646C
647      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
648      INTEGER  klab(klon,klev), kcbot(klon)
649C
650      LOGICAL llflag(klon), ldcum(klon)
651      INTEGER i, k, icall, is
652      REAL zbuo, zqold(klon)
653C----------------------------------------------------------------------
654C INITIALIZE VALUES AT LIFTING LEVEL
655C----------------------------------------------------------------------
656      DO i = 1, klon
657         klab(i,klev)=1
658         kcbot(i)=klev-1
659         ldcum(i)=.FALSE.
660      ENDDO
661C----------------------------------------------------------------------
662C DO ASCENT IN SUBCLOUD LAYER,
663C CHECK FOR EXISTENCE OF CONDENSATION LEVEL,
664C ADJUST T,Q AND L ACCORDINGLY
665C CHECK FOR BUOYANCY AND SET FLAGS
666C----------------------------------------------------------------------
667      DO 290 k = klev-1, 2, -1
668c
669      is = 0
670      DO i = 1, klon
671         IF (klab(i,k+1).EQ.1) is = is + 1
672         llflag(i) = .FALSE.
673         IF (klab(i,k+1).EQ.1) llflag(i) = .TRUE.
674      ENDDO
675      IF (is.EQ.0) GOTO 290
676c
677      DO i = 1, klon
678      IF(llflag(i)) THEN
679         pqu(i,k) = pqu(i,k+1)
680         ptu(i,k) = ptu(i,k+1)+(pgeoh(i,k+1)-pgeoh(i,k))/RCPD
681         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
682     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))+0.5
683         IF (zbuo.GT.0.) klab(i,k) = 1
684         zqold(i) = pqu(i,k)
685      ENDIF
686      ENDDO
687c
688      icall=1
689      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
690c
691      DO i = 1, klon
692      IF (llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN
693         klab(i,k) = 2
694         plu(i,k) = plu(i,k) + zqold(i)-pqu(i,k)
695         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
696     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))+0.5
697         IF (zbuo.GT.0.) kcbot(i) = k
698         IF (zbuo.GT.0.) ldcum(i) = .TRUE.
699      ENDIF
700      ENDDO
701c
702  290 CONTINUE
703c
704      RETURN
705      END
706      SUBROUTINE flxasc(pdtime, ptenh, pqenh, pten, pqen, pqsen,
707     .     pgeo, pgeoh, pap, paph, pqte, pvervel,
708     .     ldland, ldcum, ktype, klab, ptu, pqu, plu,
709     .     pmfu, pmfub, pentr, pmfus, pmfuq,
710     .     pmful, plude, pdmfup, kcbot, kctop, kctop0, kcum,
711     .     pen_u, pde_u)
712      IMPLICIT none
713C----------------------------------------------------------------------
714C THIS ROUTINE DOES THE CALCULATIONS FOR CLOUD ASCENTS
715C FOR CUMULUS PARAMETERIZATION
716C----------------------------------------------------------------------
717#include "dimensions.h"
718#include "dimphy.h"
719#include "YOMCST.h"
720#include "YOETHF.h"
721#include "YOECUMF.h"
722C
723      REAL pdtime
724      REAL pten(klon,klev), ptenh(klon,klev)
725      REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)
726      REAL pgeo(klon,klev), pgeoh(klon,klev)
727      REAL pap(klon,klev), paph(klon,klev+1)
728      REAL pqte(klon,klev)
729      REAL pvervel(klon,klev) ! vitesse verticale en Pa/s
730C
731      REAL pmfub(klon), pentr(klon)
732      REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)
733      REAL plude(klon,klev)
734      REAL pmfu(klon,klev), pmfus(klon,klev)
735      REAL pmfuq(klon,klev), pmful(klon,klev)
736      REAL pdmfup(klon,klev)
737      INTEGER ktype(klon), klab(klon,klev), kcbot(klon), kctop(klon)
738      INTEGER kctop0(klon)
739      LOGICAL ldland(klon), ldcum(klon)
740C
741      REAL pen_u(klon,klev), pde_u(klon,klev)
742      REAL zqold(klon)
743      REAL zdland(klon)
744      LOGICAL llflag(klon)
745      INTEGER k, i, is, icall, kcum
746      REAL ztglace, zdphi, zqeen, zseen, zscde, zqude
747      REAL zmfusk, zmfuqk, zmfulk, zbuo, zdnoprc, zprcon, zlnew
748c
749      REAL zpbot(klon), zptop(klon), zrho(klon)
750      REAL zdprho, zentr, zpmid, zmftest, zmfmax
751      LOGICAL llo1, llo2
752c
753      REAL zwmax(klon), zzzmb
754      INTEGER klwmin(klon) ! level of maximum vertical velocity
755C----------------------------------------------------------------------
756      ztglace = RTT - 13.
757c
758c Chercher le niveau ou la vitesse verticale est maximale:
759      DO i = 1, klon
760         klwmin(i) = klev
761         zwmax(i) = 0.0
762      ENDDO
763      DO k = klev, 3, -1
764      DO i = 1, klon
765      IF (pvervel(i,k).LT.zwmax(i)) THEN
766         zwmax(i) = pvervel(i,k)
767         klwmin(i) = k
768      ENDIF
769      ENDDO
770      ENDDO
771C----------------------------------------------------------------------
772C SET DEFAULT VALUES
773C----------------------------------------------------------------------
774      DO i = 1, klon
775         IF (.NOT.ldcum(i)) ktype(i)=0
776      ENDDO
777c
778      DO k=1,klev
779      DO i = 1, klon
780         plu(i,k)=0.
781         pmfu(i,k)=0.
782         pmfus(i,k)=0.
783         pmfuq(i,k)=0.
784         pmful(i,k)=0.
785         plude(i,k)=0.
786         pdmfup(i,k)=0.
787         IF(.NOT.ldcum(i).OR.ktype(i).EQ.3) klab(i,k)=0
788         IF(.NOT.ldcum(i).AND.paph(i,k).LT.4.E4) kctop0(i)=k
789      ENDDO
790      ENDDO
791c
792      DO i = 1, klon
793      IF (ldland(i)) THEN
794         zdland(i)=3.0E4
795         zdphi=pgeoh(i,kctop0(i))-pgeoh(i,kcbot(i))
796         IF (ptu(i,kctop0(i)).GE.ztglace) zdland(i)=zdphi
797         zdland(i)=MAX(3.0E4,zdland(i))
798         zdland(i)=MIN(5.0E4,zdland(i))
799      ENDIF
800      ENDDO
801C
802C Initialiser les valeurs au niveau d'ascendance
803C
804      DO i = 1, klon
805         kctop(i) = klev-1
806         IF (.NOT.ldcum(i)) THEN
807            kcbot(i) = klev-1
808            pmfub(i) = 0.
809            pqu(i,klev) = 0.
810         ENDIF
811         pmfu(i,klev) = pmfub(i)
812         pmfus(i,klev) = pmfub(i)*(RCPD*ptu(i,klev)+pgeoh(i,klev))
813         pmfuq(i,klev) = pmfub(i)*pqu(i,klev)
814      ENDDO
815c
816      DO i = 1, klon
817         ldcum(i) = .FALSE.
818      ENDDO
819C----------------------------------------------------------------------
820C  DO ASCENT: SUBCLOUD LAYER (klab=1) ,CLOUDS (klab=2)
821C  BY DOING FIRST DRY-ADIABATIC ASCENT AND THEN
822C  BY ADJUSTING T,Q AND L ACCORDINGLY IN *flxadjtq*,
823C  THEN CHECK FOR BUOYANCY AND SET FLAGS ACCORDINGLY
824C----------------------------------------------------------------------
825      DO 480 k = klev-1,3,-1
826c
827      IF (LMFMID .AND. k.LT.klev-1 .AND. k.GT.klev/2) THEN
828         DO i = 1, klon
829         IF (.NOT.ldcum(i) .AND. klab(i,k+1).EQ.0 .AND.
830     .       pqen(i,k).GT.0.9*pqsen(i,k)) THEN
831            ptu(i,k+1) = pten(i,k) +(pgeo(i,k)-pgeoh(i,k+1))/RCPD
832            pqu(i,k+1) = pqen(i,k)
833            plu(i,k+1) = 0.0
834            zzzmb = MAX(CMFCMIN, -pvervel(i,k)/RG)
835            zmfmax = (paph(i,k)-paph(i,k-1))/(RG*pdtime)
836            pmfub(i) = MIN(zzzmb,zmfmax)
837            pmfu(i,k+1) = pmfub(i)
838            pmfus(i,k+1) = pmfub(i)*(RCPD*ptu(i,k+1)+pgeoh(i,k+1))
839            pmfuq(i,k+1) = pmfub(i)*pqu(i,k+1)
840            pmful(i,k+1) = 0.0
841            pdmfup(i,k+1) = 0.0
842            kcbot(i) = k
843            klab(i,k+1) = 1
844            ktype(i) = 3
845            pentr(i) = ENTRMID
846         ENDIF
847         ENDDO
848      ENDIF
849c
850      is = 0
851      DO i = 1, klon
852         is = is + klab(i,k+1)
853         IF (klab(i,k+1) .EQ. 0) klab(i,k) = 0
854         llflag(i) = .FALSE.
855         IF (klab(i,k+1) .GT. 0) llflag(i) = .TRUE.
856      ENDDO
857      IF (is .EQ. 0) GOTO 480
858c
859c calculer le taux d'entrainement et de detrainement
860c
861      DO i = 1, klon
862         pen_u(i,k) = 0.0
863         pde_u(i,k) = 0.0
864         zrho(i)=paph(i,k+1)/(RD*ptenh(i,k+1))
865         zpbot(i)=paph(i,kcbot(i))
866         zptop(i)=paph(i,kctop0(i))
867      ENDDO
868c
869      DO 125 i = 1, klon
870      IF(ldcum(i)) THEN
871         zdprho=(paph(i,k+1)-paph(i,k))/(RG*zrho(i))
872         zentr=pentr(i)*pmfu(i,k+1)*zdprho
873         llo1=k.LT.kcbot(i)
874         IF(llo1) pde_u(i,k)=zentr
875         zpmid=0.5*(zpbot(i)+zptop(i))
876         llo2=llo1.AND.ktype(i).EQ.2.AND.
877     .        (zpbot(i)-paph(i,k).LT.0.2E5.OR.
878     .         paph(i,k).GT.zpmid)
879         IF(llo2) pen_u(i,k)=zentr
880         llo2=llo1.AND.(ktype(i).EQ.1.OR.ktype(i).EQ.3).AND.
881     .        (k.GE.MAX(klwmin(i),kctop0(i)+2).OR.pap(i,k).GT.zpmid)
882         IF(llo2) pen_u(i,k)=zentr
883         llo1=pen_u(i,k).GT.0..AND.(ktype(i).EQ.1.OR.ktype(i).EQ.2)
884         IF(llo1) THEN
885            zentr=zentr*(1.+3.*(1.-MIN(1.,(zpbot(i)-pap(i,k))/1.5E4)))
886            pen_u(i,k)=pen_u(i,k)*(1.+3.*(1.-MIN(1.,
887     .                 (zpbot(i)-pap(i,k))/1.5E4)))
888            pde_u(i,k)=pde_u(i,k)*(1.+3.*(1.-MIN(1.,
889     .                 (zpbot(i)-pap(i,k))/1.5E4)))
890         ENDIF
891         IF(llo2.AND.pqenh(i,k+1).GT.1.E-5)
892     .   pen_u(i,k)=zentr+MAX(pqte(i,k),0.)/pqenh(i,k+1)*
893     .              zrho(i)*zdprho
894      ENDIF
895  125 CONTINUE
896c
897C----------------------------------------------------------------------
898c DO ADIABATIC ASCENT FOR ENTRAINING/DETRAINING PLUME
899C----------------------------------------------------------------------
900c
901      DO 420 i = 1, klon
902      IF (llflag(i)) THEN
903         IF (k.LT.kcbot(i)) THEN
904            zmftest = pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)
905            zmfmax = MIN(zmftest,(paph(i,k)-paph(i,k-1))/(RG*pdtime))
906            pen_u(i,k)=MAX(pen_u(i,k)-MAX(0.0,zmftest-zmfmax),0.0)
907         ENDIF
908         pde_u(i,k)=MIN(pde_u(i,k),0.75*pmfu(i,k+1))
909c calculer le flux de masse du niveau k a partir de celui du k+1
910         pmfu(i,k)=pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)
911c calculer les valeurs Su, Qu et l du niveau k dans le panache montant
912         zqeen=pqenh(i,k+1)*pen_u(i,k)
913         zseen=(RCPD*ptenh(i,k+1)+pgeoh(i,k+1))*pen_u(i,k)
914         zscde=(RCPD*ptu(i,k+1)+pgeoh(i,k+1))*pde_u(i,k)
915         zqude=pqu(i,k+1)*pde_u(i,k)
916         plude(i,k)=plu(i,k+1)*pde_u(i,k)
917         zmfusk=pmfus(i,k+1)+zseen-zscde
918         zmfuqk=pmfuq(i,k+1)+zqeen-zqude
919         zmfulk=pmful(i,k+1)    -plude(i,k)
920         plu(i,k)=zmfulk*(1./MAX(CMFCMIN,pmfu(i,k)))
921         pqu(i,k)=zmfuqk*(1./MAX(CMFCMIN,pmfu(i,k)))
922         ptu(i,k)=(zmfusk*(1./MAX(CMFCMIN,pmfu(i,k)))-
923     1               pgeoh(i,k))/RCPD
924         ptu(i,k)=MAX(100.,ptu(i,k))
925         ptu(i,k)=MIN(400.,ptu(i,k))
926         zqold(i)=pqu(i,k)
927      ELSE
928         zqold(i)=0.0
929      ENDIF
930  420 CONTINUE
931c
932C----------------------------------------------------------------------
933c DO CORRECTIONS FOR MOIST ASCENT BY ADJUSTING T,Q AND L
934C----------------------------------------------------------------------
935c
936      icall = 1
937      CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
938C
939      DO 440 i = 1, klon
940      IF(llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN
941         klab(i,k) = 2
942         plu(i,k) = plu(i,k)+zqold(i)-pqu(i,k)
943         zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-
944     .          ptenh(i,k)*(1.+RETV*pqenh(i,k))
945         IF (klab(i,k+1).EQ.1) zbuo=zbuo+0.5
946         IF (zbuo.GT.0..AND.pmfu(i,k).GE.0.1*pmfub(i)) THEN
947            kctop(i) = k
948            ldcum(i) = .TRUE.
949            zdnoprc = 1.5E4
950            IF (ldland(i)) zdnoprc = zdland(i)
951            zprcon = CPRCON
952            IF ((zpbot(i)-paph(i,k)).LT.zdnoprc) zprcon = 0.0
953            zlnew=plu(i,k)/(1.+zprcon*(pgeoh(i,k)-pgeoh(i,k+1)))
954            pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))
955            plu(i,k)=zlnew
956         ELSE
957            klab(i,k)=0
958            pmfu(i,k)=0.
959         ENDIF
960      ENDIF
961  440 CONTINUE
962      DO 455 i = 1, klon
963      IF (llflag(i)) THEN
964         pmful(i,k)=plu(i,k)*pmfu(i,k)
965         pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)
966         pmfuq(i,k)=pqu(i,k)*pmfu(i,k)
967      ENDIF
968  455 CONTINUE
969C
970  480 CONTINUE
971C----------------------------------------------------------------------
972C DETERMINE CONVECTIVE FLUXES ABOVE NON-BUOYANCY LEVEL
973C    (NOTE: CLOUD VARIABLES LIKE T,Q AND L ARE NOT
974C           AFFECTED BY DETRAINMENT AND ARE ALREADY KNOWN
975C           FROM PREVIOUS CALCULATIONS ABOVE)
976C----------------------------------------------------------------------
977      DO i = 1, klon
978         IF (kctop(i).EQ.klev-1) ldcum(i) = .FALSE.
979         kcbot(i) = MAX(kcbot(i),kctop(i))
980      ENDDO
981c
982      is = 0
983      DO i = 1, klon
984         if (ldcum(i)) is = is + 1
985      ENDDO
986      kcum = is
987      IF (is.EQ.0) GOTO 800
988c
989      DO 530 i = 1, klon
990      IF (ldcum(i)) THEN
991         k=kctop(i)-1
992         pde_u(i,k)=(1.-CMFCTOP)*pmfu(i,k+1)
993         plude(i,k)=pde_u(i,k)*plu(i,k+1)
994         pmfu(i,k)=pmfu(i,k+1)-pde_u(i,k)
995         zlnew=plu(i,k)
996         pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))
997         plu(i,k)=zlnew
998         pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)
999         pmfuq(i,k)=pqu(i,k)*pmfu(i,k)
1000         pmful(i,k)=plu(i,k)*pmfu(i,k)
1001         plude(i,k-1)=pmful(i,k)
1002      ENDIF
1003  530 CONTINUE
1004C
1005  800 CONTINUE
1006      RETURN
1007      END
1008      SUBROUTINE flxflux(pdtime, pqen, pqsen, ptenh, pqenh, pap
1009     .  ,  paph, ldland, pgeoh, kcbot, kctop, lddraf, kdtop
1010     .  ,  ktype, ldcum, pmfu, pmfd, pmfus, pmfds
1011     .  ,  pmfuq, pmfdq, pmful, plude, pdmfup, pdmfdp
1012     .  ,  pten, prfl, psfl, pdpmel, ktopm2
1013     .  ,  pmflxr, pmflxs)
1014      IMPLICIT none
1015C----------------------------------------------------------------------
1016C THIS ROUTINE DOES THE FINAL CALCULATION OF CONVECTIVE
1017C FLUXES IN THE CLOUD LAYER AND IN THE SUBCLOUD LAYER
1018C----------------------------------------------------------------------
1019#include "dimensions.h"
1020#include "dimphy.h"
1021#include "YOMCST.h"
1022#include "YOETHF.h"
1023#include "YOECUMF.h"
1024C
1025      REAL cevapcu(klev)
1026C     -----------------------------------------------------------------
1027      REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)
1028      REAL pten(klon,klev), ptenh(klon,klev)
1029      REAL paph(klon,klev+1), pgeoh(klon,klev)
1030c
1031      REAL pap(klon,klev)
1032      REAL ztmsmlt, zdelta, zqsat
1033C
1034      REAL pmfu(klon,klev), pmfus(klon,klev)
1035      REAL pmfd(klon,klev), pmfds(klon,klev)
1036      REAL pmfuq(klon,klev), pmful(klon,klev)
1037      REAL pmfdq(klon,klev)
1038      REAL plude(klon,klev)
1039      REAL pdmfup(klon,klev), pdpmel(klon,klev)
1040      REAL pdmfdp(klon,klev)
1041      REAL prfl(klon), psfl(klon)
1042      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
1043      INTEGER  kcbot(klon), kctop(klon), ktype(klon)
1044      LOGICAL  ldland(klon), ldcum(klon)
1045      INTEGER k, i
1046      REAL zcons1, zcons2, zcucov, ztmelp2
1047      REAL pdtime, zdp, zzp, zfac, zsnmlt, zrfl, zrnew
1048      REAL zrmin, zrfln, zdrfl
1049      REAL zpds, zpdr, zdenom
1050      INTEGER ktopm2, itop, ikb
1051c
1052      LOGICAL lddraf(klon)
1053      INTEGER kdtop(klon)
1054c
1055#include "FCTTRE.h"
1056c
1057      DO 101 k=1,klev
1058      CEVAPCU(k)=1.93E-6*261.*SQRT(1.E3/(38.3*0.293)
1059     1 *SQRT(0.5*(paph(1,k)+paph(1,k+1))/paph(1,klev+1)) ) * 0.5/RG
1060 101  CONTINUE
1061c
1062c SPECIFY CONSTANTS
1063c
1064      zcons1 = RCPD/(RLMLT*RG*pdtime)
1065      zcons2 = 1./(RG*pdtime)
1066      zcucov = 0.05
1067      ztmelp2 = RTT + 2.
1068c
1069c DETERMINE FINAL CONVECTIVE FLUXES
1070c
1071      itop=klev
1072      DO 110 i = 1, klon
1073         itop=MIN(itop,kctop(i))
1074         IF (.NOT.ldcum(i) .OR. kdtop(i).LT.kctop(i)) lddraf(i)=.FALSE.
1075         IF(.NOT.ldcum(i)) ktype(i)=0
1076  110 CONTINUE
1077c
1078      ktopm2=itop-2
1079      DO 120 k=ktopm2,klev
1080      DO 115 i = 1, klon
1081      IF(ldcum(i).AND.k.GE.kctop(i)-1) THEN
1082         pmfus(i,k)=pmfus(i,k)-pmfu(i,k)*
1083     .                (RCPD*ptenh(i,k)+pgeoh(i,k))
1084         pmfuq(i,k)=pmfuq(i,k)-pmfu(i,k)*pqenh(i,k)
1085         zdp = 1.5E4
1086         IF ( ldland(i) ) zdp = 3.E4
1087c
1088c        l'eau liquide detrainee est precipitee quand certaines
1089c        conditions sont reunies (sinon, elle est consideree
1090c        evaporee dans l'environnement)
1091c
1092         IF(paph(i,kcbot(i))-paph(i,kctop(i)).GE.zdp.AND.
1093     .      pqen(i,k-1).GT.0.8*pqsen(i,k-1))
1094     .      pdmfup(i,k-1)=pdmfup(i,k-1)+plude(i,k-1)
1095c
1096         IF(lddraf(i).AND.k.GE.kdtop(i)) THEN
1097            pmfds(i,k)=pmfds(i,k)-pmfd(i,k)*
1098     .                   (RCPD*ptenh(i,k)+pgeoh(i,k))
1099            pmfdq(i,k)=pmfdq(i,k)-pmfd(i,k)*pqenh(i,k)
1100         ELSE
1101            pmfd(i,k)=0.
1102            pmfds(i,k)=0.
1103            pmfdq(i,k)=0.
1104            pdmfdp(i,k-1)=0.
1105         END IF
1106      ELSE
1107         pmfu(i,k)=0.
1108         pmfus(i,k)=0.
1109         pmfuq(i,k)=0.
1110         pmful(i,k)=0.
1111         pdmfup(i,k-1)=0.
1112         plude(i,k-1)=0.
1113         pmfd(i,k)=0.
1114         pmfds(i,k)=0.
1115         pmfdq(i,k)=0.
1116         pdmfdp(i,k-1)=0.
1117      ENDIF
1118  115 CONTINUE
1119  120 CONTINUE
1120c
1121      DO 130 k=ktopm2,klev
1122      DO 125 i = 1, klon
1123      IF(ldcum(i).AND.k.GT.kcbot(i)) THEN
1124         ikb=kcbot(i)
1125         zzp=((paph(i,klev+1)-paph(i,k))/
1126     .        (paph(i,klev+1)-paph(i,ikb)))
1127         IF (ktype(i).EQ.3) zzp = zzp**2
1128         pmfu(i,k)=pmfu(i,ikb)*zzp
1129         pmfus(i,k)=pmfus(i,ikb)*zzp
1130         pmfuq(i,k)=pmfuq(i,ikb)*zzp
1131         pmful(i,k)=pmful(i,ikb)*zzp
1132      ENDIF
1133  125 CONTINUE
1134  130 CONTINUE
1135c
1136c CALCULATE RAIN/SNOW FALL RATES
1137c CALCULATE MELTING OF SNOW
1138c CALCULATE EVAPORATION OF PRECIP
1139c
1140      DO k = 1, klev+1
1141      DO i = 1, klon
1142         pmflxr(i,k) = 0.0
1143         pmflxs(i,k) = 0.0
1144      ENDDO
1145      ENDDO
1146      DO k = ktopm2, klev
1147      DO i = 1, klon
1148      IF (ldcum(i)) THEN
1149         IF (pmflxs(i,k).GT.0.0 .AND. pten(i,k).GT.ztmelp2) THEN
1150            zfac=zcons1*(paph(i,k+1)-paph(i,k))
1151            zsnmlt=MIN(pmflxs(i,k),zfac*(pten(i,k)-ztmelp2))
1152            pdpmel(i,k)=zsnmlt
1153            ztmsmlt=pten(i,k)-zsnmlt/zfac
1154            zdelta=MAX(0.,SIGN(1.,RTT-ztmsmlt))
1155            zqsat=R2ES*FOEEW(ztmsmlt, zdelta) / pap(i,k)
1156            zqsat=MIN(0.5,zqsat)
1157            zqsat=zqsat/(1.-RETV  *zqsat)
1158            pqsen(i,k) = zqsat
1159         ENDIF
1160         IF (pten(i,k).GT.RTT) THEN
1161           pmflxr(i,k+1)=pmflxr(i,k)+pdmfup(i,k)+pdmfdp(i,k)+pdpmel(i,k)
1162         ELSE
1163           pmflxs(i,k+1)=pmflxs(i,k)+pdmfup(i,k)+pdmfdp(i,k)-pdpmel(i,k)
1164         ENDIF
1165c        si la precipitation est negative, on ajuste le plux du
1166c        panache descendant pour eliminer la negativite
1167         IF ((pmflxr(i,k+1)+pmflxs(i,k+1)).LT.0.0) THEN
1168            pdmfdp(i,k) = -pmflxr(i,k)-pmflxs(i,k)-pdmfup(i,k)
1169            pmflxr(i,k+1) = 0.0
1170            pmflxs(i,k+1) = 0.0
1171            pdpmel(i,k) = 0.0
1172         ENDIF
1173      ENDIF
1174      ENDDO
1175      ENDDO
1176c
1177      DO k = ktopm2, klev
1178      DO i = 1, klon
1179      IF (ldcum(i) .AND. k.GE.kcbot(i)) THEN
1180         zrfl = pmflxr(i,k) + pmflxs(i,k)
1181         IF (zrfl.GT.1.0E-20) THEN
1182            zrnew=(MAX(0.,SQRT(zrfl/zcucov)-
1183     .            CEVAPCU(k)*(paph(i,k+1)-paph(i,k))*
1184     .            MAX(0.,pqsen(i,k)-pqen(i,k))))**2*zcucov
1185            zrmin=zrfl-zcucov*MAX(0.,0.8*pqsen(i,k)-pqen(i,k))
1186     .            *zcons2*(paph(i,k+1)-paph(i,k))
1187            zrnew=MAX(zrnew,zrmin)
1188            zrfln=MAX(zrnew,0.)
1189            zdrfl=MIN(0.,zrfln-zrfl)
1190            zdenom=1.0/MAX(1.0E-20,pmflxr(i,k)+pmflxs(i,k))
1191            IF (pten(i,k).GT.RTT) THEN
1192               zpdr = pdmfdp(i,k)
1193               zpds = 0.0
1194            ELSE
1195               zpdr = 0.0
1196               zpds = pdmfdp(i,k)
1197            ENDIF
1198            pmflxr(i,k+1) = pmflxr(i,k) + zpdr + pdpmel(i,k)
1199     .                    + zdrfl*pmflxr(i,k)*zdenom
1200            pmflxs(i,k+1) = pmflxs(i,k) + zpds - pdpmel(i,k)
1201     .                    + zdrfl*pmflxs(i,k)*zdenom
1202            pdmfup(i,k) = pdmfup(i,k) + zdrfl
1203         ELSE
1204            pmflxr(i,k+1) = 0.0
1205            pmflxs(i,k+1) = 0.0
1206            pdmfdp(i,k) = 0.0
1207            pdpmel(i,k) = 0.0
1208         ENDIF
1209      ENDIF
1210      ENDDO
1211      ENDDO
1212c
1213      DO 210 i = 1, klon
1214         prfl(i) = pmflxr(i,klev+1)
1215         psfl(i) = pmflxs(i,klev+1)
1216  210 CONTINUE
1217c
1218      RETURN
1219      END
1220      SUBROUTINE flxdtdq(pdtime, ktopm2, paph, ldcum, pten
1221     .  ,  pmfus, pmfds, pmfuq, pmfdq, pmful, pdmfup, pdmfdp
1222     .  ,  pdpmel, dt_con, dq_con)
1223      IMPLICIT none
1224c----------------------------------------------------------------------
1225c calculer les tendances T et Q
1226c----------------------------------------------------------------------
1227#include "dimensions.h"
1228#include "dimphy.h"
1229#include "YOMCST.h"
1230#include "YOETHF.h"
1231#include "YOECUMF.h"
1232C     -----------------------------------------------------------------
1233      LOGICAL  llo1
1234C
1235      REAL pten(klon,klev), paph(klon,klev+1)
1236      REAL pmfus(klon,klev), pmfuq(klon,klev), pmful(klon,klev)
1237      REAL pmfds(klon,klev), pmfdq(klon,klev)
1238      REAL pdmfup(klon,klev)
1239      REAL pdmfdp(klon,klev)
1240      REAL pdpmel(klon,klev)
1241      LOGICAL ldcum(klon)
1242      REAL dt_con(klon,klev), dq_con(klon,klev)
1243c
1244      INTEGER ktopm2
1245      REAL pdtime
1246c
1247      INTEGER i, k
1248      REAL zalv, zdtdt, zdqdt
1249c
1250      DO 210 k=ktopm2,klev-1
1251      DO 220 i = 1, klon
1252      IF (ldcum(i)) THEN
1253         llo1 = (pten(i,k)-RTT).GT.0.
1254         zalv = RLSTT
1255         IF (llo1) zalv = RLVTT
1256         zdtdt=RG/(paph(i,k+1)-paph(i,k))/RCPD
1257     .        *(pmfus(i,k+1)-pmfus(i,k)
1258     .         +pmfds(i,k+1)-pmfds(i,k)
1259     .          -RLMLT*pdpmel(i,k)
1260     .          -zalv*(pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1261     .         )
1262         dt_con(i,k)=zdtdt
1263         zdqdt=RG/(paph(i,k+1)-paph(i,k))
1264     .        *(pmfuq(i,k+1)-pmfuq(i,k)
1265     .         +pmfdq(i,k+1)-pmfdq(i,k)
1266     .          +pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))
1267         dq_con(i,k)=zdqdt
1268      ENDIF
1269  220 CONTINUE
1270  210 CONTINUE
1271C
1272      k = klev
1273      DO 230 i = 1, klon
1274      IF (ldcum(i)) THEN
1275         llo1 = (pten(i,k)-RTT).GT.0.
1276         zalv = RLSTT
1277         IF (llo1) zalv = RLVTT
1278         zdtdt=-RG/(paph(i,k+1)-paph(i,k))/RCPD
1279     .         *(pmfus(i,k)+pmfds(i,k)+RLMLT*pdpmel(i,k)
1280     .           -zalv*(pmful(i,k)+pdmfup(i,k)+pdmfdp(i,k)))
1281         dt_con(i,k)=zdtdt
1282         zdqdt=-RG/(paph(i,k+1)-paph(i,k))
1283     .            *(pmfuq(i,k)+pmfdq(i,k)+pmful(i,k)
1284     .             +pdmfup(i,k)+pdmfdp(i,k))
1285         dq_con(i,k)=zdqdt
1286      ENDIF
1287  230 CONTINUE
1288C
1289      RETURN
1290      END
1291      SUBROUTINE flxdlfs(ptenh, pqenh, pgeoh, paph, ptu, pqu,
1292     .     ldcum, kcbot, kctop, pmfub, prfl, ptd, pqd,
1293     .     pmfd, pmfds, pmfdq, pdmfdp, kdtop, lddraf)
1294      IMPLICIT none
1295C
1296C----------------------------------------------------------------------
1297C THIS ROUTINE CALCULATES LEVEL OF FREE SINKING FOR
1298C CUMULUS DOWNDRAFTS AND SPECIFIES T,Q,U AND V VALUES
1299C
1300C TO PRODUCE LFS-VALUES FOR CUMULUS DOWNDRAFTS
1301C FOR MASSFLUX CUMULUS PARAMETERIZATION
1302C
1303C INPUT ARE ENVIRONMENTAL VALUES OF T,Q,U,V,P,PHI
1304C AND UPDRAFT VALUES T,Q,U AND V AND ALSO
1305C CLOUD BASE MASSFLUX AND CU-PRECIPITATION RATE.
1306C IT RETURNS T,Q,U AND V VALUES AND MASSFLUX AT LFS.
1307C
1308C CHECK FOR NEGATIVE BUOYANCY OF AIR OF EQUAL PARTS OF
1309C MOIST ENVIRONMENTAL AIR AND CLOUD AIR.
1310C----------------------------------------------------------------------
1311#include "dimensions.h"
1312#include "dimphy.h"
1313#include "YOMCST.h"
1314#include "YOETHF.h"
1315#include "YOECUMF.h"
1316C
1317      REAL ptenh(klon,klev)
1318      REAL pqenh(klon,klev)
1319      REAL pgeoh(klon,klev), paph(klon,klev+1)
1320      REAL ptu(klon,klev), pqu(klon,klev)
1321      REAL pmfub(klon)
1322      REAL prfl(klon)
1323C
1324      REAL ptd(klon,klev), pqd(klon,klev)
1325      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1326      REAL pdmfdp(klon,klev)
1327      INTEGER  kcbot(klon), kctop(klon), kdtop(klon)
1328      LOGICAL  ldcum(klon), lddraf(klon)
1329C
1330      REAL ztenwb(klon,klev), zqenwb(klon,klev), zcond(klon)
1331      REAL zttest, zqtest, zbuo, zmftop
1332      LOGICAL  llo2(klon)
1333      INTEGER i, k, is, icall
1334C----------------------------------------------------------------------
1335      DO i= 1, klon
1336         lddraf(i)=.FALSE.
1337         kdtop(i)=klev+1
1338      ENDDO
1339C
1340C----------------------------------------------------------------------
1341C DETERMINE LEVEL OF FREE SINKING BY
1342C DOING A SCAN FROM TOP TO BASE OF CUMULUS CLOUDS
1343C
1344C FOR EVERY POINT AND PROCEED AS FOLLOWS:
1345C     (1) DETEMINE WET BULB ENVIRONMENTAL T AND Q
1346C     (2) DO MIXING WITH CUMULUS CLOUD AIR
1347C     (3) CHECK FOR NEGATIVE BUOYANCY
1348C
1349C THE ASSUMPTION IS THAT AIR OF DOWNDRAFTS IS MIXTURE
1350C OF 50% CLOUD AIR + 50% ENVIRONMENTAL AIR AT WET BULB
1351C TEMPERATURE (I.E. WHICH BECAME SATURATED DUE TO
1352C EVAPORATION OF RAIN AND CLOUD WATER)
1353C----------------------------------------------------------------------
1354C
1355      DO 290 k = 3, klev-3
1356C
1357      is=0
1358      DO 212 i= 1, klon
1359         ztenwb(i,k)=ptenh(i,k)
1360         zqenwb(i,k)=pqenh(i,k)
1361         llo2(i) = ldcum(i).AND.prfl(i).GT.0.
1362     .             .AND..NOT.lddraf(i)
1363     .             .AND.(k.LT.kcbot(i).AND.k.GT.kctop(i))
1364         IF ( llo2(i) ) is = is + 1
1365  212 CONTINUE
1366      IF(is.EQ.0) GO TO 290
1367C
1368      icall=2
1369      CALL flxadjtq(paph(1,k), ztenwb(1,k), zqenwb(1,k), llo2, icall)
1370C
1371C----------------------------------------------------------------------
1372C DO MIXING OF CUMULUS AND ENVIRONMENTAL AIR
1373C AND CHECK FOR NEGATIVE BUOYANCY.
1374C THEN SET VALUES FOR DOWNDRAFT AT LFS.
1375C----------------------------------------------------------------------
1376      DO 222 i= 1, klon
1377      IF (llo2(i)) THEN
1378         zttest=0.5*(ptu(i,k)+ztenwb(i,k))
1379         zqtest=0.5*(pqu(i,k)+zqenwb(i,k))
1380         zbuo=zttest*(1.+RETV*zqtest)-
1381     .        ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1382         zcond(i)=pqenh(i,k)-zqenwb(i,k)
1383         zmftop=-CMFDEPS*pmfub(i)
1384         IF (zbuo.LT.0..AND.prfl(i).GT.10.*zmftop*zcond(i)) THEN
1385            kdtop(i)=k
1386            lddraf(i)=.TRUE.
1387            ptd(i,k)=zttest
1388            pqd(i,k)=zqtest
1389            pmfd(i,k)=zmftop
1390            pmfds(i,k)=pmfd(i,k)*(RCPD*ptd(i,k)+pgeoh(i,k))
1391            pmfdq(i,k)=pmfd(i,k)*pqd(i,k)
1392            pdmfdp(i,k-1)=-0.5*pmfd(i,k)*zcond(i)
1393            prfl(i)=prfl(i)+pdmfdp(i,k-1)
1394         ENDIF
1395      ENDIF
1396  222 CONTINUE
1397c
1398  290 CONTINUE
1399C
1400      RETURN
1401      END
1402      SUBROUTINE flxddraf(ptenh, pqenh, pgeoh, paph, prfl,
1403     .           ptd, pqd, pmfd, pmfds, pmfdq, pdmfdp,
1404     .           lddraf, pen_d, pde_d)
1405      IMPLICIT none
1406C
1407C----------------------------------------------------------------------
1408C          THIS ROUTINE CALCULATES CUMULUS DOWNDRAFT DESCENT
1409C
1410C          TO PRODUCE THE VERTICAL PROFILES FOR CUMULUS DOWNDRAFTS
1411C          (I.E. T,Q,U AND V AND FLUXES)
1412C
1413C          INPUT IS T,Q,P,PHI,U,V AT HALF LEVELS.
1414C          IT RETURNS FLUXES OF S,Q AND EVAPORATION RATE
1415C          AND U,V AT LEVELS WHERE DOWNDRAFT OCCURS
1416C
1417C          CALCULATE MOIST DESCENT FOR ENTRAINING/DETRAINING PLUME BY
1418C          A) MOVING AIR DRY-ADIABATICALLY TO NEXT LEVEL BELOW AND
1419C          B) CORRECTING FOR EVAPORATION TO OBTAIN SATURATED STATE.
1420C
1421C----------------------------------------------------------------------
1422#include "dimensions.h"
1423#include "dimphy.h"
1424#include "YOMCST.h"
1425#include "YOETHF.h"
1426#include "YOECUMF.h"
1427C
1428      REAL ptenh(klon,klev), pqenh(klon,klev)
1429      REAL pgeoh(klon,klev), paph(klon,klev+1)
1430C
1431      REAL ptd(klon,klev), pqd(klon,klev)
1432      REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)
1433      REAL pdmfdp(klon,klev)
1434      REAL prfl(klon)
1435      LOGICAL lddraf(klon)
1436C
1437      REAL pen_d(klon,klev), pde_d(klon,klev), zcond(klon)
1438      LOGICAL llo2(klon), llo1
1439      INTEGER i, k, is, icall, itopde
1440      REAL zentr, zseen, zqeen, zsdde, zqdde, zmfdsk, zmfdqk, zdmfdp
1441      REAL zbuo
1442C----------------------------------------------------------------------
1443C CALCULATE MOIST DESCENT FOR CUMULUS DOWNDRAFT BY
1444C       (A) CALCULATING ENTRAINMENT RATES, ASSUMING
1445C           LINEAR DECREASE OF MASSFLUX IN PBL
1446C       (B) DOING MOIST DESCENT - EVAPORATIVE COOLING
1447C           AND MOISTENING IS CALCULATED IN *flxadjtq*
1448C       (C) CHECKING FOR NEGATIVE BUOYANCY AND
1449C           SPECIFYING FINAL T,Q,U,V AND DOWNWARD FLUXES
1450C
1451      DO 180 k = 3, klev
1452c
1453      is = 0
1454      DO i = 1, klon
1455         llo2(i)=lddraf(i).AND.pmfd(i,k-1).LT.0.
1456         IF (llo2(i)) is = is + 1
1457      ENDDO
1458      IF (is.EQ.0) GOTO 180
1459c
1460      DO i = 1, klon
1461      IF (llo2(i)) THEN
1462         zentr = ENTRDD*pmfd(i,k-1)*RD*ptenh(i,k-1)/
1463     .           (RG*paph(i,k-1))*(paph(i,k)-paph(i,k-1))
1464         pen_d(i,k) = zentr
1465         pde_d(i,k) = zentr
1466      ENDIF
1467      ENDDO
1468c
1469      itopde = klev-2
1470      IF (k.GT.itopde) THEN
1471         DO i = 1, klon
1472         IF (llo2(i)) THEN
1473            pen_d(i,k)=0.
1474            pde_d(i,k)=pmfd(i,itopde)*
1475     .      (paph(i,k)-paph(i,k-1))/(paph(i,klev+1)-paph(i,itopde))
1476         ENDIF
1477         ENDDO
1478      ENDIF
1479C
1480      DO i = 1, klon
1481      IF (llo2(i)) THEN
1482         pmfd(i,k) = pmfd(i,k-1)+pen_d(i,k)-pde_d(i,k)
1483         zseen = (RCPD*ptenh(i,k-1)+pgeoh(i,k-1))*pen_d(i,k)
1484         zqeen = pqenh(i,k-1)*pen_d(i,k)
1485         zsdde = (RCPD*ptd(i,k-1)+pgeoh(i,k-1))*pde_d(i,k)
1486         zqdde = pqd(i,k-1)*pde_d(i,k)
1487         zmfdsk = pmfds(i,k-1)+zseen-zsdde
1488         zmfdqk = pmfdq(i,k-1)+zqeen-zqdde
1489         pqd(i,k) = zmfdqk*(1./MIN(-CMFCMIN,pmfd(i,k)))
1490         ptd(i,k) = (zmfdsk*(1./MIN(-CMFCMIN,pmfd(i,k)))-
1491     .               pgeoh(i,k))/RCPD
1492         ptd(i,k) = MIN(400.,ptd(i,k))
1493         ptd(i,k) = MAX(100.,ptd(i,k))
1494         zcond(i) = pqd(i,k)
1495      ENDIF
1496      ENDDO
1497C
1498      icall = 2
1499      CALL flxadjtq(paph(1,k), ptd(1,k), pqd(1,k), llo2, icall)
1500C
1501      DO i = 1, klon
1502      IF (llo2(i)) THEN
1503         zcond(i) = zcond(i)-pqd(i,k)
1504         zbuo = ptd(i,k)*(1.+RETV  *pqd(i,k))-
1505     .          ptenh(i,k)*(1.+RETV  *pqenh(i,k))
1506         llo1 = zbuo.LT.0..AND.(prfl(i)-pmfd(i,k)*zcond(i).GT.0.)
1507         IF (.not.llo1) pmfd(i,k) = 0.0
1508         pmfds(i,k) = (RCPD*ptd(i,k)+pgeoh(i,k))*pmfd(i,k)
1509         pmfdq(i,k) = pqd(i,k)*pmfd(i,k)
1510         zdmfdp = -pmfd(i,k)*zcond(i)
1511         pdmfdp(i,k-1) = zdmfdp
1512         prfl(i) = prfl(i)+zdmfdp
1513      ENDIF
1514      ENDDO
1515c
1516  180 CONTINUE
1517      RETURN
1518      END
1519      SUBROUTINE flxadjtq(pp, pt, pq, ldflag, kcall)
1520      IMPLICIT none
1521c======================================================================
1522c Objet: ajustement entre T et Q
1523c======================================================================
1524C NOTE: INPUT PARAMETER kcall DEFINES CALCULATION AS
1525C        kcall=0    ENV. T AND QS IN*CUINI*
1526C        kcall=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)
1527C        kcall=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1528C
1529#include "dimensions.h"
1530#include "dimphy.h"
1531#include "YOMCST.h"
1532C
1533      REAL pt(klon), pq(klon), pp(klon)
1534      LOGICAL ldflag(klon)
1535      INTEGER kcall
1536c
1537      REAL zcond(klon), zcond1
1538      REAL Z5alvcp, z5alscp, zalvdcp, zalsdcp
1539      REAL zdelta, zcvm5, zldcp, zqsat, zcor
1540      INTEGER is, i
1541#include "YOETHF.h"
1542#include "FCTTRE.h"
1543C
1544      z5alvcp = r5les*RLVTT/RCPD
1545      z5alscp = r5ies*RLSTT/RCPD
1546      zalvdcp = rlvtt/RCPD
1547      zalsdcp = rlstt/RCPD
1548C
1549      DO 210 i =1, klon
1550      IF (ldflag(i)) THEN
1551         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1552         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1553         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1554         zqsat = R2ES*FOEEW(pt(i),zdelta) / pp(i)
1555         zqsat = MIN(0.5,zqsat)
1556         zcor = 1./(1.-RETV*zqsat)
1557         zqsat = zqsat*zcor
1558         zcond(i) = (pq(i)-zqsat)
1559     .     / (1. + FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor))
1560         IF (kcall.EQ.1) zcond(i) = MAX(zcond(i),0.)
1561         IF (kcall.EQ.2) zcond(i) = MIN(zcond(i),0.)
1562         pt(i) = pt(i) + zldcp*zcond(i)
1563         pq(i) = pq(i) - zcond(i)
1564      ENDIF
1565  210 CONTINUE
1566C
1567      is = 0
1568      DO i =1, klon
1569         IF (zcond(i).NE.0.) is = is + 1
1570      ENDDO
1571      IF (is.EQ.0) GOTO 230
1572C
1573      DO 220 i = 1, klon
1574      IF(ldflag(i).AND.zcond(i).NE.0.) THEN
1575         zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))
1576         zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp
1577         zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp
1578         zqsat = R2ES* FOEEW(pt(i),zdelta) / pp(i)
1579         zqsat = MIN(0.5,zqsat)
1580         zcor = 1./(1.-RETV*zqsat)
1581         zqsat = zqsat*zcor
1582         zcond1 = (pq(i)-zqsat)
1583     .     / (1. + FOEDE(pt(i),zdelta,zcvm5,zqsat,zcor))
1584         pt(i) = pt(i) + zldcp*zcond1
1585         pq(i) = pq(i) - zcond1
1586      ENDIF
1587  220 CONTINUE
1588C
1589  230 CONTINUE
1590      RETURN
1591      END
1592      SUBROUTINE flxsetup
1593      IMPLICIT none
1594C
1595C     THIS ROUTINE DEFINES DISPOSABLE PARAMETERS FOR MASSFLUX SCHEME
1596C
1597#include "YOECUMF.h"
1598C
1599      ENTRPEN=1.0E-4  ! ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
1600      ENTRSCV=3.0E-4  ! ENTRAINMENT RATE FOR SHALLOW CONVECTION
1601      ENTRMID=1.0E-4  ! ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
1602      ENTRDD =2.0E-4  ! ENTRAINMENT RATE FOR DOWNDRAFTS
1603      CMFCTOP=0.33  ! RELATIVE CLOUD MASSFLUX AT LEVEL ABOVE NONBUO LEVEL
1604      CMFCMAX=1.0  ! MAXIMUM MASSFLUX VALUE ALLOWED FOR UPDRAFTS ETC
1605      CMFCMIN=1.E-10  ! MINIMUM MASSFLUX VALUE (FOR SAFETY)
1606      CMFDEPS=0.3  ! FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
1607      CPRCON =2.0E-4  ! CONVERSION FROM CLOUD WATER TO RAIN
1608      RHCDD=1.  ! RELATIVE SATURATION IN DOWNDRAFRS (NO LONGER USED)
1609c                 (FORMULATION IMPLIES SATURATION)
1610      LMFPEN = .TRUE.
1611      LMFSCV = .TRUE.
1612      LMFMID = .TRUE.
1613      LMFDD = .TRUE.
1614      LMFDUDV = .TRUE.
1615c
1616      RETURN
1617      END
Note: See TracBrowser for help on using the repository browser.