Context Navigation

cv3_routines.F90 @ 2415

Last change on this file since 2415 was 2398, checked in by jyg, 9 years ago

Introduction of two new flags (ok_conv_stop
[Def=F], ok_intermittent [Def=F]) and one new
parameter (tau_stop [Def=15000]) in
conv_param.data:

. ok_conv_stop=T => convective mass fluxes are

set to zero if there is no trigger for a time
lapse longer than tau_stop.

. ok_intermittent=T => intermittent convection is

represented; the change concerns the upper
bound on the cloud base mass flux in
cv3p2_closure.F90.

The upper bound on Alp in physiq.F90 should

also be changed: this is still to be done.

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

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

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

Download in other formats:

Original Format