Context Navigation

conlmd.F @ 3536

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

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

Download in other formats:

Original Format