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

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

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