Context Navigation

cv3_routines.F90 @ 2509

Last change on this file since 2509 was 2508, checked in by jyg, 9 years ago

Optimization of cv3_yield (in cv3_routines.F90):
suppression of the triple vertical loops. The
optimized code is activated by:

ok_optim_yield=y

which changes numerical results.

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

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

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

Download in other formats:

Original Format