Context Navigation

conlmd.f90 @ 5441

Last change on this file since 5441 was 5285, checked in by abarral, 8 weeks ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
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: 63.6 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
[5285]	7	USE yomcst_mod_h
[5284]	8	USE yoethf_mod_h
[5274]	9	IMPLICIT NONE
[1992]	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	! ======================================================================
[5274]	15
[1992]	16	! Arguments:
	17
	18	REAL dtime ! pas d'integration (s)
	19	REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
	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)
	33	PARAMETER (usekuo=.TRUE.)
	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
	49	! Manabe
	50	(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, ibas, itop)
	51
	52	IF (usekuo) THEN
	53	! cc CALL fiajc ! ancienne version de Convection Kuo
	54	CALL conkuo & ! nouvelle version de Convection
	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)
	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
	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)))
	80	zz = d_ql(i, k) ! re-evap. de l'eau liquide
	81	zb = max(0.0, zz)
	82	za = -max(0.0, zz)(zlvdcp(1.-zdelta)+zlsdcp*zdelta)
	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	RETURN
	89	END SUBROUTINE conlmd
	90	SUBROUTINE conman(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, &
	91	snow, ibas, itop)
	92	USE dimphy
[5285]	93	USE yomcst_mod_h
[5284]	94	USE yoethf_mod_h
[1992]	95	IMPLICIT NONE
	96	! ======================================================================
	97	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19970324
	98	! Objet: ajustement humide convectif avec la possibilite de faire
	99	! l'ajustement sur une fraction de la maille.
	100	! Methode: On impose une distribution uniforme pour la vapeur d'eau
	101	! au sein d'une maille. On applique la procedure d'ajustement
	102	! successivement a la totalite, 75%, 50%, 25% et 5% de la maille
	103	! jusqu'a ce que l'ajustement a lieu. J'espere que ceci augmente
	104	! les activites convectives et corrige le biais "trop froid et sec"
	105	! du modele.
	106	! ======================================================================
	107	REAL dtime ! pas d'integration (s)
	108	REAL t(klon, klev) ! temperature (K)
	109	REAL q(klon, klev) ! humidite specifique (kg/kg)
	110	REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
	111	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	112
	113	REAL d_t(klon, klev) ! incrementation temperature
	114	REAL d_q(klon, klev) ! incrementation humidite
	115	REAL d_ql(klon, klev) ! incrementation eau liquide
	116	REAL rneb(klon, klev) ! nebulosite
	117	REAL rain(klon) ! pluies (mm/s)
	118	REAL snow(klon) ! neige (mm/s)
	119	INTEGER ibas(klon) ! niveau du bas
	120	INTEGER itop(klon) ! niveau du haut
	121
	122	LOGICAL afaire(klon) ! .TRUE. implique l'ajustement
	123	LOGICAL accompli(klon) ! .TRUE. si l'ajustement est effectif
	124
	125	INTEGER nb ! nombre de sous-fractions a considere
	126	PARAMETER (nb=1)
	127	! cc PARAMETER (nb=3)
	128
	129	REAL ratqs ! largeur de la distribution pour vapeur d'eau
	130	PARAMETER (ratqs=0.05)
	131
	132	REAL w_q(klon, klev)
	133	REAL w_d_t(klon, klev), w_d_q(klon, klev), w_d_ql(klon, klev)
	134	REAL w_rneb(klon, klev)
	135	REAL w_rain(klon), w_snow(klon)
	136	INTEGER w_ibas(klon), w_itop(klon)
	137	REAL zq1, zq2
	138	INTEGER i, k, n
	139
	140	REAL t_coup
	141	PARAMETER (t_coup=234.0)
	142	REAL zdp1, zdp2
	143	REAL zqs1, zqs2, zdqs1, zdqs2
	144	REAL zgamdz
	145	REAL zflo ! flotabilite
	146	REAL zsat ! sur-saturation
	147	REAL zdelta, zcor, zcvm5
	148	LOGICAL imprim
	149
	150	INTEGER ncpt
	151	SAVE ncpt
	152	!$OMP THREADPRIVATE(ncpt)
	153	REAL frac(nb) ! valeur de la maille fractionnelle
	154	SAVE frac
	155	!$OMP THREADPRIVATE(frac)
	156	INTEGER opt_cld(nb) ! option pour le modele nuageux
	157	SAVE opt_cld
	158	!$OMP THREADPRIVATE(opt_cld)
	159	LOGICAL appel1er
	160	SAVE appel1er
	161	!$OMP THREADPRIVATE(appel1er)
	162
	163	! Fonctions thermodynamiques:
	164
	165	include "FCTTRE.h"
	166
	167	DATA frac/1.0/
	168	DATA opt_cld/4/
	169	! cc DATA frac / 1.0, 0.50, 0.25/
	170	! cc DATA opt_cld / 4, 4, 4/
	171
	172	DATA appel1er/.TRUE./
	173	DATA ncpt/0/
	174
	175	IF (appel1er) THEN
	176	PRINT *, 'conman, nb:', nb
	177	PRINT *, 'conman, frac:', frac
	178	PRINT *, 'conman, opt_cld:', opt_cld
	179	appel1er = .FALSE.
	180	END IF
	181
	182	! Initialiser les sorties a zero:
	183
	184	DO k = 1, klev
	185	DO i = 1, klon
	186	d_t(i, k) = 0.0
	187	d_q(i, k) = 0.0
	188	d_ql(i, k) = 0.0
	189	rneb(i, k) = 0.0
	190	END DO
	191	END DO
	192	DO i = 1, klon
	193	ibas(i) = klev
	194	itop(i) = 1
	195	rain(i) = 0.0
	196	snow(i) = 0.0
	197	END DO
	198
	199	! S'il n'y a pas d'instabilite conditionnelle,
	200	! pas la penne de se fatiguer:
	201
	202	DO i = 1, klon
	203	afaire(i) = .FALSE.
	204	END DO
	205	DO k = 1, klev - 1
	206	DO i = 1, klon
	207	IF (thermcep) THEN
	208	zdelta = max(0., sign(1.,rtt-t(i,k)))
	209	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	210	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
	211	zqs1 = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
	212	zqs1 = min(0.5, zqs1)
	213	zcor = 1./(1.-retv*zqs1)
	214	zqs1 = zqs1*zcor
	215	zdqs1 = foede(t(i,k), zdelta, zcvm5, zqs1, zcor)
	216
	217	zdelta = max(0., sign(1.,rtt-t(i,k+1)))
	218	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	219	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k+1))
	220	zqs2 = r2es*foeew(t(i,k+1), zdelta)/pplay(i, k+1)
	221	zqs2 = min(0.5, zqs2)
	222	zcor = 1./(1.-retv*zqs2)
	223	zqs2 = zqs2*zcor
	224	zdqs2 = foede(t(i,k+1), zdelta, zcvm5, zqs2, zcor)
	225	ELSE
	226	IF (t(i,k)<t_coup) THEN
	227	zqs1 = qsats(t(i,k))/pplay(i, k)
	228	zdqs1 = dqsats(t(i,k), zqs1)
	229
	230	zqs2 = qsats(t(i,k+1))/pplay(i, k+1)
	231	zdqs2 = dqsats(t(i,k+1), zqs2)
	232	ELSE
	233	zqs1 = qsatl(t(i,k))/pplay(i, k)
	234	zdqs1 = dqsatl(t(i,k), zqs1)
	235
	236	zqs2 = qsatl(t(i,k+1))/pplay(i, k+1)
	237	zdqs2 = dqsatl(t(i,k+1), zqs2)
	238	END IF
	239	END IF
	240	zdp1 = paprs(i, k) - paprs(i, k+1)
	241	zdp2 = paprs(i, k+1) - paprs(i, k+2)
	242	zgamdz = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd(rd(t(i, &
	243	k)zdp1+t(i,k+1)zdp2)/(zdp1+zdp2)+rlvtt(zqs1zdp1+zqs2*zdp2)/(zdp1+ &
	244	zdp2))/(1.0+(zdqs1zdp1+zdqs2zdp2)/(zdp1+zdp2))
	245	zflo = t(i, k) + zgamdz - t(i, k+1)
	246	zsat = (q(i,k)-zqs1)zdp1 + (q(i,k+1)-zqs2)zdp2
	247	IF (zflo>0.0) afaire(i) = .TRUE.
	248	! erreur IF (zflo.GT.0.0 .AND. zsat.GT.0.0) afaire(i) = .TRUE.
	249	END DO
	250	END DO
	251
	252	imprim = mod(ncpt, 48) == 0
	253	DO n = 1, nb
	254
	255	DO k = 1, klev
