Context Navigation

cv30_routines_mod.f90 @ 5428

Last change on this file since 5428 was 5346, checked in by fhourdin, 3 weeks ago

Debut de replaysation de la convection profonde.

Regroupement de cvparam, cv3param et cvthermo (récemment
passés de statut de .h à module, dans un unique module
lmdz_cv_ini.f90

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: 108.8 KB

Rev	Line
[5283]	1	MODULE cv30_routines_mod
	2	!------------------------------------------------------------
	3	! Parameters for convectL, iflag_con=30:
	4	! (includes - microphysical parameters,
	5	! - parameters that control the rate of approach
	6	! to quasi-equilibrium)
	7	! - noff & minorig (previously in input of convect1)
	8	!------------------------------------------------------------
[1992]	9
[5283]	10	IMPLICIT NONE; PRIVATE
	11	PUBLIC sigd, spfac, pbcrit, ptcrit, omtrain, dtovsh, dpbase, dttrig, dtcrit, &
	12	tau, beta, alpha, delta, betad, noff, minorig, nl, nlp, nlm, &
	13	cv30_param, cv30_prelim, cv30_feed, cv30_undilute1, cv30_trigger, &
	14	cv30_compress, cv30_undilute2, cv30_closure, cv30_mixing, cv30_unsat, &
	15	cv30_yield, cv30_tracer, cv30_uncompress, cv30_epmax_fn_cape
[879]	16
[5283]	17	INTEGER noff, minorig, nl, nlp, nlm
	18	REAL sigd, spfac
	19	REAL pbcrit, ptcrit
	20	REAL omtrain
	21	REAL dtovsh, dpbase, dttrig
	22	REAL dtcrit, tau, beta, alpha
	23	REAL delta
	24	REAL betad
[879]	25
[5283]	26	!$OMP THREADPRIVATE(sigd, spfac, pbcrit, ptcrit, omtrain, dtovsh, dpbase, dttrig, dtcrit, &
	27	!$OMP tau, beta, alpha, delta, betad, noff, minorig, nl, nlp, nlm)
	28	CONTAINS
[879]	29
[5283]	30	SUBROUTINE cv30_param(nd, delt)
	31	USE conema3_mod_h
[879]	32
[5283]	33	IMPLICIT NONE
[879]	34
[5283]	35	! ------------------------------------------------------------
	36	! Set parameters for convectL for iflag_con = 3
	37	! ------------------------------------------------------------
[879]	38
	39
[5283]	40	! * PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *
	41	! * PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *
	42	! * PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *
	43	! * EFFICIENCY IS ASSUMED TO BE UNITY *
	44	! * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
	45	! * SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
	46	! * OF CLOUD *
[879]	47
[5283]	48	! [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
	49	! * ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *
	50	! * APPROACH TO QUASI-EQUILIBRIUM *
	51	! * (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *
	52	! * (BETA MUST BE LESS THAN OR EQUAL TO 1) *
[879]	53
[5283]	54	! * DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *
	55	! * APPROACH TO QUASI-EQUILIBRIUM *
	56	! * IT MUST BE LESS THAN 0 *
[879]	57
[5283]	58	INTEGER nd
	59	REAL delt ! timestep (seconds)
[879]	60
[5283]	61	! noff: integer limit for convection (nd-noff)
	62	! minorig: First level of convection
[879]	63
[5283]	64	! -- limit levels for convection:
[879]	65
[5283]	66	noff = 1
	67	minorig = 1
	68	nl = nd - noff
	69	nlp = nl + 1
	70	nlm = nl - 1
[879]	71
[5283]	72	! -- "microphysical" parameters:
[879]	73
[5283]	74	sigd = 0.01
	75	spfac = 0.15
	76	pbcrit = 150.0
	77	ptcrit = 500.0
	78	! IM cf. FH epmax = 0.993
[879]	79
[5283]	80	omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
[879]	81
[5283]	82	! -- misc:
[879]	83
[5283]	84	dtovsh = -0.2 ! dT for overshoot
	85	dpbase = -40. ! definition cloud base (400m above LCL)
	86	dttrig = 5. ! (loose) condition for triggering
[879]	87
[5283]	88	! -- rate of approach to quasi-equilibrium:
[879]	89
[5283]	90	dtcrit = -2.0
	91	tau = 8000.
	92	beta = 1.0 - delt / tau
	93	alpha = 1.5E-3 * delt / tau
	94	! increase alpha to compensate W decrease:
	95	alpha = alpha * 1.5
[879]	96
[5283]	97	! -- interface cloud parameterization:
[879]	98
[5283]	99	delta = 0.01 ! cld
[879]	100
[5283]	101	! -- interface with boundary-layer (gust factor): (sb)
[879]	102
[5283]	103	betad = 10.0 ! original value (from convect 4.3)
[879]	104
[5283]	105	END SUBROUTINE cv30_param
[879]	106
[5283]	107	SUBROUTINE cv30_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm, &
	108	th)
[5346]	109	USE lmdz_cv_ini, ONLY : rrd,rrv,lv0,eps,cpv,cpd,clmcpv,cl,rrv
[879]	110
[5283]	111	IMPLICIT NONE
[879]	112
[5283]	113	! =====================================================================
	114	! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
	115	! "ori": from convect4.3 (vectorized)
	116	! "convect3": to be exactly consistent with convect3
	117	! =====================================================================
[879]	118
[5283]	119	! inputs:
	120	INTEGER len, nd, ndp1
	121	REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
[879]	122
[5283]	123	! outputs:
	124	REAL lv(len, nd), cpn(len, nd), tv(len, nd)
	125	REAL gz(len, nd), h(len, nd), hm(len, nd)
	126	REAL th(len, nd)
[879]	127
[5283]	128	! local variables:
	129	INTEGER k, i
	130	REAL rdcp
	131	REAL tvx, tvy ! convect3
	132	REAL cpx(len, nd)
[879]	133
[5283]	134	! ori do 110 k=1,nlp
	135	DO k = 1, nl ! convect3
	136	DO i = 1, len
	137	! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
	138	lv(i, k) = lv0 - clmcpv * (t(i, k) - 273.15)
	139	cpn(i, k) = cpd * (1.0 - q(i, k)) + cpv * q(i, k)
	140	cpx(i, k) = cpd * (1.0 - q(i, k)) + cl * q(i, k)
	141	! ori tv(i,k)=t(i,k)(1.0+q(i,k)epsim1)
	142	tv(i, k) = t(i, k) * (1.0 + q(i, k) / eps - q(i, k))
	143	rdcp = (rrd * (1. - q(i, k)) + q(i, k) * rrv) / cpn(i, k)
	144	th(i, k) = t(i, k) * (1000.0 / p(i, k))**rdcp
	145	END DO
[1992]	146	END DO
[879]	147
[5283]	148	! gz = phi at the full levels (same as p).
[879]	149
[1992]	150	DO i = 1, len
[5283]	151	gz(i, 1) = 0.0
	152	END DO
	153	! ori do 140 k=2,nlp
	154	DO k = 2, nl ! convect3
	155	DO i = 1, len
	156	tvx = t(i, k) * (1. + q(i, k) / eps - q(i, k)) !convect3
	157	tvy = t(i, k - 1) * (1. + q(i, k - 1) / eps - q(i, k - 1)) !convect3
	158	gz(i, k) = gz(i, k - 1) + 0.5 * rrd * (tvx + tvy) & !convect3
	159	* (p(i, k - 1) - p(i, k)) / ph(i, k) !convect3
[879]	160
[5283]	161	! ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
	162	! ori & *(p(i,k-1)-p(i,k))/ph(i,k)
	163	END DO
[1992]	164	END DO
[879]	165
[5283]	166	! h = phi + cpT (dry static energy).
	167	! hm = phi + cp(T-Tbase)+Lq
[879]	168
[5283]	169	! ori do 170 k=1,nlp
	170	DO k = 1, nl ! convect3
	171	DO i = 1, len
	172	h(i, k) = gz(i, k) + cpn(i, k) * t(i, k)
	173	hm(i, k) = gz(i, k) + cpx(i, k) * (t(i, k) - t(i, 1)) + lv(i, k) * q(i, k)
	174	END DO
[1992]	175	END DO
[879]	176
[5283]	177	END SUBROUTINE cv30_prelim
[879]	178
[5283]	179	SUBROUTINE cv30_feed(len, nd, t, q, qs, p, ph, hm, gz, nk, icb, icbmax, &
	180	iflag, tnk, qnk, gznk, plcl)
[879]	181
[5283]	182	IMPLICIT NONE
[879]	183
[5283]	184	! ================================================================
	185	! Purpose: CONVECTIVE FEED
[879]	186
[5283]	187	! Main differences with cv_feed:
	188	! - ph added in input
	189	! - here, nk(i)=minorig
	190	! - icb defined differently (plcl compared with ph instead of p)
[879]	191
[5283]	192	! Main differences with convect3:
	193	! - we do not compute dplcldt and dplcldr of CLIFT anymore
	194	! - values iflag different (but tests identical)
	195	! - A,B explicitely defined (!...)
	196	! ================================================================
[879]	197
[5283]	198	! inputs:
	199	INTEGER len, nd
	200	REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
	201	REAL hm(len, nd), gz(len, nd)
	202	REAL ph(len, nd + 1)
[879]	203
[5283]	204	! outputs:
	205	INTEGER iflag(len), nk(len), icb(len), icbmax
	206	REAL tnk(len), qnk(len), gznk(len), plcl(len)
[879]	207
[5283]	208	! local variables:
	209	INTEGER i, k
	210	INTEGER ihmin(len)
	211	REAL work(len)
	212	REAL pnk(len), qsnk(len), rh(len), chi(len)
	213	REAL a, b ! convect3
	214	! ym
	215	plcl = 0.0
	216	! @ !-------------------------------------------------------------------
	217	! @ ! --- Find level of minimum moist static energy
	218	! @ ! --- If level of minimum moist static energy coincides with
	219	! @ ! --- or is lower than minimum allowable parcel origin level,
	220	! @ ! --- set iflag to 6.
	221	! @ !-------------------------------------------------------------------
	222	! @
	223	! @ do 180 i=1,len
	224	! @ work(i)=1.0e12
	225	! @ ihmin(i)=nl
	226	! @ 180 continue
	227	! @ do 200 k=2,nlp
	228	! @ do 190 i=1,len
	229	! @ if((hm(i,k).lt.work(i)).AND.
	230	! @ & (hm(i,k).lt.hm(i,k-1)))THEN
	231	! @ work(i)=hm(i,k)
	232	! @ ihmin(i)=k
	233	! @ endif
	234	! @ 190 continue
	235	! @ 200 continue
	236	! @ do 210 i=1,len
	237	! @ ihmin(i)=min(ihmin(i),nlm)
	238	! @ IF(ihmin(i).le.minorig)THEN
	239	! @ iflag(i)=6
	240	! @ endif
	241	! @ 210 continue
	242	! @ c
	243	! @ !-------------------------------------------------------------------
	244	! @ ! --- Find that model level below the level of minimum moist static
	245	! @ ! --- energy that has the maximum value of moist static energy
	246	! @ !-------------------------------------------------------------------
	247	! @
	248	! @ do 220 i=1,len
	249	! @ work(i)=hm(i,minorig)
	250	! @ nk(i)=minorig
	251	! @ 220 continue
	252	! @ do 240 k=minorig+1,nl
	253	! @ do 230 i=1,len
	254	! @ if((hm(i,k).gt.work(i)).AND.(k.le.ihmin(i)))THEN
	255	! @ work(i)=hm(i,k)
	256	! @ nk(i)=k
	257	! @ endif
	258	! @ 230 continue
	259	! @ 240 continue
[879]	260
[5283]	261	! -------------------------------------------------------------------
	262	! --- Origin level of ascending parcels for convect3:
	263	! -------------------------------------------------------------------
[879]	264
[5283]	265	DO i = 1, len
	266	nk(i) = minorig
	267	END DO
[879]	268
[5283]	269	! -------------------------------------------------------------------
	270	! --- Check whether parcel level temperature and specific humidity
	271	! --- are reasonable
	272	! -------------------------------------------------------------------
	273	DO i = 1, len
	274	IF (((t(i, nk(i))<250.0) .OR. (q(i, nk(i))<=0.0)) & ! @ & .OR.(
	275	! p(i,ihmin(i)).lt.400.0
	276	! ) )
	277	.AND. (iflag(i)==0)) iflag(i) = 7
	278	END DO
	279	! -------------------------------------------------------------------
	280	! --- Calculate lifted condensation level of air at parcel origin level
	281	! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
	282	! -------------------------------------------------------------------
[879]	283
[5283]	284	a = 1669.0 ! convect3
	285	b = 122.0 ! convect3
[879]	286
[5283]	287	DO i = 1, len
[879]	288
[5283]	289	IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002
[879]	290
[5283]	291	tnk(i) = t(i, nk(i))
	292	qnk(i) = q(i, nk(i))
	293	gznk(i) = gz(i, nk(i))
	294	pnk(i) = p(i, nk(i))
	295	qsnk(i) = qs(i, nk(i))
[879]	296
[5283]	297	rh(i) = qnk(i) / qsnk(i)
	298	! ori rh(i)=min(1.0,rh(i)) ! removed for convect3
	299	! ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
	300	chi(i) = tnk(i) / (a - b * rh(i) - tnk(i)) ! convect3
	301	plcl(i) = pnk(i) * (rh(i)**chi(i))
	302	IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag &
	303	(i) = 8
[879]	304
[5283]	305	END IF ! iflag=7
[879]	306
[5283]	307	END DO
[879]	308
[5283]	309	! -------------------------------------------------------------------
	310	! --- Calculate first level above lcl (=icb)
	311	! -------------------------------------------------------------------
[879]	312
[5283]	313	! @ do 270 i=1,len
	314	! @ icb(i)=nlm
	315	! @ 270 continue
	316	! @c
	317	! @ do 290 k=minorig,nl
	318	! @ do 280 i=1,len
	319	! @ if((k.ge.(nk(i)+1)).AND.(p(i,k).lt.plcl(i)))
	320	! @ & icb(i)=min(icb(i),k)
	321	! @ 280 continue
	322	! @ 290 continue
	323	! @c
	324	! @ do 300 i=1,len
	325	! @ if((icb(i).ge.nlm).AND.(iflag(i).EQ.0))iflag(i)=9
	326	! @ 300 continue
[879]	327
[5283]	328	DO i = 1, len
	329	icb(i) = nlm
	330	END DO
[879]	331
[5283]	332	! la modification consiste a comparer plcl a ph et non a p:
	333	! icb est defini par : ph(icb)<plcl<ph(icb-1)
	334	! @ do 290 k=minorig,nl
	335	DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
	336	DO i = 1, len
	337	IF (ph(i, k)<plcl(i)) icb(i) = min(icb(i), k)
	338	END DO
	339	END DO
	340
[1992]	341	DO i = 1, len
[5283]	342	! @ if((icb(i).ge.nlm).AND.(iflag(i).EQ.0))iflag(i)=9
	343	IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
[1992]	344	END DO
[879]	345
[5283]	346	DO i = 1, len
	347	icb(i) = icb(i) - 1 ! icb sup ou egal a 2
	348	END DO
[879]	349
[5283]	350	! Compute icbmax.
[879]	351
[5283]	352	icbmax = 2
	353	DO i = 1, len
	354	! icbmax=max(icbmax,icb(i))
	355	IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
	356	END DO
[879]	357
[5283]	358	END SUBROUTINE cv30_feed
[879]	359
[5283]	360	SUBROUTINE cv30_undilute1(len, nd, t, q, qs, gz, plcl, p, nk, icb, tp, tvp, &
	361	clw, icbs)
[5346]	362	USE lmdz_cv_ini, ONLY : cl,clmcpv,cpd,cpv,eps,lv0,rrv
[879]	363
[5283]	364	IMPLICIT NONE
[879]	365
[5283]	366	! ----------------------------------------------------------------
	367	! Equivalent de TLIFT entre NK et ICB+1 inclus
[879]	368
[5283]	369	! Differences with convect4:
	370	! - specify plcl in input
	371	! - icbs is the first level above LCL (may differ from icb)
	372	! - in the iterations, used x(icbs) instead x(icb)
	373	! - many minor differences in the iterations
	374	! - tvp is computed in only one time
	375	! - icbs: first level above Plcl (IMIN de TLIFT) in output
	376	! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
	377	! ----------------------------------------------------------------
[879]	378
	379
	380
[5283]	381	! inputs:
	382	INTEGER len, nd
	383	INTEGER nk(len), icb(len)
	384	REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
	385	REAL p(len, nd)
	386	REAL plcl(len) ! convect3
[879]	387
[5283]	388	! outputs:
	389	REAL tp(len, nd), tvp(len, nd), clw(len, nd)
[879]	390
[5283]	391	! local variables:
	392	INTEGER i, k
	393	INTEGER icb1(len), icbs(len), icbsmax2 ! convect3
	394	REAL tg, qg, alv, s, ahg, tc, denom, es, rg
	395	REAL ah0(len), cpp(len)
	396	REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
	397	REAL qsicb(len) ! convect3
	398	REAL cpinv(len) ! convect3
[879]	399
[5283]	400	! -------------------------------------------------------------------
	401	! --- Calculates the lifted parcel virtual temperature at nk,
	402	! --- the actual temperature, and the adiabatic
	403	! --- liquid water content. The procedure is to solve the equation.
	404	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
	405	! -------------------------------------------------------------------
[879]	406
[5283]	407	DO i = 1, len
	408	tnk(i) = t(i, nk(i))
	409	qnk(i) = q(i, nk(i))
	410	gznk(i) = gz(i, nk(i))
	411	! ori ticb(i)=t(i,icb(i))
	412	! ori gzicb(i)=gz(i,icb(i))
	413	END DO
[879]	414
[5283]	415	! * Calculate certain parcel quantities, including static energy *
[879]	416
[5283]	417	DO i = 1, len
	418	ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
	419	273.15)) + gznk(i)
	420	cpp(i) = cpd * (1. - qnk(i)) + qnk(i) * cpv
	421	cpinv(i) = 1. / cpp(i)
	422	END DO
