Context Navigation

cv3_routines.F90 @ 3624

Last change on this file since 3624 was 3624, checked in by jyg, 5 years ago
Correction to commit 3614: A new formula is used for "phinu2p" (hopefully correct) in subroutine cv3_undilute2 in cv3_routines.F90.
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: 178.8 KB

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

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

Download in other formats:

Original Format