[524]	256	DO i = 1, klon
[1992]	257	IF (afaire(i)) THEN
	258	zq1 = q(i, k)*(1.0-ratqs)
	259	zq2 = q(i, k)*(1.0+ratqs)
	260	w_q(i, k) = zq2 - frac(n)/2.0*(zq2-zq1)
	261	END IF
	262	END DO
	263	END DO
	264
	265	CALL conmanv(dtime, paprs, pplay, t, w_q, afaire, opt_cld(n), w_d_t, &
	266	w_d_q, w_d_ql, w_rneb, w_rain, w_snow, w_ibas, w_itop, accompli, &
	267	imprim)
	268	DO k = 1, klev
[524]	269	DO i = 1, klon
[1992]	270	IF (afaire(i) .AND. accompli(i)) THEN
	271	d_t(i, k) = w_d_t(i, k)*frac(n)
	272	d_q(i, k) = w_d_q(i, k)*frac(n)
	273	d_ql(i, k) = w_d_ql(i, k)*frac(n)
	274	IF (nint(w_rneb(i,k))==1) rneb(i, k) = frac(n)
	275	END IF
	276	END DO
	277	END DO
	278	DO i = 1, klon
[524]	279	IF (afaire(i) .AND. accompli(i)) THEN
[1992]	280	rain(i) = w_rain(i)*frac(n)
	281	snow(i) = w_snow(i)*frac(n)
	282	ibas(i) = min(ibas(i), w_ibas(i))
	283	itop(i) = max(itop(i), w_itop(i))
	284	END IF
	285	END DO
	286	DO i = 1, klon
	287	IF (afaire(i) .AND. accompli(i)) afaire(i) = .FALSE.
	288	END DO
	289
	290	END DO
	291
	292	ncpt = ncpt + 1
	293
	294	RETURN
	295	END SUBROUTINE conman
	296	SUBROUTINE conmanv(dtime, paprs, pplay, t, q, afaire, opt_cld, d_t, d_q, &
	297	d_ql, rneb, rain, snow, ibas, itop, accompli, imprim)
	298	USE dimphy
[5285]	299	USE yomcst_mod_h
[5284]	300	USE yoethf_mod_h
[5274]	301	IMPLICIT NONE
[1992]	302	! ======================================================================
	303	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	304	! Objet: ajustement humide (convection proposee par Manabe).
	305	! Pour une colonne verticale, il peut avoir plusieurs blocs
	306	! necessitant l'ajustement. ibas est le bas du plus bas bloc
	307	! et itop est le haut du plus haut bloc
	308	! ======================================================================
	309
[5274]	310
[1992]	311	! Arguments:
	312
	313	REAL dtime ! pas d'integration (s)
	314	REAL t(klon, klev) ! temperature (K)
	315	REAL q(klon, klev) ! humidite specifique (kg/kg)
	316	REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
	317	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	318	INTEGER opt_cld ! comment traiter l'eau liquide
	319	LOGICAL afaire(klon) ! .TRUE. si le point est a faire (Input)
	320	LOGICAL imprim ! .T. pour imprimer quelques diagnostiques
	321
	322	REAL d_t(klon, klev) ! incrementation temperature
	323	REAL d_q(klon, klev) ! incrementation humidite
	324	REAL d_ql(klon, klev) ! incrementation eau liquide
	325	REAL rneb(klon, klev) ! nebulosite
	326	REAL rain(klon) ! pluies (mm/s)
	327	REAL snow(klon) ! neige (mm/s)
	328	INTEGER ibas(klon) ! niveau du bas
	329	INTEGER itop(klon) ! niveau du haut
	330	LOGICAL accompli(klon) ! .TRUE. si l'ajustement a eu lieu (Output)
	331
	332	! Quelques options:
	333
	334	LOGICAL new_top ! re-calculer sommet quand re-ajustement est fait
	335	PARAMETER (new_top=.FALSE.)
	336	LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
	337	PARAMETER (evap_prec=.TRUE.)
	338	REAL coef_eva
	339	PARAMETER (coef_eva=1.0E-05)
	340	REAL t_coup
	341	PARAMETER (t_coup=234.0)
	342	REAL seuil_vap
	343	PARAMETER (seuil_vap=1.0E-10)
	344	LOGICAL old_tau ! implique precip nulle, si vrai.
	345	PARAMETER (old_tau=.FALSE.)
	346	REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
	347	REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
	348	PARAMETER (dpmin=0.15, tomax=0.97)
	349	REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
	350	PARAMETER (dpmax=0.30, tomin=0.05)
	351	REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
	352	PARAMETER (deep_sig=0.50, deep_to=0.05)
	353	LOGICAL exigent ! implique un calcul supplementaire pour Qs
	354	PARAMETER (exigent=.FALSE.)
	355
	356	INTEGER kbase
	357	PARAMETER (kbase=0)
	358
	359	! Variables locales:
	360
	361	INTEGER nexpo
	362	INTEGER i, k, k1min, k1max, k2min, k2max, is
	363	REAL zgamdz(klon, klev-1)
	364	REAL zt(klon, klev), zq(klon, klev)
	365	REAL zqs(klon, klev), zdqs(klon, klev)
	366	REAL zqmqsdp(klon, klev)
	367	REAL ztnew(klon, klev), zqnew(klon, klev)
	368	REAL zcond(klon), zvapo(klon), zrapp(klon)
	369	REAL zrfl(klon), zrfln, zqev, zqevt
	370	REAL zsat(klon) ! sur-saturation
	371	REAL zflo(klon) ! flotabilite
	372	REAL za(klon), zb(klon), zc(klon)
	373	INTEGER k1(klon), k2(klon)
	374	REAL zdelta, zcor, zcvm5
	375	REAL delp(klon, klev)
	376	LOGICAL possible(klon), todo(klon), etendre(klon)
	377	LOGICAL aller(klon), todobis(klon)
	378	REAL zalfa
	379	INTEGER nbtodo, nbdone
	380
	381	! Fonctions thermodynamiques:
	382
	383	include "FCTTRE.h"
	384
	385	DO k = 1, klev
	386	DO i = 1, klon
	387	delp(i, k) = paprs(i, k) - paprs(i, k+1)
	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
	420	zdelta = max(0., sign(1.,rtt-zt(i,k)))
	421	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	422	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i,k))
	423	zqs(i, k) = r2es*foeew(zt(i,k), zdelta)/pplay(i, k)
	424	zqs(i, k) = min(0.5, zqs(i,k))
	425	zcor = 1./(1.-retv*zqs(i,k))
	426	zqs(i, k) = zqs(i, k)*zcor
	427	zdqs(i, k) = foede(zt(i,k), zdelta, zcvm5, zqs(i,k), zcor)
	428	ELSE
	429	IF (zt(i,k)<t_coup) THEN
	430	zqs(i, k) = qsats(zt(i,k))/pplay(i, k)
	431	zdqs(i, k) = dqsats(zt(i,k), zqs(i,k))
	432	ELSE
	433	zqs(i, k) = qsatl(zt(i,k))/pplay(i, k)
	434	zdqs(i, k) = dqsatl(zt(i,k), zqs(i,k))
	435	END IF
	436	END IF
	437
	438	! Calculer (q-qs)*dp
	439	zqmqsdp(i, k) = (zq(i,k)-zqs(i,k))*delp(i, k)
	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
[1992]	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)))
	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
[1992]	470	zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
	471	zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k-1)
	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
	484	810 CONTINUE ! chercher le bas de la colonne a ajuster
	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
	496	zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k+1)
	497	zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k+1)
	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
	532	820 CONTINUE ! chercher le haut de la colonne
	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)
	544	zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
	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
	581	830 CONTINUE ! faire l'ajustement en sachant k1 et k2
	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.
[1992]	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))
	613	END IF
	614	END DO
	615
	616	DO k = k1min, k2max
	617	DO i = 1, klon
	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))
	623	END IF
	624	END DO
	625	END DO
	626
	627	DO i = 1, klon
	628	IF (todo(i)) THEN
	629	k = k1(i)
	630	ztnew(i, k) = zb(i)/zc(i)
	631	zqnew(i, k) = zqs(i, k) + (ztnew(i,k)-zt(i,k))rcpd/rlvttzdqs(i, k)
	632	END IF
	633	END DO
	634
	635	DO k = k1min, k2max
	636	DO i = 1, klon
	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/rlvttzdqs(i, k)
	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
	653	zcond(i) = zcond(i) + (zq(i,k)-zqnew(i,k))*delp(i, k)/rg
	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
	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
	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
	708	IF (todo(i)) zrfl(i) = zcond(i)/dtime
	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) + &
	719	zqnew(i, k)*delp(i, k)/rg
	720	END DO
	721	END DO
	722	DO i = 1, klon
