/[lmdze]/trunk/libf/phylmd/conflx.f
ViewVC logotype

Annotation of /trunk/libf/phylmd/conflx.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 38 - (hide annotations)
Thu Jan 6 17:52:19 2011 UTC (13 years, 4 months ago) by guez
File size: 56599 byte(s) 
Extracted ASCII art from "inigeom" into a separate text file in the
documentation.

"test_disvert" now creates a separate file for layer thicknesses.

Moved variables from module "yomcst" to module "suphec_m" because this
is where those variables are defined. Kept in "yomcst" only parameters
of Earth orbit. Gave the attribute "parameter" to some variables of
module "suphec_m".

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