Context Navigation

conlmd.F @ 1071

Last change on this file since 1071 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 70.0 KB

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

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

Download in other formats:

Original Format