[524]	723	IF (todo(i)) THEN
[1992]	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
	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
	732	d_ql(i, k) = d_ql(i, k) + zrapp(i)*zqnew(i, k)
	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) + &
	745	delp(i, k)/rg
	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
	751	d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)
	752	END IF
	753	END DO
	754	END DO
	755	DO i = 1, klon
	756	IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	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) + &
	767	max(0.0, zq(i,k)-zqnew(i,k))*delp(i, k)/rg
	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, &
	773	k) = d_ql(i, k) + toliq(i)zcond(i)/zvapo(i)max(0.0, zq(i,k)-zqnew &
	774	(i,k))
	775	END DO
	776	END DO
	777	DO i = 1, klon
	778	IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	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
[1992]	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
	793	END DO
	794	END DO
	795	DO k = k1min, k2max
[524]	796	DO i = 1, klon
[1992]	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
	799	END DO
	800	END DO
	801	DO i = 1, klon
	802	IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	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
	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)))rgdtime/delp(i, k)
	823	zqev = min(zqev, zqevt)
	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)))dtimerlvtt/rcpd/(1.0+rvtmp2*zq(i,k))
	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
	838	IF (zt(i,1)>rtt) THEN
	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
	864	zdelta = max(0., sign(1.,rtt-zt(i,k)))
	865	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	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)
	872	ELSE
	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))
	876	ELSE
	877	zqs(i, k) = qsatl(zt(i,k))/pplay(i, k)
	878	zdqs(i, k) = dqsatl(zt(i,k), zqs(i,k))
	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
	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)))
	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
[1992]	905	zqmqsdp(i, k) = (zq(i,k)-zqs(i,k))*delp(i, k)
	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
	918	IF (todo(i) .AND. k1(i)>(kbase+1)) THEN
	919	k = k1(i)
	920	zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
	921	zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k-1)
	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
	934	IF (k1max>(kbase+1)) THEN
	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)
	939	zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k+1)
	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
	1001	860 CONTINUE ! Calculer les tendances et diagnostiques
	1002	! cc print*, "Apres 860"
	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)
	1008	zq(i, k) = max(zq(i,k), seuil_vap)
	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
[1992]	1017	IF (rneb(i,k)>0.0) THEN
	1018	ibas(i) = k
	1019	GO TO 807
	1020	END IF
	1021	END DO
	1022	807 CONTINUE
	1023	DO k = klev, 1, -1
	1024	IF (rneb(i,k)>0.0) THEN
	1025	itop(i) = k
	1026	GO TO 808
	1027	END IF
	1028	END DO
	1029	808 CONTINUE
	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	RETURN
	1044	END SUBROUTINE conmanv
	1045	SUBROUTINE conkuo(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
	1046	rain, snow, ibas, itop)
	1047	USE dimphy
[5285]	1048	USE yomcst_mod_h
[5284]	1049	USE yoethf_mod_h
[5274]	1050	IMPLICIT NONE
[1992]	1051	! ======================================================================
	1052	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	1053	! Objet: Schema de convection de type Kuo (1965).
	1054	! Cette version du code peut calculer le niveau de depart
	1055	! N.B. version vectorielle (le 6 oct. 1997)
	1056	! ======================================================================
	1057
[5274]	1058
[1992]	1059	! Arguments:
	1060
	1061	REAL dtime ! intervalle du temps (s)
	1062	REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
	1063	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	1064	REAL t(klon, klev) ! temperature (K)
	1065	REAL q(klon, klev) ! humidite specifique
	1066	REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)
	1067
	1068	REAL d_t(klon, klev) ! incrementation temperature
	1069	REAL d_q(klon, klev) ! incrementation humidite
	1070	REAL d_ql(klon, klev) ! incrementation eau liquide
	1071	REAL rneb(klon, klev) ! nebulosite
	1072	REAL rain(klon) ! pluies (mm/s)
	1073	REAL snow(klon) ! neige (mm/s)
	1074	INTEGER itop(klon) ! niveau du sommet
	1075	INTEGER ibas(klon) ! niveau du bas
	1076
	1077	LOGICAL ldcum(klon) ! convection existe
	1078	LOGICAL todo(klon)
	1079
	1080	! Quelsques options:
	1081
	1082	LOGICAL calcfcl ! calculer le niveau de convection libre
	1083	PARAMETER (calcfcl=.TRUE.)
	1084	INTEGER ldepar ! niveau fixe de convection libre
	1085	PARAMETER (ldepar=4)
	1086	INTEGER opt_cld ! comment traiter l'eau liquide
	1087	PARAMETER (opt_cld=4) ! valeur possible: 0, 1, 2, 3 ou 4
	1088	LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
	1089	PARAMETER (evap_prec=.TRUE.)
	1090	REAL coef_eva
	1091	PARAMETER (coef_eva=1.0E-05)
	1092	LOGICAL new_deh ! nouvelle facon de calculer dH
	1093	PARAMETER (new_deh=.FALSE.)
	1094	REAL t_coup
	1095	PARAMETER (t_coup=234.0)
	1096	LOGICAL old_tau ! implique precipitation nulle
	1097	PARAMETER (old_tau=.FALSE.)
	1098	REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
	1099	REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
	1100	PARAMETER (dpmin=0.15, tomax=0.97)
	1101	REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
	1102	PARAMETER (dpmax=0.30, tomin=0.05)
	1103	REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
	1104	PARAMETER (deep_sig=0.50, deep_to=0.05)
	1105
	1106	! Variables locales:
	1107
	1108	INTEGER nexpo
	1109	LOGICAL nuage(klon)
	1110	INTEGER i, k, kbmin, kbmax, khmax
	1111	REAL ztotal(klon, klev), zdeh(klon, klev)
	1112	REAL zgz(klon, klev)
	1113	REAL zqs(klon, klev)
	1114	REAL zdqs(klon, klev)
	1115	REAL ztemp(klon, klev)
	1116	REAL zpres(klon, klev)
	1117	REAL zconv(klon) ! convergence d'humidite
	1118	REAL zvirt(klon) ! convergence virtuelle d'humidite
	1119	REAL zfrac(klon) ! fraction convective
	1120	INTEGER kb(klon), kh(klon)
	1121
	1122	REAL zcond(klon), zvapo(klon), zrapp(klon)
	1123	REAL zrfl(klon), zrfln, zqev, zqevt
	1124	REAL zdelta, zcvm5, zcor
	1125	REAL zvar
	1126
	1127	LOGICAL appel1er
	1128	SAVE appel1er
	1129	!$OMP THREADPRIVATE(appel1er)
	1130
	1131	! Fonctions thermodynamiques
	1132
	1133	include "FCTTRE.h"
	1134
	1135	DATA appel1er/.TRUE./
	1136
	1137	IF (appel1er) THEN
	1138	PRINT *, 'conkuo, calcfcl:', calcfcl
	1139	IF (.NOT. calcfcl) PRINT *, 'conkuo, ldepar:', ldepar
	1140	PRINT *, 'conkuo, opt_cld:', opt_cld
	1141	PRINT *, 'conkuo, evap_prec:', evap_prec
	1142	PRINT *, 'conkuo, new_deh:', new_deh
	1143	appel1er = .FALSE.
	1144	END IF
	1145
	1146	! Initialiser les sorties a zero
	1147
	1148	DO k = 1, klev
	1149	DO i = 1, klon
	1150	d_q(i, k) = 0.0
	1151	d_t(i, k) = 0.0
	1152	d_ql(i, k) = 0.0
	1153	rneb(i, k) = 0.0
	1154	END DO
	1155	END DO
	1156	DO i = 1, klon
	1157	rain(i) = 0.0
	1158	snow(i) = 0.0
	1159	ibas(i) = 0
	1160	itop(i) = 0
	1161	END DO
	1162
	1163	! Calculer la vapeur d'eau saturante Qs et sa derive L/Cp * dQs/dT
	1164
	1165	DO k = 1, klev
	1166	DO i = 1, klon
	1167	IF (thermcep) THEN
	1168	zdelta = max(0., sign(1.,rtt-t(i,k)))
	1169	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1170	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
	1171	zqs(i, k) = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
	1172	zqs(i, k) = min(0.5, zqs(i,k))
	1173	zcor = 1./(1.-retv*zqs(i,k))
	1174	zqs(i, k) = zqs(i, k)*zcor
	1175	zdqs(i, k) = foede(t(i,k), zdelta, zcvm5, zqs(i,k), zcor)
	1176	ELSE
	1177	IF (t(i,k)<t_coup) THEN
	1178	zqs(i, k) = qsats(t(i,k))/pplay(i, k)
	1179	zdqs(i, k) = dqsats(t(i,k), zqs(i,k))
	1180	ELSE
	1181	zqs(i, k) = qsatl(t(i,k))/pplay(i, k)
	1182	zdqs(i, k) = dqsatl(t(i,k), zqs(i,k))
	1183	END IF
	1184	END IF
	1185	END DO
	1186	END DO
	1187
	1188	! Calculer gz (energie potentielle)
	1189
	1190	DO i = 1, klon
	1191	zgz(i, 1) = rdt(i, 1)/(0.5(paprs(i,1)+pplay(i, &
	1192	1)))*(paprs(i,1)-pplay(i,1))
	1193	END DO
	1194	DO k = 2, klev
	1195	DO i = 1, klon
	1196	zgz(i, k) = zgz(i, k-1) + rd0.5(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i &
	1197	,k-1)-pplay(i,k))
	1198	END DO
	1199	END DO
	1200
	1201	! Calculer l'energie statique humide saturee (CpT + gz + LQs)
	1202
	1203	DO k = 1, klev
	1204	DO i = 1, klon
	1205	ztotal(i, k) = rcpdt(i, k) + rlvttzqs(i, k) + zgz(i, k)
	1206	END DO
	1207	END DO
	1208
	1209	! Determiner le niveau de depart et calculer la difference de
	1210	! l'energie statique humide saturee (ztotal) entre la couche
	1211	! de depart et chaque couche au-dessus.
	1212
	1213	IF (calcfcl) THEN
	1214	DO k = 1, klev