[879]	423
[5283]	424	! * Calculate lifted parcel quantities below cloud base *
[879]	425
[5283]	426	DO i = 1, len !convect3
	427	icb1(i) = min(max(icb(i), 2), nl)
	428	! if icb is below LCL, start loop at ICB+1:
	429	! (icbs est le premier niveau au-dessus du LCL)
	430	icbs(i) = icb1(i) !convect3
	431	IF (plcl(i)<p(i, icb1(i))) THEN
	432	icbs(i) = min(icbs(i) + 1, nl) !convect3
	433	END IF
	434	END DO !convect3
[879]	435
[5283]	436	DO i = 1, len !convect3
	437	ticb(i) = t(i, icbs(i)) !convect3
	438	gzicb(i) = gz(i, icbs(i)) !convect3
	439	qsicb(i) = qs(i, icbs(i)) !convect3
	440	END DO !convect3
[879]	441
	442
[5283]	443	! Re-compute icbsmax (icbsmax2): !convect3
	444	!convect3
	445	icbsmax2 = 2 !convect3
	446	DO i = 1, len !convect3
	447	icbsmax2 = max(icbsmax2, icbs(i)) !convect3
	448	END DO !convect3
	449
	450	! initialization outputs:
	451
	452	DO k = 1, icbsmax2 ! convect3
	453	DO i = 1, len ! convect3
	454	tp(i, k) = 0.0 ! convect3
	455	tvp(i, k) = 0.0 ! convect3
	456	clw(i, k) = 0.0 ! convect3
	457	END DO ! convect3
[1992]	458	END DO ! convect3
[879]	459
[5283]	460	! tp and tvp below cloud base:
[879]	461
[5283]	462	DO k = minorig, icbsmax2 - 1
	463	DO i = 1, len
	464	tp(i, k) = tnk(i) - (gz(i, k) - gznk(i)) * cpinv(i)
	465	tvp(i, k) = tp(i, k) * (1. + qnk(i) / eps - qnk(i)) !whole thing (convect3)
	466	END DO
[1992]	467	END DO
[879]	468
[5283]	469	! * Find lifted parcel quantities above cloud base *
[879]	470
[5283]	471	DO i = 1, len
	472	tg = ticb(i)
	473	! ori qg=qs(i,icb(i))
	474	qg = qsicb(i) ! convect3
	475	! debug alv=lv0-clmcpv*(ticb(i)-t0)
	476	alv = lv0 - clmcpv * (ticb(i) - 273.15)
[879]	477
[5283]	478	! First iteration.
[879]	479
[5283]	480	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	481	s = cpd * (1. - qnk(i)) + cl * qnk(i) & ! convect3
	482	+ alv * alv * qg / (rrv * ticb(i) * ticb(i)) ! convect3
	483	s = 1. / s
	484	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	485	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gzicb(i) ! convect3
	486	tg = tg + s * (ah0(i) - ahg)
	487	! ori tg=max(tg,35.0)
	488	! debug tc=tg-t0
	489	tc = tg - 273.15
	490	denom = 243.5 + tc
	491	denom = max(denom, 1.0) ! convect3
	492	! ori IF(tc.ge.0.0)THEN
	493	es = 6.112 * exp(17.67 * tc / denom)
	494	! ori else
	495	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	496	! ori endif
	497	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	498	qg = eps * es / (p(i, icbs(i)) - es * (1. - eps))
[879]	499
[5283]	500	! Second iteration.
[879]	501
	502
[5283]	503	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	504	! ori s=1./s
	505	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	506	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gzicb(i) ! convect3
	507	tg = tg + s * (ah0(i) - ahg)
	508	! ori tg=max(tg,35.0)
	509	! debug tc=tg-t0
	510	tc = tg - 273.15
	511	denom = 243.5 + tc
	512	denom = max(denom, 1.0) ! convect3
	513	! ori IF(tc.ge.0.0)THEN
	514	es = 6.112 * exp(17.67 * tc / denom)
	515	! ori else
	516	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	517	! ori end if
	518	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	519	qg = eps * es / (p(i, icbs(i)) - es * (1. - eps))
[879]	520
[5283]	521	alv = lv0 - clmcpv * (ticb(i) - 273.15)
[879]	522
[5283]	523	! ori c approximation here:
	524	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	525	! ori & -gz(i,icb(i))-alv*qg)/cpd
[879]	526
[5283]	527	! convect3: no approximation:
	528	tp(i, icbs(i)) = (ah0(i) - gz(i, icbs(i)) - alv * qg) / (cpd + (cl - cpd) * qnk(i))
[879]	529
[5283]	530	! ori clw(i,icb(i))=qnk(i)-qg
	531	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	532	clw(i, icbs(i)) = qnk(i) - qg
	533	clw(i, icbs(i)) = max(0.0, clw(i, icbs(i)))
[879]	534
[5283]	535	rg = qg / (1. - qnk(i))
	536	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	537	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
	538	tvp(i, icbs(i)) = tp(i, icbs(i)) * (1. + qg / eps - qnk(i)) !whole thing
[879]	539
[5283]	540	END DO
[879]	541
[5283]	542	! ori do 380 k=minorig,icbsmax2
	543	! ori do 370 i=1,len
	544	! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
	545	! ori 370 continue
	546	! ori 380 continue
[879]	547
	548
[5283]	549	! -- The following is only for convect3:
[879]	550
[5283]	551	! * icbs is the first level above the LCL:
	552	! if plcl<p(icb), then icbs=icb+1
	553	! if plcl>p(icb), then icbs=icb
[1403]	554
[5283]	555	! * the routine above computes tvp from minorig to icbs (included).
[879]	556
[5283]	557	! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
	558	! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
[879]	559
[5283]	560	! * therefore, in the case icbs=icb, we compute tvp at level icb+1
	561	! (tvp at other levels will be computed in cv3_undilute2.F)
[879]	562
[5283]	563	DO i = 1, len
	564	ticb(i) = t(i, icb(i) + 1)
	565	gzicb(i) = gz(i, icb(i) + 1)
	566	qsicb(i) = qs(i, icb(i) + 1)
	567	END DO
[879]	568
[5283]	569	DO i = 1, len
	570	tg = ticb(i)
	571	qg = qsicb(i) ! convect3
	572	! debug alv=lv0-clmcpv*(ticb(i)-t0)
	573	alv = lv0 - clmcpv * (ticb(i) - 273.15)
[879]	574
[5283]	575	! First iteration.
[879]	576
[5283]	577	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	578	s = cpd * (1. - qnk(i)) + cl * qnk(i) & ! convect3
	579	+ alv * alv * qg / (rrv * ticb(i) * ticb(i)) ! convect3
	580	s = 1. / s
	581	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	582	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gzicb(i) ! convect3
	583	tg = tg + s * (ah0(i) - ahg)
	584	! ori tg=max(tg,35.0)
	585	! debug tc=tg-t0
	586	tc = tg - 273.15
	587	denom = 243.5 + tc
	588	denom = max(denom, 1.0) ! convect3
	589	! ori IF(tc.ge.0.0)THEN
	590	es = 6.112 * exp(17.67 * tc / denom)
	591	! ori else
	592	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	593	! ori endif
	594	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	595	qg = eps * es / (p(i, icb(i) + 1) - es * (1. - eps))
[879]	596
[5283]	597	! Second iteration.
[879]	598
	599
[5283]	600	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	601	! ori s=1./s
	602	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	603	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gzicb(i) ! convect3
	604	tg = tg + s * (ah0(i) - ahg)
	605	! ori tg=max(tg,35.0)
	606	! debug tc=tg-t0
	607	tc = tg - 273.15
	608	denom = 243.5 + tc
	609	denom = max(denom, 1.0) ! convect3
	610	! ori IF(tc.ge.0.0)THEN
	611	es = 6.112 * exp(17.67 * tc / denom)
	612	! ori else
	613	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	614	! ori end if
	615	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	616	qg = eps * es / (p(i, icb(i) + 1) - es * (1. - eps))
[879]	617
[5283]	618	alv = lv0 - clmcpv * (ticb(i) - 273.15)
[879]	619
[5283]	620	! ori c approximation here:
	621	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	622	! ori & -gz(i,icb(i))-alv*qg)/cpd
[879]	623
[5283]	624	! convect3: no approximation:
	625	tp(i, icb(i) + 1) = (ah0(i) - gz(i, icb(i) + 1) - alv * qg) / (cpd + (cl - cpd) * qnk(i))
[879]	626
[5283]	627	! ori clw(i,icb(i))=qnk(i)-qg
	628	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	629	clw(i, icb(i) + 1) = qnk(i) - qg
	630	clw(i, icb(i) + 1) = max(0.0, clw(i, icb(i) + 1))
[879]	631
[5283]	632	rg = qg / (1. - qnk(i))
	633	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	634	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
	635	tvp(i, icb(i) + 1) = tp(i, icb(i) + 1) * (1. + qg / eps - qnk(i)) !whole thing
[879]	636
[5283]	637	END DO
[879]	638
[5283]	639	END SUBROUTINE cv30_undilute1
[879]	640
[5283]	641	SUBROUTINE cv30_trigger(len, nd, icb, plcl, p, th, tv, tvp, pbase, buoybase, &
	642	iflag, sig, w0)
	643	IMPLICIT NONE
[879]	644
[5283]	645	! -------------------------------------------------------------------
	646	! --- TRIGGERING
[879]	647
[5283]	648	! - computes the cloud base
	649	! - triggering (crude in this version)
	650	! - relaxation of sig and w0 when no convection
[879]	651
[5283]	652	! Caution1: if no convection, we set iflag=4
	653	! (it used to be 0 in convect3)
[879]	654
[5283]	655	! Caution2: at this stage, tvp (and thus buoy) are know up
	656	! through icb only!
	657	! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
	658	! -------------------------------------------------------------------
[879]	659
	660
	661
[5283]	662	! input:
	663	INTEGER len, nd
	664	INTEGER icb(len)
	665	REAL plcl(len), p(len, nd)
	666	REAL th(len, nd), tv(len, nd), tvp(len, nd)
[879]	667
[5283]	668	! output:
	669	REAL pbase(len), buoybase(len)
[879]	670
[5283]	671	! input AND output:
	672	INTEGER iflag(len)
	673	REAL sig(len, nd), w0(len, nd)
[879]	674
[5283]	675	! local variables:
	676	INTEGER i, k
	677	REAL tvpbase, tvbase, tdif, ath, ath1
[879]	678
	679
[5283]	680	! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
[879]	681
[5283]	682	DO i = 1, len
	683	pbase(i) = plcl(i) + dpbase
	684	tvpbase = tvp(i, icb(i)) * (pbase(i) - p(i, icb(i) + 1)) / &
	685	(p(i, icb(i)) - p(i, icb(i) + 1)) + tvp(i, icb(i) + 1) * (p(i, icb(i)) - pbase(i)) / (&
	686	p(i, icb(i)) - p(i, icb(i) + 1))
	687	tvbase = tv(i, icb(i)) * (pbase(i) - p(i, icb(i) + 1)) / &
	688	(p(i, icb(i)) - p(i, icb(i) + 1)) + tv(i, icb(i) + 1) * (p(i, icb(i)) - pbase(i)) / (p &
	689	(i, icb(i)) - p(i, icb(i) + 1))
	690	buoybase(i) = tvpbase - tvbase
	691	END DO
[879]	692
	693
[5283]	694	! * make sure that column is dry adiabatic between the surface *
	695	! * and cloud base, and that lifted air is positively buoyant *
	696	! * at cloud base *
	697	! * if not, return to calling program after resetting *
	698	! * sig(i) and w0(i) *
[879]	699
	700
[5283]	701	! oct3 do 200 i=1,len
	702	! oct3
	703	! oct3 tdif = buoybase(i)
	704	! oct3 ath1 = th(i,1)
	705	! oct3 ath = th(i,icb(i)-1) - dttrig
	706	! oct3
	707	! oct3 if (tdif.lt.dtcrit .OR. ath.gt.ath1) THEN
	708	! oct3 do 60 k=1,nl
	709	! oct3 sig(i,k) = betasig(i,k) - 2.alphatdiftdif
	710	! oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
	711	! oct3 w0(i,k) = beta*w0(i,k)
	712	! oct3 60 continue
	713	! oct3 iflag(i)=4 ! pour version vectorisee
	714	! oct3c convect3 iflag(i)=0
	715	! oct3cccc RETURN
	716	! oct3 endif
	717	! oct3
	718	! oct3200 continue
[879]	719
[5283]	720	! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
[879]	721
[5283]	722	DO k = 1, nl
	723	DO i = 1, len
[879]	724
[5283]	725	tdif = buoybase(i)
	726	ath1 = th(i, 1)
	727	ath = th(i, icb(i) - 1) - dttrig
[879]	728
[5283]	729	IF (tdif<dtcrit .OR. ath>ath1) THEN
	730	sig(i, k) = beta * sig(i, k) - 2. * alpha * tdif * tdif
	731	sig(i, k) = amax1(sig(i, k), 0.0)
	732	w0(i, k) = beta * w0(i, k)
	733	iflag(i) = 4 ! pour version vectorisee
	734	! convect3 iflag(i)=0
	735	END IF
[879]	736
[5283]	737	END DO
[1992]	738	END DO
[879]	739
[5283]	740	! fin oct3 --
[879]	741
[5283]	742	END SUBROUTINE cv30_trigger
[879]	743
[5283]	744	SUBROUTINE cv30_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, &
	745	plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, &
	746	th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, &
	747	iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, &
	748	v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0)
	749	USE print_control_mod, ONLY: lunout
	750	IMPLICIT NONE
[879]	751
	752
	753
[5283]	754	! inputs:
	755	INTEGER len, ncum, nd, ntra, nloc
	756	INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
	757	REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
	758	REAL pbase1(len), buoybase1(len)
	759	REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
	760	REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
	761	REAL p1(len, nd), ph1(len, nd + 1), tv1(len, nd), tp1(len, nd)
	762	REAL tvp1(len, nd), clw1(len, nd)
	763	REAL th1(len, nd)
	764	REAL sig1(len, nd), w01(len, nd)
	765	REAL tra1(len, nd, ntra)
[879]	766
[5283]	767	! outputs:
	768	! en fait, on a nloc=len pour l'instant (cf cv_driver)
	769	INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
	770	REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
	771	REAL pbase(nloc), buoybase(nloc)
	772	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
	773	REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
	774	REAL p(nloc, nd), ph(nloc, nd + 1), tv(nloc, nd), tp(nloc, nd)
	775	REAL tvp(nloc, nd), clw(nloc, nd)
	776	REAL th(nloc, nd)
	777	REAL sig(nloc, nd), w0(nloc, nd)
	778	REAL tra(nloc, nd, ntra)
[879]	779
[5283]	780	! local variables:
	781	INTEGER i, k, nn, j
[879]	782
[5283]	783	CHARACTER (LEN = 20) :: modname = 'cv30_compress'
	784	CHARACTER (LEN = 80) :: abort_message
[879]	785
[5283]	786	DO k = 1, nl + 1
	787	nn = 0
	788	DO i = 1, len
	789	IF (iflag1(i)==0) THEN
	790	nn = nn + 1
	791	sig(nn, k) = sig1(i, k)
	792	w0(nn, k) = w01(i, k)
	793	t(nn, k) = t1(i, k)
	794	q(nn, k) = q1(i, k)
	795	qs(nn, k) = qs1(i, k)
	796	u(nn, k) = u1(i, k)
	797	v(nn, k) = v1(i, k)
	798	gz(nn, k) = gz1(i, k)
	799	h(nn, k) = h1(i, k)
	800	lv(nn, k) = lv1(i, k)
	801	cpn(nn, k) = cpn1(i, k)
	802	p(nn, k) = p1(i, k)
	803	ph(nn, k) = ph1(i, k)
	804	tv(nn, k) = tv1(i, k)
	805	tp(nn, k) = tp1(i, k)
	806	tvp(nn, k) = tvp1(i, k)
	807	clw(nn, k) = clw1(i, k)
	808	th(nn, k) = th1(i, k)
	809	END IF
	810	END DO
	811	END DO
	812
	813	! do 121 j=1,ntra
	814	! do 111 k=1,nd
	815	! nn=0
	816	! do 101 i=1,len
	817	! IF(iflag1(i).EQ.0)THEN
	818	! nn=nn+1
	819	! tra(nn,k,j)=tra1(i,k,j)
	820	! END IF
	821	! 101 continue
	822	! 111 continue
	823	! 121 continue
	824
	825	IF (nn/=ncum) THEN
	826	WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
	827	abort_message = ''
	828	CALL abort_physic(modname, abort_message, 1)
	829	END IF
	830
[1992]	831	nn = 0
	832	DO i = 1, len
	833	IF (iflag1(i)==0) THEN
	834	nn = nn + 1
[5283]	835	pbase(nn) = pbase1(i)
	836	buoybase(nn) = buoybase1(i)
	837	plcl(nn) = plcl1(i)
	838	tnk(nn) = tnk1(i)
	839	qnk(nn) = qnk1(i)
	840	gznk(nn) = gznk1(i)
	841	nk(nn) = nk1(i)
	842	icb(nn) = icb1(i)
	843	icbs(nn) = icbs1(i)
	844	iflag(nn) = iflag1(i)
[1992]	845	END IF
	846	END DO
[879]	847
[5283]	848	END SUBROUTINE cv30_compress
[879]	849
[5283]	850	SUBROUTINE cv30_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
	851	q, qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
	852	ep, sigp, buoy)
	853	! epmax_cape: ajout arguments
	854	USE conema3_mod_h
[5346]	855	USE lmdz_cv_ini, ONLY : eps,lv0,rrv,cl,clmcpv,cpd,cpv
[879]	856
[5283]	857	IMPLICIT NONE
[879]	858
[5283]	859	! ---------------------------------------------------------------------
	860	! Purpose:
	861	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	862	! &
	863	! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
	864	! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
	865	! &
	866	! FIND THE LEVEL OF NEUTRAL BUOYANCY
[879]	867
[5283]	868	! Main differences convect3/convect4:
	869	! - icbs (input) is the first level above LCL (may differ from icb)
	870	! - many minor differences in the iterations
	871	! - condensed water not removed from tvp in convect3
	872	! - vertical profile of buoyancy computed here (use of buoybase)
	873	! - the determination of inb is different
	874	! - no inb1, ONLY inb in output
	875	! ---------------------------------------------------------------------
[879]	876
	877
	878
