Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

conlmd.F @ 4249

Last change on this file since 4249 was 701, checked in by (none), 19 years ago
This commit was manufactured by cvs2svn to create branch 'pre_V3'.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 69.7 KB

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

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

Download in other formats:

Original Format