[524]	1215	DO i = 1, klon
[1992]	1216	zpres(i, k) = pplay(i, k)
	1217	ztemp(i, k) = t(i, k)
	1218	END DO
	1219	END DO
	1220	CALL kuofcl(ztemp, q, zgz, zpres, ldcum, kb)
	1221	DO i = 1, klon
	1222	IF (ldcum(i)) THEN
	1223	k = kb(i)
	1224	IF (new_deh) THEN
	1225	zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
	1226	ELSE
	1227	zdeh(i, k) = rcpd(t(i,k-1)-t(i,k)) - rd0.5*(t(i,k-1)+t(i,k))/ &
	1228	paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
	1229	rlvtt*(zqs(i,k-1)-zqs(i,k))
	1230	END IF
	1231	zdeh(i, k) = zdeh(i, k)*0.5
	1232	END IF
	1233	END DO
	1234	DO k = 1, klev
[524]	1235	DO i = 1, klon
[1992]	1236	IF (ldcum(i) .AND. k>=(kb(i)+1)) THEN
	1237	IF (new_deh) THEN
	1238	zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
	1239	ELSE
	1240	zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
	1241	rd0.5(t(i,k-1)+t(i,k))/paprs(i, k)* &
	1242	(pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
	1243	END IF
	1244	END IF
	1245	END DO
	1246	END DO
	1247	ELSE
	1248	DO i = 1, klon
	1249	k = ldepar
	1250	kb(i) = ldepar
	1251	ldcum(i) = .TRUE.
	1252	IF (new_deh) THEN
	1253	zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
[524]	1254	ELSE
[1992]	1255	zdeh(i, k) = rcpd(t(i,k-1)-t(i,k)) - rd0.5*(t(i,k-1)+t(i,k))/paprs( &
	1256	i, k)(pplay(i,k-1)-pplay(i,k)) + rlvtt(zqs(i,k-1)-zqs(i,k))
	1257	END IF
	1258	zdeh(i, k) = zdeh(i, k)*0.5
	1259	END DO
	1260	DO k = ldepar + 1, klev
[524]	1261	DO i = 1, klon
[1992]	1262	IF (new_deh) THEN
	1263	zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
	1264	ELSE
	1265	zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
	1266	rd0.5(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
	1267	rlvtt*(zqs(i,k-1)-zqs(i,k))
	1268	END IF
	1269	END DO
	1270	END DO
	1271	END IF
	1272
	1273	! -----Chercher le sommet du nuage
	1274	! -----Calculer la convergence de l'humidite (en kg/m**2 a un facteur
	1275	! -----psolpa/RG pres) du bas jusqu'au sommet du nuage.
	1276	! -----Calculer la convergence virtuelle pour que toute la maille
	1277	! -----deviennt nuageuse (du bas jusqu'au sommet du nuage)
	1278
	1279	DO i = 1, klon
	1280	nuage(i) = .TRUE.
	1281	zconv(i) = 0.0
	1282	zvirt(i) = 0.0
	1283	kh(i) = -999
	1284	END DO
	1285	DO k = 1, klev
	1286	DO i = 1, klon
	1287	IF (k>=kb(i) .AND. ldcum(i)) THEN
	1288	nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
	1289	IF (nuage(i)) THEN
	1290	kh(i) = k
	1291	zconv(i) = zconv(i) + conv_q(i, k)dtime(paprs(i,k)-paprs(i,k+1))
	1292	zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k) &
	1293	-paprs(i,k+1))
	1294	END IF
	1295	END IF
	1296	END DO
	1297	END DO
	1298
	1299	DO i = 1, klon
	1300	todo(i) = ldcum(i) .AND. kh(i) > kb(i) .AND. zconv(i) > 0.0
	1301	END DO
	1302
	1303	kbmin = klev
	1304	kbmax = 0
	1305	khmax = 0
	1306	DO i = 1, klon
	1307	IF (todo(i)) THEN
	1308	kbmin = min(kbmin, kb(i))
	1309	kbmax = max(kbmax, kb(i))
	1310	khmax = max(khmax, kh(i))
	1311	END IF
	1312	END DO
	1313
	1314	! -----Calculer la surface couverte par le nuage
	1315
	1316	DO i = 1, klon
	1317	IF (todo(i)) THEN
	1318	zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
	1319	END IF
	1320	END DO
	1321
	1322	! -----Calculs essentiels:
	1323
	1324	DO i = 1, klon
	1325	IF (todo(i)) THEN
	1326	zcond(i) = 0.0
	1327	END IF
	1328	END DO
	1329	DO k = kbmin, khmax
	1330	DO i = 1, klon
	1331	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1332	zvar = zdeh(i, k)/(1.+zdqs(i,k))
	1333	d_t(i, k) = zvar*zfrac(i)/rcpd
	1334	d_q(i, k) = (zvarzdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))zfrac(i) - &
	1335	conv_q(i, k)*dtime
	1336	zcond(i) = zcond(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
	1337	rneb(i, k) = zfrac(i)
	1338	END IF
	1339	END DO
	1340	END DO
	1341
	1342	DO i = 1, klon
	1343	IF (todo(i) .AND. zcond(i)<0.0) THEN
	1344	PRINT *, 'WARNING: cond. negative (Kuo) ', i, kb(i), kh(i), zcond(i)
	1345	zcond(i) = 0.0
	1346	DO k = kb(i), kh(i)
	1347	d_t(i, k) = 0.0
	1348	d_q(i, k) = 0.0
	1349	END DO
	1350	todo(i) = .FALSE. ! effort totalement perdu
	1351	END IF
	1352	END DO
	1353
	1354	! =====
	1355	! Une fois que la condensation a lieu, on doit construire un
	1356	! "modele nuageux" pour partager la condensation entre l'eau
	1357	! liquide nuageuse et la precipitation (leur rapport toliq
	1358	! est calcule selon l'epaisseur nuageuse). Je suppose que
	1359	! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
	1360	! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
	1361	! lineaire entre dpmin et dpmax).
	1362	! =====
	1363	DO i = 1, klon
	1364	IF (todo(i)) THEN
	1365	toliq(i) = tomax - ((paprs(i,kb(i))-paprs(i,kh(i)+1))/paprs(i,1)-dpmin) &
	1366	*(tomax-tomin)/(dpmax-dpmin)
	1367	toliq(i) = max(tomin, min(tomax,toliq(i)))
	1368	IF (pplay(i,kh(i))/paprs(i,1)<=deep_sig) toliq(i) = deep_to
	1369	IF (old_tau) toliq(i) = 1.0
	1370	END IF
	1371	END DO
	1372	! =====
	1373	! On doit aussi determiner la distribution verticale de
	1374	! l'eau nuageuse. Plusieurs options sont proposees:
	1375
	1376	! (0) La condensation precipite integralement (toliq ne sera
	1377	! pas utilise).
	1378	! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
	1379	! a la vapeur d'eau locale.
	1380	! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
	1381	! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
	1382	! est effectivement diminuee pendant le processus d'ajustement.
	1383	! (4) Elle est en fonction (lineaire ou exponentielle) de la
	1384	! distance (epaisseur en pression) avec le niveau k1 (la couche
	1385	! k1 n'aura donc pas d'eau liquide).
	1386	! =====
	1387
	1388	IF (opt_cld==0) THEN
	1389
	1390	DO i = 1, klon
	1391	IF (todo(i)) zrfl(i) = zcond(i)/dtime
	1392	END DO
	1393
	1394	ELSE IF (opt_cld==1) THEN
	1395
	1396	DO i = 1, klon
	1397	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
	1398	END DO
	1399	DO k = kbmin, khmax
[524]	1400	DO i = 1, klon
[1992]	1401	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1402	zvapo(i) = zvapo(i) + (q(i,k)+d_q(i,k))*(paprs(i,k)-paprs(i,k+1))/ &
	1403	rg
	1404	END IF
	1405	END DO
	1406	END DO
	1407	DO i = 1, klon
[524]	1408	IF (todo(i)) THEN
[1992]	1409	zrapp(i) = toliq(i)*zcond(i)/zvapo(i)
	1410	zrapp(i) = max(0., min(1.,zrapp(i)))
	1411	END IF
	1412	END DO
	1413	DO k = kbmin, khmax
[524]	1414	DO i = 1, klon
[1992]	1415	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1416	d_ql(i, k) = zrapp(i)*(q(i,k)+d_q(i,k))
	1417	END IF
	1418	END DO
	1419	END DO
	1420	DO i = 1, klon
[524]	1421	IF (todo(i)) THEN
[1992]	1422	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	1423	END IF
	1424	END DO
	1425
	1426	ELSE IF (opt_cld==2) THEN
	1427
	1428	DO i = 1, klon
	1429	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
	1430	END DO
	1431	DO k = kbmin, khmax
[524]	1432	DO i = 1, klon
[1992]	1433	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1434	zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg
	1435	END IF
	1436	END DO
	1437	END DO
	1438	DO k = kbmin, khmax
	1439	DO i = 1, klon
	1440	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1441	d_ql(i, k) = toliq(i)*zcond(i)/zvapo(i)
	1442	END IF
	1443	END DO
	1444	END DO
	1445	DO i = 1, klon
[524]	1446	IF (todo(i)) THEN
[1992]	1447	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	1448	END IF
	1449	END DO
	1450
	1451	ELSE IF (opt_cld==3) THEN
	1452
	1453	DO i = 1, klon
	1454	IF (todo(i)) THEN
	1455	zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
	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	zvapo(i) = zvapo(i) + max(0.0, -d_q(i,k))*(paprs(i,k)-paprs(i,k+1)) &
	1462	/rg
	1463	END IF
	1464	END DO
	1465	END DO
	1466	DO k = kbmin, khmax
[524]	1467	DO i = 1, klon
[1992]	1468	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i) .AND. zvapo(i)>0.0) THEN
	1469	d_ql(i, k) = d_ql(i, k) + toliq(i)zcond(i)/zvapo(i)max(0.0, -d_q( &
	1470	i,k))
	1471	END IF
	1472	END DO
	1473	END DO
	1474	DO i = 1, klon
