Context Navigation

conlmd.F @ 634

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

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

Download in other formats:

Original Format