Context Navigation

cv3_routines.F90 @ 2300

Last change on this file since 2300 was 2298, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes -r2237:2291 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 130.0 KB

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

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

Download in other formats:

Original Format