Context Navigation

conlmd.f90 @ 5276

Last change on this file since 5276 was 5274, checked in by abarral, 9 months ago
Replace yomcst.h by existing module
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 69.8 KB

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

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

Download in other formats:

Original Format