[524]	1475	IF (todo(i)) THEN
[1992]	1476	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	1477	END IF
	1478	END DO
	1479
	1480	ELSE IF (opt_cld==4) THEN
	1481
	1482	nexpo = 3
	1483	! cc nexpo = 1 ! distribution lineaire
	1484
	1485	DO i = 1, klon
	1486	IF (todo(i)) THEN
	1487	zvapo(i) = 0.0 ! quantite integrale de masse (avec ponderation)
	1488	END IF
	1489	END DO
	1490	DO k = kbmin, khmax
[524]	1491	DO i = 1, klon
[1992]	1492	IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
	1493	zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg*(pplay(i,kb(i))- &
	1494	pplay(i,k))**nexpo
	1495	END IF
	1496	END DO
	1497	END DO
	1498	DO k = kbmin, khmax
[524]	1499	DO i = 1, klon
[1992]	1500	IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
	1501	d_ql(i, k) = d_ql(i, k) + toliq(i)zcond(i)/zvapo(i)(pplay(i,kb(i) &
	1502	)-pplay(i,k))**nexpo
	1503	END IF
	1504	END DO
	1505	END DO
	1506	DO i = 1, klon
[524]	1507	IF (todo(i)) THEN
[1992]	1508	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	1509	END IF
	1510	END DO
	1511
	1512	ELSE ! valeur non-prevue pour opt_cld
	1513
	1514	PRINT *, 'opt_cld est faux:', opt_cld
	1515	CALL abort
	1516
	1517	END IF ! fin de opt_cld
	1518
	1519	! L'eau precipitante peut etre re-evaporee:
	1520
	1521	IF (evap_prec .AND. kbmax>=2) THEN
	1522	DO k = kbmax, 1, -1
	1523	DO i = 1, klon
	1524	IF (todo(i) .AND. k<=(kb(i)-1) .AND. zrfl(i)>0.0) THEN
	1525	zqev = max(0.0, (zqs(i,k)-q(i,k))*zfrac(i))
	1526	zqevt = coef_eva(1.0-q(i,k)/zqs(i,k))sqrt(zrfl(i))* &
	1527	(paprs(i,k)-paprs(i,k+1))/pplay(i, k)t(i, k)rd/rg
	1528	zqevt = max(0.0, min(zqevt,zrfl(i)))rgdtime/ &
	1529	(paprs(i,k)-paprs(i,k+1))
	1530	zqev = min(zqev, zqevt)
	1531	zrfln = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime
	1532	d_q(i, k) = -(zrfln-zrfl(i))(rg/(paprs(i,k)-paprs(i,k+1)))dtime
	1533	d_t(i, k) = (zrfln-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
	1534	k+1)))dtimerlvtt/rcpd
	1535	zrfl(i) = zrfln
	1536	END IF
	1537	END DO
	1538	END DO
	1539	END IF
	1540
	1541	! La temperature de la premiere couche determine la pluie ou la neige:
	1542
	1543	DO i = 1, klon
	1544	IF (todo(i)) THEN
	1545	IF (t(i,1)>rtt) THEN
	1546	rain(i) = rain(i) + zrfl(i)
[524]	1547	ELSE
[1992]	1548	snow(i) = snow(i) + zrfl(i)
	1549	END IF
	1550	END IF
	1551	END DO
	1552
	1553	RETURN
	1554	END SUBROUTINE conkuo
	1555	SUBROUTINE kuofcl(pt, pq, pg, pp, ldcum, kcbot)
[5285]	1556	USE yomcst_mod_h
[1992]	1557	USE dimphy
[5284]	1558	USE yoethf_mod_h
[1992]	1559	IMPLICIT NONE
	1560	! ======================================================================
	1561	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
	1562	! adaptation du code de Tiedtke du ECMWF
	1563	! Objet: calculer le niveau de convection libre
	1564	! (FCL: Free Convection Level)
	1565	! ======================================================================
	1566	! Arguments:
	1567	! pt---input-R- temperature (K)
	1568	! pq---input-R- vapeur d'eau (kg/kg)
	1569	! pg---input-R- geopotentiel (g*z ou z est en metre)
	1570	! pp---input-R- pression (Pa)
	1571
	1572	! LDCUM---output-L- Y-t-il la convection
	1573	! kcbot---output-I- Niveau du bas de la convection
	1574	! ======================================================================
	1575
	1576	REAL pt(klon, klev), pq(klon, klev), pg(klon, klev), pp(klon, klev)
	1577	INTEGER kcbot(klon)
	1578	LOGICAL ldcum(klon)
	1579
	1580	REAL ztu(klon, klev), zqu(klon, klev), zlu(klon, klev)
	1581	REAL zqold(klon), zbuo
	1582	INTEGER is, i, k
	1583
	1584	! klab=1: on est sous le nuage convectif
	1585	! klab=2: le bas du nuage convectif
	1586	! klab=0: autres couches
	1587	INTEGER klab(klon, klev)
	1588
	1589	! quand lflag=.true., on est sous le nuage, il faut donc appliquer
	1590	! le processus d'elevation.
	1591	LOGICAL lflag(klon)
	1592
	1593	DO k = 1, klev
	1594	DO i = 1, klon
	1595	ztu(i, k) = pt(i, k)
	1596	zqu(i, k) = pq(i, k)
	1597	zlu(i, k) = 0.0
	1598	klab(i, k) = 0
	1599	END DO
	1600	END DO
	1601	! ----------------------------------------------------------------------
	1602	DO i = 1, klon
	1603	klab(i, 1) = 1
	1604	kcbot(i) = 2
	1605	ldcum(i) = .FALSE.
	1606	END DO
	1607
	1608	DO k = 2, klev - 1
	1609
	1610	is = 0
	1611	DO i = 1, klon
	1612	IF (klab(i,k-1)==1) is = is + 1
	1613	lflag(i) = .FALSE.
	1614	IF (klab(i,k-1)==1) lflag(i) = .TRUE.
	1615	END DO
	1616	IF (is==0) GO TO 290
	1617
	1618	! on eleve le parcel d'air selon l'adiabatique sec
	1619
	1620	DO i = 1, klon
[524]	1621	IF (lflag(i)) THEN
[1992]	1622	zqu(i, k) = zqu(i, k-1)
	1623	ztu(i, k) = ztu(i, k-1) + (pg(i,k-1)-pg(i,k))/rcpd
	1624	zbuo = ztu(i, k)(1.+retvzqu(i,k)) - pt(i, k)(1.+retvpq(i,k)) + &
	1625	0.5
	1626	IF (zbuo>0.) klab(i, k) = 1
	1627	zqold(i) = zqu(i, k)
