Context Navigation

conlmd.F90 @ 5501

Last change on this file since 5501 was 5153, checked in by abarral, 6 months ago
Revert FCTTRE to INCLUDE to assess impact of inlining
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: 66.2 KB

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

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

Download in other formats:

Original Format