Context Navigation

cv3_routines.F90 @ 1999

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1920:1997 into testing branch
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: 124.8 KB

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

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

Download in other formats:

Original Format