[524]	1628	END IF
[1992]	1629	END DO
[524]	1630
[1992]	1631	! on calcule la condensation eventuelle
	1632
	1633	CALL adjtq(pp(1,k), ztu(1,k), zqu(1,k), lflag, 1)
	1634
	1635	! s'il y a la condensation et la "buoyancy" force est positive
	1636	! c'est bien le bas de la tour de convection
	1637
	1638	DO i = 1, klon
	1639	IF (lflag(i) .AND. zqu(i,k)/=zqold(i)) THEN
	1640	klab(i, k) = 2
	1641	zlu(i, k) = zlu(i, k) + zqold(i) - zqu(i, k)
	1642	zbuo = ztu(i, k)(1.+retvzqu(i,k)) - pt(i, k)(1.+retvpq(i,k)) + &
	1643	0.5
	1644	IF (zbuo>0.) THEN
	1645	kcbot(i) = k
	1646	ldcum(i) = .TRUE.
	1647	END IF
	1648	END IF
	1649	END DO
	1650
	1651	290 END DO
	1652
	1653	RETURN
	1654	END SUBROUTINE kuofcl
	1655	SUBROUTINE adjtq(pp, pt, pq, ldflag, kcall)
[5285]	1656	USE yomcst_mod_h
[1992]	1657	USE dimphy
[5284]	1658	USE yoethf_mod_h
[1992]	1659	IMPLICIT NONE
	1660	! ======================================================================
	1661	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
	1662	! adaptation du code de Tiedtke du ECMWF
	1663	! Objet: ajustement entre T et Q
	1664	! ======================================================================
	1665	! Arguments:
	1666	! pp---input-R- pression (Pa)
	1667	! pt---input/output-R- temperature (K)
	1668	! pq---input/output-R- vapeur d'eau (kg/kg)
	1669	! ======================================================================
	1670	! TO PRODUCE T,Q AND L VALUES FOR CLOUD ASCENT
	1671
	1672	! NOTE: INPUT PARAMETER KCALL DEFINES CALCULATION AS
	1673	! KCALL=0 ENV. T AND QS INCUINI
	1674	! KCALL=1 CONDENSATION IN UPDRAFTS (E.G. CUBASE, CUASC)
	1675	! KCALL=2 EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
	1676
	1677
	1678	REAL pt(klon), pq(klon), pp(klon)
	1679	LOGICAL ldflag(klon)
	1680	INTEGER kcall
	1681
	1682	REAL t_coup
	1683	PARAMETER (t_coup=234.0)
	1684
	1685	REAL zcond(klon), zcond1
	1686	REAL zdelta, zcvm5, zldcp, zqsat, zcor, zdqsat
	1687	INTEGER is, i
	1688	include "FCTTRE.h"
	1689
	1690	DO i = 1, klon
	1691	zcond(i) = 0.0
	1692	END DO
	1693
	1694	DO i = 1, klon
	1695	IF (ldflag(i)) THEN
	1696	zdelta = max(0., sign(1.,rtt-pt(i)))
	1697	zldcp = rlvtt(1.-zdelta) + zdeltarlstt
	1698	zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
	1699	IF (thermcep) THEN
	1700	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1701	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
	1702	zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
	1703	zqsat = min(0.5, zqsat)
	1704	zcor = 1./(1.-retv*zqsat)
	1705	zqsat = zqsat*zcor
	1706	zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
[524]	1707	ELSE
[1992]	1708	IF (pt(i)<t_coup) THEN
	1709	zqsat = qsats(pt(i))/pp(i)
	1710	zdqsat = dqsats(pt(i), zqsat)
	1711	ELSE
	1712	zqsat = qsatl(pt(i))/pp(i)
	1713	zdqsat = dqsatl(pt(i), zqsat)
	1714	END IF
	1715	END IF
	1716	zcond(i) = (pq(i)-zqsat)/(1.+zdqsat)
	1717	IF (kcall==1) zcond(i) = max(zcond(i), 0.)
	1718	IF (kcall==2) zcond(i) = min(zcond(i), 0.)
	1719	pt(i) = pt(i) + zldcp*zcond(i)
	1720	pq(i) = pq(i) - zcond(i)
	1721	END IF
	1722	END DO
	1723
	1724	is = 0
	1725	DO i = 1, klon
	1726	IF (zcond(i)/=0.) is = is + 1
	1727	END DO
	1728	IF (is==0) GO TO 230
	1729
	1730	DO i = 1, klon
	1731	IF (ldflag(i) .AND. zcond(i)/=0.) THEN
	1732	zdelta = max(0., sign(1.,rtt-pt(i)))
	1733	zldcp = rlvtt(1.-zdelta) + zdeltarlstt
	1734	zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
	1735	IF (thermcep) THEN
	1736	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1737	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
	1738	zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
	1739	zqsat = min(0.5, zqsat)
	1740	zcor = 1./(1.-retv*zqsat)
	1741	zqsat = zqsat*zcor
	1742	zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
[524]	1743	ELSE
[1992]	1744	IF (pt(i)<t_coup) THEN
	1745	zqsat = qsats(pt(i))/pp(i)
	1746	zdqsat = dqsats(pt(i), zqsat)
	1747	ELSE
	1748	zqsat = qsatl(pt(i))/pp(i)
	1749	zdqsat = dqsatl(pt(i), zqsat)
	1750	END IF
	1751	END IF
	1752	zcond1 = (pq(i)-zqsat)/(1.+zdqsat)
	1753	pt(i) = pt(i) + zldcp*zcond1
	1754	pq(i) = pq(i) - zcond1
	1755	END IF
	1756	END DO
	1757
	1758	230 CONTINUE
	1759	RETURN
	1760	END SUBROUTINE adjtq
	1761	SUBROUTINE fiajh(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, &
	1762	ibas, itop)
	1763	USE dimphy
[5285]	1764	USE yomcst_mod_h
[5284]	1765	USE yoethf_mod_h
[5274]	1766	IMPLICIT NONE
[1992]	1767
	1768	! Ajustement humide (Schema de convection de Manabe)
	1769	! .
	1770
