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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 10 - (hide annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years, 1 month ago) by guez
File size: 56517 byte(s) 
Added NetCDF directory "/home/guez/include" in "g95.mk" and
"nag_tools.mk".

Added some "intent" attributes in "PVtheta", "advtrac", "caladvtrac",
"calfis", "diagedyn", "dissip", "vlspltqs", "aeropt", "ajsec",
"calltherm", "clmain", "cltrac", "cltracrn", "concvl", "conema3",
"conflx", "fisrtilp", "newmicro", "nuage", "diagcld1", "diagcld2",
"drag_noro", "lift_noro", "SUGWD", "physiq", "phytrac", "radlwsw", "thermcell".

Removed the case "ierr == 0" in "abort_gcm"; moved call to "histclo"
and messages for end of run from "abort_gcm" to "gcm"; replaced call
to "abort_gcm" in "leapfrog" by exit from outer loop.

In "calfis": removed argument "pp" and variable "unskap"; changed
"pksurcp" from scalar to rank 2; use "pressure_var"; rewrote
computation of "zplev", "zplay", "ztfi", "pcvgt" using "dyn_phy";
added computation of "pls".

Removed unused variable in "dynredem0".

In "exner_hyb": changed "dellta" from scalar to rank 1; replaced call
to "ssum" by call to "sum"; removed variables "xpn" and "xps";
replaced some loops by array expressions.

In "leapfrog": use "pressure_var"; deleted variables "p", "longcles".

Converted common blocks "YOECUMF", "YOEGWD" to modules.

Removed argument "pplay" in "cvltr", "diagetpq", "nflxtr".

Created module "raddimlw" from include file "raddimlw.h".

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