[5283]	879	! inputs:
	880	INTEGER ncum, nd, nloc
	881	INTEGER icb(nloc), icbs(nloc), nk(nloc)
	882	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
	883	REAL p(nloc, nd)
	884	REAL tnk(nloc), qnk(nloc), gznk(nloc)
	885	REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
	886	REAL pbase(nloc), buoybase(nloc), plcl(nloc)
[879]	887
[5283]	888	! outputs:
	889	INTEGER inb(nloc)
	890	REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
	891	REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
	892	REAL buoy(nloc, nd)
[879]	893
[5283]	894	! local variables:
	895	INTEGER i, k
	896	REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
	897	REAL by, defrac, pden
	898	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
	899	LOGICAL lcape(nloc)
[879]	900
[5283]	901	! =====================================================================
	902	! --- SOME INITIALIZATIONS
	903	! =====================================================================
[879]	904
[5283]	905	DO k = 1, nl
	906	DO i = 1, ncum
	907	ep(i, k) = 0.0
	908	sigp(i, k) = spfac
	909	END DO
[1992]	910	END DO
[879]	911
[5283]	912	! =====================================================================
	913	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	914	! =====================================================================
[879]	915
[5283]	916	! --- The procedure is to solve the equation.
	917	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
[879]	918
[5283]	919	! * Calculate certain parcel quantities, including static energy *
[879]	920
[5283]	921	DO i = 1, ncum
	922	ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) & ! debug &
	923	! +qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
	924	+ qnk(i) * (lv0 - clmcpv * (tnk(i) - 273.15)) + gznk(i)
	925	END DO
[879]	926
	927
[5283]	928	! * Find lifted parcel quantities above cloud base *
[879]	929
[5283]	930	DO k = minorig + 1, nl
	931	DO i = 1, ncum
	932	! ori IF(k.ge.(icb(i)+1))THEN
	933	IF (k>=(icbs(i) + 1)) THEN ! convect3
	934	tg = t(i, k)
	935	qg = qs(i, k)
	936	! debug alv=lv0-clmcpv*(t(i,k)-t0)
	937	alv = lv0 - clmcpv * (t(i, k) - 273.15)
[879]	938
[5283]	939	! First iteration.
[879]	940
[5283]	941	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	942	s = cpd * (1. - qnk(i)) + cl * qnk(i) & ! convect3
	943	+ alv * alv * qg / (rrv * t(i, k) * t(i, k)) ! convect3
	944	s = 1. / s
	945	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	946	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k) ! convect3
	947	tg = tg + s * (ah0(i) - ahg)
	948	! ori tg=max(tg,35.0)
	949	! debug tc=tg-t0
	950	tc = tg - 273.15
	951	denom = 243.5 + tc
	952	denom = max(denom, 1.0) ! convect3
	953	! ori IF(tc.ge.0.0)THEN
	954	es = 6.112 * exp(17.67 * tc / denom)
	955	! ori else
	956	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	957	! ori endif
	958	qg = eps * es / (p(i, k) - es * (1. - eps))
[879]	959
[5283]	960	! Second iteration.
[879]	961
[5283]	962	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	963	! ori s=1./s
	964	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	965	ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k) ! convect3
	966	tg = tg + s * (ah0(i) - ahg)
	967	! ori tg=max(tg,35.0)
	968	! debug tc=tg-t0
	969	tc = tg - 273.15
	970	denom = 243.5 + tc
	971	denom = max(denom, 1.0) ! convect3
	972	! ori IF(tc.ge.0.0)THEN
	973	es = 6.112 * exp(17.67 * tc / denom)
	974	! ori else
	975	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	976	! ori endif
	977	qg = eps * es / (p(i, k) - es * (1. - eps))
[879]	978
[5283]	979	! debug alv=lv0-clmcpv*(t(i,k)-t0)
	980	alv = lv0 - clmcpv * (t(i, k) - 273.15)
	981	! PRINT*,'cpd dans convect2 ',cpd
	982	! PRINT*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
	983	! PRINT*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
[879]	984
[5283]	985	! ori c approximation here:
	986	! ori
	987	! tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
[879]	988
[5283]	989	! convect3: no approximation:
	990	tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) / (cpd + (cl - cpd) * qnk(i))
[879]	991
[5283]	992	clw(i, k) = qnk(i) - qg
	993	clw(i, k) = max(0.0, clw(i, k))
	994	rg = qg / (1. - qnk(i))
	995	! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
	996	! convect3: (qg utilise au lieu du vrai mixing ratio rg):
	997	tvp(i, k) = tp(i, k) * (1. + qg / eps - qnk(i)) ! whole thing
	998	END IF
	999	END DO
[1992]	1000	END DO
[879]	1001
[5283]	1002	! =====================================================================
	1003	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
	1004	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
	1005	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
	1006	! =====================================================================
[879]	1007
[5283]	1008	! ori do 320 k=minorig+1,nl
	1009	DO k = 1, nl ! convect3
	1010	DO i = 1, ncum
	1011	pden = ptcrit - pbcrit
	1012	ep(i, k) = (plcl(i) - p(i, k) - pbcrit) / pden * epmax
	1013	ep(i, k) = amax1(ep(i, k), 0.0)
	1014	ep(i, k) = amin1(ep(i, k), epmax)
	1015	sigp(i, k) = spfac
	1016	! ori IF(k.ge.(nk(i)+1))THEN
	1017	! ori tca=tp(i,k)-t0
	1018	! ori IF(tca.ge.0.0)THEN
	1019	! ori elacrit=elcrit
	1020	! ori else
	1021	! ori elacrit=elcrit*(1.0-tca/tlcrit)
	1022	! ori endif
	1023	! ori elacrit=max(elacrit,0.0)
	1024	! ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
	1025	! ori ep(i,k)=max(ep(i,k),0.0 )
	1026	! ori ep(i,k)=min(ep(i,k),1.0 )
	1027	! ori sigp(i,k)=sigs
	1028	! ori endif
	1029	END DO
[1992]	1030	END DO
[879]	1031
[5283]	1032	! =====================================================================
	1033	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
	1034	! --- VIRTUAL TEMPERATURE
	1035	! =====================================================================
[879]	1036
[5283]	1037	! dans convect3, tvp est calcule en une seule fois, et sans retirer
	1038	! l'eau condensee (~> reversible CAPE)
[879]	1039
[5283]	1040	! ori do 340 k=minorig+1,nl
	1041	! ori do 330 i=1,ncum
	1042	! ori IF(k.ge.(icb(i)+1))THEN
	1043	! ori tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
	1044	! oric PRINT*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
	1045	! oric PRINT*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
	1046	! ori endif
	1047	! ori 330 continue
	1048	! ori 340 continue
[879]	1049
[5283]	1050	! ori do 350 i=1,ncum
	1051	! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
	1052	! ori 350 continue
[879]	1053
[5283]	1054	DO i = 1, ncum ! convect3
	1055	tp(i, nlp) = tp(i, nl) ! convect3
	1056	END DO ! convect3
[879]	1057
[5283]	1058	! =====================================================================
	1059	! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
	1060	! =====================================================================
[879]	1061
[5283]	1062	! -- this is for convect3 only:
[879]	1063
[5283]	1064	! first estimate of buoyancy:
[879]	1065
[5283]	1066	DO i = 1, ncum
	1067	DO k = 1, nl
	1068	buoy(i, k) = tvp(i, k) - tv(i, k)
	1069	END DO
[1992]	1070	END DO
[879]	1071
[5283]	1072	! set buoyancy=buoybase for all levels below base
	1073	! for safety, set buoy(icb)=buoybase
[879]	1074
[5283]	1075	DO i = 1, ncum
	1076	DO k = 1, nl
	1077	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i, k)>=pbase(i))) THEN
	1078	buoy(i, k) = buoybase(i)
	1079	END IF
	1080	END DO
	1081	! IM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i)
	1082	buoy(i, icb(i)) = buoybase(i)
[1992]	1083	END DO
[879]	1084
[5283]	1085	! -- end convect3
[879]	1086
[5283]	1087	! =====================================================================
	1088	! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
	1089	! --- LEVEL OF NEUTRAL BUOYANCY
	1090	! =====================================================================
[879]	1091
[5283]	1092	! -- this is for convect3 only:
[879]	1093
[5283]	1094	DO i = 1, ncum
	1095	inb(i) = nl - 1
	1096	END DO
[879]	1097
[5283]	1098	DO i = 1, ncum
	1099	DO k = 1, nl - 1
	1100	IF ((k>=icb(i)) .AND. (buoy(i, k)<dtovsh)) THEN
	1101	inb(i) = min(inb(i), k)
	1102	END IF
	1103	END DO
[1992]	1104	END DO
[879]	1105
[5283]	1106	! -- end convect3
[879]	1107
[5283]	1108	! ori do 510 i=1,ncum
	1109	! ori cape(i)=0.0
	1110	! ori capem(i)=0.0
	1111	! ori inb(i)=icb(i)+1
	1112	! ori inb1(i)=inb(i)
	1113	! ori 510 continue
[879]	1114
[5283]	1115	! Originial Code
[879]	1116
[5283]	1117	! do 530 k=minorig+1,nl-1
	1118	! do 520 i=1,ncum
	1119	! IF(k.ge.(icb(i)+1))THEN
	1120	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1121	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1122	! cape(i)=cape(i)+by
	1123	! IF(by.ge.0.0)inb1(i)=k+1
	1124	! IF(cape(i).gt.0.0)THEN
	1125	! inb(i)=k+1
	1126	! capem(i)=cape(i)
	1127	! END IF
	1128	! END IF
	1129	! 520 continue
	1130	! 530 continue
	1131	! do 540 i=1,ncum
	1132	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
	1133	! cape(i)=capem(i)+byp
	1134	! defrac=capem(i)-cape(i)
	1135	! defrac=max(defrac,0.001)
	1136	! frac(i)=-cape(i)/defrac
	1137	! frac(i)=min(frac(i),1.0)
	1138	! frac(i)=max(frac(i),0.0)
	1139	! 540 continue
[879]	1140
[5283]	1141	! K Emanuel fix
[879]	1142
[5283]	1143	! CALL zilch(byp,ncum)
	1144	! do 530 k=minorig+1,nl-1
	1145	! do 520 i=1,ncum
	1146	! IF(k.ge.(icb(i)+1))THEN
	1147	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1148	! cape(i)=cape(i)+by
	1149	! IF(by.ge.0.0)inb1(i)=k+1
	1150	! IF(cape(i).gt.0.0)THEN
	1151	! inb(i)=k+1
	1152	! capem(i)=cape(i)
	1153	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1154	! END IF
	1155	! END IF
	1156	! 520 continue
	1157	! 530 continue
	1158	! do 540 i=1,ncum
	1159	! inb(i)=max(inb(i),inb1(i))
	1160	! cape(i)=capem(i)+byp(i)
	1161	! defrac=capem(i)-cape(i)
	1162	! defrac=max(defrac,0.001)
	1163	! frac(i)=-cape(i)/defrac
	1164	! frac(i)=min(frac(i),1.0)
	1165	! frac(i)=max(frac(i),0.0)
	1166	! 540 continue
[879]	1167
[5283]	1168	! J Teixeira fix
[879]	1169
[5283]	1170	! ori CALL zilch(byp,ncum)
	1171	! ori do 515 i=1,ncum
	1172	! ori lcape(i)=.TRUE.
	1173	! ori 515 continue
	1174	! ori do 530 k=minorig+1,nl-1
	1175	! ori do 520 i=1,ncum
	1176	! ori IF(cape(i).lt.0.0)lcape(i)=.FALSE.
	1177	! ori if((k.ge.(icb(i)+1)).AND.lcape(i))THEN
	1178	! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1179	! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1180	! ori cape(i)=cape(i)+by
	1181	! ori IF(by.ge.0.0)inb1(i)=k+1
	1182	! ori IF(cape(i).gt.0.0)THEN
	1183	! ori inb(i)=k+1
	1184	! ori capem(i)=cape(i)
	1185	! ori endif
	1186	! ori endif
	1187	! ori 520 continue
	1188	! ori 530 continue
	1189	! ori do 540 i=1,ncum
	1190	! ori cape(i)=capem(i)+byp(i)
	1191	! ori defrac=capem(i)-cape(i)
	1192	! ori defrac=max(defrac,0.001)
	1193	! ori frac(i)=-cape(i)/defrac
	1194	! ori frac(i)=min(frac(i),1.0)
	1195	! ori frac(i)=max(frac(i),0.0)
	1196	! ori 540 continue
[879]	1197
[5283]	1198	! =====================================================================
	1199	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
	1200	! =====================================================================
[879]	1201
[5283]	1202	! ym do i=1,ncum*nlp
	1203	! ym hp(i,1)=h(i,1)
	1204	! ym enddo
[879]	1205
[5283]	1206	DO k = 1, nlp
	1207	DO i = 1, ncum
	1208	hp(i, k) = h(i, k)
	1209	END DO
[1992]	1210	END DO
[879]	1211
[5283]	1212	DO k = minorig + 1, nl
	1213	DO i = 1, ncum
	1214	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
	1215	hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k &
	1216	)
	1217	END IF
	1218	END DO
[1992]	1219	END DO
[879]	1220
[5283]	1221	END SUBROUTINE cv30_undilute2
[879]	1222
[5283]	1223	SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
	1224	sig, w0, cape, m)
[5346]	1225	USE lmdz_cv_ini, ONLY : rrd
[879]	1226
[5283]	1227	IMPLICIT NONE
[879]	1228
[5283]	1229	! ===================================================================
	1230	! --- CLOSURE OF CONVECT3
[879]	1231
[5283]	1232	! vectorization: S. Bony
	1233	! ===================================================================
[879]	1234
	1235
	1236
[5283]	1237	! input:
	1238	INTEGER ncum, nd, nloc
	1239	INTEGER icb(nloc), inb(nloc)
	1240	REAL pbase(nloc)
	1241	REAL p(nloc, nd), ph(nloc, nd + 1)
	1242	REAL tv(nloc, nd), buoy(nloc, nd)
[879]	1243
[5283]	1244	! input/output:
	1245	REAL sig(nloc, nd), w0(nloc, nd)
[879]	1246
[5283]	1247	! output:
	1248	REAL cape(nloc)
	1249	REAL m(nloc, nd)
[879]	1250
[5283]	1251	! local variables:
	1252	INTEGER i, j, k, icbmax
	1253	REAL deltap, fac, w, amu
	1254	REAL dtmin(nloc, nd), sigold(nloc, nd)
	1255
	1256	! -------------------------------------------------------
	1257	! -- Initialization
	1258	! -------------------------------------------------------
	1259
	1260	DO k = 1, nl
	1261	DO i = 1, ncum
	1262	m(i, k) = 0.0
	1263	END DO
[1992]	1264	END DO
[879]	1265
[5283]	1266	! -------------------------------------------------------
	1267	! -- Reset sig(i) and w0(i) for i>inb and i<icb
	1268	! -------------------------------------------------------
[879]	1269
[5283]	1270	! update sig and w0 above LNB:
[879]	1271
[5283]	1272	DO k = 1, nl - 1
	1273	DO i = 1, ncum
	1274	IF ((inb(i)<(nl - 1)) .AND. (k>=(inb(i) + 1))) THEN
	1275	sig(i, k) = beta * sig(i, k) + 2. * alpha * buoy(i, inb(i)) * abs(buoy(i, inb(&
	1276	i)))
	1277	sig(i, k) = amax1(sig(i, k), 0.0)
	1278	w0(i, k) = beta * w0(i, k)
	1279	END IF
	1280	END DO
[1992]	1281	END DO
[879]	1282
[5283]	1283	! compute icbmax:
[879]	1284
[5283]	1285	icbmax = 2
	1286	DO i = 1, ncum
	1287	icbmax = max(icbmax, icb(i))
	1288	END DO
[879]	1289
[5283]	1290	! update sig and w0 below cloud base:
[879]	1291
[5283]	1292	DO k = 1, icbmax
	1293	DO i = 1, ncum
	1294	IF (k<=icb(i)) THEN
	1295	sig(i, k) = beta * sig(i, k) - 2. * alpha * buoy(i, icb(i)) * buoy(i, icb(i))
	1296	sig(i, k) = amax1(sig(i, k), 0.0)
	1297	w0(i, k) = beta * w0(i, k)
	1298	END IF
	1299	END DO
[1992]	1300	END DO
[879]	1301
[5283]	1302	! IF(inb.lt.(nl-1))THEN
	1303	! do 85 i=inb+1,nl-1
	1304	! sig(i)=betasig(i)+2.alphabuoy(inb)
	1305	! 1 abs(buoy(inb))
	1306	! sig(i)=amax1(sig(i),0.0)
	1307	! w0(i)=beta*w0(i)
	1308	! 85 continue
	1309	! end if
[879]	1310
[5283]	1311	! do 87 i=1,icb
	1312	! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
	1313	! sig(i)=amax1(sig(i),0.0)
	1314	! w0(i)=beta*w0(i)
	1315	! 87 continue
[1742]	1316
[5283]	1317	! -------------------------------------------------------------
	1318	! -- Reset fractional areas of updrafts and w0 at initial time
	1319	! -- and after 10 time steps of no convection
	1320	! -------------------------------------------------------------
[1742]	1321
[5283]	1322	DO k = 1, nl - 1
	1323	DO i = 1, ncum
	1324	IF (sig(i, nd)<1.5 .OR. sig(i, nd)>12.0) THEN
	1325	sig(i, k) = 0.0
	1326	w0(i, k) = 0.0
	1327	END IF
	1328	END DO
[1992]	1329	END DO
[879]	1330
[5283]	1331	! -------------------------------------------------------------
	1332	! -- Calculate convective available potential energy (cape),
	1333	! -- vertical velocity (w), fractional area covered by
	1334	! -- undilute updraft (sig), and updraft mass flux (m)
	1335	! -------------------------------------------------------------
[879]	1336
[5283]	1337	DO i = 1, ncum
	1338	cape(i) = 0.0
	1339	END DO
[879]	1340
[5283]	1341	! compute dtmin (minimum buoyancy between ICB and given level k):
[879]	1342
[5283]	1343	DO i = 1, ncum
	1344	DO k = 1, nl
	1345	dtmin(i, k) = 100.0
	1346	END DO