[5274]	1771
[1992]	1772	! Arguments:
	1773
	1774	REAL dtime ! intervalle du temps (s)
	1775	REAL t(klon, klev) ! temperature (K)
	1776	REAL q(klon, klev) ! humidite specifique (kg/kg)
	1777	REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
	1778	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	1779
	1780	REAL d_t(klon, klev) ! incrementation pour la temperature
	1781	REAL d_q(klon, klev) ! incrementation pour vapeur d'eau
	1782	REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
	1783	REAL rneb(klon, klev) ! fraction nuageuse
	1784
	1785	REAL rain(klon) ! variable non utilisee
	1786	REAL snow(klon) ! variable non utilisee
	1787	INTEGER ibas(klon) ! variable non utilisee
	1788	INTEGER itop(klon) ! variable non utilisee
	1789
	1790	REAL t_coup
	1791	PARAMETER (t_coup=234.0)
	1792	REAL seuil_vap
	1793	PARAMETER (seuil_vap=1.0E-10)
	1794
	1795	! Variables locales:
	1796
	1797	INTEGER i, k
	1798	INTEGER k1, k1p, k2, k2p
	1799	LOGICAL itest(klon)
	1800	REAL delta_q(klon, klev)
	1801	REAL cp_new_t(klev)
	1802	REAL cp_delta_t(klev)
	1803	REAL new_qb(klev)
	1804	REAL v_cptj(klev), v_cptjk1, v_ssig
	1805	REAL v_cptt(klon, klev), v_p, v_t
	1806	REAL v_qs(klon, klev), v_qsd(klon, klev)
	1807	REAL zq1(klon), zq2(klon)
	1808	REAL gamcpdz(klon, 2:klev)
	1809	REAL zdp, zdpm
	1810
	1811	REAL zsat ! sur-saturation
	1812	REAL zflo ! flotabilite
	1813
	1814	REAL local_q(klon, klev), local_t(klon, klev)
	1815
	1816	REAL zdelta, zcor, zcvm5
	1817
	1818	include "FCTTRE.h"
	1819
	1820	DO k = 1, klev
	1821	DO i = 1, klon
	1822	local_q(i, k) = q(i, k)
	1823	local_t(i, k) = t(i, k)
	1824	rneb(i, k) = 0.0
	1825	d_ql(i, k) = 0.0
	1826	d_t(i, k) = 0.0
	1827	d_q(i, k) = 0.0
	1828	END DO
	1829	END DO
	1830	DO i = 1, klon
	1831	rain(i) = 0.0
	1832	snow(i) = 0.0
	1833	ibas(i) = 0
	1834	itop(i) = 0
	1835	END DO
	1836
	1837	! Calculer v_qs et v_qsd:
	1838
	1839	DO k = 1, klev
	1840	DO i = 1, klon
	1841	v_cptt(i, k) = rcpd*local_t(i, k)
	1842	v_t = local_t(i, k)
	1843	v_p = pplay(i, k)
	1844
	1845	IF (thermcep) THEN
	1846	zdelta = max(0., sign(1.,rtt-v_t))
	1847	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1848	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
	1849	v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
	1850	v_qs(i, k) = min(0.5, v_qs(i,k))
	1851	zcor = 1./(1.-retv*v_qs(i,k))
	1852	v_qs(i, k) = v_qs(i, k)*zcor
	1853	v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
	1854	ELSE
	1855	IF (v_t<t_coup) THEN
	1856	v_qs(i, k) = qsats(v_t)/v_p
	1857	v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
	1858	ELSE
	1859	v_qs(i, k) = qsatl(v_t)/v_p
	1860	v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
	1861	END IF
	1862	END IF
	1863	END DO
	1864	END DO
	1865
	1866	! Calculer Gamma * Cp * dz: (gamm est le gradient critique)
	1867
	1868	DO k = 2, klev
	1869	DO i = 1, klon
	1870	zdp = paprs(i, k) - paprs(i, k+1)
	1871	zdpm = paprs(i, k-1) - paprs(i, k)
	1872	gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)(v_cptt(i,k-1)zdpm+ &
	1873	v_cptt(i,k)zdp)+rlvtt/(zdpm+zdp)(v_qs(i,k-1)*zdpm+ &
	1874	v_qs(i,k)zdp))(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
	1875	k-1)zdpm+v_qsd(i,k)zdp)/(zdpm+zdp))
	1876	END DO
	1877	END DO
	1878
	1879	! ------------------------------------ modification des profils instables
	1880	DO i = 1, klon
	1881	itest(i) = .FALSE.
	1882
	1883	k1 = 0
	1884	k2 = 1
	1885
	1886	810 CONTINUE ! chercher k1, le bas de la colonne
	1887	k2 = k2 + 1
	1888	IF (k2>klev) GO TO 9999
	1889	zflo = v_cptt(i, k2-1) - v_cptt(i, k2) - gamcpdz(i, k2)
	1890	zsat = (local_q(i,k2-1)-v_qs(i,k2-1))*(paprs(i,k2-1)-paprs(i,k2)) + &
	1891	(local_q(i,k2)-v_qs(i,k2))*(paprs(i,k2)-paprs(i,k2+1))
	1892	IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 810
	1893	k1 = k2 - 1
	1894	itest(i) = .TRUE.
	1895
	1896	820 CONTINUE ! chercher k2, le haut de la colonne
	1897	IF (k2==klev) GO TO 821
	1898	k2p = k2 + 1
	1899	zsat = zsat + (paprs(i,k2p)-paprs(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
	1900	zflo = v_cptt(i, k2p-1) - v_cptt(i, k2p) - gamcpdz(i, k2p)
	1901	IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 821
	1902	k2 = k2p
	1903	GO TO 820
	1904	821 CONTINUE
	1905
	1906	! ------------------------------------------------------ ajustement local
	1907	830 CONTINUE ! ajustement proprement dit
	1908	v_cptj(k1) = 0.0
	1909	zdp = paprs(i, k1) - paprs(i, k1+1)
	1910	v_cptjk1 = ((1.0+v_qsd(i,k1))(v_cptt(i,k1)+v_cptj(k1))+rlvtt(local_q(i, &
	1911	k1)-v_qs(i,k1)))*zdp
	1912	v_ssig = zdp*(1.0+v_qsd(i,k1))
	1913
	1914	k1p = k1 + 1
	1915	DO k = k1p, k2
	1916	zdp = paprs(i, k) - paprs(i, k+1)
	1917	v_cptj(k) = v_cptj(k-1) + gamcpdz(i, k)
	1918	v_cptjk1 = v_cptjk1 + zdp((1.0+v_qsd(i,k))(v_cptt(i, &
	1919	k)+v_cptj(k))+rlvtt*(local_q(i,k)-v_qs(i,k)))
	1920	v_ssig = v_ssig + zdp*(1.0+v_qsd(i,k))
	1921	END DO
	1922
	1923	DO k = k1, k2
	1924	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
	1925	cp_delta_t(k) = cp_new_t(k) - v_cptt(i, k)
	1926	new_qb(k) = v_qs(i, k) + v_qsd(i, k)*cp_delta_t(k)/rlvtt
	1927	local_q(i, k) = new_qb(k)
	1928	local_t(i, k) = cp_new_t(k)/rcpd
	1929	END DO
	1930
	1931	! --------------------------------------------------- sondage vers le bas
	1932	! -- on redefinit les variables prognostiques dans
	1933	! -- la colonne qui vient d'etre ajustee
	1934
	1935	DO k = k1, k2
	1936	v_cptt(i, k) = rcpd*local_t(i, k)
	1937	v_t = local_t(i, k)
	1938	v_p = pplay(i, k)
	1939
	1940	IF (thermcep) THEN
	1941	zdelta = max(0., sign(1.,rtt-v_t))
	1942	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1943	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
	1944	v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
	1945	v_qs(i, k) = min(0.5, v_qs(i,k))
	1946	zcor = 1./(1.-retv*v_qs(i,k))
	1947	v_qs(i, k) = v_qs(i, k)*zcor
	1948	v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
	1949	ELSE
	1950	IF (v_t<t_coup) THEN
	1951	v_qs(i, k) = qsats(v_t)/v_p
	1952	v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
	1953	ELSE
	1954	v_qs(i, k) = qsatl(v_t)/v_p
	1955	v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
	1956	END IF
	1957	END IF
	1958	END DO
	1959	DO k = 2, klev
	1960	zdpm = paprs(i, k-1) - paprs(i, k)
	1961	zdp = paprs(i, k) - paprs(i, k+1)
	1962	gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)(v_cptt(i,k-1)zdpm+ &
	1963	v_cptt(i,k)zdp)+rlvtt/(zdpm+zdp)(v_qs(i,k-1)*zdpm+ &
	1964	v_qs(i,k)zdp))(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
	1965	k-1)zdpm+v_qsd(i,k)zdp)/(zdpm+zdp))
	1966	END DO
	1967
	1968	! Verifier si l'on peut etendre la colonne vers le bas
	1969
	1970	IF (k1==1) GO TO 841 ! extension echouee
	1971	zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
	1972	zsat = (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1)) + &
	1973	(local_q(i,k1)-v_qs(i,k1))*(paprs(i,k1)-paprs(i,k1+1))
	1974	IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 841 ! extension echouee
	1975
	1976	840 CONTINUE
	1977	k1 = k1 - 1
	1978	IF (k1==1) GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	1979	zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1))
	1980	zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
	1981	IF (zflo>0.0 .AND. zsat>0.0) THEN
	1982	GO TO 840
	1983	ELSE
	1984	GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	1985	END IF
	1986	841 CONTINUE
	1987
	1988	GO TO 810 ! chercher d'autres blocks en haut
	1989
	1990	9999 END DO ! boucle sur tous les points
	1991	! -----------------------------------------------------------------------
	1992
	1993	! Determiner la fraction nuageuse (hypothese: la nebulosite a lieu
	1994	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
	1995
	1996	DO k = 1, klev
	1997	DO i = 1, klon
	1998	IF (itest(i)) THEN
	1999	delta_q(i, k) = local_q(i, k) - q(i, k)
	2000	IF (delta_q(i,k)<0.) rneb(i, k) = 1.0
	2001	END IF
	2002	END DO
	2003	END DO
	2004
	2005	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
	2006	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
	2007	! diminue et d'une maniere proportionnelle a cet diminution):
	2008
	2009	DO i = 1, klon
	2010	IF (itest(i)) THEN
	2011	zq1(i) = 0.0
	2012	zq2(i) = 0.0
	2013	END IF
	2014	END DO
	2015	DO k = 1, klev
	2016	DO i = 1, klon
	2017	IF (itest(i)) THEN
	2018	zdp = paprs(i, k) - paprs(i, k+1)
	2019	zq1(i) = zq1(i) - delta_q(i, k)*zdp
	2020	zq2(i) = zq2(i) - min(0.0, delta_q(i,k))*zdp
	2021	END IF
	2022	END DO
	2023	END DO
	2024	DO k = 1, klev
	2025	DO i = 1, klon
	2026	IF (itest(i)) THEN
	2027	IF (zq2(i)/=0.0) d_ql(i, k) = -min(0.0, delta_q(i,k))*zq1(i)/zq2(i)
	2028	END IF
	2029	END DO
	2030	END DO
	2031
	2032	DO k = 1, klev
	2033	DO i = 1, klon
	2034	local_q(i, k) = max(local_q(i,k), seuil_vap)
	2035	END DO
	2036	END DO
	2037
	2038	DO k = 1, klev
	2039	DO i = 1, klon
	2040	d_t(i, k) = local_t(i, k) - t(i, k)
	2041	d_q(i, k) = local_q(i, k) - q(i, k)
	2042	END DO
	2043	END DO
	2044
	2045	RETURN
	2046	END SUBROUTINE fiajh
	2047	SUBROUTINE fiajc(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
	2048	rain, snow, ibas, itop)
	2049	USE dimphy
[5285]	2050	USE yomcst_mod_h
[5284]	2051	USE yoethf_mod_h
[5274]	2052	IMPLICIT NONE
[1992]	2053
	2054
[5274]	2055
[1992]	2056	! Options:
	2057
	2058	INTEGER plb ! niveau de depart pour la convection
	2059	PARAMETER (plb=4)
	2060
	2061	! Mystere: cette option n'est pas innocente pour les resultats !
	2062	! Qui peut resoudre ce mystere ? (Z.X.Li mars 1995)
	2063	LOGICAL vector ! calcul vectorise
	2064	PARAMETER (vector=.FALSE.)
	2065
	2066	REAL t_coup
	2067	PARAMETER (t_coup=234.0)
	2068
	2069	! Arguments:
	2070
	2071	REAL q(klon, klev) ! humidite specifique (kg/kg)
	2072	REAL t(klon, klev) ! temperature (K)
	2073	REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
	2074	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	2075	REAL dtime ! intervalle du temps (s)
	2076	REAL conv_q(klon, klev) ! taux de convergence de l'humidite
	2077	REAL rneb(klon, klev) ! fraction nuageuse
	2078	REAL d_q(klon, klev) ! incrementaion pour la vapeur d'eau
	2079	REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
	2080	REAL d_t(klon, klev) ! incrementation pour la temperature
	2081	REAL rain(klon) ! variable non-utilisee
	2082	REAL snow(klon) ! variable non-utilisee
	2083	INTEGER itop(klon) ! variable non-utilisee
	2084	INTEGER ibas(klon) ! variable non-utilisee
	2085
	2086	INTEGER kh(klon), i, k
	2087	LOGICAL nuage(klon), test(klon, klev)
	2088	REAL zconv(klon), zdeh(klon, klev), zvirt(klon)
	2089	REAL zdqs(klon, klev), zqs(klon, klev)
	2090	REAL ztt, zvar, zfrac(klon)
	2091	REAL zq1(klon), zq2(klon)
	2092	REAL zdelta, zcor, zcvm5
	2093
	2094	include "FCTTRE.h"
	2095
	2096	! Initialiser les sorties:
	2097
	2098	DO k = 1, klev
	2099	DO i = 1, klon
	2100	rneb(i, k) = 0.0
	2101	d_ql(i, k) = 0.0
	2102	d_t(i, k) = 0.0
	2103	d_q(i, k) = 0.0
	2104	END DO
	2105	END DO
	2106	DO i = 1, klon
	2107	itop(i) = 0
	2108	ibas(i) = 0
	2109	rain(i) = 0.0
	2110	snow(i) = 0.0
	2111	END DO
	2112
	2113	! Calculer Qs et L/Cp * dQs/dT:
	2114
	2115	DO k = 1, klev
	2116	DO i = 1, klon
	2117	ztt = t(i, k)
	2118	IF (thermcep) THEN
	2119	zdelta = max(0., sign(1.,rtt-ztt))
	2120	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	2121	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
	2122	zqs(i, k) = r2es*foeew(ztt, zdelta)/pplay(i, k)
	2123	zqs(i, k) = min(0.5, zqs(i,k))
	2124	zcor = 1./(1.-retv*zqs(i,k))
	2125	zqs(i, k) = zqs(i, k)*zcor
	2126	zdqs(i, k) = foede(ztt, zdelta, zcvm5, zqs(i,k), zcor)
	2127	ELSE
	2128	IF (ztt<t_coup) THEN
	2129	zqs(i, k) = qsats(ztt)/pplay(i, k)
	2130	zdqs(i, k) = dqsats(ztt, zqs(i,k))
	2131	ELSE
	2132	zqs(i, k) = qsatl(ztt)/pplay(i, k)
	2133	zdqs(i, k) = dqsatl(ztt, zqs(i,k))
	2134	END IF
	2135	END IF
	2136	END DO
	2137	END DO
	2138
	2139	! Determiner la difference de l'energie totale saturee:
	2140
	2141	DO i = 1, klon
	2142	k = plb
	2143	zdeh(i, k) = rcpd(t(i,k-1)-t(i,k)) - rd0.5*(t(i,k-1)+t(i,k))/paprs(i, k &
	2144	)(pplay(i,k-1)-pplay(i,k)) + rlvtt(zqs(i,k-1)-zqs(i,k))
	2145	zdeh(i, k) = zdeh(i, k)*0.5 ! on prend la moitie
	2146	END DO
	2147	DO k = plb + 1, klev
	2148	DO i = 1, klon
	2149	zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
	2150	rd0.5(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
	2151	rlvtt*(zqs(i,k-1)-zqs(i,k))
	2152	END DO
	2153	END DO
	2154
	2155	! Determiner le sommet du nuage selon l'instabilite
	2156	! Calculer les convergences d'humidite (reelle et virtuelle)
	2157
	2158	DO i = 1, klon
	2159	nuage(i) = .TRUE.
	2160	zconv(i) = 0.0
	2161	zvirt(i) = 0.0
	2162	kh(i) = -999
	2163	END DO
	2164	DO k = plb, klev
	2165	DO i = 1, klon
	2166	nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
	2167	IF (nuage(i)) THEN
	2168	kh(i) = k
	2169	zconv(i) = zconv(i) + conv_q(i, k)dtime(paprs(i,k)-paprs(i,k+1))
	2170	zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k)- &
	2171	paprs(i,k+1))
	2172	END IF
	2173	END DO
	2174	END DO
	2175
	2176	IF (vector) THEN
	2177
	2178
	2179	DO k = plb, klev
