Context Navigation

conlmd.F @ 34

Last change on this file since 34 was 2, checked in by lmdz, 25 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 69.1 KB

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

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

Download in other formats:

Original Format