[1992]	1347	END DO
[879]	1348
[5283]	1349	DO i = 1, ncum
	1350	DO k = 1, nl
	1351	DO j = minorig, nl
	1352	IF ((k>=(icb(i) + 1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k - &
	1353	1))) THEN
	1354	dtmin(i, k) = amin1(dtmin(i, k), buoy(i, j))
	1355	END IF
	1356	END DO
[1992]	1357	END DO
	1358	END DO
[879]	1359
[5283]	1360	! the interval on which cape is computed starts at pbase :
	1361	DO k = 1, nl
	1362	DO i = 1, ncum
[879]	1363
[5283]	1364	IF ((k>=(icb(i) + 1)) .AND. (k<=inb(i))) THEN
[879]	1365
[5283]	1366	deltap = min(pbase(i), ph(i, k - 1)) - min(pbase(i), ph(i, k))
	1367	cape(i) = cape(i) + rrd * buoy(i, k - 1) * deltap / p(i, k - 1)
	1368	cape(i) = amax1(0.0, cape(i))
	1369	sigold(i, k) = sig(i, k)
[879]	1370
[5283]	1371	! dtmin(i,k)=100.0
	1372	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
	1373	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
	1374	! 97 continue
[879]	1375
[5283]	1376	sig(i, k) = beta * sig(i, k) + alpha * dtmin(i, k) * abs(dtmin(i, k))
	1377	sig(i, k) = amax1(sig(i, k), 0.0)
	1378	sig(i, k) = amin1(sig(i, k), 0.01)
	1379	fac = amin1(((dtcrit - dtmin(i, k)) / dtcrit), 1.0)
	1380	w = (1. - beta) * fac * sqrt(cape(i)) + beta * w0(i, k)
	1381	amu = 0.5 * (sig(i, k) + sigold(i, k)) * w
	1382	m(i, k) = amu * 0.007 * p(i, k) * (ph(i, k) - ph(i, k + 1)) / tv(i, k)
	1383	w0(i, k) = w
	1384	END IF
[879]	1385
[5283]	1386	END DO
[1992]	1387	END DO
[879]	1388
[5283]	1389	DO i = 1, ncum
	1390	w0(i, icb(i)) = 0.5 * w0(i, icb(i) + 1)
	1391	m(i, icb(i)) = 0.5 * m(i, icb(i) + 1) * (ph(i, icb(i)) - ph(i, icb(i) + 1)) / &
	1392	(ph(i, icb(i) + 1) - ph(i, icb(i) + 2))
	1393	sig(i, icb(i)) = sig(i, icb(i) + 1)
	1394	sig(i, icb(i) - 1) = sig(i, icb(i))
	1395	END DO
[879]	1396
	1397
[5283]	1398	! cape=0.0
	1399	! do 98 i=icb+1,inb
	1400	! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
	1401	! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
	1402	! dcape=rrdbuoy(i-1)deltap/p(i-1)
	1403	! dlnp=deltap/p(i-1)
	1404	! cape=amax1(0.0,cape)
	1405	! sigold=sig(i)
[879]	1406
[5283]	1407	! dtmin=100.0
	1408	! do 97 j=icb,i-1
	1409	! dtmin=amin1(dtmin,buoy(j))
	1410	! 97 continue
[879]	1411
[5283]	1412	! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
	1413	! sig(i)=amax1(sig(i),0.0)
	1414	! sig(i)=amin1(sig(i),0.01)
	1415	! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
	1416	! w=(1.-beta)facsqrt(cape)+beta*w0(i)
	1417	! amu=0.5(sig(i)+sigold)w
	1418	! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
	1419	! w0(i)=w
	1420	! 98 continue
	1421	! w0(icb)=0.5*w0(icb+1)
	1422	! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
	1423	! sig(icb)=sig(icb+1)
	1424	! sig(icb-1)=sig(icb)
[879]	1425
[5283]	1426	END SUBROUTINE cv30_closure
[879]	1427
[5283]	1428	SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
	1429	u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, &
	1430	vent, sij, elij, ments, qents, traent)
[5346]	1431	USE lmdz_cv_ini, ONLY : cpd,cpv,rrv
[879]	1432
[5283]	1433	IMPLICIT NONE
[879]	1434
[5283]	1435	! ---------------------------------------------------------------------
	1436	! a faire:
	1437	! - changer rr(il,1) -> qnk(il)
	1438	! - vectorisation de la partie normalisation des flux (do 789...)
	1439	! ---------------------------------------------------------------------
[879]	1440
	1441
	1442
[5283]	1443	! inputs:
	1444	INTEGER ncum, nd, na, ntra, nloc
	1445	INTEGER icb(nloc), inb(nloc), nk(nloc)
	1446	REAL sig(nloc, nd)
	1447	REAL qnk(nloc)
	1448	REAL ph(nloc, nd + 1)
	1449	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
	1450	REAL u(nloc, nd), v(nloc, nd)
	1451	REAL tra(nloc, nd, ntra) ! input of convect3
	1452	REAL lv(nloc, na), h(nloc, na), hp(nloc, na)
	1453	REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
	1454	REAL m(nloc, na) ! input of convect3
[879]	1455
[5283]	1456	! outputs:
	1457	REAL ment(nloc, na, na), qent(nloc, na, na)
	1458	REAL uent(nloc, na, na), vent(nloc, na, na)
	1459	REAL sij(nloc, na, na), elij(nloc, na, na)
	1460	REAL traent(nloc, nd, nd, ntra)
	1461	REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
	1462	REAL sigij(nloc, nd, nd)
[879]	1463
[5283]	1464	! local variables:
	1465	INTEGER i, j, k, il, im, jm
	1466	INTEGER num1, num2
	1467	INTEGER nent(nloc, na)
	1468	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
	1469	REAL alt, smid, sjmin, sjmax, delp, delm
	1470	REAL asij(nloc), smax(nloc), scrit(nloc)
	1471	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
	1472	REAL wgh
	1473	REAL zm(nloc, na)
	1474	LOGICAL lwork(nloc)
[879]	1475
[5283]	1476	! =====================================================================
	1477	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	1478	! =====================================================================
	1479
	1480	! ori do 360 i=1,ncum*nlp
	1481	DO j = 1, nl
[1992]	1482	DO i = 1, ncum
[5283]	1483	nent(i, j) = 0
	1484	! in convect3, m is computed in cv3_closure
	1485	! ori m(i,1)=0.0
[1992]	1486	END DO
	1487	END DO
[879]	1488
[5283]	1489	! ori do 400 k=1,nlp
	1490	! ori do 390 j=1,nlp
	1491	DO j = 1, nl
	1492	DO k = 1, nl
	1493	DO i = 1, ncum
	1494	qent(i, k, j) = rr(i, j)
	1495	uent(i, k, j) = u(i, j)
	1496	vent(i, k, j) = v(i, j)
	1497	elij(i, k, j) = 0.0
	1498	! ym ment(i,k,j)=0.0
	1499	! ym sij(i,k,j)=0.0
	1500	END DO
[1992]	1501	END DO
	1502	END DO
[879]	1503
[5283]	1504	! ym
	1505	ment(1:ncum, 1:nd, 1:nd) = 0.0
	1506	sij(1:ncum, 1:nd, 1:nd) = 0.0
	1507
[1992]	1508	! do k=1,ntra
[5283]	1509	! do j=1,nd ! instead nlp
	1510	! do i=1,nd ! instead nlp
[1992]	1511	! do il=1,ncum
[5283]	1512	! traent(il,i,j,k)=tra(il,j,k)
[1992]	1513	! enddo
	1514	! enddo
	1515	! enddo
[5283]	1516	! enddo
	1517	zm(:, :) = 0.
[879]	1518
[5283]	1519	! =====================================================================
	1520	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	1521	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	1522	! --- FRACTION (sij)
	1523	! =====================================================================
[879]	1524
[5283]	1525	DO i = minorig + 1, nl
[879]	1526
[5283]	1527	DO j = minorig, nl
	1528	DO il = 1, ncum
	1529	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il) - &
	1530	1)) .AND. (j<=inb(il))) THEN
[879]	1531
[5283]	1532	rti = rr(il, 1) - ep(il, i) * clw(il, i)
	1533	bf2 = 1. + lv(il, j) * lv(il, j) * rs(il, j) / (rrv * t(il, j) * t(il, j) * cpd)
	1534	anum = h(il, j) - hp(il, i) + (cpv - cpd) * t(il, j) * (rti - rr(il, j))
	1535	denom = h(il, i) - hp(il, i) + (cpd - cpv) * (rr(il, i) - rti) * t(il, j)
	1536	dei = denom
	1537	IF (abs(dei)<0.01) dei = 0.01
	1538	sij(il, i, j) = anum / dei
	1539	sij(il, i, i) = 1.0
	1540	altem = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) * rti - rs(il, j)
	1541	altem = altem / bf2
	1542	cwat = clw(il, j) * (1. - ep(il, j))
	1543	stemp = sij(il, i, j)
	1544	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
	1545	anum = anum - lv(il, j) * (rti - rs(il, j) - cwat * bf2)
	1546	denom = denom + lv(il, j) * (rr(il, i) - rti)
	1547	IF (abs(denom)<0.01) denom = 0.01
	1548	sij(il, i, j) = anum / denom
	1549	altem = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) * rti - &
	1550	rs(il, j)
	1551	altem = altem - (bf2 - 1.) * cwat
	1552	END IF
	1553	IF (sij(il, i, j)>0.0 .AND. sij(il, i, j)<0.95) THEN
	1554	qent(il, i, j) = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) * rti
	1555	uent(il, i, j) = sij(il, i, j) * u(il, i) + &
	1556	(1. - sij(il, i, j)) * u(il, nk(il))
	1557	vent(il, i, j) = sij(il, i, j) * v(il, i) + &
	1558	(1. - sij(il, i, j)) * v(il, nk(il))
	1559	! !!! do k=1,ntra
	1560	! !!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1561	! !!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1562	! !!! END DO
	1563	elij(il, i, j) = altem
	1564	elij(il, i, j) = amax1(0.0, elij(il, i, j))
	1565	ment(il, i, j) = m(il, i) / (1. - sij(il, i, j))
	1566	nent(il, i) = nent(il, i) + 1
	1567	END IF
	1568	sij(il, i, j) = amax1(0.0, sij(il, i, j))
	1569	sij(il, i, j) = amin1(1.0, sij(il, i, j))
	1570	END IF ! new
	1571	END DO
	1572	END DO
[879]	1573
[5283]	1574	! do k=1,ntra
	1575	! do j=minorig,nl
	1576	! do il=1,ncum
	1577	! IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
	1578	! : (j.ge.(icb(il)-1)).AND.(j.le.inb(il)))THEN
	1579	! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1580	! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1581	! END IF
	1582	! enddo
	1583	! enddo
	1584	! enddo
[879]	1585
	1586
[5283]	1587	! *** if no air can entrain at level i assume that updraft detrains
	1588	! ***
	1589	! *** at that level and calculate detrained air flux and properties
	1590	! ***
	1591
	1592
	1593	! @ do 170 i=icb(il),inb(il)
	1594
[1992]	1595	DO il = 1, ncum
[5283]	1596	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il, i)==0)) THEN
	1597	! @ IF(nent(il,i).EQ.0)THEN
	1598	ment(il, i, i) = m(il, i)
	1599	qent(il, i, i) = rr(il, nk(il)) - ep(il, i) * clw(il, i)
	1600	uent(il, i, i) = u(il, nk(il))
	1601	vent(il, i, i) = v(il, nk(il))
	1602	elij(il, i, i) = clw(il, i)
	1603	! MAF sij(il,i,i)=1.0
	1604	sij(il, i, i) = 0.0
[1992]	1605	END IF
	1606	END DO
	1607	END DO
[879]	1608
[5283]	1609	! do j=1,ntra
	1610	! do i=minorig+1,nl
	1611	! do il=1,ncum
	1612	! if (i.ge.icb(il) .AND. i.le.inb(il) .AND. nent(il,i).EQ.0) THEN
	1613	! traent(il,i,i,j)=tra(il,nk(il),j)
	1614	! END IF
	1615	! enddo
	1616	! enddo
	1617	! enddo
[879]	1618
[5283]	1619	DO j = minorig, nl
	1620	DO i = minorig, nl
	1621	DO il = 1, ncum
	1622	IF ((j>=(icb(il) - 1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
	1623	inb(il))) THEN
	1624	sigij(il, i, j) = sij(il, i, j)
	1625	END IF
	1626	END DO
	1627	END DO
	1628	END DO
	1629	! @ enddo
[879]	1630
[5283]	1631	! @170 continue
[879]	1632
[5283]	1633	! =====================================================================
	1634	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	1635	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	1636	! =====================================================================
[879]	1637
[5283]	1638	! ym CALL zilch(asum,ncum*nd)
	1639	! ym CALL zilch(bsum,ncum*nd)
	1640	! ym CALL zilch(csum,ncum*nd)
	1641	CALL zilch(asum, nloc * nd)
	1642	CALL zilch(csum, nloc * nd)
	1643	CALL zilch(csum, nloc * nd)
[879]	1644
[1992]	1645	DO il = 1, ncum
[5283]	1646	lwork(il) = .FALSE.
[1992]	1647	END DO
[879]	1648
[5283]	1649	DO i = minorig + 1, nl
[879]	1650
[5283]	1651	num1 = 0
[1992]	1652	DO il = 1, ncum
[5283]	1653	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
[1992]	1654	END DO
[5283]	1655	IF (num1<=0) GO TO 789
[879]	1656
[1992]	1657	DO il = 1, ncum
[5283]	1658	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	1659	lwork(il) = (nent(il, i)/=0)
	1660	qp = rr(il, 1) - ep(il, i) * clw(il, i)
	1661	anum = h(il, i) - hp(il, i) - lv(il, i) * (qp - rs(il, i)) + &
	1662	(cpv - cpd) * t(il, i) * (qp - rr(il, i))
	1663	denom = h(il, i) - hp(il, i) + lv(il, i) * (rr(il, i) - qp) + &
	1664	(cpd - cpv) * t(il, i) * (rr(il, i) - qp)
	1665	IF (abs(denom)<0.01) denom = 0.01
	1666	scrit(il) = anum / denom
	1667	alt = qp - rs(il, i) + scrit(il) * (rr(il, i) - qp)
	1668	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
	1669	smax(il) = 0.0
	1670	asij(il) = 0.0
	1671	END IF
	1672	END DO
[879]	1673
[5283]	1674	DO j = nl, minorig, -1
	1675
	1676	num2 = 0
	1677	DO il = 1, ncum
	1678	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb(&
	1679	il) - 1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
	1680	END DO
	1681	IF (num2<=0) GO TO 175
	1682
	1683	DO il = 1, ncum
	1684	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb(&
	1685	il) - 1) .AND. j<=inb(il) .AND. lwork(il)) THEN
	1686
	1687	IF (sij(il, i, j)>1.0E-16 .AND. sij(il, i, j)<0.95) THEN
	1688	wgh = 1.0
	1689	IF (j>i) THEN
	1690	sjmax = amax1(sij(il, i, j + 1), smax(il))
	1691	sjmax = amin1(sjmax, scrit(il))
	1692	smax(il) = amax1(sij(il, i, j), smax(il))
	1693	sjmin = amax1(sij(il, i, j - 1), smax(il))
	1694	sjmin = amin1(sjmin, scrit(il))
	1695	IF (sij(il, i, j)<(smax(il) - 1.0E-16)) wgh = 0.0
	1696	smid = amin1(sij(il, i, j), scrit(il))
	1697	ELSE
	1698	sjmax = amax1(sij(il, i, j + 1), scrit(il))
	1699	smid = amax1(sij(il, i, j), scrit(il))
	1700	sjmin = 0.0
	1701	IF (j>1) sjmin = sij(il, i, j - 1)
	1702	sjmin = amax1(sjmin, scrit(il))
	1703	END IF
	1704	delp = abs(sjmax - smid)
	1705	delm = abs(sjmin - smid)
	1706	asij(il) = asij(il) + wgh * (delp + delm)
	1707	ment(il, i, j) = ment(il, i, j) * (delp + delm) * wgh
[1992]	1708	END IF
	1709	END IF
[5283]	1710	END DO
[879]	1711
[5283]	1712	175 END DO
[879]	1713
[1992]	1714	DO il = 1, ncum
[5283]	1715	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	1716	asij(il) = amax1(1.0E-16, asij(il))
	1717	asij(il) = 1.0 / asij(il)
	1718	asum(il, i) = 0.0
	1719	bsum(il, i) = 0.0
	1720	csum(il, i) = 0.0
[1992]	1721	END IF
	1722	END DO
[879]	1723
[5283]	1724	DO j = minorig, nl
	1725	DO il = 1, ncum
	1726	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(&
	1727	il) - 1) .AND. j<=inb(il)) THEN
	1728	ment(il, i, j) = ment(il, i, j) * asij(il)
	1729	END IF
	1730	END DO
[1992]	1731	END DO
[879]	1732
[5283]	1733	DO j = minorig, nl
	1734	DO il = 1, ncum
	1735	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(&
	1736	il) - 1) .AND. j<=inb(il)) THEN
	1737	asum(il, i) = asum(il, i) + ment(il, i, j)
	1738	ment(il, i, j) = ment(il, i, j) * sig(il, j)
	1739	bsum(il, i) = bsum(il, i) + ment(il, i, j)
	1740	END IF
	1741	END DO
	1742	END DO
[879]	1743
[1992]	1744	DO il = 1, ncum
[5283]	1745	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	1746	bsum(il, i) = amax1(bsum(il, i), 1.0E-16)
	1747	bsum(il, i) = 1.0 / bsum(il, i)
[1992]	1748	END IF
	1749	END DO
[879]	1750
[5283]	1751	DO j = minorig, nl
	1752	DO il = 1, ncum
	1753	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(&
	1754	il) - 1) .AND. j<=inb(il)) THEN
	1755	ment(il, i, j) = ment(il, i, j) * asum(il, i) * bsum(il, i)
	1756	END IF
	1757	END DO
	1758	END DO
	1759
	1760	DO j = minorig, nl
	1761	DO il = 1, ncum
	1762	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(&
	1763	il) - 1) .AND. j<=inb(il)) THEN
	1764	csum(il, i) = csum(il, i) + ment(il, i, j)
	1765	END IF
	1766	END DO
	1767	END DO
	1768
