Context Navigation

conlmd.F @ 1095

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

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

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

Download in other formats:

Original Format