Context Navigation

conlmd.F @ 1068

Last change on this file since 1068 was 814, checked in by Laurent Fairhead, 17 years ago
Rajout de la physique utilisant les thermiques FH LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 70.0 KB

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

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

Download in other formats:

Original Format