Context Navigation

lmdz_cv30.f90 @ 5442

Last change on this file since 5442 was 5159, checked in by abarral, 5 months ago
Put dimensions.h and paramet.h into modules
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 108.7 KB

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

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

Download in other formats:

Original Format