Context Navigation

cv3_routines.F90 @ 2489

Last change on this file since 2489 was 2481, checked in by fhourdin, 8 years ago
Introduction d'une dependance epmax=f(Cape) sur proposition de Camille Risi
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: 140.6 KB

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

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

Download in other formats:

Original Format