[524]	2180	DO i = 1, klon
[1992]	2181	IF (k<=kh(i) .AND. kh(i)>plb .AND. zconv(i)>0.0) THEN
	2182	test(i, k) = .TRUE.
	2183	zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
	2184	ELSE
	2185	test(i, k) = .FALSE.
	2186	END IF
	2187	END DO
	2188	END DO
	2189
	2190	DO k = plb, klev
[524]	2191	DO i = 1, klon
[1992]	2192	IF (test(i,k)) THEN
	2193	zvar = zdeh(i, k)/(1.0+zdqs(i,k))
	2194	d_q(i, k) = (zvarzdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))zfrac(i) - &
	2195	conv_q(i, k)*dtime
	2196	d_t(i, k) = zvar*zfrac(i)/rcpd
	2197	END IF
	2198	END DO
	2199	END DO
	2200
	2201	DO i = 1, klon
	2202	zq1(i) = 0.0
	2203	zq2(i) = 0.0
	2204	END DO
	2205	DO k = plb, klev
[524]	2206	DO i = 1, klon
[1992]	2207	IF (test(i,k)) THEN
	2208	IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
	2209	zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
	2210	zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
	2211	END IF
	2212	END DO
	2213	END DO
	2214
	2215	DO k = plb, klev
[524]	2216	DO i = 1, klon
[1992]	2217	IF (test(i,k)) THEN
	2218	IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
	2219	END IF
	2220	END DO
	2221	END DO
	2222
	2223	ELSE ! (.NOT. vector)
	2224
	2225	DO i = 1, klon
	2226	IF (kh(i)>plb .AND. zconv(i)>0.0) THEN
	2227	! cc IF (kh(i).LE.plb) GOTO 999 ! il n'y a pas d'instabilite
	2228	! cc IF (zconv(i).LE.0.0) GOTO 999 ! convergence insuffisante
	2229	zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
	2230	DO k = plb, kh(i)
	2231	zvar = zdeh(i, k)/(1.0+zdqs(i,k))
	2232	d_q(i, k) = (zvarzdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))zfrac(i) - &
	2233	conv_q(i, k)*dtime
	2234	d_t(i, k) = zvar*zfrac(i)/rcpd
	2235	END DO
	2236
	2237	zq1(i) = 0.0
	2238	zq2(i) = 0.0
	2239	DO k = plb, kh(i)
	2240	IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
	2241	zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
	2242	zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
	2243	END DO
	2244	DO k = plb, kh(i)
	2245	IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
	2246	END DO
	2247	END IF
	2248	END DO
	2249
	2250	END IF ! fin de teste sur vector
	2251
	2252	RETURN
	2253	END SUBROUTINE fiajc

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

Download in other formats:

Original Format