[1992]	1769	DO il = 1, ncum
[5283]	1770	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	1771	csum(il, i)<m(il, i)) THEN
	1772	nent(il, i) = 0
	1773	ment(il, i, i) = m(il, i)
	1774	qent(il, i, i) = rr(il, 1) - ep(il, i) * clw(il, i)
	1775	uent(il, i, i) = u(il, nk(il))
	1776	vent(il, i, i) = v(il, nk(il))
	1777	elij(il, i, i) = clw(il, i)
	1778	! MAF sij(il,i,i)=1.0
	1779	sij(il, i, i) = 0.0
[1992]	1780	END IF
[5283]	1781	END DO ! il
[879]	1782
[5283]	1783	! do j=1,ntra
	1784	! do il=1,ncum
	1785	! if ( i.ge.icb(il) .AND. i.le.inb(il) .AND. lwork(il)
	1786	! : .AND. csum(il,i).lt.m(il,i) ) THEN
	1787	! traent(il,i,i,j)=tra(il,nk(il),j)
	1788	! END IF
	1789	! enddo
	1790	! enddo
	1791	789 END DO
[879]	1792
[5283]	1793	! MAF: renormalisation de MENT
	1794	DO jm = 1, nd
	1795	DO im = 1, nd
	1796	DO il = 1, ncum
	1797	zm(il, im) = zm(il, im) + (1. - sij(il, im, jm)) * ment(il, im, jm)
	1798	END DO
[1992]	1799	END DO
	1800	END DO
[879]	1801
[5283]	1802	DO jm = 1, nd
	1803	DO im = 1, nd
	1804	DO il = 1, ncum
	1805	IF (zm(il, im)/=0.) THEN
	1806	ment(il, im, jm) = ment(il, im, jm) * m(il, im) / zm(il, im)
	1807	END IF
	1808	END DO
[1992]	1809	END DO
	1810	END DO
[879]	1811
[5283]	1812	DO jm = 1, nd
	1813	DO im = 1, nd
	1814	DO il = 1, ncum
	1815	qents(il, im, jm) = qent(il, im, jm)
	1816	ments(il, im, jm) = ment(il, im, jm)
	1817	END DO
[1992]	1818	END DO
	1819	END DO
[879]	1820
[5283]	1821	END SUBROUTINE cv30_mixing
[879]	1822
	1823
[5283]	1824	SUBROUTINE cv30_unsat(nloc, ncum, nd, na, ntra, icb, inb, t, rr, rs, gz, u, &
	1825	v, tra, p, ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, &
	1826	mp, rp, up, vp, trap, wt, water, evap, b & ! RomP-jyg
	1827	, wdtraina, wdtrainm) ! 26/08/10 RomP-jyg
	1828	USE cvflag_mod_h
[5346]	1829	USE lmdz_cv_ini, ONLY : cpd,grav,ginv
[879]	1830
[5283]	1831	IMPLICIT NONE
[879]	1832
	1833
	1834
[5283]	1835	! inputs:
	1836	INTEGER ncum, nd, na, ntra, nloc
	1837	INTEGER icb(nloc), inb(nloc)
	1838	REAL delt, plcl(nloc)
	1839	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
	1840	REAL u(nloc, nd), v(nloc, nd)
	1841	REAL tra(nloc, nd, ntra)
	1842	REAL p(nloc, nd), ph(nloc, nd + 1)
	1843	REAL th(nloc, na), gz(nloc, na)
	1844	REAL lv(nloc, na), ep(nloc, na), sigp(nloc, na), clw(nloc, na)
	1845	REAL cpn(nloc, na), tv(nloc, na)
	1846	REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
[879]	1847
[5283]	1848	! outputs:
	1849	REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
	1850	REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
	1851	REAL trap(nloc, na, ntra)
	1852	REAL b(nloc, na)
	1853	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
	1854	! lascendance adiabatique et des flux melanges Pa et Pm.
	1855	! Distinction des wdtrain
	1856	! Pa = wdtrainA Pm = wdtrainM
	1857	REAL wdtraina(nloc, na), wdtrainm(nloc, na)
[879]	1858
[5283]	1859	! local variables
	1860	INTEGER i, j, k, il, num1
	1861	REAL tinv, delti
	1862	REAL awat, afac, afac1, afac2, bfac
	1863	REAL pr1, pr2, sigt, b6, c6, revap, tevap, delth
	1864	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
	1865	REAL ampmax
	1866	REAL lvcp(nloc, na)
	1867	REAL wdtrain(nloc)
	1868	LOGICAL lwork(nloc)
[879]	1869
	1870
[5283]	1871	! ------------------------------------------------------
[879]	1872
[5283]	1873	delti = 1. / delt
	1874	tinv = 1. / 3.
[879]	1875
[5283]	1876	mp(:, :) = 0.
[879]	1877
[5283]	1878	DO i = 1, nl
	1879	DO il = 1, ncum
	1880	mp(il, i) = 0.0
	1881	rp(il, i) = rr(il, i)
	1882	up(il, i) = u(il, i)
	1883	vp(il, i) = v(il, i)
	1884	wt(il, i) = 0.001
	1885	water(il, i) = 0.0
	1886	evap(il, i) = 0.0
	1887	b(il, i) = 0.0
	1888	lvcp(il, i) = lv(il, i) / cpn(il, i)
	1889	END DO
	1890	END DO
[879]	1891
[5283]	1892	! do k=1,ntra
	1893	! do i=1,nd
	1894	! do il=1,ncum
	1895	! trap(il,i,k)=tra(il,i,k)
	1896	! enddo
	1897	! enddo
	1898	! enddo
	1899	! RomP >>>
	1900	DO i = 1, nd
	1901	DO il = 1, ncum
	1902	wdtraina(il, i) = 0.0
	1903	wdtrainm(il, i) = 0.0
	1904	END DO
	1905	END DO
	1906	! RomP <<<
[879]	1907
[5283]	1908	! * check whether ep(inb)=0, if so, skip precipitating *
	1909	! * downdraft calculation *
[879]	1910
[1992]	1911	DO il = 1, ncum
[5283]	1912	lwork(il) = .TRUE.
	1913	IF (ep(il, inb(il))<0.0001) lwork(il) = .FALSE.
[1992]	1914	END DO
[879]	1915
[5283]	1916	CALL zilch(wdtrain, ncum)
[879]	1917
[5283]	1918	DO i = nl + 1, 1, -1
[879]	1919
[5283]	1920	num1 = 0
	1921	DO il = 1, ncum
	1922	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
	1923	END DO
	1924	IF (num1<=0) GO TO 400
[879]	1925
	1926
[5283]	1927	! * integrate liquid water equation to find condensed water *
	1928	! * and condensed water flux *
[879]	1929
	1930
	1931
[5283]	1932	! * begin downdraft loop *
[879]	1933
	1934
	1935
[5283]	1936	! * calculate detrained precipitation *
	1937
	1938	DO il = 1, ncum
	1939	IF (i<=inb(il) .AND. lwork(il)) THEN
	1940	IF (cvflag_grav) THEN
	1941	wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)
	1942	wdtraina(il, i) = wdtrain(il) / grav ! Pa 26/08/10 RomP
	1943	ELSE
	1944	wdtrain(il) = 10.0 * ep(il, i) * m(il, i) * clw(il, i)
	1945	wdtraina(il, i) = wdtrain(il) / 10. ! Pa 26/08/10 RomP
[1992]	1946	END IF
	1947	END IF
	1948	END DO
[879]	1949
[5283]	1950	IF (i>1) THEN
[879]	1951
[5283]	1952	DO j = 1, i - 1
	1953	DO il = 1, ncum
	1954	IF (i<=inb(il) .AND. lwork(il)) THEN
	1955	awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
	1956	awat = amax1(awat, 0.0)
	1957	IF (cvflag_grav) THEN
	1958	wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i)
	1959	ELSE
	1960	wdtrain(il) = wdtrain(il) + 10.0 * awat * ment(il, j, i)
	1961	END IF
	1962	END IF
	1963	END DO
	1964	END DO
	1965	DO il = 1, ncum
	1966	IF (cvflag_grav) THEN
	1967	wdtrainm(il, i) = wdtrain(il) / grav - wdtraina(il, i) ! Pm 26/08/10 RomP
	1968	ELSE
	1969	wdtrainm(il, i) = wdtrain(il) / 10. - wdtraina(il, i) ! Pm 26/08/10 RomP
	1970	END IF
	1971	END DO
[879]	1972
[5283]	1973	END IF
[879]	1974
	1975
[5283]	1976	! * find rain water and evaporation using provisional *
	1977	! * estimates of rp(i)and rp(i-1) *
[879]	1978
[5283]	1979	DO il = 1, ncum
[879]	1980
[5283]	1981	IF (i<=inb(il) .AND. lwork(il)) THEN
[879]	1982
[5283]	1983	wt(il, i) = 45.0
[879]	1984
[5283]	1985	IF (i<inb(il)) THEN
	1986	rp(il, i) = rp(il, i + 1) + (cpd * (t(il, i + 1) - t(il, &
	1987	i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i)
	1988	rp(il, i) = 0.5 * (rp(il, i) + rr(il, i))
	1989	END IF
	1990	rp(il, i) = amax1(rp(il, i), 0.0)
	1991	rp(il, i) = amin1(rp(il, i), rs(il, i))
	1992	rp(il, inb(il)) = rr(il, inb(il))
[879]	1993
[5283]	1994	IF (i==1) THEN
	1995	afac = p(il, 1) * (rs(il, 1) - rp(il, 1)) / (1.0E4 + 2000.0 * p(il, 1) * rs(il, 1))
	1996	ELSE
	1997	rp(il, i - 1) = rp(il, i) + (cpd * (t(il, i) - t(il, &
	1998	i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i)
	1999	rp(il, i - 1) = 0.5 * (rp(il, i - 1) + rr(il, i - 1))
	2000	rp(il, i - 1) = amin1(rp(il, i - 1), rs(il, i - 1))
	2001	rp(il, i - 1) = amax1(rp(il, i - 1), 0.0)
	2002	afac1 = p(il, i) * (rs(il, i) - rp(il, i)) / (1.0E4 + 2000.0 * p(il, i) * rs(il, i) &
	2003	)
	2004	afac2 = p(il, i - 1) * (rs(il, i - 1) - rp(il, i - 1)) / &
	2005	(1.0E4 + 2000.0 * p(il, i - 1) * rs(il, i - 1))
	2006	afac = 0.5 * (afac1 + afac2)
	2007	END IF
	2008	IF (i==inb(il)) afac = 0.0
	2009	afac = amax1(afac, 0.0)
	2010	bfac = 1. / (sigd * wt(il, i))
[879]	2011
[5283]	2012	! jyg1
	2013	! cc sigt=1.0
	2014	! cc IF(i.ge.icb)sigt=sigp(i)
	2015	! prise en compte de la variation progressive de sigt dans
	2016	! les couches icb et icb-1:
	2017	! pour plcl<ph(i+1), pr1=0 & pr2=1
	2018	! pour plcl>ph(i), pr1=1 & pr2=0
	2019	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
	2020	! sur le nuage, et pr2 est la proportion sous la base du
	2021	! nuage.
	2022	pr1 = (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))
	2023	pr1 = max(0., min(1., pr1))
	2024	pr2 = (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))
	2025	pr2 = max(0., min(1., pr2))
	2026	sigt = sigp(il, i) * pr1 + pr2
	2027	! jyg2
[879]	2028
[5283]	2029	b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt * afac
	2030	c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac * (ph(il, i) - ph(&
	2031	il, i + 1)) * evap(il, i + 1)
	2032	IF (c6>0.0) THEN
	2033	revap = 0.5 * (-b6 + sqrt(b6 * b6 + 4. * c6))
	2034	evap(il, i) = sigt * afac * revap
	2035	water(il, i) = revap * revap
[1992]	2036	ELSE
[5283]	2037	evap(il, i) = -evap(il, i + 1) + 0.02 * (wdtrain(il) + sigd * wt(il, i) * &
	2038	water(il, i + 1)) / (sigd * (ph(il, i) - ph(il, i + 1)))
[1992]	2039	END IF
[879]	2040
[5283]	2041	! * calculate precipitating downdraft mass flux under *
	2042	! * hydrostatic approximation *
[879]	2043
[5283]	2044	IF (i/=1) THEN
[879]	2045
[5283]	2046	tevap = amax1(0.0, evap(il, i))
	2047	delth = amax1(0.001, (th(il, i) - th(il, i - 1)))
[1992]	2048	IF (cvflag_grav) THEN
[5283]	2049	mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap * (p(il, i - 1) - p(il, i)) / &
	2050	delth
[1992]	2051	ELSE
[5283]	2052	mp(il, i) = 10. * lvcp(il, i) * sigd * tevap * (p(il, i - 1) - p(il, i)) / delth
[1992]	2053	END IF
[879]	2054
[5283]	2055	! * if hydrostatic assumption fails, *
	2056	! * solve cubic difference equation for downdraft theta *
	2057	! * and mass flux from two simultaneous differential eqns *
[879]	2058
[5283]	2059	amfac = sigd * sigd * 70.0 * ph(il, i) * (p(il, i - 1) - p(il, i)) * &
	2060	(th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))
	2061	amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) * mp(il, i))
	2062	IF (amp2>(0.1 * amfac)) THEN
	2063	xf = 100.0 * sigd * sigd * sigd * (ph(il, i) - ph(il, i + 1))
	2064	tf = b(il, i) - 5.0 * (th(il, i) - th(il, i - 1)) * t(il, i) / (lvcp(il, i) * &
	2065	sigd * th(il, i))
	2066	af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv
	2067	bf = 2. * (tinv * mp(il, i + 1))*3 + tinv mp(il, i + 1) * xf * tf + &
	2068	50. * (p(il, i - 1) - p(il, i)) * xf * tevap
	2069	fac2 = 1.0
	2070	IF (bf<0.0) fac2 = -1.0
	2071	bf = abs(bf)
	2072	ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
	2073	IF (ur>=0.0) THEN
	2074	sru = sqrt(ur)
	2075	fac = 1.0
	2076	IF ((0.5 * bf - sru)<0.0) fac = -1.0
	2077	mp(il, i) = mp(il, i + 1) * tinv + (0.5 * bf + sru)**tinv + &
	2078	fac * (abs(0.5 * bf - sru))**tinv
	2079	ELSE
	2080	d = atan(2. * sqrt(-ur) / (bf + 1.0E-28))
	2081	IF (fac2<0.0) d = 3.14159 - d
	2082	mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) * cos(d * tinv)
	2083	END IF
	2084	mp(il, i) = amax1(0.0, mp(il, i))
[879]	2085
[5283]	2086	IF (cvflag_grav) THEN
	2087	! jyg : il y a vraisemblablement une erreur dans la ligne 2
	2088	! suivante:
	2089	! il faut diviser par (mp(il,i)sigdgrav) et non par
	2090	! (mp(il,i)+sigd*0.1).
	2091	! Et il faut bien revoir les facteurs 100.
	2092	b(il, i - 1) = b(il, i) + 100.0 * (p(il, i - 1) - p(il, i)) * tevap / (mp(il, &
	2093	i) + sigd * 0.1) - 10.0 * (th(il, i) - th(il, i - 1)) * t(il, i) / (lvcp(il, i &
	2094	) * sigd * th(il, i))
	2095	ELSE
	2096	b(il, i - 1) = b(il, i) + 100.0 * (p(il, i - 1) - p(il, i)) * tevap / (mp(il, &
	2097	i) + sigd * 0.1) - 10.0 * (th(il, i) - th(il, i - 1)) * t(il, i) / (lvcp(il, i &
	2098	) * sigd * th(il, i))
	2099	END IF
	2100	b(il, i - 1) = amax1(b(il, i - 1), 0.0)
	2101	END IF
[879]	2102
[5283]	2103	! *** limit magnitude of mp(i) to meet cfl condition
	2104	! ***
[879]	2105
[5283]	2106	ampmax = 2.0 * (ph(il, i) - ph(il, i + 1)) * delti
	2107	amp2 = 2.0 * (ph(il, i - 1) - ph(il, i)) * delti
	2108	ampmax = amin1(ampmax, amp2)
	2109	mp(il, i) = amin1(mp(il, i), ampmax)
[879]	2110
[5283]	2111	! *** force mp to decrease linearly to zero
	2112	! ***
	2113	! *** between cloud base and the surface
	2114	! ***
[879]	2115
[5283]	2116	IF (p(il, i)>p(il, icb(il))) THEN
	2117	mp(il, i) = mp(il, icb(il)) * (p(il, 1) - p(il, i)) / &
	2118	(p(il, 1) - p(il, icb(il)))
	2119	END IF
[879]	2120
[5283]	2121	END IF ! i.EQ.1
[879]	2122
[5283]	2123	! * find mixing ratio of precipitating downdraft *
[879]	2124
[5283]	2125	IF (i/=inb(il)) THEN
[879]	2126
[5283]	2127	rp(il, i) = rr(il, i)
[879]	2128
[5283]	2129	IF (mp(il, i)>mp(il, i + 1)) THEN
[879]	2130
[1992]	2131	IF (cvflag_grav) THEN
[5283]	2132	rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + &
	2133	rr(il, i) * (mp(il, i) - mp(il, i + 1)) + 100. * ginv * 0.5 * sigd * (ph(il, i &
	2134	) - ph(il, i + 1)) * (evap(il, i + 1) + evap(il, i))
[1992]	2135	ELSE
[5283]	2136	rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + &
	2137	rr(il, i) * (mp(il, i) - mp(il, i + 1)) + 5. * sigd * (ph(il, i) - ph(il, i + 1 &
	2138	)) * (evap(il, i + 1) + evap(il, i))
[1992]	2139	END IF
[5283]	2140	rp(il, i) = rp(il, i) / mp(il, i)
	2141	up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) * (mp(il, i) - mp(il, i + &
	2142	1))
	2143	up(il, i) = up(il, i) / mp(il, i)
	2144	vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) * (mp(il, i) - mp(il, i + &
	2145	1))
	2146	vp(il, i) = vp(il, i) / mp(il, i)
[1742]	2147
[1992]	2148	! do j=1,ntra
[5283]	2149	! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
	2150	! testmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
	2151	! : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
	2152	! trap(il,i,j)=trap(il,i,j)/mp(il,i)
	2153	! END DO
