Context Navigation

cv3_routines.F90 @ 2126

Last change on this file since 2126 was 2110, checked in by lguez, 11 years ago

Abort if surface pressure becomes negative. The call to abort_gcm
already was in the sequential version but not in dyn3dpar nor in
dyn3dmem. In dyn3dmem, there is a variable checksum_all, which is
always true. The call to MPI_ALLREDUCE which should update
checksum_all is commented out. I do not know why (performance?).

Non-ASCII characters in comments are not always rendered properly and
they risk being lost. See revision [1740].

Bug fix in procedure convect2: the dimension len must be declared
before the array idcum which has this dimension. Bug fix in procedure
icefrac: the dimensions nl and len must be declared before the arrays
which have these dimensions.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 124.1 KB

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

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

Download in other formats:

Original Format