Context Navigation

cv3_routines.F90 @ 3946

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

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

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

Download in other formats:

Original Format