[1742]	2154
[5283]	2155	ELSE
	2156
	2157	IF (mp(il, i + 1)>1.0E-16) THEN
	2158	IF (cvflag_grav) THEN
	2159	rp(il, i) = rp(il, i + 1) + 100. * ginv * 0.5 * sigd * (ph(il, i) - ph(il, &
	2160	i + 1)) * (evap(il, i + 1) + evap(il, i)) / mp(il, i + 1)
	2161	ELSE
	2162	rp(il, i) = rp(il, i + 1) + 5. * sigd * (ph(il, i) - ph(il, i + 1)) * (evap &
	2163	(il, i + 1) + evap(il, i)) / mp(il, i + 1)
	2164	END IF
	2165	up(il, i) = up(il, i + 1)
	2166	vp(il, i) = vp(il, i + 1)
	2167
	2168	! do j=1,ntra
	2169	! trap(il,i,j)=trap(il,i+1,j)
	2170	! END DO
	2171
	2172	END IF
[1992]	2173	END IF
[5283]	2174	rp(il, i) = amin1(rp(il, i), rs(il, i))
	2175	rp(il, i) = amax1(rp(il, i), 0.0)
	2176
[1992]	2177	END IF
	2178	END IF
[5283]	2179	END DO
[1742]	2180
[5283]	2181	400 END DO
[879]	2182
[5283]	2183	END SUBROUTINE cv30_unsat
[879]	2184
[5283]	2185	SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, &
	2186	tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, &
	2187	wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, &
	2188	tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, &
	2189	mike, tls, tps, qcondc, wd)
	2190	USE conema3_mod_h
	2191	USE cvflag_mod_h
[5346]	2192	USE lmdz_cv_ini, ONLY : cl,cpd,rrv,rrd,rowl,grav,cpv
[879]	2193
[5283]	2194	IMPLICIT NONE
[879]	2195
[5283]	2196	! inputs:
	2197	INTEGER ncum, nd, na, ntra, nloc
	2198	INTEGER icb(nloc), inb(nloc)
	2199	REAL delt
	2200	REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
	2201	REAL tra(nloc, nd, ntra), sig(nloc, nd)
	2202	REAL gz(nloc, na), ph(nloc, nd + 1), h(nloc, na), hp(nloc, na)
	2203	REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
	2204	REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
	2205	REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
	2206	REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
	2207	REAL water(nloc, na), evap(nloc, na), b(nloc, na)
	2208	REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
	2209	! ym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
	2210	REAL vent(nloc, na, na), elij(nloc, na, na)
	2211	INTEGER nent(nloc, na)
	2212	REAL traent(nloc, na, na, ntra)
	2213	REAL tv(nloc, nd), tvp(nloc, nd)
[879]	2214
[5283]	2215	! input/output:
	2216	INTEGER iflag(nloc)
[879]	2217
[5283]	2218	! outputs:
	2219	REAL precip(nloc)
	2220	REAL vprecip(nloc, nd + 1)
	2221	REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	2222	REAL ftra(nloc, nd, ntra)
	2223	REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
	2224	REAL dnwd0(nloc, nd), mike(nloc, nd)
	2225	REAL tls(nloc, nd), tps(nloc, nd)
	2226	REAL qcondc(nloc, nd) ! cld
	2227	REAL wd(nloc) ! gust
[879]	2228
[5283]	2229	! local variables:
	2230	INTEGER i, k, il, n, j, num1
	2231	REAL rat, awat, delti
	2232	REAL ax, bx, cx, dx, ex
	2233	REAL cpinv, rdcp, dpinv
	2234	REAL lvcp(nloc, na), mke(nloc, na)
	2235	REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
	2236	! !! real up1(nloc), dn1(nloc)
	2237	REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
	2238	REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
	2239	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
	2240	REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
[879]	2241
	2242
[5283]	2243	! -------------------------------------------------------------
[879]	2244
[5283]	2245	! initialization:
[879]	2246
[5283]	2247	delti = 1.0 / delt
[879]	2248
[1992]	2249	DO il = 1, ncum
[5283]	2250	precip(il) = 0.0
	2251	wd(il) = 0.0 ! gust
	2252	vprecip(il, nd + 1) = 0.
[1992]	2253	END DO
	2254
[5283]	2255	DO i = 1, nd
	2256	DO il = 1, ncum
	2257	vprecip(il, i) = 0.0
	2258	ft(il, i) = 0.0
	2259	fr(il, i) = 0.0
	2260	fu(il, i) = 0.0
	2261	fv(il, i) = 0.0
	2262	qcondc(il, i) = 0.0 ! cld
	2263	qcond(il, i) = 0.0 ! cld
	2264	nqcond(il, i) = 0.0 ! cld
	2265	END DO
[1992]	2266	END DO
	2267
[5283]	2268	! do j=1,ntra
	2269	! do i=1,nd
	2270	! do il=1,ncum
	2271	! ftra(il,i,j)=0.0
	2272	! enddo
	2273	! enddo
	2274	! enddo
[1992]	2275
[5283]	2276	DO i = 1, nl
	2277	DO il = 1, ncum
	2278	lvcp(il, i) = lv(il, i) / cpn(il, i)
	2279	END DO
	2280	END DO
[1992]	2281
	2282
	2283
[5283]	2284	! * calculate surface precipitation in mm/day *
[1992]	2285
	2286	DO il = 1, ncum
[5283]	2287	IF (ep(il, inb(il))>=0.0001) THEN
[1992]	2288	IF (cvflag_grav) THEN
[5283]	2289	precip(il) = wt(il, 1) * sigd * water(il, 1) * 86400. * 1000. / (rowl * grav)
[1992]	2290	ELSE
[5283]	2291	precip(il) = wt(il, 1) * sigd * water(il, 1) * 8640.
[1992]	2292	END IF
	2293	END IF
	2294	END DO
	2295
[5283]	2296	! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
[1992]	2297
[5283]	2298	! MAF rajout pour lessivage
	2299	DO k = 1, nl
	2300	DO il = 1, ncum
	2301	IF (k<=inb(il)) THEN
	2302	IF (cvflag_grav) THEN
	2303	vprecip(il, k) = wt(il, k) * sigd * water(il, k) / grav
	2304	ELSE
	2305	vprecip(il, k) = wt(il, k) * sigd * water(il, k) / 10.
	2306	END IF
	2307	END IF
	2308	END DO
	2309	END DO
[1992]	2310
	2311
[5283]	2312	! * Calculate downdraft velocity scale *
	2313	! * NE PAS UTILISER POUR L'INSTANT *
[1992]	2314
[5283]	2315	! do il=1,ncum
	2316	! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
	2317	! : /(sigd*p(il,icb(il)))
	2318	! enddo
[1992]	2319
	2320
[5283]	2321	! * calculate tendencies of lowest level potential temperature *
	2322	! * and mixing ratio *
	2323
[1992]	2324	DO il = 1, ncum
[5283]	2325	work(il) = 1.0 / (ph(il, 1) - ph(il, 2))
	2326	am(il) = 0.0
[1992]	2327	END DO
	2328
[5283]	2329	DO k = 2, nl
	2330	DO il = 1, ncum
	2331	IF (k<=inb(il)) THEN
	2332	am(il) = am(il) + m(il, k)
	2333	END IF
	2334	END DO
	2335	END DO
[1992]	2336
[5283]	2337	DO il = 1, ncum
[1992]	2338
[5283]	2339	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
	2340	IF (cvflag_grav) THEN
	2341	IF ((0.01 * grav * work(il) * am(il))>=delti) iflag(il) = 1 !consist vect
	2342	ft(il, 1) = 0.01 * grav * work(il) * am(il) * (t(il, 2) - t(il, 1) + (gz(il, 2) - gz(il, &
	2343	1)) / cpn(il, 1))
	2344	ELSE
	2345	IF ((0.1 * work(il) * am(il))>=delti) iflag(il) = 1 !consistency vect
	2346	ft(il, 1) = 0.1 * work(il) * am(il) * (t(il, 2) - t(il, 1) + (gz(il, 2) - gz(il, &
	2347	1)) / cpn(il, 1))
	2348	END IF
[1992]	2349
[5283]	2350	ft(il, 1) = ft(il, 1) - 0.5 * lvcp(il, 1) * sigd * (evap(il, 1) + evap(il, 2))
[1992]	2351
[5283]	2352	IF (cvflag_grav) THEN
	2353	ft(il, 1) = ft(il, 1) - 0.009 * grav * sigd * mp(il, 2) * t(il, 1) * b(il, 1) * &
	2354	work(il)
	2355	ELSE
	2356	ft(il, 1) = ft(il, 1) - 0.09 * sigd * mp(il, 2) * t(il, 1) * b(il, 1) * work(il)
	2357	END IF
[1992]	2358
[5283]	2359	ft(il, 1) = ft(il, 1) + 0.01 * sigd * wt(il, 1) * (cl - cpd) * water(il, 2) * (t(il, 2 &
	2360	) - t(il, 1)) * work(il) / cpn(il, 1)
[1992]	2361
[5283]	2362	IF (cvflag_grav) THEN
	2363	! jyg1 Correction pour mieux conserver l'eau (conformite avec
	2364	! CONVECT4.3)
	2365	! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas
	2366	! evap)
	2367	fr(il, 1) = 0.01 * grav * mp(il, 2) * (rp(il, 2) - rr(il, 1)) * work(il) + &
	2368	sigd * 0.5 * (evap(il, 1) + evap(il, 2))
	2369	! +tard : +sigd*evap(il,1)
[1992]	2370
[5283]	2371	fr(il, 1) = fr(il, 1) + 0.01 * grav * am(il) * (rr(il, 2) - rr(il, 1)) * work(il)
[1992]	2372
[5283]	2373	fu(il, 1) = fu(il, 1) + 0.01 * grav * work(il) * (mp(il, 2) * (up(il, 2) - u(il, &
	2374	1)) + am(il) * (u(il, 2) - u(il, 1)))
	2375	fv(il, 1) = fv(il, 1) + 0.01 * grav * work(il) * (mp(il, 2) * (vp(il, 2) - v(il, &
	2376	1)) + am(il) * (v(il, 2) - v(il, 1)))
	2377	ELSE ! cvflag_grav
	2378	fr(il, 1) = 0.1 * mp(il, 2) * (rp(il, 2) - rr(il, 1)) * work(il) + &
	2379	sigd * 0.5 * (evap(il, 1) + evap(il, 2))
	2380	fr(il, 1) = fr(il, 1) + 0.1 * am(il) * (rr(il, 2) - rr(il, 1)) * work(il)
	2381	fu(il, 1) = fu(il, 1) + 0.1 * work(il) * (mp(il, 2) * (up(il, 2) - u(il, &
	2382	1)) + am(il) * (u(il, 2) - u(il, 1)))
	2383	fv(il, 1) = fv(il, 1) + 0.1 * work(il) * (mp(il, 2) * (vp(il, 2) - v(il, &
	2384	1)) + am(il) * (v(il, 2) - v(il, 1)))
	2385	END IF ! cvflag_grav
[1992]	2386
[5283]	2387	END DO ! il
[1992]	2388
[5283]	2389	! do j=1,ntra
	2390	! do il=1,ncum
	2391	! if (cvflag_grav) THEN
	2392	! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
	2393	! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2394	! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2395	! else
	2396	! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
	2397	! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2398	! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2399	! END IF
	2400	! enddo
	2401	! enddo
	2402
	2403	DO j = 2, nl
	2404	DO il = 1, ncum
	2405	IF (j<=inb(il)) THEN
	2406	IF (cvflag_grav) THEN
	2407	fr(il, 1) = fr(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) * (qent(il, &
	2408	j, 1) - rr(il, 1))
	2409	fu(il, 1) = fu(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) * (uent(il, &
	2410	j, 1) - u(il, 1))
	2411	fv(il, 1) = fv(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) * (vent(il, &
	2412	j, 1) - v(il, 1))
	2413	ELSE ! cvflag_grav
	2414	fr(il, 1) = fr(il, 1) + 0.1 * work(il) * ment(il, j, 1) * (qent(il, j, 1) - &
	2415	rr(il, 1))
	2416	fu(il, 1) = fu(il, 1) + 0.1 * work(il) * ment(il, j, 1) * (uent(il, j, 1) - u &
	2417	(il, 1))
	2418	fv(il, 1) = fv(il, 1) + 0.1 * work(il) * ment(il, j, 1) * (vent(il, j, 1) - v &
	2419	(il, 1))
	2420	END IF ! cvflag_grav
	2421	END IF ! j
	2422	END DO
[1992]	2423	END DO
	2424
[5283]	2425	! do k=1,ntra
	2426	! do j=2,nl
	2427	! do il=1,ncum
	2428	! if (j.le.inb(il)) THEN
	2429	! if (cvflag_grav) THEN
	2430	! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
	2431	! : *(traent(il,j,1,k)-tra(il,1,k))
	2432	! else
	2433	! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
	2434	! : *(traent(il,j,1,k)-tra(il,1,k))
	2435	! END IF
[1992]	2436
[5283]	2437	! END IF
	2438	! enddo
	2439	! enddo
	2440	! enddo
[1992]	2441
	2442
[5283]	2443	! * calculate tendencies of potential temperature and mixing ratio *
	2444	! * at levels above the lowest level *
[1992]	2445
[5283]	2446	! * first find the net saturated updraft and downdraft mass fluxes *
	2447	! * through each level *
[1992]	2448
[5283]	2449	DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
[1992]	2450
[5283]	2451	num1 = 0
[1992]	2452	DO il = 1, ncum
[5283]	2453	IF (i<=inb(il)) num1 = num1 + 1
[1992]	2454	END DO
[5283]	2455	IF (num1<=0) GO TO 500
[1992]	2456
[5283]	2457	CALL zilch(amp1, ncum)
	2458	CALL zilch(ad, ncum)
	2459
	2460	DO k = i + 1, nl + 1
[1992]	2461	DO il = 1, ncum
[5283]	2462	IF (i<=inb(il) .AND. k<=(inb(il) + 1)) THEN
	2463	amp1(il) = amp1(il) + m(il, k)
[1992]	2464	END IF
	2465	END DO
	2466	END DO
	2467
[5283]	2468	DO k = 1, i
	2469	DO j = i + 1, nl + 1
	2470	DO il = 1, ncum
	2471	IF (i<=inb(il) .AND. j<=(inb(il) + 1)) THEN
	2472	amp1(il) = amp1(il) + ment(il, k, j)
	2473	END IF
	2474	END DO
[1992]	2475	END DO
	2476	END DO
	2477
[5283]	2478	DO k = 1, i - 1
	2479	DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
	2480	DO il = 1, ncum
	2481	IF (i<=inb(il) .AND. j<=inb(il)) THEN
	2482	ad(il) = ad(il) + ment(il, j, k)
	2483	END IF
	2484	END DO
	2485	END DO
	2486	END DO
[1992]	2487
	2488	DO il = 1, ncum
	2489	IF (i<=inb(il)) THEN
[5283]	2490	dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
	2491	cpinv = 1.0 / cpn(il, i)
[1992]	2492
[5283]	2493	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
	2494	IF (cvflag_grav) THEN
	2495	IF ((0.01 * grav * dpinv * amp1(il))>=delti) iflag(il) = 1 ! vecto
	2496	ELSE
	2497	IF ((0.1 * dpinv * amp1(il))>=delti) iflag(il) = 1 ! vecto
	2498	END IF
[1992]	2499
	2500	IF (cvflag_grav) THEN
