Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

conflx.F90 @ 25

Last change on this file since 25 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 60.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format