Context Navigation

conlmd.F @ 1024

Last change on this file since 1024 was 42, checked in by lmdz, 25 years ago
Ajustement de parametres et introduction de modifs de S. Cheinet. L.Li Lf
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 69.7 KB

Line
1	SUBROUTINE conlmd (dtime, paprs, pplay, t, q, conv_q,
2	s d_t, d_q, rain, snow, ibas, itop)
3	IMPLICIT none
4	c======================================================================
5	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
6	c Objet: Schema de convection utilis'e dans le modele du LMD
7	c Ajustement humide (Manabe) + Ajustement convectif (Kuo)
8	c======================================================================
9	#include "dimensions.h"
10	#include "dimphy.h"
11	#include "YOMCST.h"
12	#include "YOETHF.h"
13	c
14	c Arguments:
15	c
16	REAL dtime ! pas d'integration (s)
17	REAL paprs(klon,klev+1) ! pression inter-couche (Pa)
18	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
19	REAL t(klon,klev) ! temperature (K)
20	REAL q(klon,klev) ! humidite specifique (kg/kg)
21	REAL conv_q(klon,klev) ! taux de convergence humidite (g/g/s)
22	c
23	REAL d_t(klon,klev) ! incrementation temperature
24	REAL d_q(klon,klev) ! incrementation humidite
25	REAL rain(klon) ! pluies (mm/s)
26	REAL snow(klon) ! neige (mm/s)
27	INTEGER ibas(klon) ! niveau du bas
28	INTEGER itop(klon) ! niveau du haut
29	c
30	LOGICAL usekuo ! utiliser convection profonde (schema Kuo)
31	PARAMETER (usekuo=.TRUE.)
32	c
33	REAL d_t_bis(klon,klev)
34	REAL d_q_bis(klon,klev)
35	REAL rain_bis(klon)
36	REAL snow_bis(klon)
37	INTEGER ibas_bis(klon)
38	INTEGER itop_bis(klon)
39	REAL d_ql(klon,klev), d_ql_bis(klon,klev)
40	REAL rneb(klon,klev), rneb_bis(klon,klev)
41	c
42	INTEGER i, k
43	REAL zlvdcp, zlsdcp, zdelta, zz, za, zb
44	c
45	ccc CALL fiajh ! ancienne version de Convection Manabe
46	CALL conman ! nouvelle version de Convection Manabe
47	e (dtime, paprs, pplay, t, q,
48	s d_t, d_q, d_ql, rneb,
49	s rain, snow, ibas, itop)
50	c
51	IF (usekuo) THEN
52	ccc CALL fiajc ! ancienne version de Convection Kuo
53	CALL conkuo ! nouvelle version de Convection Kuo
54	e (dtime, paprs, pplay, t, q, conv_q,
55	s d_t_bis, d_q_bis, d_ql_bis, rneb_bis,
56	s rain_bis, snow_bis, ibas_bis, itop_bis)
57	DO k = 1, klev
58	DO i = 1, klon
59	d_t(i,k) = d_t(i,k) + d_t_bis(i,k)
60	d_q(i,k) = d_q(i,k) + d_q_bis(i,k)
61	d_ql(i,k) = d_ql(i,k) + d_ql_bis(i,k)
62	ENDDO
63	ENDDO
64	DO i = 1, klon
65	rain(i) = rain(i) + rain_bis(i)
66	snow(i) = snow(i) + snow_bis(i)
67	ibas(i) = MIN(ibas(i),ibas_bis(i))
68	itop(i) = MAX(itop(i),itop_bis(i))
69	ENDDO
70	ENDIF
71	c
72	c L'eau liquide convective est dispersee dans l'air:
73	c
74	DO k = 1, klev
75	DO i = 1, klon
76	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q(i,k))
77	zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q(i,k))
78	zdelta = MAX(0.,SIGN(1.,RTT-t(i,k)))
79	zz = d_ql(i,k) ! re-evap. de l'eau liquide
80	zb = MAX(0.0,zz)
81	za = - MAX(0.0,zz) * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
82	d_t(i,k) = d_t(i,k) + za
83	d_q(i,k) = d_q(i,k) + zb
84	ENDDO
85	ENDDO
86	c
87	RETURN
88	END
89	SUBROUTINE conman (dtime, paprs, pplay, t, q,
90	s d_t, d_q, d_ql, rneb,
91	s rain, snow, ibas, itop)
92	IMPLICIT none
93	c======================================================================
94	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19970324
95	c Objet: ajustement humide convectif avec la possibilite de faire
96	c l'ajustement sur une fraction de la maille.
97	c Methode: On impose une distribution uniforme pour la vapeur d'eau
98	c au sein d'une maille. On applique la procedure d'ajustement
99	c successivement a la totalite, 75%, 50%, 25% et 5% de la maille
100	c jusqu'a ce que l'ajustement a lieu. J'espere que ceci augmente
101	c les activites convectives et corrige le biais "trop froid et sec"
102	c du modele.
103	c======================================================================
104	#include "dimensions.h"
105	#include "dimphy.h"
106	#include "YOMCST.h"
107	c
108	REAL dtime ! pas d'integration (s)
109	REAL t(klon,klev) ! temperature (K)
110	REAL q(klon,klev) ! humidite specifique (kg/kg)
111	REAL paprs(klon,klev+1) ! pression inter-couche (Pa)
112	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
113	c
114	REAL d_t(klon,klev) ! incrementation temperature
115	REAL d_q(klon,klev) ! incrementation humidite
116	REAL d_ql(klon,klev) ! incrementation eau liquide
117	REAL rneb(klon,klev) ! nebulosite
118	REAL rain(klon) ! pluies (mm/s)
119	REAL snow(klon) ! neige (mm/s)
120	INTEGER ibas(klon) ! niveau du bas
121	INTEGER itop(klon) ! niveau du haut
122	c
123	LOGICAL afaire(klon) ! .TRUE. implique l'ajustement
124	LOGICAL accompli(klon) ! .TRUE. si l'ajustement est effectif
125	c
126	INTEGER nb ! nombre de sous-fractions a considere
127	PARAMETER (nb=1)
128	ccc PARAMETER (nb=3)
129	c
130	REAL ratqs ! largeur de la distribution pour vapeur d'eau
131	PARAMETER (ratqs=0.05)
132	c
133	REAL w_q(klon,klev)
134	REAL w_d_t(klon,klev), w_d_q(klon,klev), w_d_ql(klon,klev)
135	REAL w_rneb(klon,klev)
136	REAL w_rain(klon), w_snow(klon)
137	INTEGER w_ibas(klon), w_itop(klon)
138	REAL zq1, zq2
139	INTEGER i, k, n
140	c
141	REAL t_coup
142	PARAMETER (t_coup=234.0)
143	REAL zdp1, zdp2
144	REAL zqs1, zqs2, zdqs1, zdqs2
145	REAL zgamdz
146	REAL zflo ! flotabilite
147	REAL zsat ! sur-saturation
148	REAL zdelta, zcor, zcvm5
149	LOGICAL imprim
150	c
151	INTEGER ncpt
152	SAVE ncpt
153	REAL frac(nb) ! valeur de la maille fractionnelle
154	SAVE frac
155	INTEGER opt_cld(nb) ! option pour le modele nuageux
156	SAVE opt_cld
157	LOGICAL appel1er
158	SAVE appel1er
159	c
160	c Fonctions thermodynamiques:
161	c
162	#include "YOETHF.h"
163	#include "FCTTRE.h"
164	c
165	DATA frac / 1.0 /
166	DATA opt_cld / 4 /
167	ccc DATA frac / 1.0, 0.50, 0.25/
168	ccc DATA opt_cld / 4, 4, 4/
169	c
170	DATA appel1er /.TRUE./
171	DATA ncpt /0/
172	c
173	IF (appel1er) THEN
174	PRINT*, 'conman, nb:', nb
175	PRINT*, 'conman, frac:', frac
176	PRINT*, 'conman, opt_cld:', opt_cld
177	appel1er = .FALSE.
178	ENDIF
179	c
180	c Initialiser les sorties a zero:
181	c
182	DO k = 1, klev
183	DO i = 1, klon
184	d_t(i,k) = 0.0
185	d_q(i,k) = 0.0
186	d_ql(i,k) = 0.0
187	rneb(i,k) = 0.0
188	ENDDO
189	ENDDO
190	DO i = 1, klon
191	ibas(i) = klev
192	itop(i) = 1
193	rain(i) = 0.0
194	snow(i) = 0.0
195	ENDDO
196	c
197	c S'il n'y a pas d'instabilite conditionnelle,
198	c pas la penne de se fatiguer:
199	c
200	DO i = 1, klon
201	afaire(i) = .FALSE.
202	ENDDO
203	DO k = 1, klev-1
204	DO i = 1, klon
205	IF (thermcep) THEN
206	zdelta=MAX(0.,SIGN(1.,RTT-t(i,k)))
207	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
208	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*q(i,k))
209	zqs1= R2ES*FOEEW(t(i,k),zdelta)/pplay(i,k)
210	zqs1=MIN(0.5,zqs1)
211	zcor=1./(1.-RETV*zqs1)
212	zqs1=zqs1*zcor
213	zdqs1 =FOEDE(t(i,k),zdelta,zcvm5,zqs1,zcor)
214	c
215	zdelta=MAX(0.,SIGN(1.,RTT-t(i,k+1)))
216	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
217	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*q(i,k+1))
218	zqs2= R2ES*FOEEW(t(i,k+1),zdelta)/pplay(i,k+1)
219	zqs2=MIN(0.5,zqs2)
220	zcor=1./(1.-RETV*zqs2)
221	zqs2=zqs2*zcor
222	zdqs2 =FOEDE(t(i,k+1),zdelta,zcvm5,zqs2,zcor)
223	ELSE
224	IF (t(i,k) .LT. t_coup) THEN
225	zqs1= qsats(t(i,k)) / pplay(i,k)
226	zdqs1= dqsats(t(i,k),zqs1)
227	c
228	zqs2= qsats(t(i,k+1)) / pplay(i,k+1)
229	zdqs2= dqsats(t(i,k+1),zqs2)
230	ELSE
231	zqs1= qsatl(t(i,k)) / pplay(i,k)
232	zdqs1= dqsatl(t(i,k),zqs1)
233	c
234	zqs2= qsatl(t(i,k+1)) / pplay(i,k+1)
235	zdqs2= dqsatl(t(i,k+1),zqs2)
236	ENDIF
237	ENDIF
238	zdp1 = paprs(i,k) - paprs(i,k+1)
239	zdp2 = paprs(i,k+1) - paprs(i,k+2)
240	zgamdz = - (pplay(i,k)-pplay(i,k+1))/paprs(i,k+1)/RCPD
241	. ( RD(t(i,k)zdp1+t(i,k+1)zdp2)/(zdp1+zdp2)
242	. +RLVTT(zqs1zdp1+zqs2*zdp2)/(zdp1+zdp2)
243	. ) / (1.0+(zdqs1zdp1+zdqs2zdp2)/(zdp1+zdp2) )
244	zflo = t(i,k) + zgamdz - t(i,k+1)
245	zsat = (q(i,k)-zqs1)zdp1 + (q(i,k+1)-zqs2)zdp2
246	IF (zflo.GT.0.0) afaire(i) = .TRUE.
247	c erreur IF (zflo.GT.0.0 .AND. zsat.GT.0.0) afaire(i) = .TRUE.
248	ENDDO
249	ENDDO
250	c
251	imprim = MOD(ncpt,48).EQ.0
252	DO 99999 n = 1, nb
253	c
254	DO k = 1, klev
255	DO i = 1, klon
256	IF (afaire(i)) THEN
257	zq1 = q(i,k) * (1.0-ratqs)
258	zq2 = q(i,k) * (1.0+ratqs)
259	w_q(i,k) = zq2 - frac(n)/2.0 * (zq2-zq1)
260	ENDIF
261	ENDDO
262	ENDDO
263	c
264	CALL conmanv (dtime, paprs, pplay, t, w_q,
265	e afaire, opt_cld(n),
266	s w_d_t, w_d_q, w_d_ql, w_rneb,
267	s w_rain, w_snow, w_ibas, w_itop,accompli,imprim)
268	DO k = 1, klev
269	DO i = 1, klon
270	IF (afaire(i) .AND. accompli(i)) THEN
271	d_t(i,k) = w_d_t(i,k) * frac(n)
272	d_q(i,k) = w_d_q(i,k) * frac(n)
273	d_ql(i,k) = w_d_ql(i,k) * frac(n)
274	IF (NINT(w_rneb(i,k)).EQ.1) rneb(i,k) = frac(n)
275	ENDIF
276	ENDDO
277	ENDDO
278	DO i = 1, klon
279	IF (afaire(i) .AND. accompli(i)) THEN
280	rain(i) = w_rain(i) * frac(n)
281	snow(i) = w_snow(i) * frac(n)
282	ibas(i) = MIN(ibas(i),w_ibas(i))
283	itop(i) = MAX(itop(i),w_itop(i))
284	ENDIF
285	ENDDO
286	DO i = 1, klon
287	IF(afaire(i) .AND. accompli(i)) afaire(i) = .FALSE.
288	ENDDO
289	c
290	99999 CONTINUE
291	c
292	ncpt = ncpt + 1
293	c
294	RETURN
295	END
296	SUBROUTINE conmanv (dtime, paprs, pplay, t, q,
297	e afaire, opt_cld,
298	s d_t, d_q, d_ql, rneb,
299	s rain, snow, ibas, itop,accompli,imprim)
300	IMPLICIT none
301	c======================================================================
302	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
303	c Objet: ajustement humide (convection proposee par Manabe).
304	c Pour une colonne verticale, il peut avoir plusieurs blocs
305	c necessitant l'ajustement. ibas est le bas du plus bas bloc
306	c et itop est le haut du plus haut bloc
307	c======================================================================
308	#include "dimensions.h"
309	#include "dimphy.h"
310	#include "YOMCST.h"
311	c
312	c Arguments:
313	c
314	REAL dtime ! pas d'integration (s)
315	REAL t(klon,klev) ! temperature (K)
316	REAL q(klon,klev) ! humidite specifique (kg/kg)
317	REAL paprs(klon,klev+1) ! pression inter-couche (Pa)
318	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
319	INTEGER opt_cld ! comment traiter l'eau liquide
320	LOGICAL afaire(klon) ! .TRUE. si le point est a faire (Input)
321	LOGICAL imprim ! .T. pour imprimer quelques diagnostiques
322	c
323	REAL d_t(klon,klev) ! incrementation temperature
324	REAL d_q(klon,klev) ! incrementation humidite
325	REAL d_ql(klon,klev) ! incrementation eau liquide
326	REAL rneb(klon,klev) ! nebulosite
327	REAL rain(klon) ! pluies (mm/s)
328	REAL snow(klon) ! neige (mm/s)
329	INTEGER ibas(klon) ! niveau du bas
330	INTEGER itop(klon) ! niveau du haut
331	LOGICAL accompli(klon) ! .TRUE. si l'ajustement a eu lieu (Output)
332	c
333	c Quelques options:
334	c
335	LOGICAL new_top ! re-calculer sommet quand re-ajustement est fait
336	PARAMETER (new_top=.FALSE.)
337	LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
338	PARAMETER (evap_prec=.TRUE.)
339	REAL coef_eva
340	PARAMETER (coef_eva=1.0E-05)
341	REAL t_coup
342	PARAMETER (t_coup=234.0)
343	REAL seuil_vap
344	PARAMETER (seuil_vap=1.0E-10)
345	LOGICAL old_tau ! implique precip nulle, si vrai.
346	PARAMETER (old_tau=.FALSE.)
347	REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
348	REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
349	PARAMETER (dpmin=0.15, tomax=0.97)
350	REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
351	PARAMETER (dpmax=0.30, tomin=0.05)
352	REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
353	PARAMETER (deep_sig=0.50, deep_to=0.05)
354	LOGICAL exigent ! implique un calcul supplementaire pour Qs
355	PARAMETER (exigent=.FALSE.)
356	c
357	INTEGER kbase
358	PARAMETER (kbase=0)
359	c
360	c Variables locales:
361	c
362	INTEGER nexpo
363	INTEGER i, k, k1min, k1max, k2min, k2max, is
364	REAL zgamdz(klon,klev-1)
365	REAL zt(klon,klev), zq(klon,klev)
366	REAL zqs(klon,klev), zdqs(klon,klev)
367	REAL zqmqsdp(klon,klev)
368	REAL ztnew(klon,klev), zqnew(klon,klev)
369	REAL zcond(klon), zvapo(klon), zrapp(klon)
370	REAL zrfl(klon), zrfln, zqev, zqevt
371	REAL zsat(klon) ! sur-saturation
372	REAL zflo(klon) ! flotabilite
373	REAL za(klon), zb(klon), zc(klon)
374	INTEGER k1(klon), k2(klon)
375	REAL zdelta, zcor, zcvm5
376	REAL delp(klon,klev)
377	LOGICAL possible(klon), todo(klon), etendre(klon)
378	LOGICAL aller(klon), todobis(klon)
379	REAL zalfa
380	INTEGER nbtodo, nbdone
381	c
382	c Fonctions thermodynamiques:
383	c
384	#include "YOETHF.h"
385	#include "FCTTRE.h"
386	c
387	DO k = 1, klev
388	DO i = 1, klon
389	delp(i,k) = paprs(i,k) - paprs(i,k+1)
390	ENDDO
391	ENDDO
392	c
393	c Initialiser les sorties a zero
394	c
395	DO k = 1, klev
396	DO i = 1, klon
397	d_t(i,k) = 0.0
398	d_q(i,k) = 0.0
399	d_ql(i,k) = 0.0
400	rneb(i,k) = 0.0
401	ENDDO
402	ENDDO
403	DO i = 1, klon
404	ibas(i) = klev
405	itop(i) = 1
406	rain(i) = 0.0
407	snow(i) = 0.0
408	accompli(i) = .FALSE.
409	ENDDO
410	c
411	c Preparations
412	c
413	DO k = 1, klev
414	DO i = 1, klon
415	IF (afaire(i)) THEN
416	zt(i,k) = t(i,k)
417	zq(i,k) = q(i,k)
418	c
419	c Calculer Qs et L/Cp*dQs/dT
420	c
421	IF (thermcep) THEN
422	zdelta=MAX(0.,SIGN(1.,RTT-zt(i,k)))
423	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
424	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i,k))
425	zqs(i,k)= R2ES*FOEEW(zt(i,k),zdelta)/pplay(i,k)
426	zqs(i,k)=MIN(0.5,zqs(i,k))
427	zcor=1./(1.-RETV*zqs(i,k))
428	zqs(i,k)=zqs(i,k)*zcor
429	zdqs(i,k) =FOEDE(zt(i,k),zdelta,zcvm5,zqs(i,k),zcor)
430	ELSE
431	IF (zt(i,k) .LT. t_coup) THEN
432	zqs(i,k)= qsats(zt(i,k)) / pplay(i,k)
433	zdqs(i,k)= dqsats(zt(i,k),zqs(i,k))
434	ELSE
435	zqs(i,k)= qsatl(zt(i,k)) / pplay(i,k)
436	zdqs(i,k)= dqsatl(zt(i,k),zqs(i,k))
437	ENDIF
438	ENDIF
439	c
440	c Calculer (q-qs)*dp
441	zqmqsdp(i,k) = (zq(i,k)-zqs(i,k)) * delp(i,k)
442	ENDIF
443	ENDDO
444	ENDDO
445	c
446	c-----zgama is the moist convective lapse rate (-dT/dz).
447	c-----zgamdz(*,k) est la difference minimale autorisee des temperatures
448	c-----entre deux couches (k et k+1), c.a.d. si T(k+1)-T(k) est inferieur
449	c-----a zgamdz(*,k), alors ces 2 couches sont instables conditionnellement
450	c
451	DO k = 1, klev-1
452	DO i = 1, klon
453	IF (afaire(i)) THEN
454	zgamdz(i,k) = - (pplay(i,k)-pplay(i,k+1))/paprs(i,k+1)/RCPD
455	. ( RD(zt(i,k)delp(i,k)+zt(i,k+1)delp(i,k+1))
456	. /(delp(i,k)+delp(i,k+1))
457	. +RLVTT(zqs(i,k)delp(i,k)+zqs(i,k+1)*delp(i,k+1))
458	. /(delp(i,k)+delp(i,k+1))
459	. ) / (1.0+(zdqs(i,k)delp(i,k)+zdqs(i,k+1)delp(i,k+1))
460	. /(delp(i,k)+delp(i,k+1)) )
461	ENDIF
462	ENDDO
463	ENDDO
464	c
465	c On cherche la presence simultanee d'instabilite conditionnelle
466	c et de sur-saturation. Sinon, pas la penne de se fatiguer:
467	c
468	DO i = 1, klon
469	possible(i) = .FALSE.
470	ENDDO
471	DO k = 2, klev
472	DO i = 1, klon
473	IF (afaire(i)) THEN
474	zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
475	zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1)
476	IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) possible(i) = .TRUE.
477	ENDIF
478	ENDDO
479	ENDDO
480	c
481	DO i = 1, klon
482	IF (possible(i)) THEN
483	k1(i) = kbase
484	k2(i) = k1(i) + 1
485	ENDIF
486	ENDDO
487	c
488	810 CONTINUE ! chercher le bas de la colonne a ajuster
489	c
490	k2min = klev
491	DO i = 1, klon
492	todo(i) = .FALSE.
493	aller(i) = .TRUE.
494	IF (possible(i)) k2min = MIN(k2min,k2(i))
495	ENDDO
496	IF (k2min.EQ.klev) GOTO 860
497	DO k = k2min, klev-1
498	DO i = 1, klon
499	IF (possible(i) .AND. k.GE.k2(i) .AND. aller(i)) THEN
500	zflo(i) = zt(i,k) + zgamdz(i,k) - zt(i,k+1)
501	zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k+1)
502	IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN
503	k1(i) = k
504	k2(i) = k+1
505	todo(i) = .TRUE.
506	aller(i) = .FALSE.
507	ENDIF
508	ENDIF
509	ENDDO
510	ENDDO
511	DO i = 1, klon
512	IF (possible(i).AND.aller(i)) THEN
513	todo(i) = .FALSE.
514	k1(i) = klev
515	k2(i) = klev
516	ENDIF
517	ENDDO
518	c
519	CCC DO i = 1, klon
520	CCC IF (possible(i)) THEN
521	CCC 811 k2(i) = k2(i) + 1
522	CCC IF (k2(i) .GT. klev) THEN
523	CCC todo(i) = .FALSE.
524	CCC GOTO 812
525	CCC ENDIF
526	CCC k = k2(i)
527	CCC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
528	CCC zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1)
529	CCC IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 811
530	CCC k1(i) = k2(i) - 1
531	CCC todo(i) = .TRUE.
532	CCC ENDIF
533	CCC 812 CONTINUE
534	CCC ENDDO
535	c
536	820 CONTINUE ! chercher le haut de la colonne
537	c
538	k2min = klev
539	DO i = 1, klon
540	aller(i) = .TRUE.
541	IF (todo(i)) k2min = MIN(k2min,k2(i))
542	ENDDO
543	IF (k2min.LT.klev) THEN
544	DO k = k2min, klev
545	DO i = 1, klon
546	IF (todo(i) .AND. k.GT.k2(i) .AND. aller(i)) THEN
547	zsat(i) = zsat(i) + zqmqsdp(i,k)
548	zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
549	IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) THEN
550	aller(i) = .FALSE.
551	ELSE
552	k2(i) = k
553	ENDIF
554	ENDIF
555	ENDDO
556	ENDDO
557	c error is = 0
558	c error DO i = 1, klon
559	c error IF(todo(i).AND.aller(i)) THEN
560	c error is = is + 1
561	c error todo(i) = .FALSE.
562	c error k2(i) = klev
563	c error ENDIF
564	c error ENDDO
565	c error IF (is.GT.0) THEN
566	c error PRINT*, "Bizard. je pourrais continuer mais j arrete"
567	c error CALL abort
568	c error ENDIF
569	ENDIF
570	c
571	CCC DO i = 1, klon
572	CCC IF (todo(i)) THEN
573	CCC 821 CONTINUE
574	CCC IF (k2(i) .EQ. klev) GOTO 822
575	CCC k = k2(i) + 1
576	CCC zsat(i) = zsat(i) + zqmqsdp(i,k)
577	CCC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
578	CCC IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 822
579	CCC k2(i) = k
580	CCC GOTO 821
581	CCC ENDIF
582	CCC 822 CONTINUE
583	CCC ENDDO
584	c
585	830 CONTINUE ! faire l'ajustement en sachant k1 et k2
586	c
587	is = 0
588	DO i = 1, klon
589	IF (todo(i)) THEN
590	IF (k2(i).LE.k1(i)) is = is + 1
591	ENDIF
592	ENDDO
593	IF (is.GT.0) THEN
594	PRINT*, "Impossible: k1 trop grand ou k2 trop petit"
595	PRINT*, "is=", is
596	CALL abort
597	ENDIF
598	c
599	k1min = klev
600	k1max = 1
601	k2max = 1
602	DO i = 1, klon
603	IF (todo(i)) THEN
604	k1min = MIN(k1min,k1(i))
605	k1max = MAX(k1max,k1(i))
606	k2max = MAX(k2max,k2(i))
607	ENDIF
608	ENDDO
609	c
610	DO i = 1, klon
611	IF (todo(i)) THEN
612	k = k1(i)
613	za(i) = 0.
614	zb(i) = ( RCPD(1.+zdqs(i,k))(zt(i,k)-za(i))
615	. -RLVTT(zqs(i,k)-zq(i,k)) )delp(i,k)
616	zc(i) = delp(i,k) * RCPD*(1.+zdqs(i,k))
617	ENDIF
618	ENDDO
619	c
620	DO k = k1min, k2max
621	DO i = 1, klon
622	IF (todo(i) .AND. k.GE.(k1(i)+1) .AND. k.LE.k2(i)) THEN
623	za(i) = za(i) + zgamdz(i,k-1)
624	zb(i) = zb(i)+(RCPD(1.+zdqs(i,k))(zt(i,k)-za(i))
625	. -RLVTT(zqs(i,k)-zq(i,k)) ) delp(i,k)
626	zc(i) = zc(i) + delp(i,k)RCPD(1.+zdqs(i,k))
627	ENDIF
628	ENDDO
629	ENDDO
630	c
631	DO i = 1, klon
632	IF (todo(i)) THEN
633	k = k1(i)
634	ztnew(i,k) = zb(i)/zc(i)
635	zqnew(i,k) = zqs(i,k) + (ztnew(i,k)-zt(i,k))
636	. RCPD/RLVTTzdqs(i,k)
637	ENDIF
638	ENDDO
639	c
640	DO k = k1min, k2max
641	DO i = 1, klon
642	IF (todo(i) .AND. k.GE.(k1(i)+1) .AND. k.LE.k2(i)) THEN
643	ztnew(i,k) = ztnew(i,k-1) + zgamdz(i,k-1)
644	zqnew(i,k) = zqs(i,k) + (ztnew(i,k)-zt(i,k))
645	. RCPD/RLVTTzdqs(i,k)
646	ENDIF
647	ENDDO
648	ENDDO
649	c
650	c Quantite de condensation produite pendant l'ajustement:
651	c
652	DO i = 1, klon
653	zcond(i) = 0.0
654	ENDDO
655	DO k = k1min, k2max
656	DO i = 1, klon
657	IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) THEN
658	rneb(i,k) = 1.0
659	zcond(i) = zcond(i) + (zq(i,k)-zqnew(i,k)) *delp(i,k)/RG
660	ENDIF
661	ENDDO
662	ENDDO
663	c
664	c Si condensation negative, effort completement perdu:
665	c
666	DO i = 1, klon
667	IF (todo(i).AND.zcond(i).LE.0.) todo(i) = .FALSE.
668	ENDDO
669	c
670	c L'ajustement a ete accompli, meme les calculs accessoires
671	c ne sont pas encore faits:
672	c
673	DO i = 1, klon
674	IF (todo(i)) accompli(i) = .TRUE.
675	ENDDO
676	c
677	c=====
678	c Une fois que la condensation a lieu, on doit construire un
679	c "modele nuageux" pour partager la condensation entre l'eau
680	c liquide nuageuse et la precipitation (leur rapport toliq
681	c est calcule selon l'epaisseur nuageuse). Je suppose que
682	c toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
683	c et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
684	c lineaire entre dpmin et dpmax).
685	c=====
686	DO i = 1, klon
687	IF (todo(i)) THEN
688	toliq(i) = tomax-((paprs(i,k1(i))-paprs(i,k2(i)+1))
689	. /paprs(i,1)-dpmin)
690	. *(tomax-tomin)/(dpmax-dpmin)
691	toliq(i) = MAX(tomin,MIN(tomax,toliq(i)))
692	IF (pplay(i,k2(i))/paprs(i,1) .LE. deep_sig) toliq(i) = deep_to
693	IF (old_tau) toliq(i) = 1.0
694	ENDIF
695	ENDDO
696	c=====
697	c On doit aussi determiner la distribution verticale de
698	c l'eau nuageuse. Plusieurs options sont proposees:
699	c
700	c (0) La condensation precipite integralement (toliq ne sera
701	c pas utilise).
702	c (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
703	c a la vapeur d'eau locale.
704	c (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
705	c (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
706	c est effectivement diminuee pendant le processus d'ajustement.
707	c (4) Elle est en fonction (lineaire ou exponentielle) de la
708	c distance (epaisseur en pression) avec le niveau k1 (la couche
709	c k1 n'aura donc pas d'eau liquide).
710	c=====
711	c
712	IF (opt_cld.EQ.0) THEN
713	c
714	DO i = 1, klon
715	IF (todo(i)) zrfl(i) = zcond(i) / dtime
716	ENDDO
717	c
718	ELSE IF (opt_cld.EQ.1) THEN
719	c
720	DO i = 1, klon
721	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
722	ENDDO
723	DO k = k1min, k2max
724	DO i = 1, klon
725	IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i))
726	. zvapo(i) = zvapo(i) + zqnew(i,k)*delp(i,k)/RG
727	ENDDO
728	ENDDO
729	DO i = 1, klon
730	IF (todo(i)) THEN
731	zrapp(i) = toliq(i) * zcond(i) / zvapo(i)
732	zrapp(i) = MAX(0.,MIN(1.,zrapp(i)))
733	zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime
734	ENDIF
735	ENDDO
736	DO k = k1min, k2max
737	DO i = 1, klon
738	IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) THEN
739	d_ql(i,k) = d_ql(i,k) + zrapp(i) * zqnew(i,k)
740	ENDIF
741	ENDDO
742	ENDDO
743	c
744	ELSE IF (opt_cld.EQ.2) THEN
745	c
746	DO i = 1, klon
747	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
748	ENDDO
749	DO k = k1min, k2max
750	DO i = 1, klon
751	IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i))
752	. zvapo(i) = zvapo(i) + delp(i,k)/RG
753	ENDDO
754	ENDDO
755	DO k = k1min, k2max
756	DO i = 1, klon
757	IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) THEN
758	d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i)
759	ENDIF
760	ENDDO
761	ENDDO
762	DO i = 1, klon
763	IF (todo(i)) zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime
764	ENDDO
765	c
766	ELSE IF (opt_cld.EQ.3) THEN
767	c
768	DO i = 1, klon
769	IF (todo(i)) zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
770	ENDDO
771	DO k = k1min, k2max
772	DO i = 1, klon
773	IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i))
774	. zvapo(i) = zvapo(i) + MAX(0.0,zq(i,k)-zqnew(i,k))
775	. * delp(i,k)/RG
776	ENDDO
777	ENDDO
778	DO k = k1min, k2max
779	DO i = 1, klon
780	IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i) .AND.
781	. zvapo(i).GT.0.0)
782	. d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i)
783	. * MAX(0.0,zq(i,k)-zqnew(i,k))
784	ENDDO
785	ENDDO
786	DO i = 1, klon
787	IF (todo(i)) zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime
788	ENDDO
789	c
790	ELSE IF (opt_cld.EQ.4) THEN
791	c
792	nexpo = 3
793	ccc nexpo = 1 ! distribution lineaire
794	c
795	DO i = 1, klon
796	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
797	ENDDO ! (avec ponderation)
798	DO k = k1min, k2max
799	DO i = 1, klon
800	IF (todo(i) .AND. k.GE.(k1(i)+1) .AND. k.LE.k2(i))
801	. zvapo(i) = zvapo(i) + delp(i,k) / RG
802	. * (pplay(i,k1(i))-pplay(i,k))**nexpo
803	ENDDO
804	ENDDO
805	DO k = k1min, k2max
806	DO i = 1, klon
807	IF (todo(i) .AND. k.GE.(k1(i)+1) .AND. k.LE.k2(i))
808	. d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i)
809	. * (pplay(i,k1(i))-pplay(i,k))**nexpo
810	ENDDO
811	ENDDO
812	DO i = 1, klon
813	IF (todo(i)) zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime
814	ENDDO
815	c
816	ELSE ! valeur non-prevue pour opt_cld
817	c
818	PRINT*, "opt_cld est faux:", opt_cld
819	CALL abort
820	c
821	ENDIF ! fin de opt_cld
822	c
823	c L'eau precipitante peut etre evaporee:
824	c
825	zalfa = 0.05
826	IF (evap_prec .AND. (k1max.GE.2)) THEN
827	DO k = k1max-1, 1, -1
828	DO i = 1, klon
829	IF (todo(i) .AND. k.LT.k1(i) .AND. zrfl(i).GT.0.0) THEN
830	zqev = MAX (0.0, (zqs(i,k)-zq(i,k))*zalfa )
831	zqevt = coef_eva * (1.0-zq(i,k)/zqs(i,k))*SQRT(zrfl(i))
832	. * delp(i,k)/pplay(i,k)zt(i,k)RD/RG
833	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * RG*dtime/delp(i,k)
834	zqev = MIN (zqev, zqevt)
835	zrfln = zrfl(i) - zqev*(delp(i,k))/RG/dtime
836	zq(i,k) = zq(i,k) - (zrfln-zrfl(i))
837	. * (RG/(delp(i,k)))*dtime
838	zt(i,k) = zt(i,k) + (zrfln-zrfl(i))
839	. * (RG/(delp(i,k)))*dtime
840	. * RLVTT/RCPD/(1.0+RVTMP2*zq(i,k))
841	zrfl(i) = zrfln
842	ENDIF
843	ENDDO
844	ENDDO
845	ENDIF
846	c
847	c La temperature de la premiere couche determine la pluie ou la neige:
848	c
849	DO i = 1, klon
850	IF (todo(i)) THEN
851	IF (zt(i,1) .GT. RTT) THEN
852	rain(i) = rain(i) + zrfl(i)
853	ELSE
854	snow(i) = snow(i) + zrfl(i)
855	ENDIF
856	ENDIF
857	ENDDO
858	c
859	c Mise a jour de la temperature et de l'humidite
860	c
861	DO k = k1min, k2max
862	DO i = 1, klon
863	IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) THEN
864	zt(i,k) = ztnew(i,k)
865	zq(i,k) = zqnew(i,k)
866	ENDIF
867	ENDDO
868	ENDDO
869	c
870	c Re-calculer certaines variables pour etendre et re-ajuster la colonne
871	c
872	IF (exigent) THEN
873	DO k = 1, klev
874	DO i = 1, klon
875	IF (todo(i)) THEN
876	IF (thermcep) THEN
877	zdelta=MAX(0.,SIGN(1.,RTT-zt(i,k)))
878	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
879	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i,k))
880	zqs(i,k)= R2ES*FOEEW(zt(i,k),zdelta)/pplay(i,k)
881	zqs(i,k)=MIN(0.5,zqs(i,k))
882	zcor=1./(1.-RETV*zqs(i,k))
883	zqs(i,k)=zqs(i,k)*zcor
884	zdqs(i,k) =FOEDE(zt(i,k),zdelta,zcvm5,zqs(i,k),zcor)
885	ELSE
886	IF (zt(i,k) .LT. t_coup) THEN
887	zqs(i,k)= qsats(zt(i,k)) / pplay(i,k)
888	zdqs(i,k)= dqsats(zt(i,k),zqs(i,k))
889	ELSE
890	zqs(i,k)= qsatl(zt(i,k)) / pplay(i,k)
891	zdqs(i,k)= dqsatl(zt(i,k),zqs(i,k))
892	ENDIF
893	ENDIF
894	ENDIF
895	ENDDO
896	ENDDO
897	ENDIF
898	c
899	IF (exigent) THEN
900	DO k = 1, klev-1
901	DO i = 1, klon
902	IF (todo(i)) THEN
903	zgamdz(i,k) = - (pplay(i,k)-pplay(i,k+1))/paprs(i,k+1)/RCPD
904	. ( RD(zt(i,k)delp(i,k)+zt(i,k+1)delp(i,k+1))
905	. /(delp(i,k)+delp(i,k+1))
906	. +RLVTT(zqs(i,k)delp(i,k)+zqs(i,k+1)*delp(i,k+1))
907	. /(delp(i,k)+delp(i,k+1))
908	. ) / (1.0+(zdqs(i,k)delp(i,k)+zdqs(i,k+1)delp(i,k+1))
909	. /(delp(i,k)+delp(i,k+1)) )
910	ENDIF
911	ENDDO
912	ENDDO
913	ENDIF
914	c
915	c Puisque l'humidite a ete modifiee, on re-fait (q-qs)*dp
916	c
917	DO k = 1, klev
918	DO i = 1, klon
919	IF (todo(i)) THEN
920	zqmqsdp(i,k) = (zq(i,k)-zqs(i,k))*delp(i,k)
921	ENDIF
922	ENDDO
923	ENDDO
924	c
925	c Verifier si l'on peut etendre le bas de la colonne
926	c
927	DO i = 1, klon
928	etendre(i) = .FALSE.
929	ENDDO
930	c
931	k1max = 1
932	DO i = 1, klon
933	IF (todo(i) .AND. k1(i).GT.(kbase+1)) THEN
934	k = k1(i)
935	zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
936	zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1)
937	csc voici l'ancienne ligne:
938	csc IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) THEN
939	csc sylvain: il faut RESPECTER les 2 criteres:
940	IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN
941	etendre(i) = .TRUE.
942	k1(i) = k1(i) - 1
943	k1max = MAX(k1max,k1(i))
944	aller(i) = .TRUE.
945	ENDIF
946	ENDIF
947	ENDDO
948	c
949	IF (k1max.GT.(kbase+1)) THEN
950	DO k = k1max, kbase+1, -1
951	DO i = 1, klon
952	IF (etendre(i) .AND. k.LT.k1(i) .AND. aller(i)) THEN
953	zsat(i) = zsat(i) + zqmqsdp(i,k)
954	zflo(i) = zt(i,k) + zgamdz(i,k) - zt(i,k+1)
955	IF (zsat(i).LE.0.0 .OR. zflo(i).LE.0.0) THEN
956	aller(i) = .FALSE.
957	ELSE
958	k1(i) = k
959	ENDIF
960	ENDIF
961	ENDDO
962	ENDDO
963	DO i = 1, klon
964	IF (etendre(i).AND.aller(i)) THEN
965	k1(i) = 1
966	ENDIF
967	ENDDO
968	ENDIF
969	c
970	CCC DO i = 1, klon
971	CCC IF (etendre(i)) THEN
972	CCC 840 k = k1(i)
973	CCC IF (k.GT.1) THEN
974	CCC zsat(i) = zsat(i) + zqmqsdp(i,k-1)
975	CCC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
976	CCC IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN
977	CCC k1(i) = k - 1
978	CCC GOTO 840
979	CCC ENDIF
980	CCC ENDIF
981	CCC ENDIF
982	CCC ENDDO
983	c
984	DO i = 1, klon
985	todobis(i) = todo(i)
986	todo(i) = .FALSE.
987	ENDDO
988	is = 0
989	DO i = 1, klon
990	IF (etendre(i)) THEN
991	todo(i) = .TRUE.
992	is = is + 1
993	ENDIF
994	ENDDO
995	IF (is.GT.0) THEN
996	IF (new_top) THEN
997	GOTO 820 ! chercher de nouveau le sommet k2
998	ELSE
999	GOTO 830 ! supposer que le sommet est celui deja trouve
1000	ENDIF
1001	ENDIF
1002	c
1003	DO i = 1, klon
1004	possible(i) = .FALSE.
1005	ENDDO
1006	is = 0
1007	DO i = 1, klon
1008	IF (todobis(i) .AND. k2(i).LT.klev) THEN
1009	is = is + 1
1010	possible(i) = .TRUE.
1011	ENDIF
1012	ENDDO
1013	IF (is.GT.0) GOTO 810 !on cherche en haut d'autres blocks
1014	c a ajuster a partir du sommet de la colonne precedente
1015	c
1016	860 CONTINUE ! Calculer les tendances et diagnostiques
1017	ccc print*, "Apres 860"
1018	c
1019	DO k = 1, klev
1020	DO i = 1, klon
1021	IF (accompli(i)) THEN
1022	d_t(i,k) = zt(i,k) - t(i,k)
1023	zq(i,k) = MAX(zq(i,k),seuil_vap)
1024	d_q(i,k) = zq(i,k) - q(i,k)
1025	ENDIF
1026	ENDDO
1027	ENDDO
1028	c
1029	DO 888 i = 1, klon
1030	IF (accompli(i)) THEN
1031	DO k = 1, klev
1032	IF (rneb(i,k).GT.0.0) THEN
1033	ibas(i) = k
1034	GOTO 807
1035	ENDIF
1036	ENDDO
1037	807 CONTINUE
1038	DO k = klev, 1, -1
1039	IF (rneb(i,k).GT.0.0) THEN
1040	itop(i) = k
1041	GOTO 808
1042	ENDIF
1043	ENDDO
1044	808 CONTINUE
1045	ENDIF
1046	888 CONTINUE
1047	c
1048	IF (imprim) THEN
1049	nbtodo = 0
1050	nbdone = 0
1051	DO i = 1, klon
1052	IF (afaire(i)) nbtodo = nbtodo + 1
1053	IF (accompli(i)) nbdone = nbdone + 1
1054	ENDDO
1055	PRINT*, "nbTodo, nbDone=", nbtodo, nbdone
1056	ENDIF
1057	c
1058	RETURN
1059	END
1060	SUBROUTINE conkuo(dtime, paprs, pplay, t, q, conv_q,
1061	s d_t, d_q, d_ql, rneb,
1062	s rain, snow, ibas, itop)
1063	IMPLICIT none
1064	c======================================================================
1065	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1066	c Objet: Schema de convection de type Kuo (1965).
1067	c Cette version du code peut calculer le niveau de depart
1068	c N.B. version vectorielle (le 6 oct. 1997)
1069	c======================================================================
1070	#include "dimensions.h"
1071	#include "dimphy.h"
1072	#include "YOMCST.h"
1073	c
1074	c Arguments:
1075	c
1076	REAL dtime ! intervalle du temps (s)
1077	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
1078	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
1079	REAL t(klon,klev) ! temperature (K)
1080	REAL q(klon,klev) ! humidite specifique
1081	REAL conv_q(klon,klev) ! taux de convergence humidite (g/g/s)
1082	c
1083	REAL d_t(klon,klev) ! incrementation temperature
1084	REAL d_q(klon,klev) ! incrementation humidite
1085	REAL d_ql(klon,klev) ! incrementation eau liquide
1086	REAL rneb(klon,klev) ! nebulosite
1087	REAL rain(klon) ! pluies (mm/s)
1088	REAL snow(klon) ! neige (mm/s)
1089	INTEGER itop(klon) ! niveau du sommet
1090	INTEGER ibas(klon) ! niveau du bas
1091	c
1092	LOGICAL ldcum(klon) ! convection existe
1093	LOGICAL todo(klon)
1094	c
1095	c Quelsques options:
1096	c
1097	LOGICAL calcfcl ! calculer le niveau de convection libre
1098	PARAMETER (calcfcl=.TRUE.)
1099	INTEGER ldepar ! niveau fixe de convection libre
1100	PARAMETER (ldepar=4)
1101	INTEGER opt_cld ! comment traiter l'eau liquide
1102	PARAMETER (opt_cld=4) ! valeur possible: 0, 1, 2, 3 ou 4
1103	LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
1104	PARAMETER (evap_prec=.TRUE.)
1105	REAL coef_eva
1106	PARAMETER (coef_eva=1.0E-05)
1107	LOGICAL new_deh ! nouvelle facon de calculer dH
1108	PARAMETER (new_deh=.FALSE.)
1109	REAL t_coup
1110	PARAMETER (t_coup=234.0)
1111	LOGICAL old_tau ! implique precipitation nulle
1112	PARAMETER (old_tau=.FALSE.)
1113	REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
1114	REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
1115	PARAMETER (dpmin=0.15, tomax=0.97)
1116	REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
1117	PARAMETER (dpmax=0.30, tomin=0.05)
1118	REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
1119	PARAMETER (deep_sig=0.50, deep_to=0.05)
1120	c
1121	c Variables locales:
1122	c
1123	INTEGER nexpo
1124	LOGICAL nuage(klon)
1125	INTEGER i, k, kbmin, kbmax, khmax
1126	REAL ztotal(klon,klev), zdeh(klon,klev)
1127	REAL zgz(klon,klev)
1128	REAL zqs(klon,klev)
1129	REAL zdqs(klon,klev)
1130	REAL ztemp(klon,klev)
1131	REAL zpres(klon,klev)
1132	REAL zconv(klon) ! convergence d'humidite
1133	REAL zvirt(klon) ! convergence virtuelle d'humidite
1134	REAL zfrac(klon) ! fraction convective
1135	INTEGER kb(klon), kh(klon)
1136	c
1137	REAL zcond(klon), zvapo(klon), zrapp(klon)
1138	REAL zrfl(klon), zrfln, zqev, zqevt
1139	REAL zdelta, zcvm5, zcor
1140	REAL zvar
1141	c
1142	LOGICAL appel1er
1143	SAVE appel1er
1144	c
1145	c Fonctions thermodynamiques
1146	c
1147	#include "YOETHF.h"
1148	#include "FCTTRE.h"
1149	c
1150	DATA appel1er /.TRUE./
1151	c
1152	IF (appel1er) THEN
1153	PRINT*, 'conkuo, calcfcl:', calcfcl
1154	IF (.NOT.calcfcl) PRINT*, 'conkuo, ldepar:', ldepar
1155	PRINT*, 'conkuo, opt_cld:', opt_cld
1156	PRINT*, 'conkuo, evap_prec:', evap_prec
1157	PRINT*, 'conkuo, new_deh:', new_deh
1158	appel1er = .FALSE.
1159	ENDIF
1160	c
1161	c Initialiser les sorties a zero
1162	c
1163	DO k = 1, klev
1164	DO i = 1, klon
1165	d_q(i,k) = 0.0
1166	d_t(i,k) = 0.0
1167	d_ql(i,k) = 0.0
1168	rneb(i,k) = 0.0
1169	ENDDO
1170	ENDDO
1171	DO i = 1, klon
1172	rain(i) = 0.0
1173	snow(i) = 0.0
1174	ibas(i) = 0
1175	itop(i) = 0
1176	ENDDO
1177	c
1178	c Calculer la vapeur d'eau saturante Qs et sa derive L/Cp * dQs/dT
1179	c
1180	DO k = 1, klev
1181	DO i = 1, klon
1182	IF (thermcep) THEN
1183	zdelta=MAX(0.,SIGN(1.,RTT-t(i,k)))
1184	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
1185	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*q(i,k))
1186	zqs(i,k)=R2ES*FOEEW(t(i,k),zdelta)/pplay(i,k)
1187	zqs(i,k)=MIN(0.5,zqs(i,k))
1188	zcor=1./(1.-RETV*zqs(i,k))
1189	zqs(i,k)=zqs(i,k)*zcor
1190	zdqs(i,k) =FOEDE(t(i,k),zdelta,zcvm5,zqs(i,k),zcor)
1191	ELSE
1192	IF (t(i,k).LT.t_coup) THEN
1193	zqs(i,k) = qsats(t(i,k))/pplay(i,k)
1194	zdqs(i,k) = dqsats(t(i,k),zqs(i,k))
1195	ELSE
1196	zqs(i,k) = qsatl(t(i,k))/pplay(i,k)
1197	zdqs(i,k) = dqsatl(t(i,k),zqs(i,k))
1198	ENDIF
1199	ENDIF
1200	ENDDO
1201	ENDDO
1202	c
1203	c Calculer gz (energie potentielle)
1204	c
1205	DO i = 1, klon
1206	zgz(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1207	. * (paprs(i,1)-pplay(i,1))
1208	ENDDO
1209	DO k = 2, klev
1210	DO i = 1, klon
1211	zgz(i,k) = zgz(i,k-1)
1212	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1213	. * (pplay(i,k-1)-pplay(i,k))
1214	ENDDO
1215	ENDDO
1216	c
1217	c Calculer l'energie statique humide saturee (CpT + gz + LQs)
1218	c
1219	DO k = 1, klev
1220	DO i = 1, klon
1221	ztotal(i,k) = RCPDt(i,k) + RLVTTzqs(i,k) + zgz(i,k)
1222	ENDDO
1223	ENDDO
1224	c
1225	c Determiner le niveau de depart et calculer la difference de
1226	c l'energie statique humide saturee (ztotal) entre la couche
1227	c de depart et chaque couche au-dessus.
1228	c
1229	IF (calcfcl) THEN
1230	DO k = 1, klev
1231	DO i = 1, klon
1232	zpres(i,k) = pplay(i,k)
1233	ztemp(i,k) = t(i,k)
1234	ENDDO
1235	ENDDO
1236	CALL kuofcl(ztemp, q, zgz, zpres, ldcum, kb)
1237	DO i = 1, klon
1238	IF (ldcum(i)) THEN
1239	k = kb(i)
1240	IF (new_deh) THEN
1241	zdeh(i,k) = ztotal(i,k-1) - ztotal(i,k)
1242	ELSE
1243	zdeh(i,k) = RCPD * (t(i,k-1)-t(i,k))
1244	. - RD 0.5(t(i,k-1)+t(i,k))/paprs(i,k)
1245	. *(pplay(i,k-1)-pplay(i,k))
1246	. + RLVTT*(zqs(i,k-1)-zqs(i,k))
1247	ENDIF
1248	zdeh(i,k) = zdeh(i,k) * 0.5
1249	ENDIF
1250	ENDDO
1251	DO k = 1, klev
1252	DO i = 1, klon
1253	IF (ldcum(i) .AND. k.GE.(kb(i)+1)) THEN
1254	IF (new_deh) THEN
1255	zdeh(i,k) = zdeh(i,k-1) + (ztotal(i,k-1)-ztotal(i,k))
1256	ELSE
1257	zdeh(i,k) = zdeh(i,k-1)
1258	. + RCPD * (t(i,k-1)-t(i,k))
1259	. - RD 0.5(t(i,k-1)+t(i,k))/paprs(i,k)
1260	. *(pplay(i,k-1)-pplay(i,k))
1261	. + RLVTT*(zqs(i,k-1)-zqs(i,k))
1262	ENDIF
1263	ENDIF
1264	ENDDO
1265	ENDDO
1266	ELSE
1267	DO i = 1, klon
1268	k = ldepar
1269	kb(i) = ldepar
1270	ldcum(i) = .TRUE.
1271	IF (new_deh) THEN
1272	zdeh(i,k) = ztotal(i,k-1) - ztotal(i,k)
1273	ELSE
1274	zdeh(i,k) = RCPD * (t(i,k-1)-t(i,k))
1275	. - RD 0.5(t(i,k-1)+t(i,k))/paprs(i,k)
1276	. *(pplay(i,k-1)-pplay(i,k))
1277	. + RLVTT*(zqs(i,k-1)-zqs(i,k))
1278	ENDIF
1279	zdeh(i,k) = zdeh(i,k) * 0.5
1280	ENDDO
1281	DO k = ldepar+1, klev
1282	DO i = 1, klon
1283	IF (new_deh) THEN
1284	zdeh(i,k) = zdeh(i,k-1) + (ztotal(i,k-1)-ztotal(i,k))
1285	ELSE
1286	zdeh(i,k) = zdeh(i,k-1)
1287	. + RCPD * (t(i,k-1)-t(i,k))
1288	. - RD 0.5(t(i,k-1)+t(i,k))/paprs(i,k)
1289	. *(pplay(i,k-1)-pplay(i,k))
1290	. + RLVTT*(zqs(i,k-1)-zqs(i,k))
1291	ENDIF
1292	ENDDO
1293	ENDDO
1294	ENDIF
1295	c
1296	c-----Chercher le sommet du nuage
1297	c-----Calculer la convergence de l'humidite (en kg/m**2 a un facteur
1298	c-----psolpa/RG pres) du bas jusqu'au sommet du nuage.
1299	c-----Calculer la convergence virtuelle pour que toute la maille
1300	c-----deviennt nuageuse (du bas jusqu'au sommet du nuage)
1301	c
1302	DO i = 1, klon
1303	nuage(i) = .TRUE.
1304	zconv(i) = 0.0
1305	zvirt(i) = 0.0
1306	kh(i) = -999
1307	ENDDO
1308	DO k = 1, klev
1309	DO i = 1, klon
1310	IF (k.GE.kb(i) .AND. ldcum(i)) THEN
1311	nuage(i) = nuage(i) .AND. zdeh(i,k).GT.0.0
1312	IF (nuage(i)) THEN
1313	kh(i) = k
1314	zconv(i)=zconv(i)+conv_q(i,k)*dtime
1315	. *(paprs(i,k)-paprs(i,k+1))
1316	zvirt(i)=zvirt(i)+(zdeh(i,k)/RLVTT+zqs(i,k)-q(i,k))
1317	. *(paprs(i,k)-paprs(i,k+1))
1318	ENDIF
1319	ENDIF
1320	ENDDO
1321	ENDDO
1322	c
1323	DO i = 1, klon
1324	todo(i) = ldcum(i) .AND. kh(i).GT.kb(i) .AND. zconv(i).GT.0.0
1325	ENDDO
1326	c
1327	kbmin = klev
1328	kbmax = 0
1329	khmax = 0
1330	DO i = 1, klon
1331	IF (todo(i)) THEN
1332	kbmin = MIN(kbmin,kb(i))
1333	kbmax = MAX(kbmax,kb(i))
1334	khmax = MAX(khmax,kh(i))
1335	ENDIF
1336	ENDDO
1337	c
1338	c-----Calculer la surface couverte par le nuage
1339	c
1340	DO i = 1, klon
1341	IF (todo(i)) THEN
1342	zfrac(i) = MAX(0.0,MIN(zconv(i)/zvirt(i), 1.0))
1343	ENDIF
1344	ENDDO
1345	c
1346	c-----Calculs essentiels:
1347	c
1348	DO i = 1, klon
1349	IF (todo(i)) THEN
1350	zcond(i) = 0.0
1351	ENDIF
1352	ENDDO
1353	DO k = kbmin, khmax
1354	DO i = 1, klon
1355	IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN
1356	zvar = zdeh(i,k)/(1.+zdqs(i,k))
1357	d_t(i,k) = zvar * zfrac(i) / RCPD
1358	d_q(i,k) = (zvarzdqs(i,k)/RLVTT+zqs(i,k)-q(i,k))zfrac(i)
1359	. - conv_q(i,k)*dtime
1360	zcond(i) = zcond(i) - d_q(i,k) *(paprs(i,k)-paprs(i,k+1))/RG
1361	rneb(i,k) = zfrac(i)
1362	ENDIF
1363	ENDDO
1364	ENDDO
1365	c
1366	DO i = 1, klon
1367	IF (todo(i) .AND. zcond(i).LT.0.0) THEN
1368	PRINT*, 'WARNING: cond. negative (Kuo) ',
1369	. i,kb(i),kh(i), zcond(i)
1370	zcond(i) = 0.0
1371	DO k = kb(i), kh(i)
1372	d_t(i,k) = 0.0
1373	d_q(i,k) = 0.0
1374	ENDDO
1375	todo(i) = .FALSE. ! effort totalement perdu
1376	ENDIF
1377	ENDDO
1378	c
1379	c=====
1380	c Une fois que la condensation a lieu, on doit construire un
1381	c "modele nuageux" pour partager la condensation entre l'eau
1382	c liquide nuageuse et la precipitation (leur rapport toliq
1383	c est calcule selon l'epaisseur nuageuse). Je suppose que
1384	c toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
1385	c et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
1386	c lineaire entre dpmin et dpmax).
1387	c=====
1388	DO i = 1, klon
1389	IF (todo(i)) THEN
1390	toliq(i) = tomax-((paprs(i,kb(i))-paprs(i,kh(i)+1))
1391	. /paprs(i,1)-dpmin)
1392	. *(tomax-tomin)/(dpmax-dpmin)
1393	toliq(i) = MAX(tomin,MIN(tomax,toliq(i)))
1394	IF (pplay(i,kh(i))/paprs(i,1) .LE. deep_sig) toliq(i) = deep_to
1395	IF (old_tau) toliq(i) = 1.0
1396	ENDIF
1397	ENDDO
1398	c=====
1399	c On doit aussi determiner la distribution verticale de
1400	c l'eau nuageuse. Plusieurs options sont proposees:
1401	c
1402	c (0) La condensation precipite integralement (toliq ne sera
1403	c pas utilise).
1404	c (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
1405	c a la vapeur d'eau locale.
1406	c (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
1407	c (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
1408	c est effectivement diminuee pendant le processus d'ajustement.
1409	c (4) Elle est en fonction (lineaire ou exponentielle) de la
1410	c distance (epaisseur en pression) avec le niveau k1 (la couche
1411	c k1 n'aura donc pas d'eau liquide).
1412	c=====
1413	c
1414	IF (opt_cld.EQ.0) THEN
1415	c
1416	DO i = 1, klon
1417	IF (todo(i)) zrfl(i) = zcond(i) / dtime
1418	ENDDO
1419	c
1420	ELSE IF (opt_cld.EQ.1) THEN
1421	c
1422	DO i = 1, klon
1423	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
1424	ENDDO
1425	DO k = kbmin, khmax
1426	DO i = 1, klon
1427	IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN
1428	zvapo(i) = zvapo(i) + (q(i,k)+d_q(i,k))
1429	. *(paprs(i,k)-paprs(i,k+1))/RG
1430	ENDIF
1431	ENDDO
1432	ENDDO
1433	DO i = 1, klon
1434	IF (todo(i)) THEN
1435	zrapp(i) = toliq(i) * zcond(i) / zvapo(i)
1436	zrapp(i) = MAX(0.,MIN(1.,zrapp(i)))
1437	ENDIF
1438	ENDDO
1439	DO k = kbmin, khmax
1440	DO i = 1, klon
1441	IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN
1442	d_ql(i,k) = zrapp(i) * (q(i,k)+d_q(i,k))
1443	ENDIF
1444	ENDDO
1445	ENDDO
1446	DO i = 1, klon
1447	IF (todo(i)) THEN
1448	zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime
1449	ENDIF
1450	ENDDO
1451	c
1452	ELSE IF (opt_cld.EQ.2) THEN
1453	c
1454	DO i = 1, klon
1455	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
1456	ENDDO
1457	DO k = kbmin, khmax
1458	DO i = 1, klon
1459	IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN
1460	zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/RG
1461	ENDIF
1462	ENDDO
1463	ENDDO
1464	DO k = kbmin, khmax
1465	DO i = 1, klon
1466	IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN
1467	d_ql(i,k) = toliq(i) * zcond(i) / zvapo(i)
1468	ENDIF
1469	ENDDO
1470	ENDDO
1471	DO i = 1, klon
1472	IF (todo(i)) THEN
1473	zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime
1474	ENDIF
1475	ENDDO
1476	c
1477	ELSE IF (opt_cld.EQ.3) THEN
1478	c
1479	DO i = 1, klon
1480	IF (todo(i)) THEN
1481	zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
1482	ENDIF
1483	ENDDO
1484	DO k = kbmin, khmax
1485	DO i = 1, klon
1486	IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN
1487	zvapo(i) = zvapo(i) + MAX(0.0,-d_q(i,k))
1488	. * (paprs(i,k)-paprs(i,k+1))/RG
1489	ENDIF
1490	ENDDO
1491	ENDDO
1492	DO k = kbmin, khmax
1493	DO i = 1, klon
1494	IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i) .AND.
1495	. zvapo(i).GT.0.0) THEN
1496	d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i)
1497	. * MAX(0.0,-d_q(i,k))
1498	ENDIF
1499	ENDDO
1500	ENDDO
1501	DO i = 1, klon
1502	IF (todo(i)) THEN
1503	zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime
1504	ENDIF
1505	ENDDO
1506	c
1507	ELSE IF (opt_cld.EQ.4) THEN
1508	c
1509	nexpo = 3
1510	ccc nexpo = 1 ! distribution lineaire
1511	c
1512	DO i = 1, klon
1513	IF (todo(i)) THEN
1514	zvapo(i) = 0.0 ! quantite integrale de masse (avec ponderation)
1515	ENDIF
1516	ENDDO
1517	DO k = kbmin, khmax
1518	DO i = 1, klon
1519	IF (todo(i) .AND. k.GE.(kb(i)+1) .AND. k.LE.kh(i)) THEN
1520	zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1)) / RG
1521	. * (pplay(i,kb(i))-pplay(i,k))**nexpo
1522	ENDIF
1523	ENDDO
1524	ENDDO
1525	DO k = kbmin, khmax
1526	DO i = 1, klon
1527	IF (todo(i) .AND. k.GE.(kb(i)+1) .AND. k.LE.kh(i)) THEN
1528	d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i)
1529	. * (pplay(i,kb(i))-pplay(i,k))**nexpo
1530	ENDIF
1531	ENDDO
1532	ENDDO
1533	DO i = 1, klon
1534	IF (todo(i)) THEN
1535	zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime
1536	ENDIF
1537	ENDDO
1538	c
1539	ELSE ! valeur non-prevue pour opt_cld
1540	c
1541	PRINT*, "opt_cld est faux:", opt_cld
1542	CALL abort
1543	c
1544	ENDIF ! fin de opt_cld
1545	c
1546	c L'eau precipitante peut etre re-evaporee:
1547	c
1548	IF (evap_prec .AND. kbmax.GE.2) THEN
1549	DO k = kbmax, 1, -1
1550	DO i = 1, klon
1551	IF (todo(i) .AND. k.LE.(kb(i)-1) .AND. zrfl(i).GT.0.0) THEN
1552	zqev = MAX (0.0, (zqs(i,k)-q(i,k))*zfrac(i) )
1553	zqevt = coef_eva * (1.0-q(i,k)/zqs(i,k))*SQRT(zrfl(i))
1554	. * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)t(i,k)RD/RG
1555	zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
1556	. * RG*dtime/(paprs(i,k)-paprs(i,k+1))
1557	zqev = MIN (zqev, zqevt)
1558	zrfln = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
1559	. /RG/dtime
1560	d_q(i,k) = - (zrfln-zrfl(i))
1561	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
1562	d_t(i,k) = (zrfln-zrfl(i))
1563	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
1564	. * RLVTT/RCPD
1565	zrfl(i) = zrfln
1566	ENDIF
1567	ENDDO
1568	ENDDO
1569	ENDIF
1570	c
1571	c La temperature de la premiere couche determine la pluie ou la neige:
1572	c
1573	DO i = 1, klon
1574	IF (todo(i)) THEN
1575	IF (t(i,1) .GT. RTT) THEN
1576	rain(i) = rain(i) + zrfl(i)
1577	ELSE
1578	snow(i) = snow(i) + zrfl(i)
1579	ENDIF
1580	ENDIF
1581	ENDDO
1582	c
1583	RETURN
1584	END
1585	SUBROUTINE kuofcl(pt, pq, pg, pp, LDCUM, kcbot)
1586	IMPLICIT none
1587	c======================================================================
1588	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
1589	c adaptation du code de Tiedtke du ECMWF
1590	c Objet: calculer le niveau de convection libre
1591	c (FCL: Free Convection Level)
1592	c======================================================================
1593	c Arguments:
1594	c pt---input-R- temperature (K)
1595	c pq---input-R- vapeur d'eau (kg/kg)
1596	c pg---input-R- geopotentiel (g*z ou z est en metre)
1597	c pp---input-R- pression (Pa)
1598	c
1599	c LDCUM---output-L- Y-t-il la convection
1600	c kcbot---output-I- Niveau du bas de la convection
1601	c======================================================================
1602	#include "dimensions.h"
1603	#include "dimphy.h"
1604	#include "YOMCST.h"
1605	#include "YOETHF.h"
1606	C
1607	REAL pt(klon,klev), pq(klon,klev), pg(klon,klev), pp(klon,klev)
1608	INTEGER kcbot(klon)
1609	LOGICAL LDCUM(klon)
1610	C
1611	REAL ztu(klon,klev), zqu(klon,klev), zlu(klon,klev)
1612	REAL zqold(klon), zbuo
1613	INTEGER is, i, k
1614	c
1615	c klab=1: on est sous le nuage convectif
1616	c klab=2: le bas du nuage convectif
1617	c klab=0: autres couches
1618	INTEGER klab(klon,klev)
1619	c
1620	c quand lflag=.true., on est sous le nuage, il faut donc appliquer
1621	c le processus d'elevation.
1622	LOGICAL lflag(klon)
1623	C
1624	DO k = 1, klev
1625	DO i = 1, klon
1626	ztu(i,k) = pt(i,k)
1627	zqu(i,k) = pq(i,k)
1628	zlu(i,k) = 0.0
1629	klab(i,k) = 0
1630	ENDDO
1631	ENDDO
1632	C----------------------------------------------------------------------
1633	DO i = 1, klon
1634	klab(i,1)=1
1635	kcbot(i)=2
1636	LDCUM(i)=.FALSE.
1637	ENDDO
1638	C
1639	DO 290 k = 2, klev-1
1640	c
1641	is=0
1642	DO i = 1, klon
1643	if (klab(i,k-1).EQ.1) is = is + 1
1644	lflag(i) = .FALSE.
1645	if (klab(i,k-1).EQ.1) lflag(i) = .TRUE.
1646	ENDDO
1647	IF (is.EQ.0) GOTO 290
1648	c
1649	c on eleve le parcel d'air selon l'adiabatique sec
1650	c
1651	DO i = 1, klon
1652	IF (lflag(i)) THEN
1653	zqu(i,k) = zqu(i,k-1)
1654	ztu(i,k) = ztu(i,k-1) + (pg(i,k-1)-pg(i,k))/RCPD
1655	zbuo = ztu(i,k)(1.+RETVzqu(i,k))-
1656	. pt(i,k)(1.+RETVpq(i,k))+0.5
1657	IF (zbuo.GT.0.) klab(i,k)=1
1658	zqold(i) = zqu(i,k)
1659	ENDIF
1660	ENDDO
1661	c
1662	c on calcule la condensation eventuelle
1663	c
1664	CALL adjtq(pp(1,k), ztu(1,k), zqu(1,k), lflag, 1)
1665	c
1666	c s'il y a la condensation et la "buoyancy" force est positive
1667	c c'est bien le bas de la tour de convection
1668	c
1669	DO i=1, klon
1670	IF(lflag(i).AND.zqu(i,k).NE.zqold(i)) THEN
1671	klab(i,k) = 2
1672	zlu(i,k) = zlu(i,k)+zqold(i)-zqu(i,k)
1673	zbuo = ztu(i,k)(1.+RETVzqu(i,k))-
1674	. pt(i,k)(1.+RETVpq(i,k))+0.5
1675	IF (zbuo.GT.0.) THEN
1676	kcbot(i) = k
1677	LDCUM(i) = .TRUE.
1678	ENDIF
1679	ENDIF
1680	ENDDO
1681	C
1682	290 CONTINUE
1683	C
1684	RETURN
1685	END
1686	SUBROUTINE adjtq(pp, pt, pq, LDFLAG, KCALL)
1687	IMPLICIT none
1688	c======================================================================
1689	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
1690	c adaptation du code de Tiedtke du ECMWF
1691	c Objet: ajustement entre T et Q
1692	c======================================================================
1693	c Arguments:
1694	c pp---input-R- pression (Pa)
1695	c pt---input/output-R- temperature (K)
1696	c pq---input/output-R- vapeur d'eau (kg/kg)
1697	c======================================================================
1698	C TO PRODUCE T,Q AND L VALUES FOR CLOUD ASCENT
1699	C
1700	C NOTE: INPUT PARAMETER KCALL DEFINES CALCULATION AS
1701	C KCALL=0 ENV. T AND QS INCUINI
1702	C KCALL=1 CONDENSATION IN UPDRAFTS (E.G. CUBASE, CUASC)
1703	C KCALL=2 EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
1704	C
1705	#include "dimensions.h"
1706	#include "dimphy.h"
1707	#include "YOMCST.h"
1708	C
1709	REAL pt(klon), pq(klon), pp(klon)
1710	LOGICAL ldflag(klon)
1711	INTEGER KCALL
1712	c
1713	REAL t_coup
1714	PARAMETER (t_coup=234.0)
1715	c
1716	REAL zcond(klon), zcond1
1717	REAL zdelta, zcvm5, zldcp, zqsat, zcor, zdqsat
1718	INTEGER is, i
1719	#include "YOETHF.h"
1720	#include "FCTTRE.h"
1721	c
1722	DO i = 1, klon
1723	zcond(i) = 0.0
1724	ENDDO
1725	C
1726	DO 210 i=1, klon
1727	IF (LDFLAG(i)) THEN
1728	zdelta=MAX(0.,SIGN(1.,RTT-pt(i)))
1729	zldcp = RLVTT(1.-zdelta) + zdeltaRLSTT
1730	zldcp = zldcp / RCPD/(1.0+RVTMP2*pq(i))
1731	IF (thermcep) THEN
1732	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
1733	zcvm5 = zcvm5 / RCPD/(1.0+RVTMP2*pq(i))
1734	zqsat=R2ES*FOEEW (pt(i), zdelta) / pp(i)
1735	zqsat=MIN(0.5,zqsat)
1736	zcor=1./(1.-RETV *zqsat)
1737	zqsat=zqsat*zcor
1738	zdqsat = FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor)
1739	ELSE
1740	IF (pt(i).LT.t_coup) THEN
1741	zqsat = qsats(pt(i))/pp(i)
1742	zdqsat = dqsats(pt(i),zqsat)
1743	ELSE
1744	zqsat = qsatl(pt(i))/pp(i)
1745	zdqsat = dqsatl(pt(i),zqsat)
1746	ENDIF
1747	ENDIF
1748	zcond(i)=(pq(i)-zqsat) / (1. + zdqsat)
1749	IF(KCALL.EQ.1) zcond(i)=MAX(zcond(i),0.)
1750	IF(KCALL.EQ.2) zcond(i)=MIN(zcond(i),0.)
1751	pt(i)=pt(i)+zldcp*zcond(i)
1752	pq(i)=pq(i)-zcond(i)
1753	ENDIF
1754	210 CONTINUE
1755	C
1756	is = 0
1757	DO i=1, klon
1758	if (zcond(i).NE.0.) is = is + 1
1759	ENDDO
1760	IF(is.EQ.0) GOTO 230
1761	C
1762	DO 220 i = 1, klon
1763	IF(LDFLAG(i).AND.zcond(i).NE.0.) THEN
1764	zdelta=MAX(0.,SIGN(1.,RTT-pt(i)))
1765	zldcp = RLVTT(1.-zdelta) + zdeltaRLSTT
1766	zldcp = zldcp / RCPD/(1.0+RVTMP2*pq(i))
1767	IF (thermcep) THEN
1768	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
1769	zcvm5 = zcvm5 / RCPD/(1.0+RVTMP2*pq(i))
1770	zqsat=R2ES*FOEEW (pt(i), zdelta) / pp(i)
1771	zqsat=MIN(0.5,zqsat)
1772	zcor=1./(1.-RETV *zqsat)
1773	zqsat=zqsat*zcor
1774	zdqsat = FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor)
1775	ELSE
1776	IF (pt(i).LT.t_coup) THEN
1777	zqsat = qsats(pt(i))/pp(i)
1778	zdqsat = dqsats(pt(i),zqsat)
1779	ELSE
1780	zqsat = qsatl(pt(i))/pp(i)
1781	zdqsat = dqsatl(pt(i),zqsat)
1782	ENDIF
1783	ENDIF
1784	zcond1=(pq(i)-zqsat) / (1.+zdqsat)
1785	pt(i)=pt(i)+zldcp*zcond1
1786	pq(i)=pq(i)-zcond1
1787	END IF
1788	220 CONTINUE
1789	C
1790	230 CONTINUE
1791	RETURN
1792	END
1793	SUBROUTINE fiajh(dtime, paprs, pplay, t, q,
1794	. d_t, d_q, d_ql, rneb,
1795	. rain, snow, ibas, itop)
1796	IMPLICIT NONE
1797	c
1798	c Ajustement humide (Schema de convection de Manabe)
1799	C.
1800	#include "dimensions.h"
1801	#include "dimphy.h"
1802	#include "YOMCST.h"
1803	c
1804	c Arguments:
1805	c
1806	REAL dtime ! intervalle du temps (s)
1807	REAL t(klon,klev) ! temperature (K)
1808	REAL q(klon,klev) ! humidite specifique (kg/kg)
1809	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
1810	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
1811	c
1812	REAL d_t(klon,klev) ! incrementation pour la temperature
1813	REAL d_q(klon,klev) ! incrementation pour vapeur d'eau
1814	REAL d_ql(klon,klev) ! incrementation pour l'eau liquide
1815	REAL rneb(klon,klev) ! fraction nuageuse
1816	c
1817	REAL rain(klon) ! variable non utilisee
1818	REAL snow(klon) ! variable non utilisee
1819	INTEGER ibas(klon) ! variable non utilisee
1820	INTEGER itop(klon) ! variable non utilisee
1821
1822	REAL t_coup
1823	PARAMETER (t_coup=234.0)
1824	REAL seuil_vap
1825	PARAMETER (seuil_vap=1.0E-10)
1826	c
1827	c Variables locales:
1828	c
1829	INTEGER i, k
1830	INTEGER k1, k1p, k2, k2p
1831	LOGICAL itest(klon)
1832	REAL delta_q(klon, klev)
1833	REAL cp_new_t(klev)
1834	REAL cp_delta_t(klev)
1835	REAL new_qb(klev)
1836	REAL v_cptj(klev), v_cptjk1, v_ssig
1837	REAL v_cptt(klon,klev), v_p, v_t
1838	REAL v_qs(klon,klev), v_qsd(klon,klev)
1839	REAL zq1(klon), zq2(klon)
1840	REAL gamcpdz(klon,2:klev)
1841	REAL zdp, zdpm
1842	c
1843	REAL zsat ! sur-saturation
1844	REAL zflo ! flotabilite
1845	c
1846	REAL local_q(klon,klev),local_t(klon,klev)
1847	c
1848	REAL zdelta, zcor, zcvm5
1849	C
1850	#include "YOETHF.h"
1851	#include "FCTTRE.h"
1852	C
1853	DO k = 1, klev
1854	DO i = 1, klon
1855	local_q(i,k) = q(i,k)
1856	local_t(i,k) = t(i,k)
1857	rneb(i,k) = 0.0
1858	d_ql(i,k) = 0.0
1859	d_t(i,k) = 0.0
1860	d_q(i,k) = 0.0
1861	ENDDO
1862	ENDDO
1863	DO i = 1, klon
1864	rain(i) = 0.0
1865	snow(i) = 0.0
1866	ibas(i) = 0
1867	itop(i) = 0
1868	ENDDO
1869	c
1870	c Calculer v_qs et v_qsd:
1871	c
1872	DO k = 1, klev
1873	DO i = 1, klon
1874	v_cptt(i,k) = RCPD * local_t(i,k)
1875	v_t = local_t(i,k)
1876	v_p = pplay(i,k)
1877	c
1878	IF (thermcep) THEN
1879	zdelta=MAX(0.,SIGN(1.,RTT-v_t))
1880	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
1881	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*local_q(i,k))
1882	v_qs(i,k)= R2ES * FOEEW(v_t,zdelta)/v_p
1883	v_qs(i,k)=MIN(0.5,v_qs(i,k))
1884	zcor=1./(1.-RETV*v_qs(i,k))
1885	v_qs(i,k)=v_qs(i,k)*zcor
1886	v_qsd(i,k) =FOEDE(v_t,zdelta,zcvm5,v_qs(i,k),zcor)
1887	ELSE
1888	IF (v_t.LT.t_coup) THEN
1889	v_qs(i,k) = qsats(v_t) / v_p
1890	v_qsd(i,k) = dqsats(v_t,v_qs(i,k))
1891	ELSE
1892	v_qs(i,k) = qsatl(v_t) / v_p
1893	v_qsd(i,k) = dqsatl(v_t,v_qs(i,k))
1894	ENDIF
1895	ENDIF
1896	ENDDO
1897	ENDDO
1898	c
1899	c Calculer Gamma * Cp * dz: (gamm est le gradient critique)
1900	c
1901	DO k = 2, klev
1902	DO i = 1, klon
1903	zdp = paprs(i,k)-paprs(i,k+1)
1904	zdpm = paprs(i,k-1)-paprs(i,k)
1905	gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
1906	. (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp)
1907	. +RLVTT /(zdpm+zdp) *
1908	. (v_qs(i,k-1)zdpm + v_qs(i,k)zdp)
1909	. )* (pplay(i,k-1)-pplay(i,k)) / paprs(i,k) )
1910	. / (1.0+(v_qsd(i,k-1)*zdpm+
1911	. v_qsd(i,k)*zdp)/(zdpm+zdp) )
1912	ENDDO
1913	ENDDO
1914	C
1915	C------------------------------------ modification des profils instables
1916	DO 9999 i = 1, klon
1917	itest(i) = .FALSE.
1918	C
1919	k1 = 0
1920	k2 = 1
1921	C
1922	810 CONTINUE ! chercher k1, le bas de la colonne
1923	k2 = k2 + 1
1924	IF (k2 .GT. klev) GOTO 9999
1925	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
1926	zsat=(local_q(i,k2-1)-v_qs(i,k2-1))*(paprs(i,k2-1)-paprs(i,k2))
1927	. +(local_q(i,k2)-v_qs(i,k2))*(paprs(i,k2)-paprs(i,k2+1))
1928	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
1929	k1 = k2 - 1
1930	itest(i) = .TRUE.
1931	C
1932	820 CONTINUE ! chercher k2, le haut de la colonne
1933	IF (k2 .EQ. klev) GOTO 821
1934	k2p = k2 + 1
1935	zsat=zsat +(paprs(i,k2p)-paprs(i,k2p+1))
1936	. *(local_q(i,k2p)-v_qs(i,k2p))
1937	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
1938	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821
1939	k2 = k2p
1940	GOTO 820
1941	821 CONTINUE
1942	C
1943	C------------------------------------------------------ ajustement local
1944	830 CONTINUE ! ajustement proprement dit
1945	v_cptj(k1) = 0.0
1946	zdp = paprs(i,k1)-paprs(i,k1+1)
1947	v_cptjk1 = ( (1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1))
1948	. + RLVTT(local_q(i,k1)-v_qs(i,k1)) ) zdp
1949	v_ssig = zdp * (1.0+v_qsd(i,k1))
1950	C
1951	k1p = k1 + 1
1952	DO k = k1p, k2
1953	zdp = paprs(i,k)-paprs(i,k+1)
1954	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
1955	v_cptjk1 = v_cptjk1 + zdp
1956	. * ( (1.0+v_qsd(i, k))*(v_cptt(i,k)+v_cptj(k))
1957	. + RLVTT*(local_q(i,k)-v_qs(i,k)) )
1958	v_ssig = v_ssig + zdp *(1.0+v_qsd(i,k))
1959	ENDDO
1960	C
1961	DO k = k1, k2
1962	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
1963	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
1964	new_qb(k) = v_qs(i,k) + v_qsd(i,k)*cp_delta_t(k)/RLVTT
1965	local_q(i,k) = new_qb(k)
1966	local_t(i,k) = cp_new_t(k) / RCPD
1967	ENDDO
1968	C
1969	C--------------------------------------------------- sondage vers le bas
1970	C -- on redefinit les variables prognostiques dans
1971	C -- la colonne qui vient d'etre ajustee
1972	C
1973	DO k = k1, k2
1974	v_cptt(i,k) = RCPD * local_t(i,k)
1975	v_t = local_t(i,k)
1976	v_p = pplay(i,k)
1977	C
1978	IF (thermcep) THEN
1979	zdelta=MAX(0.,SIGN(1.,RTT-v_t))
1980	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
1981	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*local_q(i,k))
1982	v_qs(i,k)= R2ES * FOEEW(v_t,zdelta)/v_p
1983	v_qs(i,k)=MIN(0.5,v_qs(i,k))
1984	zcor=1./(1.-RETV*v_qs(i,k))
1985	v_qs(i,k)=v_qs(i,k)*zcor
1986	v_qsd(i,k) =FOEDE(v_t,zdelta,zcvm5,v_qs(i,k),zcor)
1987	ELSE
1988	IF (v_t.LT.t_coup) THEN
1989	v_qs(i,k) = qsats(v_t) / v_p
1990	v_qsd(i,k) = dqsats(v_t,v_qs(i,k))
1991	ELSE
1992	v_qs(i,k) = qsatl(v_t) / v_p
1993	v_qsd(i,k) = dqsatl(v_t,v_qs(i,k))
1994	ENDIF
1995	ENDIF
1996	ENDDO
1997	DO k = 2, klev
1998	zdpm = paprs(i,k-1) - paprs(i,k)
1999	zdp = paprs(i,k) - paprs(i,k+1)
2000	gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
2001	. (v_cptt(i,k-1)zdpm+v_cptt(i,k)zdp)
2002	. +RLVTT /(zdpm+zdp) *
2003	. (v_qs(i,k-1)zdpm+v_qs(i,k)zdp)
2004	. )* (pplay(i,k-1)-pplay(i,k)) / paprs(i,k) )
2005	. / (1.0+(v_qsd(i,k-1)zdpm+v_qsd(i,k)zdp)
2006	. /(zdpm+zdp) )
2007	ENDDO
2008	C
2009	C Verifier si l'on peut etendre la colonne vers le bas
2010	C
2011	IF (k1 .EQ. 1) GOTO 841 ! extension echouee
2012	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
2013	zsat=(local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1))
2014	. + (local_q(i,k1)-v_qs(i,k1))*(paprs(i,k1)-paprs(i,k1+1))
2015	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! extension echouee
2016	C
2017	840 CONTINUE
2018	k1 = k1 - 1
2019	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
2020	zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1))
2021	. *(paprs(i,k1-1)-paprs(i,k1))
2022	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
2023	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
2024	GOTO 840
2025	ELSE
2026	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
2027	ENDIF
2028	841 CONTINUE
2029	C
2030	GOTO 810 ! chercher d'autres blocks en haut
2031	C
2032	9999 CONTINUE ! boucle sur tous les points
2033	C-----------------------------------------------------------------------
2034	c
2035	c Determiner la fraction nuageuse (hypothese: la nebulosite a lieu
2036	c a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
2037	c
2038	DO k = 1, klev
2039	DO i = 1, klon
2040	IF (itest(i)) THEN
2041	delta_q(i,k) = local_q(i,k) - q(i,k)
2042	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
2043	ENDIF
2044	ENDDO
2045	ENDDO
2046	c
2047	c Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
2048	c l'eau liquide est distribuee aux endroits ou la vapeur d'eau
2049	c diminue et d'une maniere proportionnelle a cet diminution):
2050	c
2051	DO i = 1, klon
2052	IF (itest(i)) THEN
2053	zq1(i) = 0.0
2054	zq2(i) = 0.0
2055	ENDIF
2056	ENDDO
2057	DO k = 1, klev
2058	DO i = 1, klon
2059	IF (itest(i)) THEN
2060	zdp = paprs(i,k)-paprs(i,k+1)
2061	zq1(i) = zq1(i) - delta_q(i,k) * zdp
2062	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
2063	ENDIF
2064	ENDDO
2065	ENDDO
2066	DO k = 1, klev
2067	DO i = 1, klon
2068	IF (itest(i)) THEN
2069	IF (zq2(i).NE.0.0)
2070	. d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
2071	ENDIF
2072	ENDDO
2073	ENDDO
2074	C
2075	DO k = 1, klev
2076	DO i = 1, klon
2077	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
2078	ENDDO
2079	ENDDO
2080	C
2081	DO k = 1, klev
2082	DO i = 1, klon
2083	d_t(i,k) = local_t(i,k) - t(i,k)
2084	d_q(i,k) = local_q(i,k) - q(i,k)
2085	ENDDO
2086	ENDDO
2087	c
2088	RETURN
2089	END
2090	SUBROUTINE fiajc(dtime,paprs,pplay,
2091	. t, q,conv_q,
2092	. d_t, d_q, d_ql,rneb,
2093	. rain, snow, ibas, itop)
2094	IMPLICIT NONE
2095	c
2096	#include "dimensions.h"
2097	#include "dimphy.h"
2098	#include "YOMCST.h"
2099	c
2100	c Options:
2101	c
2102	INTEGER plb ! niveau de depart pour la convection
2103	PARAMETER (plb=4)
2104	c
2105	c Mystere: cette option n'est pas innocente pour les resultats !
2106	c Qui peut resoudre ce mystere ? (Z.X.Li mars 1995)
2107	LOGICAL vector ! calcul vectorise
2108	PARAMETER (vector=.FALSE.)
2109	c
2110	REAL t_coup
2111	PARAMETER (t_coup=234.0)
2112	c
2113	c Arguments:
2114	c
2115	REAL q(klon,klev) ! humidite specifique (kg/kg)
2116	REAL t(klon,klev) ! temperature (K)
2117	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
2118	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
2119	REAL dtime ! intervalle du temps (s)
2120	REAL conv_q(klon,klev) ! taux de convergence de l'humidite
2121	REAL rneb(klon,klev) ! fraction nuageuse
2122	REAL d_q(klon,klev) ! incrementaion pour la vapeur d'eau
2123	REAL d_ql(klon,klev) ! incrementation pour l'eau liquide
2124	REAL d_t(klon,klev) ! incrementation pour la temperature
2125	REAL rain(klon) ! variable non-utilisee
2126	REAL snow(klon) ! variable non-utilisee
2127	INTEGER itop(klon) ! variable non-utilisee
2128	INTEGER ibas(klon) ! variable non-utilisee
2129	c
2130	INTEGER kh(klon), i, k
2131	LOGICAL nuage(klon), test(klon,klev)
2132	REAL zconv(klon), zdeh(klon,klev), zvirt(klon)
2133	REAL zdqs(klon,klev), zqs(klon,klev)
2134	REAL ztt, zvar, zfrac(klon)
2135	REAL zq1(klon), zq2(klon)
2136	REAL zdelta, zcor, zcvm5
2137	C
2138	#include "YOETHF.h"
2139	#include "FCTTRE.h"
2140	c
2141	c Initialiser les sorties:
2142	c
2143	DO k = 1, klev
2144	DO i = 1, klon
2145	rneb(i,k) = 0.0
2146	d_ql(i,k) = 0.0
2147	d_t(i,k) = 0.0
2148	d_q(i,k) = 0.0
2149	ENDDO
2150	ENDDO
2151	DO i = 1, klon
2152	itop(i) = 0
2153	ibas(i) = 0
2154	rain(i) = 0.0
2155	snow(i) = 0.0
2156	ENDDO
2157	c
2158	c Calculer Qs et L/Cp * dQs/dT:
2159	c
2160	DO k = 1, klev
2161	DO i = 1, klon
2162	ztt = t(i,k)
2163	IF (thermcep) THEN
2164	zdelta=MAX(0.,SIGN(1.,RTT-ztt))
2165	zcvm5=R5LESRLVTT(1.-zdelta) + zdeltaR5IESRLSTT
2166	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*q(i,k))
2167	zqs(i,k)= R2ES*FOEEW(ztt,zdelta)/pplay(i,k)
2168	zqs(i,k)=MIN(0.5,zqs(i,k))
2169	zcor=1./(1.-RETV*zqs(i,k))
2170	zqs(i,k)=zqs(i,k)*zcor
2171	zdqs(i,k) =FOEDE(ztt,zdelta,zcvm5,zqs(i,k),zcor)
2172	ELSE
2173	IF (ztt .LT. t_coup) THEN
2174	zqs(i,k) = qsats(ztt) / pplay(i,k)
2175	zdqs(i,k) = dqsats(ztt,zqs(i,k))
2176	ELSE
2177	zqs(i,k) = qsatl(ztt) / pplay(i,k)
2178	zdqs(i,k) = dqsatl(ztt,zqs(i,k))
2179	ENDIF
2180	ENDIF
2181	ENDDO
2182	ENDDO
2183	c
2184	c Determiner la difference de l'energie totale saturee:
2185	c
2186	DO i = 1, klon
2187	k = plb
2188	zdeh(i,k) = RCPD * (t(i,k-1)-t(i,k))
2189	. - RD 0.5(t(i,k-1)+t(i,k))/paprs(i,k)
2190	. *(pplay(i,k-1)-pplay(i,k))
2191	. + RLVTT*(zqs(i,k-1)-zqs(i,k))
2192	zdeh(i,k) = zdeh(i,k) * 0.5 ! on prend la moitie
2193	ENDDO
2194	DO k = plb+1, klev
2195	DO i = 1, klon
2196	zdeh(i,k) = zdeh(i,k-1)
2197	. + RCPD * (t(i,k-1)-t(i,k))
2198	. - RD 0.5(t(i,k-1)+t(i,k))/paprs(i,k)
2199	. *(pplay(i,k-1)-pplay(i,k))
2200	. + RLVTT*(zqs(i,k-1)-zqs(i,k))
2201	ENDDO
2202	ENDDO
2203	c
2204	c Determiner le sommet du nuage selon l'instabilite
2205	c Calculer les convergences d'humidite (reelle et virtuelle)
2206	c
2207	DO i = 1, klon
2208	nuage(i) = .TRUE.
2209	zconv(i) = 0.0
2210	zvirt(i) = 0.0
2211	kh(i) = -999
2212	ENDDO
2213	DO k = plb, klev
2214	DO i = 1, klon
2215	nuage(i) = nuage(i) .AND. zdeh(i,k).GT.0.0
2216	IF (nuage(i)) THEN
2217	kh(i) = k
2218	zconv(i) = zconv(i)+conv_q(i,k)*dtime
2219	. *(paprs(i,k)-paprs(i,k+1))
2220	zvirt(i)=zvirt(i)+(zdeh(i,k)/RLVTT+zqs(i,k)-q(i,k))
2221	. *(paprs(i,k)-paprs(i,k+1))
2222	ENDIF
2223	ENDDO
2224	ENDDO
2225	c
2226	IF (vector) THEN
2227	c
2228	c
2229	DO k = plb, klev
2230	DO i = 1, klon
2231	IF (k.LE.kh(i) .AND. kh(i).GT.plb .AND. zconv(i).GT.0.0) THEN
2232	test(i,k) = .TRUE.
2233	zfrac(i) = MAX(0.0,MIN(zconv(i)/zvirt(i),1.0))
2234	ELSE
2235	test(i,k) = .FALSE.
2236	ENDIF
2237	ENDDO
2238	ENDDO
2239	c
2240	DO k = plb, klev
2241	DO i = 1, klon
2242	IF (test(i,k)) THEN
2243	zvar = zdeh(i,k)/(1.0+zdqs(i,k))
2244	d_q(i,k) = (zvarzdqs(i,k)/RLVTT+zqs(i,k)-q(i,k))zfrac(i)
2245	. - conv_q(i,k)*dtime
2246	d_t(i,k) = zvar * zfrac(i) / RCPD
2247	ENDIF
2248	ENDDO
2249	ENDDO
2250	c
2251	DO i = 1, klon
2252	zq1(i) = 0.0
2253	zq2(i) = 0.0
2254	ENDDO
2255	DO k = plb, klev
2256	DO i = 1, klon
2257	IF (test(i,k)) THEN
2258	IF (d_q(i,k).LT.0.0) rneb(i,k) = zfrac(i)
2259	zq1(i) = zq1(i) - d_q(i,k) * (paprs(i,k)-paprs(i,k+1))
2260	zq2(i) = zq2(i) - MIN(0.0, d_q(i,k))
2261	. * (paprs(i,k)-paprs(i,k+1))
2262	ENDIF
2263	ENDDO
2264	ENDDO
2265	c
2266	DO k = plb, klev
2267	DO i = 1, klon
2268	IF (test(i,k)) THEN
2269	IF(zq2(i).NE.0.)d_ql(i,k)=-MIN(0.0,d_q(i,k))*zq1(i)/zq2(i)
2270	ENDIF
2271	ENDDO
2272	ENDDO
2273	c
2274	ELSE ! (.NOT. vector)
2275	c
2276	DO 999 i = 1, klon
2277	IF (kh(i).GT.plb .AND. zconv(i).GT.0.0) THEN
2278	ccc IF (kh(i).LE.plb) GOTO 999 ! il n'y a pas d'instabilite
2279	ccc IF (zconv(i).LE.0.0) GOTO 999 ! convergence insuffisante
2280	zfrac(i) = MAX(0.0,MIN(zconv(i)/zvirt(i),1.0))
2281	DO k = plb, kh(i)
2282	zvar = zdeh(i,k)/(1.0+zdqs(i,k))
2283	d_q(i,k) = (zvarzdqs(i,k)/RLVTT+zqs(i,k)-q(i,k))zfrac(i)
2284	. - conv_q(i,k)*dtime
2285	d_t(i,k) = zvar * zfrac(i) / RCPD
2286	ENDDO
2287	c
2288	zq1(i) = 0.0
2289	zq2(i) = 0.0
2290	DO k = plb, kh(i)
2291	IF (d_q(i,k).LT.0.0) rneb(i,k) = zfrac(i)
2292	zq1(i) = zq1(i) - d_q(i,k) * (paprs(i,k)-paprs(i,k+1))
2293	zq2(i) = zq2(i) - MIN(0.0, d_q(i,k))
2294	. * (paprs(i,k)-paprs(i,k+1))
2295	ENDDO
2296	DO k = plb, kh(i)
2297	IF(zq2(i).NE.0.)d_ql(i,k)=-MIN(0.0,d_q(i,k))*zq1(i)/zq2(i)
2298	ENDDO
2299	ENDIF
2300	999 CONTINUE
2301	c
2302	ENDIF ! fin de teste sur vector
2303	c
2304	RETURN
2305	END

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

Download in other formats:

Original Format