[5283]	2501	ft(il, i) = 0.01 * grav * dpinv * (amp1(il) * (t(il, i + 1) - t(il, &
	2502	i) + (gz(il, i + 1) - gz(il, i)) * cpinv) - ad(il) * (t(il, i) - t(il, &
	2503	i - 1) + (gz(il, i) - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) * (evap(&
	2504	il, i) + evap(il, i + 1))
	2505	rat = cpn(il, i - 1) * cpinv
	2506	ft(il, i) = ft(il, i) - 0.009 * grav * sigd * (mp(il, i + 1) * t(il, i) * b(il, i) &
	2507	- mp(il, i) * t(il, i - 1) * rat * b(il, i - 1)) * dpinv
	2508	ft(il, i) = ft(il, i) + 0.01 * grav * dpinv * ment(il, i, i) * (hp(il, i) - h(&
	2509	il, i) + t(il, i) * (cpv - cpd) * (rr(il, i) - qent(il, i, i))) * cpinv
[1992]	2510	ELSE ! cvflag_grav
[5283]	2511	ft(il, i) = 0.1 * dpinv * (amp1(il) * (t(il, i + 1) - t(il, &
	2512	i) + (gz(il, i + 1) - gz(il, i)) * cpinv) - ad(il) * (t(il, i) - t(il, &
	2513	i - 1) + (gz(il, i) - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) * (evap(&
	2514	il, i) + evap(il, i + 1))
	2515	rat = cpn(il, i - 1) * cpinv
	2516	ft(il, i) = ft(il, i) - 0.09 * sigd * (mp(il, i + 1) * t(il, i) * b(il, i) - mp(il &
	2517	, i) * t(il, i - 1) * rat * b(il, i - 1)) * dpinv
	2518	ft(il, i) = ft(il, i) + 0.1 * dpinv * ment(il, i, i) * (hp(il, i) - h(il, i) + &
	2519	t(il, i) * (cpv - cpd) * (rr(il, i) - qent(il, i, i))) * cpinv
[1992]	2520	END IF ! cvflag_grav
	2521
[5283]	2522	ft(il, i) = ft(il, i) + 0.01 * sigd * wt(il, i) * (cl - cpd) * water(il, i + 1) * (&
	2523	t(il, i + 1) - t(il, i)) * dpinv * cpinv
[1992]	2524
	2525	IF (cvflag_grav) THEN
[5283]	2526	fr(il, i) = 0.01 * grav * dpinv * (amp1(il) * (rr(il, i + 1) - rr(il, &
	2527	i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
	2528	fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * (amp1(il) * (u(il, i + 1) - u(il, &
	2529	i)) - ad(il) * (u(il, i) - u(il, i - 1)))
	2530	fv(il, i) = fv(il, i) + 0.01 * grav * dpinv * (amp1(il) * (v(il, i + 1) - v(il, &
	2531	i)) - ad(il) * (v(il, i) - v(il, i - 1)))
[1992]	2532	ELSE ! cvflag_grav
[5283]	2533	fr(il, i) = 0.1 * dpinv * (amp1(il) * (rr(il, i + 1) - rr(il, &
	2534	i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
	2535	fu(il, i) = fu(il, i) + 0.1 * dpinv * (amp1(il) * (u(il, i + 1) - u(il, &
	2536	i)) - ad(il) * (u(il, i) - u(il, i - 1)))
	2537	fv(il, i) = fv(il, i) + 0.1 * dpinv * (amp1(il) * (v(il, i + 1) - v(il, &
	2538	i)) - ad(il) * (v(il, i) - v(il, i - 1)))
[1992]	2539	END IF ! cvflag_grav
[5283]	2540
	2541	END IF ! i
[1992]	2542	END DO
	2543
[5283]	2544	! do k=1,ntra
	2545	! do il=1,ncum
	2546	! if (i.le.inb(il)) THEN
	2547	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2548	! cpinv=1.0/cpn(il,i)
	2549	! if (cvflag_grav) THEN
	2550	! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
	2551	! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2552	! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2553	! else
	2554	! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
	2555	! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2556	! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2557	! END IF
	2558	! END IF
	2559	! enddo
	2560	! enddo
[1992]	2561
[5283]	2562	DO k = 1, i - 1
	2563	DO il = 1, ncum
	2564	IF (i<=inb(il)) THEN
	2565	dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
	2566	cpinv = 1.0 / cpn(il, i)
[1992]	2567
[5283]	2568	awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i)
	2569	awat = amax1(awat, 0.0)
[1992]	2570
[5283]	2571	IF (cvflag_grav) THEN
	2572	fr(il, i) = fr(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (qent(il, k &
	2573	, i) - awat - rr(il, i))
	2574	fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (uent(il, k &
	2575	, i) - u(il, i))
	2576	fv(il, i) = fv(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (vent(il, k &
	2577	, i) - v(il, i))
	2578	ELSE ! cvflag_grav
	2579	fr(il, i) = fr(il, i) + 0.1 * dpinv * ment(il, k, i) * (qent(il, k, i) - &
	2580	awat - rr(il, i))
	2581	fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (uent(il, k &
	2582	, i) - u(il, i))
	2583	fv(il, i) = fv(il, i) + 0.1 * dpinv * ment(il, k, i) * (vent(il, k, i) - v(&
	2584	il, i))
	2585	END IF ! cvflag_grav
[1992]	2586
[5283]	2587	! (saturated updrafts resulting from mixing) ! cld
	2588	qcond(il, i) = qcond(il, i) + (elij(il, k, i) - awat) ! cld
	2589	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2590	END IF ! i
	2591	END DO
	2592	END DO
[1992]	2593
[5283]	2594	! do j=1,ntra
	2595	! do k=1,i-1
	2596	! do il=1,ncum
	2597	! if (i.le.inb(il)) THEN
	2598	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2599	! cpinv=1.0/cpn(il,i)
	2600	! if (cvflag_grav) THEN
	2601	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2602	! : *(traent(il,k,i,j)-tra(il,i,j))
	2603	! else
	2604	! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2605	! : *(traent(il,k,i,j)-tra(il,i,j))
	2606	! END IF
	2607	! END IF
	2608	! enddo
	2609	! enddo
	2610	! enddo
[1992]	2611
[5283]	2612	DO k = i, nl + 1
	2613	DO il = 1, ncum
	2614	IF (i<=inb(il) .AND. k<=inb(il)) THEN
	2615	dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
	2616	cpinv = 1.0 / cpn(il, i)
[1992]	2617
[5283]	2618	IF (cvflag_grav) THEN
	2619	fr(il, i) = fr(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (qent(il, k &
	2620	, i) - rr(il, i))
	2621	fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (uent(il, k &
	2622	, i) - u(il, i))
	2623	fv(il, i) = fv(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (vent(il, k &
	2624	, i) - v(il, i))
	2625	ELSE ! cvflag_grav
	2626	fr(il, i) = fr(il, i) + 0.1 * dpinv * ment(il, k, i) * (qent(il, k, i) - rr &
	2627	(il, i))
	2628	fu(il, i) = fu(il, i) + 0.1 * dpinv * ment(il, k, i) * (uent(il, k, i) - u(&
	2629	il, i))
	2630	fv(il, i) = fv(il, i) + 0.1 * dpinv * ment(il, k, i) * (vent(il, k, i) - v(&
	2631	il, i))
	2632	END IF ! cvflag_grav
	2633	END IF ! i and k
	2634	END DO
	2635	END DO
[1992]	2636
[5283]	2637	! do j=1,ntra
	2638	! do k=i,nl+1
	2639	! do il=1,ncum
	2640	! if (i.le.inb(il) .AND. k.le.inb(il)) THEN
	2641	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2642	! cpinv=1.0/cpn(il,i)
	2643	! if (cvflag_grav) THEN
	2644	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2645	! : *(traent(il,k,i,j)-tra(il,i,j))
	2646	! else
	2647	! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2648	! : *(traent(il,k,i,j)-tra(il,i,j))
	2649	! END IF
	2650	! END IF ! i and k
	2651	! enddo
	2652	! enddo
	2653	! enddo
[1992]	2654
[5283]	2655	DO il = 1, ncum
	2656	IF (i<=inb(il)) THEN
	2657	dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
	2658	cpinv = 1.0 / cpn(il, i)
[1992]	2659
[5283]	2660	IF (cvflag_grav) THEN
	2661	! sb: on ne fait pas encore la correction permettant de mieux
	2662	! conserver l'eau:
	2663	fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) + evap(il, i + 1)) + &
	2664	0.01 * grav * (mp(il, i + 1) * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, &
	2665	i) - rr(il, i - 1))) * dpinv
[1992]	2666
[5283]	2667	fu(il, i) = fu(il, i) + 0.01 * grav * (mp(il, i + 1) * (up(il, i + 1) - u(il, &
	2668	i)) - mp(il, i) * (up(il, i) - u(il, i - 1))) * dpinv
	2669	fv(il, i) = fv(il, i) + 0.01 * grav * (mp(il, i + 1) * (vp(il, i + 1) - v(il, &
	2670	i)) - mp(il, i) * (vp(il, i) - v(il, i - 1))) * dpinv
	2671	ELSE ! cvflag_grav
	2672	fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) + evap(il, i + 1)) + &
	2673	0.1 * (mp(il, i + 1) * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, i) - rr(il, &
	2674	i - 1))) * dpinv
	2675	fu(il, i) = fu(il, i) + 0.1 * (mp(il, i + 1) * (up(il, i + 1) - u(il, &
	2676	i)) - mp(il, i) * (up(il, i) - u(il, i - 1))) * dpinv
	2677	fv(il, i) = fv(il, i) + 0.1 * (mp(il, i + 1) * (vp(il, i + 1) - v(il, &
	2678	i)) - mp(il, i) * (vp(il, i) - v(il, i - 1))) * dpinv
	2679	END IF ! cvflag_grav
[1992]	2680
[5283]	2681	END IF ! i
	2682	END DO
[1992]	2683
[5283]	2684	! sb: interface with the cloud parameterization: ! cld
[1992]	2685
[5283]	2686	DO k = i + 1, nl
	2687	DO il = 1, ncum
	2688	IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld
	2689	! (saturated downdrafts resulting from mixing) ! cld
	2690	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
	2691	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2692	END IF ! cld
	2693	END DO ! cld
	2694	END DO ! cld
[1992]	2695
[5283]	2696	! (particular case: no detraining level is found) ! cld
	2697	DO il = 1, ncum ! cld
	2698	IF (i<=inb(il) .AND. nent(il, i)==0) THEN ! cld
	2699	qcond(il, i) = qcond(il, i) + (1. - ep(il, i)) * clw(il, i) ! cld
	2700	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2701	END IF ! cld
	2702	END DO ! cld
[1992]	2703
[5283]	2704	DO il = 1, ncum ! cld
	2705	IF (i<=inb(il) .AND. nqcond(il, i)/=0.) THEN ! cld
	2706	qcond(il, i) = qcond(il, i) / nqcond(il, i) ! cld
	2707	END IF ! cld
	2708	END DO
[1992]	2709
[5283]	2710	! do j=1,ntra
	2711	! do il=1,ncum
	2712	! if (i.le.inb(il)) THEN
	2713	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2714	! cpinv=1.0/cpn(il,i)
[1992]	2715
[5283]	2716	! if (cvflag_grav) THEN
	2717	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
	2718	! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2719	! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
	2720	! else
	2721	! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
	2722	! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2723	! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
	2724	! END IF
	2725	! END IF ! i
	2726	! enddo
	2727	! enddo
[1992]	2728
[5283]	2729	500 END DO
[1992]	2730
	2731
[5283]	2732	! * move the detrainment at level inb down to level inb-1 *
	2733	! * in such a way as to preserve the vertically *
	2734	! * integrated enthalpy and water tendencies *
[1992]	2735
[5283]	2736	DO il = 1, ncum
[1992]	2737
[5283]	2738	ax = 0.1 * ment(il, inb(il), inb(il)) * (hp(il, inb(il)) - h(il, inb(il)) + t(il, &
	2739	inb(il)) * (cpv - cpd) * (rr(il, inb(il)) - qent(il, inb(il), &
	2740	inb(il)))) / (cpn(il, inb(il)) * (ph(il, inb(il)) - ph(il, inb(il) + 1)))
	2741	ft(il, inb(il)) = ft(il, inb(il)) - ax
	2742	ft(il, inb(il) - 1) = ft(il, inb(il) - 1) + ax * cpn(il, inb(il)) * (ph(il, inb(il &
	2743	)) - ph(il, inb(il) + 1)) / (cpn(il, inb(il) - 1) * (ph(il, inb(il) - 1) - ph(il, &
	2744	inb(il))))
[1992]	2745
[5283]	2746	bx = 0.1 * ment(il, inb(il), inb(il)) * (qent(il, inb(il), inb(il)) - rr(il, inb(&
	2747	il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
	2748	fr(il, inb(il)) = fr(il, inb(il)) - bx
	2749	fr(il, inb(il) - 1) = fr(il, inb(il) - 1) + bx * (ph(il, inb(il)) - ph(il, inb(il) + &
	2750	1)) / (ph(il, inb(il) - 1) - ph(il, inb(il)))
[1992]	2751
[5283]	2752	cx = 0.1 * ment(il, inb(il), inb(il)) * (uent(il, inb(il), inb(il)) - u(il, inb(il &
	2753	))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
	2754	fu(il, inb(il)) = fu(il, inb(il)) - cx
	2755	fu(il, inb(il) - 1) = fu(il, inb(il) - 1) + cx * (ph(il, inb(il)) - ph(il, inb(il) + &
	2756	1)) / (ph(il, inb(il) - 1) - ph(il, inb(il)))
[1992]	2757
[5283]	2758	dx = 0.1 * ment(il, inb(il), inb(il)) * (vent(il, inb(il), inb(il)) - v(il, inb(il &
	2759	))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
	2760	fv(il, inb(il)) = fv(il, inb(il)) - dx
	2761	fv(il, inb(il) - 1) = fv(il, inb(il) - 1) + dx * (ph(il, inb(il)) - ph(il, inb(il) + &
	2762	1)) / (ph(il, inb(il) - 1) - ph(il, inb(il)))
	2763
[1992]	2764	END DO
	2765
[5283]	2766	! do j=1,ntra
	2767	! do il=1,ncum
	2768	! ex=0.1*ment(il,inb(il),inb(il))
	2769	! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
	2770	! : /(ph(il,inb(il))-ph(il,inb(il)+1))
	2771	! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
	2772	! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
	2773	! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
	2774	! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2775	! enddo
	2776	! enddo
[1992]	2777
	2778
[5283]	2779	! * homoginize tendencies below cloud base *
[1992]	2780
	2781	DO il = 1, ncum
[5283]	2782	asum(il) = 0.0
	2783	bsum(il) = 0.0
	2784	csum(il) = 0.0
	2785	dsum(il) = 0.0
[1992]	2786	END DO
	2787
[5283]	2788	DO i = 1, nl
	2789	DO il = 1, ncum
	2790	IF (i<=(icb(il) - 1)) THEN
	2791	asum(il) = asum(il) + ft(il, i) * (ph(il, i) - ph(il, i + 1))
	2792	bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (cl - cpd) * (t(il, i) - t(il, &
	2793	1))) * (ph(il, i) - ph(il, i + 1))
	2794	csum(il) = csum(il) + (lv(il, i) + (cl - cpd) * (t(il, i) - t(il, &
	2795	1))) * (ph(il, i) - ph(il, i + 1))
	2796	dsum(il) = dsum(il) + t(il, i) * (ph(il, i) - ph(il, i + 1)) / th(il, i)
	2797	END IF
	2798	END DO
[1992]	2799	END DO
[5283]	2800
	2801	! !!! do 700 i=1,icb(il)-1
	2802	DO i = 1, nl
	2803	DO il = 1, ncum
	2804	IF (i<=(icb(il) - 1)) THEN
	2805	ft(il, i) = asum(il) * t(il, i) / (th(il, i) * dsum(il))
	2806	fr(il, i) = bsum(il) / csum(il)
	2807	END IF
	2808	END DO
[1992]	2809	END DO
	2810
	2811
[5283]	2812	! * reset counter and return *
	2813
[1992]	2814	DO il = 1, ncum
[5283]	2815	sig(il, nd) = 2.0
	2816	END DO
	2817
	2818	DO i = 1, nd
	2819	DO il = 1, ncum
[1992]	2820	upwd(il, i) = 0.0
	2821	dnwd(il, i) = 0.0
[5283]	2822	END DO
[1992]	2823	END DO
	2824
[5283]	2825	DO i = 1, nl
[1992]	2826	DO il = 1, ncum
[5283]	2827	dnwd0(il, i) = -mp(il, i)
[1992]	2828	END DO
	2829	END DO
[5283]	2830	DO i = nl + 1, nd
	2831	DO il = 1, ncum
	2832	dnwd0(il, i) = 0.
	2833	END DO
	2834	END DO
[1992]	2835
[5283]	2836	DO i = 1, nl
	2837	DO il = 1, ncum
	2838	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	2839	upwd(il, i) = 0.0
	2840	dnwd(il, i) = 0.0
	2841	END IF
	2842	END DO
	2843	END DO
	2844
	2845	DO i = 1, nl
	2846	DO k = 1, nl
[1992]	2847	DO il = 1, ncum
[5283]	2848	up1(il, k, i) = 0.0
	2849	dn1(il, k, i) = 0.0
[1992]	2850	END DO
	2851	END DO
	2852	END DO
	2853
[5283]	2854	DO i = 1, nl
	2855	DO k = i, nl
	2856	DO n = 1, i - 1
	2857	DO il = 1, ncum
	2858	IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
	2859	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
	2860	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
	2861	END IF
	2862	END DO
	2863	END DO
[1992]	2864	END DO
	2865	END DO
	2866
[5283]	2867	DO i = 2, nl
	2868	DO k = i, nl
	2869	DO il = 1, ncum
	2870	! test if (i.ge.icb(il).AND.i.le.inb(il).AND.k.le.inb(il))
	2871	! THEN
	2872	IF (i<=inb(il) .AND. k<=inb(il)) THEN
	2873	upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
	2874	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
	2875	END IF
	2876	END DO
	2877	END DO
	2878	END DO
[1992]	2879
	2880
[5283]	2881	! !!! DO il=1,ncum
	2882	! !!! do i=icb(il),inb(il)
	2883	! !!!
	2884	! !!! upwd(il,i)=0.0
	2885	! !!! dnwd(il,i)=0.0
	2886	! !!! do k=i,inb(il)
	2887	! !!! up1=0.0
	2888	! !!! dn1=0.0
	2889	! !!! do n=1,i-1
	2890	! !!! up1=up1+ment(il,n,k)
	2891	! !!! dn1=dn1-ment(il,k,n)
	2892	! !!! enddo
	2893	! !!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
	2894	! !!! dnwd(il,i)=dnwd(il,i)+dn1
	2895	! !!! enddo
	2896	! !!! enddo
	2897	! !!!
	2898	! !!! ENDDO
[1992]	2899
[5283]	2900	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2901	! determination de la variation de flux ascendant entre
	2902	! deux niveau non dilue mike
	2903	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1992]	2904
[5283]	2905	DO i = 1, nl
	2906	DO il = 1, ncum
	2907	mike(il, i) = m(il, i)
	2908	END DO
[1992]	2909	END DO
	2910
[5283]	2911	DO i = nl + 1, nd
	2912	DO il = 1, ncum
	2913	mike(il, i) = 0.
	2914	END DO
[1992]	2915	END DO
	2916
[5283]	2917	DO i = 1, nd
[1992]	2918	DO il = 1, ncum
[5283]	2919	ma(il, i) = 0
[1992]	2920	END DO
	2921	END DO
	2922
[5283]	2923	DO i = 1, nl
	2924	DO j = i, nl
	2925	DO il = 1, ncum
	2926	ma(il, i) = ma(il, i) + m(il, j)
	2927	END DO
	2928	END DO
[1992]	2929	END DO
	2930
[5283]	2931	DO i = nl + 1, nd
	2932	DO il = 1, ncum
	2933	ma(il, i) = 0.
	2934	END DO
[1992]	2935	END DO
	2936
[5283]	2937	DO i = 1, nl
	2938	DO il = 1, ncum
	2939	IF (i<=(icb(il) - 1)) THEN
	2940	ma(il, i) = 0
	2941	END IF
	2942	END DO
	2943	END DO
[1992]	2944
[5283]	2945	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2946	! icb represente de niveau ou se trouve la
	2947	! base du nuage , et inb le top du nuage
	2948	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2949
	2950	DO i = 1, nd
	2951	DO il = 1, ncum
	2952	mke(il, i) = upwd(il, i) + dnwd(il, i)
	2953	END DO
[1992]	2954	END DO
	2955
[5283]	2956	DO i = 1, nd
	2957	DO il = 1, ncum
	2958	rdcp = (rrd * (1. - rr(il, i)) - rr(il, i) * rrv) / (cpd * (1. - rr(il, &
	2959	i)) + rr(il, i) * cpv)
	2960	tls(il, i) = t(il, i) * (1000.0 / p(il, i))**rdcp
	2961	tps(il, i) = tp(il, i)
	2962	END DO
[1992]	2963	END DO
	2964
	2965
[5283]	2966	! * diagnose the in-cloud mixing ratio * ! cld
	2967	! * of condensed water * ! cld
	2968	! cld
[1992]	2969
[5283]	2970	DO i = 1, nd ! cld
	2971	DO il = 1, ncum ! cld
	2972	mac(il, i) = 0.0 ! cld
	2973	wa(il, i) = 0.0 ! cld
	2974	siga(il, i) = 0.0 ! cld
	2975	sax(il, i) = 0.0 ! cld
	2976	END DO ! cld
[1992]	2977	END DO ! cld
	2978
[5283]	2979	DO i = minorig, nl ! cld
	2980	DO k = i + 1, nl + 1 ! cld
	2981	DO il = 1, ncum ! cld
	2982	IF (i<=inb(il) .AND. k<=(inb(il) + 1)) THEN ! cld
	2983	mac(il, i) = mac(il, i) + m(il, k) ! cld
	2984	END IF ! cld
	2985	END DO ! cld
[1992]	2986	END DO ! cld
	2987	END DO ! cld
	2988
[5283]	2989	DO i = 1, nl ! cld
	2990	DO j = 1, i ! cld
	2991	DO il = 1, ncum ! cld
	2992	IF (i>=icb(il) .AND. i<=(inb(il) - 1) & ! cld
	2993	.AND. j>=icb(il)) THEN ! cld
	2994	sax(il, i) = sax(il, i) + rrd * (tvp(il, j) - tv(il, j)) & ! cld
	2995	* (ph(il, j) - ph(il, j + 1)) / p(il, j) ! cld
	2996	END IF ! cld
	2997	END DO ! cld
	2998	END DO ! cld
	2999	END DO ! cld
	3000
	3001	DO i = 1, nl ! cld
[1992]	3002	DO il = 1, ncum ! cld
[5283]	3003	IF (i>=icb(il) .AND. i<=(inb(il) - 1) & ! cld
	3004	.AND. sax(il, i)>0.0) THEN ! cld
	3005	wa(il, i) = sqrt(2. * sax(il, i)) ! cld
[1992]	3006	END IF ! cld
	3007	END DO ! cld
	3008	END DO ! cld
	3009
[5283]	3010	DO i = 1, nl ! cld
	3011	DO il = 1, ncum ! cld
	3012	IF (wa(il, i)>0.0) & ! cld
	3013	siga(il, i) = mac(il, i) / wa(il, i) & ! cld
	3014	* rrd * tvp(il, i) / p(il, i) / 100. / delta ! cld
	3015	siga(il, i) = min(siga(il, i), 1.0) ! cld
	3016	! IM cf. FH
	3017	IF (iflag_clw==0) THEN
	3018	qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) & ! cld
	3019	+ (1. - siga(il, i)) * qcond(il, i) ! cld
	3020	ELSE IF (iflag_clw==1) THEN
	3021	qcondc(il, i) = qcond(il, i) ! cld
	3022	END IF
	3023
	3024	END DO ! cld
[1992]	3025	END DO ! cld
	3026
[5283]	3027	END SUBROUTINE cv30_yield
[1992]	3028
[5283]	3029	!RomP >>>
	3030	SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, &
	3031	d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
	3032	IMPLICIT NONE
[1992]	3033
	3034
	3035
[5283]	3036	! inputs:
	3037	INTEGER ncum, nd, na, nloc, len
	3038	REAL ment(nloc, na, na), sij(nloc, na, na)
	3039	REAL clw(nloc, nd), elij(nloc, na, na)
	3040	REAL ep(nloc, na)
	3041	INTEGER icb(nloc), inb(nloc)
	3042	REAL vprecip(nloc, nd + 1)
	3043	! ouputs:
	3044	REAL da(nloc, na), phi(nloc, na, na)
	3045	REAL phi2(nloc, na, na)
	3046	REAL d1a(nloc, na), dam(nloc, na)
	3047	REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
	3048	! variables pour tracer dans precip de l'AA et des mel
	3049	! local variables:
	3050	INTEGER i, j, k, nam1
	3051	REAL epm(nloc, na, na)
[1992]	3052
[5283]	3053	nam1 = na - 1 ! Introduced because ep is not defined for j=na
	3054	! variables d'Emanuel : du second indice au troisieme
	3055	! ---> tab(i,k,j) -> de l origine k a l arrivee j
	3056	! ment, sij, elij
	3057	! variables personnelles : du troisieme au second indice
	3058	! ---> tab(i,j,k) -> de k a j
	3059	! phi, phi2
[1992]	3060
[5283]	3061	! initialisations
[1992]	3062	DO j = 1, na
	3063	DO i = 1, ncum
[5283]	3064	da(i, j) = 0.
	3065	d1a(i, j) = 0.
	3066	dam(i, j) = 0.
	3067	eplamm(i, j) = 0.
[1992]	3068	END DO
	3069	END DO
[5283]	3070	DO k = 1, na
	3071	DO j = 1, na
	3072	DO i = 1, ncum
	3073	epm(i, j, k) = 0.
	3074	epmlmmm(i, j, k) = 0.
	3075	phi(i, j, k) = 0.
	3076	phi2(i, j, k) = 0.
	3077	END DO
	3078	END DO
	3079	END DO
[1992]	3080
[5283]	3081	! fraction deau condensee dans les melanges convertie en precip : epm
	3082	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
	3083	DO j = 1, nam1
	3084	DO k = 1, j - 1
	3085	DO i = 1, ncum
	3086	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
	3087	!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
	3088	epm(i, j, k) = 1. - (1. - ep(i, j)) * clw(i, j) / max(elij(i, k, j), 1.E-16)
	3089
	3090	epm(i, j, k) = max(epm(i, j, k), 0.0)
	3091	END IF
	3092	END DO
[1992]	3093	END DO
	3094	END DO
	3095
[5283]	3096	DO j = 1, nam1
	3097	DO k = 1, nam1
	3098	DO i = 1, ncum
	3099	IF (k>=icb(i) .AND. k<=inb(i)) THEN
	3100	eplamm(i, j) = eplamm(i, j) + ep(i, j) * clw(i, j) * ment(i, j, k) * (1. - &
	3101	sij(i, j, k))
	3102	END IF
	3103	END DO
[1992]	3104	END DO
	3105	END DO
	3106
[5283]	3107	DO j = 1, nam1
	3108	DO k = 1, j - 1
	3109	DO i = 1, ncum
	3110	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
	3111	epmlmmm(i, j, k) = epm(i, j, k) * elij(i, k, j) * ment(i, k, j)
	3112	END IF
	3113	END DO
[1992]	3114	END DO
	3115	END DO
	3116
[5283]	3117	! matrices pour calculer la tendance des concentrations dans cvltr.F90
	3118	DO j = 1, nam1
	3119	DO k = 1, nam1
	3120	DO i = 1, ncum
	3121	da(i, j) = da(i, j) + (1. - sij(i, k, j)) * ment(i, k, j)
	3122	phi(i, j, k) = sij(i, k, j) * ment(i, k, j)
	3123	d1a(i, j) = d1a(i, j) + ment(i, k, j) * ep(i, k) * (1. - sij(i, k, j))
	3124	END DO
[1992]	3125	END DO
	3126	END DO
	3127
[5283]	3128	DO j = 1, nam1
	3129	DO k = 1, j - 1
	3130	DO i = 1, ncum
	3131	dam(i, j) = dam(i, j) + ment(i, k, j) * epm(i, j, k) * (1. - ep(i, k)) * (1. - &
	3132	sij(i, k, j))
	3133	phi2(i, j, k) = phi(i, j, k) * epm(i, j, k)
	3134	END DO
[1992]	3135	END DO
	3136	END DO
	3137
[5283]	3138	END SUBROUTINE cv30_tracer
	3139	! RomP <<<
[1992]	3140
[5283]	3141	SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
	3142	vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
	3143	dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
	3144	epmlmmm, eplamm, wdtraina, wdtrainm, epmax_diag, iflag1, precip1, vprecip1, evap1, &
	3145	ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
	3146	dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
	3147	elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1, epmax_diag1) ! epmax_cape
	3148	IMPLICIT NONE
[1992]	3149
	3150
	3151
[5283]	3152	! inputs:
	3153	INTEGER len, ncum, nd, ntra, nloc
	3154	INTEGER idcum(nloc)
	3155	INTEGER iflag(nloc)
	3156	INTEGER inb(nloc)
	3157	REAL precip(nloc)
	3158	REAL vprecip(nloc, nd + 1), evap(nloc, nd)
	3159	REAL ep(nloc, nd)
	3160	REAL sig(nloc, nd), w0(nloc, nd)
	3161	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	3162	REAL ftra(nloc, nd, ntra)
	3163	REAL ma(nloc, nd)
	3164	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
	3165	REAL qcondc(nloc, nd)
	3166	REAL wd(nloc), cape(nloc)
	3167	REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
	3168	REAL epmax_diag(nloc) ! epmax_cape
	3169	! RomP >>>
	3170	REAL phi2(nloc, nd, nd)
	3171	REAL d1a(nloc, nd), dam(nloc, nd)
	3172	REAL wdtraina(nloc, nd), wdtrainm(nloc, nd)
	3173	REAL sij(nloc, nd, nd)
	3174	REAL elij(nloc, nd, nd), clw(nloc, nd)
	3175	REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd)
	3176	! RomP <<<
[1992]	3177
[5283]	3178	! outputs:
	3179	INTEGER iflag1(len)
	3180	INTEGER inb1(len)
	3181	REAL precip1(len)
	3182	REAL vprecip1(len, nd + 1), evap1(len, nd) !<<< RomP
	3183	REAL ep1(len, nd) !<<< RomP
	3184	REAL sig1(len, nd), w01(len, nd)
	3185	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
	3186	REAL ftra1(len, nd, ntra)
	3187	REAL ma1(len, nd)
	3188	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
	3189	REAL qcondc1(nloc, nd)
	3190	REAL wd1(nloc), cape1(nloc)
	3191	REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd)
	3192	REAL epmax_diag1(len) ! epmax_cape
	3193	! RomP >>>
	3194	REAL phi21(len, nd, nd)
	3195	REAL d1a1(len, nd), dam1(len, nd)
	3196	REAL wdtraina1(len, nd), wdtrainm1(len, nd)
	3197	REAL sij1(len, nd, nd)
	3198	REAL elij1(len, nd, nd), clw1(len, nd)
	3199	REAL epmlmmm1(len, nd, nd), eplamm1(len, nd)
	3200	! RomP <<<
[1992]	3201
[5283]	3202	! local variables:
	3203	INTEGER i, k, j
[1992]	3204
	3205	DO i = 1, ncum
[5283]	3206	precip1(idcum(i)) = precip(i)
	3207	iflag1(idcum(i)) = iflag(i)
	3208	wd1(idcum(i)) = wd(i)
	3209	inb1(idcum(i)) = inb(i)
	3210	cape1(idcum(i)) = cape(i)
	3211	epmax_diag1(idcum(i)) = epmax_diag(i) ! epmax_cape
[1992]	3212	END DO
	3213
[5283]	3214	DO k = 1, nl
	3215	DO i = 1, ncum
	3216	vprecip1(idcum(i), k) = vprecip(i, k)
	3217	evap1(idcum(i), k) = evap(i, k) !<<< RomP
	3218	sig1(idcum(i), k) = sig(i, k)
	3219	w01(idcum(i), k) = w0(i, k)
	3220	ft1(idcum(i), k) = ft(i, k)
	3221	fq1(idcum(i), k) = fq(i, k)
	3222	fu1(idcum(i), k) = fu(i, k)
	3223	fv1(idcum(i), k) = fv(i, k)
	3224	ma1(idcum(i), k) = ma(i, k)
	3225	upwd1(idcum(i), k) = upwd(i, k)
	3226	dnwd1(idcum(i), k) = dnwd(i, k)
	3227	dnwd01(idcum(i), k) = dnwd0(i, k)
	3228	qcondc1(idcum(i), k) = qcondc(i, k)
	3229	da1(idcum(i), k) = da(i, k)
	3230	mp1(idcum(i), k) = mp(i, k)
	3231	! RomP >>>
	3232	ep1(idcum(i), k) = ep(i, k)
	3233	d1a1(idcum(i), k) = d1a(i, k)
	3234	dam1(idcum(i), k) = dam(i, k)
	3235	clw1(idcum(i), k) = clw(i, k)
	3236	eplamm1(idcum(i), k) = eplamm(i, k)
	3237	wdtraina1(idcum(i), k) = wdtraina(i, k)
	3238	wdtrainm1(idcum(i), k) = wdtrainm(i, k)
	3239	! RomP <<<
	3240	END DO
	3241	END DO
[1992]	3242
[5283]	3243	DO i = 1, ncum
	3244	sig1(idcum(i), nd) = sig(i, nd)
	3245	END DO
[1992]	3246
[5283]	3247
	3248	! do 2100 j=1,ntra
	3249	! do 2110 k=1,nd ! oct3
	3250	! do 2120 i=1,ncum
	3251	! ftra1(idcum(i),k,j)=ftra(i,k,j)
	3252	! 2120 continue
	3253	! 2110 continue
	3254	! 2100 continue
	3255	DO j = 1, nd
	3256	DO k = 1, nd
	3257	DO i = 1, ncum
	3258	sij1(idcum(i), k, j) = sij(i, k, j)
	3259	phi1(idcum(i), k, j) = phi(i, k, j)
	3260	phi21(idcum(i), k, j) = phi2(i, k, j)
	3261	elij1(idcum(i), k, j) = elij(i, k, j)
	3262	epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
	3263	END DO
[1992]	3264	END DO
	3265	END DO
	3266
[5283]	3267	END SUBROUTINE cv30_uncompress
[1992]	3268
[5283]	3269	SUBROUTINE cv30_epmax_fn_cape(nloc, ncum, nd &
	3270	, cape, ep, hp, icb, inb, clw, nk, t, h, lv &
	3271	, epmax_diag)
	3272	USE conema3_mod_h
[5346]	3273	USE lmdz_cv_ini, ONLY : cpd, cpv
[2481]	3274
[5283]	3275	IMPLICIT NONE
[2481]	3276
[5283]	3277	! On fait varier epmax en fn de la cape
	3278	! Il faut donc recalculer ep, et hp qui a déjà été calculé et
	3279	! qui en dépend
	3280	! Toutes les autres variables fn de ep sont calculées plus bas.
[2481]	3281
	3282
[5283]	3283
	3284	! inputs:
	3285	INTEGER ncum, nd, nloc
	3286	INTEGER icb(nloc), inb(nloc)
	3287	REAL cape(nloc)
	3288	REAL clw(nloc, nd), lv(nloc, nd), t(nloc, nd), h(nloc, nd)
	3289	INTEGER nk(nloc)
	3290	! inouts:
	3291	REAL ep(nloc, nd)
	3292	REAL hp(nloc, nd)
	3293	! outputs ou local
	3294	REAL epmax_diag(nloc)
	3295	! locals
	3296	INTEGER i, k
	3297	REAL hp_bak(nloc, nd)
	3298	CHARACTER (LEN = 20) :: modname = 'cv30_epmax_fn_cape'
	3299	CHARACTER (LEN = 80) :: abort_message
	3300
	3301	! on recalcule ep et hp
	3302
	3303	IF (coef_epmax_cape>1e-12) THEN
	3304	DO i = 1, ncum
	3305	epmax_diag(i) = epmax - coef_epmax_cape * sqrt(cape(i))
	3306	DO k = 1, nl
	3307	ep(i, k) = ep(i, k) / epmax * epmax_diag(i)
	3308	ep(i, k) = amax1(ep(i, k), 0.0)
	3309	ep(i, k) = amin1(ep(i, k), epmax_diag(i))
[2481]	3310	enddo
[5283]	3311	enddo
[2481]	3312
[5283]	3313	! On recalcule hp:
	3314	DO k = 1, nl
	3315	DO i = 1, ncum
	3316	hp_bak(i, k) = hp(i, k)
	3317	enddo
[2481]	3318	enddo
[5283]	3319	DO k = 1, nlp
	3320	DO i = 1, ncum
	3321	hp(i, k) = h(i, k)
	3322	enddo
[2481]	3323	enddo
[5283]	3324	DO k = minorig + 1, nl
	3325	DO i = 1, ncum
	3326	IF((k>=icb(i)).AND.(k<=inb(i)))THEN
	3327	hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k)
	3328	endif
	3329	enddo
[2481]	3330	enddo !do k=minorig+1,n
[5283]	3331	! WRITE(,) 'cv30_routines 6218: hp(1,20)=',hp(1,20)
	3332	DO i = 1, ncum
	3333	DO k = 1, nl
	3334	IF (abs(hp_bak(i, k) - hp(i, k))>0.01) THEN
	3335	WRITE(, ) 'i,k=', i, k
	3336	WRITE(, ) 'coef_epmax_cape=', coef_epmax_cape
	3337	WRITE(, ) 'epmax_diag(i)=', epmax_diag(i)
	3338	WRITE(, ) 'ep(i,k)=', ep(i, k)
	3339	WRITE(, ) 'hp(i,k)=', hp(i, k)
	3340	WRITE(, ) 'hp_bak(i,k)=', hp_bak(i, k)
	3341	WRITE(, ) 'h(i,k)=', h(i, k)
	3342	WRITE(, ) 'nk(i)=', nk(i)
	3343	WRITE(, ) 'h(i,nk(i))=', h(i, nk(i))
	3344	WRITE(, ) 'lv(i,k)=', lv(i, k)
	3345	WRITE(, ) 't(i,k)=', t(i, k)
	3346	WRITE(, ) 'clw(i,k)=', clw(i, k)
	3347	WRITE(, ) 'cpd,cpv=', cpd, cpv
	3348	CALL abort_physic(modname, abort_message, 1)
	3349	endif
	3350	enddo !do k=1,nl
	3351	enddo !do i=1,ncum
	3352	ENDIF !if (coef_epmax_cape.gt.1e-12) THEN
	3353	END SUBROUTINE cv30_epmax_fn_cape
[2481]	3354
	3355
[5283]	3356	END MODULE cv30_routines_mod
[2481]	3357
[5283]	3358

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

Download in other formats:

Original Format