Context Navigation

cv3_routines.F90 @ 3627

Last change on this file since 3627 was 3617, checked in by lguez, 4 years ago
Sync latest trunk changes to branch `Ocean_skin`
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.6 KB

Rev	Line
[1992]	1
[1403]	2	! $Id: cv3_routines.F90 3617 2019-12-17 14:44:15Z lguez $
[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"
[3605]	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)
[3605]	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)
	134	qsat_depends_on_qt = .FALSE.
	135	CALL getin_p('qsat_depends_on_qt',qsat_depends_on_qt)
	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)
	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
[3605]	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"
[3605]	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)
[3605]	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
[3605]	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
[3605]	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, &
[3605]	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
[3605]	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
[3605]	1189	REAL smallestreal
	1190	REAL tg, qg, dqgdT, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
	1191	REAL :: phinu2p
[1992]	1192	REAL als
[3605]	1193	REAL :: qsat_new, snew
	1194	REAL, DIMENSION (nloc,nd) :: qi
	1195	REAL, DIMENSION (nloc,nd) :: ha ! moist static energy of adiabatic ascents
	1196	! taking into account precip ejection
	1197	REAL, DIMENSION (nloc,nd) :: hla ! liquid water static energy of adiabatic ascents
	1198	! taking into account precip ejection
	1199	REAL, DIMENSION (nloc,nd) :: qcld ! specific cloud water
	1200	REAL, DIMENSION (nloc,nd) :: qhsat ! specific humidity at saturation
	1201	REAL, DIMENSION (nloc,nd) :: dqhsatdT ! dqhsat/dT
	1202	REAL, DIMENSION (nloc,nd) :: frac ! ice fraction function of envt temperature
	1203	REAL, DIMENSION (nloc,nd) :: qps ! specific solid precipitation
	1204	REAL, DIMENSION (nloc,nd) :: qpl ! specific liquid precipitation
	1205	REAL, DIMENSION (nloc) :: ah0, cape, capem, byp
	1206	LOGICAL, DIMENSION (nloc) :: lcape
	1207	INTEGER, DIMENSION (nloc) :: iposit
	1208	REAL :: denomm1
	1209	REAL :: by, defrac, pden, tbis
	1210	REAL :: fracg
	1211	REAL :: deltap
	1212	REAL, SAVE :: Tx, Tm
	1213	DATA Tx/263.15/, Tm/243.15/
	1214	!$OMP THREADPRIVATE(Tx, Tm)
	1215	REAL :: aa, bb, dd, ddelta, discr
	1216	REAL :: ff, fp
	1217	REAL :: coefx, coefm, Zx, Zm, Ux, U, Um
[524]	1218
[2638]	1219	IF (prt_level >= 10) THEN
[3605]	1220	print *,'cv3_undilute2.0. icvflag_Tpa, t(1,k), q(1,k), qs(1,k) ', &
	1221	icvflag_Tpa, (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
[2638]	1222	ENDIF
[3605]	1223	smallestreal=tiny(smallestreal)
[2638]	1224
[2007]	1225	! =====================================================================
	1226	! --- SOME INITIALIZATIONS
	1227	! =====================================================================
[524]	1228
[1992]	1229	DO k = 1, nl
	1230	DO i = 1, ncum
	1231	qi(i, k) = 0.
	1232	END DO
	1233	END DO
[524]	1234
[3605]	1235
[2007]	1236	! =====================================================================
	1237	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	1238	! =====================================================================
[524]	1239
[2007]	1240	! --- The procedure is to solve the equation.
	1241	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
[524]	1242
[2007]	1243	! * Calculate certain parcel quantities, including static energy *
[524]	1244
	1245
[1992]	1246	DO i = 1, ncum
[2007]	1247	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))*tnk(i)+ &
	1248	! debug qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
	1249	qnk(i)(lv0-clmcpv(tnk(i)-273.15)) + gznk(i)
[1992]	1250	END DO
[3605]	1251	!
	1252	! Ice fraction
	1253	!
	1254	IF (cvflag_ice) THEN
	1255	DO k = minorig, nl
	1256	DO i = 1, ncum
	1257	frac(i, k) = (Tx - t(i,k))/(Tx - Tm)
	1258	frac(i, k) = min(max(frac(i,k),0.0), 1.0)
	1259	END DO
	1260	END DO
	1261	! Below cloud base, set ice fraction to cloud base value
	1262	DO k = 1, nl
	1263	DO i = 1, ncum
	1264	IF (k<icb(i)) THEN
	1265	frac(i,k) = frac(i,icb(i))
	1266	END IF
	1267	END DO
	1268	END DO
	1269	ELSE
	1270	DO k = 1, nl
	1271	DO i = 1, ncum
	1272	frac(i,k) = 0.
	1273	END DO
	1274	END DO
	1275	ENDIF ! (cvflag_ice)
[524]	1276
	1277
[3605]	1278	DO k = minorig, nl
	1279	DO i = 1,ncum
	1280	ha(i,k) = ah0(i)
	1281	hla(i,k) = hnk(i)
	1282	qta(i,k) = qnk(i)
	1283	qpreca(i,k) = 0.
	1284	frac_a(i,k) = 0.
	1285	frac_s(i,k) = frac(i,k)
	1286	qpl(i,k) = 0.
	1287	qps(i,k) = 0.
	1288	qhsat(i,k) = qs(i,k)
	1289	qcld(i,k) = max(qta(i,k)-qhsat(i,k),0.)
	1290	IF (k <= icb(i)+1) THEN
	1291	qhsat(i,k) = qnk(i)-clw(i,k)
	1292	qcld(i,k) = clw(i,k)
	1293	ENDIF
	1294	ENDDO
	1295	ENDDO
	1296
	1297	!jyg<
	1298	! =====================================================================
	1299	! --- SET THE THE FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
	1300	! =====================================================================
	1301	DO k = 1, nl
	1302	DO i = 1, ncum
	1303	ep(i, k) = 0.0
	1304	sigp(i, k) = spfac
	1305	END DO
	1306	END DO
	1307	!>jyg
	1308	!
	1309
[2007]	1310	! * Find lifted parcel quantities above cloud base *
[524]	1311
[3605]	1312	!----------------------------------------------------------------------------
	1313	!
	1314	IF (icvflag_Tpa == 2) THEN
	1315	!
	1316	!----------------------------------------------------------------------------
	1317	!
	1318	DO k = minorig + 1, nl
	1319	DO i = 1,ncum
	1320	tp(i,k) = t(i,k)
	1321	ENDDO
	1322	!! alv = lv0 - clmcpv*(t(i,k)-273.15)
	1323	!! alf = lf0 + clmci*(t(i,k)-273.15)
	1324	!! als = alf + alv
	1325	DO j = 1,4
	1326	DO i = 1, ncum
	1327	! ori if(k.ge.(icb(i)+1))then
	1328	IF (k>=(icbs(i)+1)) THEN ! convect3
	1329	tg = tp(i, k)
	1330	IF (tg .gt. Tx) THEN
	1331	es = 6.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
	1332	qg = epses/(p(i,k)-es(1.-eps))
	1333	ELSE
	1334	esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
	1335	qg = epsesi/(p(i,k)-esi(1.-eps))
	1336	ENDIF
	1337	! Ice fraction
	1338	ff = 0.
	1339	fp = 1./(Tx - Tm)
	1340	IF (tg < Tx) THEN
	1341	IF (tg > Tm) THEN
	1342	ff = (Tx - tg)*fp
	1343	ELSE
	1344	ff = 1.
	1345	ENDIF ! (tg > Tm)
	1346	ENDIF ! (tg < Tx)
	1347	! Intermediate variables
	1348	aa = cpd + (cl-cpd)qnk(i) + lv(i,k)lv(i,k)qg/(rrvtg*tg)
	1349	ahg = (cpd + (cl-cpd)qnk(i))tg + lv(i,k)*qg - &
	1350	lf(i,k)ff(qnk(i) - qg) + gz(i,k)
	1351	dd = lf(i,k)lv(i,k)qg/(rrvtgtg)
	1352	ddelta = lf(i,k)*(qnk(i) - qg)
	1353	bb = aa + ddeltafp + ddfp*(Tx-tg)
	1354	! Compute Zx and Zm
	1355	coefx = aa
	1356	coefm = aa + dd
	1357	IF (tg .gt. Tx) THEN
	1358	Zx = ahg + coefx*(Tx - tg)
	1359	Zm = ahg - ddelta + coefm*(Tm - tg)
	1360	ELSE
	1361	IF (tg .gt. Tm) THEN
	1362	Zx = ahg + (coefx +fpddelta)(Tx - Tg)
	1363	Zm = ahg + (coefm +fpddelta)(Tm - Tg)
	1364	ELSE
	1365	Zx = ahg + ddelta + coefx*(Tx - tg)
	1366	Zm = ahg + coefm*(Tm - tg)
	1367	ENDIF ! (tg .gt. Tm)
	1368	ENDIF ! (tg .gt. Tx)
	1369	! Compute the masks Um, U, Ux
	1370	Um = (sign(1., Zm-ah0(i))+1.)/2.
	1371	Ux = (sign(1., ah0(i)-Zx)+1.)/2.
	1372	U = (1. - Um)*(1. - Ux)
	1373	! Compute the updated parcell temperature Tp : 3 cases depending on tg value
	1374	IF (tg .gt. Tx) THEN
	1375	discr = bbbb - 4ddfp(ah0(i) - ahg + ddeltafp(Tx-tg))
	1376	Tp(i,k) = tg + &
	1377	Um* (ah0(i) - ahg + ddelta) /(aa + dd) + &
	1378	U 2(ah0(i) - ahg + ddeltafp(Tx-tg))/(bb + sqrt(discr)) + &
	1379	Ux* (ah0(i) - ahg) /aa
	1380	ELSEIF (tg .gt. Tm) THEN
	1381	discr = bbbb - 4ddfp(ah0(i) - ahg)
	1382	Tp(i,k) = tg + &
	1383	Um* (ah0(i) - ahg + ddeltafp(tg-Tm))/(aa + dd) + &
	1384	U 2(ah0(i) - ahg) /(bb + sqrt(discr)) + &
	1385	Ux* (ah0(i) - ahg + ddeltafp(tg-Tx))/aa
	1386	ELSE
	1387	discr = bbbb - 4ddfp(ah0(i) - ahg + ddeltafp(Tm-tg))
	1388	Tp(i,k) = tg + &
	1389	Um* (ah0(i) - ahg) /(aa + dd) + &
	1390	U 2(ah0(i) - ahg + ddeltafp(Tm-tg))/(bb + sqrt(discr)) + &
	1391	Ux* (ah0(i) - ahg - ddelta) /aa
	1392	ENDIF ! (tg .gt. Tx)
	1393	!
	1394	!! print *,' j, k, Um, U, Ux, aa, bb, discr, dd, ddelta ', j, k, Um, U, Ux, aa, bb, discr, dd, ddelta
	1395	!! print *,' j, k, ah0(i), ahg, tg, qg, tp(i,k), ff ', j, k, ah0(i), ahg, tg, qg, tp(i,k), ff
	1396	END IF ! (k>=(icbs(i)+1))
	1397	END DO ! i = 1, ncum
	1398	END DO ! j = 1,4
	1399	DO i = 1, ncum
	1400	IF (k>=(icbs(i)+1)) THEN ! convect3
	1401	tg = tp(i, k)
	1402	IF (tg .gt. Tx) THEN
	1403	es = 6.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
	1404	qg = epses/(p(i,k)-es(1.-eps))
	1405	ELSE
	1406	esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
	1407	qg = epsesi/(p(i,k)-esi(1.-eps))
	1408	ENDIF
	1409	clw(i, k) = qnk(i) - qg
	1410	clw(i, k) = max(0.0, clw(i,k))
	1411	tvp(i, k) = max(0., tp(i,k)*(1.+qg/eps-qnk(i)))
	1412	! print*,tvp(i,k),'tvp'
	1413	IF (clw(i,k)<1.E-11) THEN
	1414	tp(i, k) = tv(i, k)
	1415	tvp(i, k) = tv(i, k)
	1416	END IF ! (clw(i,k)<1.E-11)
	1417	END IF ! (k>=(icbs(i)+1))
	1418	END DO ! i = 1, ncum
	1419	END DO ! k = minorig + 1, nl
	1420	!----------------------------------------------------------------------------
	1421	!
	1422	ELSE IF (icvflag_Tpa == 1) THEN ! (icvflag_Tpa == 2)
	1423	!
	1424	!----------------------------------------------------------------------------
	1425	!
	1426	DO k = minorig + 1, nl
	1427	DO i = 1,ncum
	1428	tp(i,k) = t(i,k)
	1429	ENDDO
	1430	!! alv = lv0 - clmcpv*(t(i,k)-273.15)
	1431	!! alf = lf0 + clmci*(t(i,k)-273.15)
	1432	!! als = alf + alv
	1433	DO j = 1,4
	1434	DO i = 1, ncum
	1435	! ori if(k.ge.(icb(i)+1))then
	1436	IF (k>=(icbs(i)+1)) THEN ! convect3
	1437	tg = tp(i, k)
	1438	IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
	1439	es = 6.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
	1440	qg = epses/(p(i,k)-es(1.-eps))
	1441	dqgdT = lv(i,k)qg/(rrvtg*tg)
	1442	ELSE
	1443	esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
	1444	qg = epsesi/(p(i,k)-esi(1.-eps))
	1445	dqgdT = (lv(i,k)+lf(i,k))qg/(rrvtg*tg)
	1446	ENDIF
	1447	IF (qsat_depends_on_qt) THEN
	1448	dqgdT = dqgdT(1.-qta(i,k-1))/(1.-qg)*2
	1449	qg = qg*(1.-qta(i,k-1))/(1.-qg)
	1450	ENDIF
	1451	ahg = (cpd + (cl-cpd)qta(i,k-1))tg + lv(i,k)*qg - &
	1452	lf(i,k)frac(i,k)(qta(i,k-1) - qg) + gz(i,k)
	1453	Tp(i,k) = tg + (ah0(i) - ahg)/ &
	1454	(cpd + (cl-cpd)qta(i,k-1) + (lv(i,k)+frac(i,k)lf(i,k))*dqgdT)
	1455	!! print *,'undilute2 iterations k, Tp(i,k), ah0(i), ahg ', &
	1456	!! k, Tp(i,k), ah0(i), ahg
	1457	END IF ! (k>=(icbs(i)+1))
	1458	END DO ! i = 1, ncum
	1459	END DO ! j = 1,4
	1460	DO i = 1, ncum
	1461	IF (k>=(icbs(i)+1)) THEN ! convect3
	1462	tg = tp(i, k)
	1463	IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
	1464	es = 6.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
	1465	qg = epses/(p(i,k)-es(1.-eps))
	1466	ELSE
	1467	esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
	1468	qg = epsesi/(p(i,k)-esi(1.-eps))
	1469	ENDIF
	1470	IF (qsat_depends_on_qt) THEN
	1471	qg = qg*(1.-qta(i,k-1))/(1.-qg)
	1472	ENDIF
	1473	qhsat(i,k) = qg
	1474	END IF ! (k>=(icbs(i)+1))
	1475	END DO ! i = 1, ncum
	1476	DO i = 1, ncum
	1477	IF (k>=(icbs(i)+1)) THEN ! convect3
	1478	clw(i, k) = qta(i,k-1) - qhsat(i,k)
	1479	clw(i, k) = max(0.0, clw(i,k))
	1480	tvp(i, k) = max(0., tp(i,k)*(1.+qhsat(i,k)/eps-qta(i,k-1)))
	1481	! print*,tvp(i,k),'tvp'
	1482	IF (clw(i,k)<1.E-11) THEN
	1483	tp(i, k) = tv(i, k)
	1484	tvp(i, k) = tv(i, k)
	1485	END IF ! (clw(i,k)<1.E-11)
	1486	END IF ! (k>=(icbs(i)+1))
	1487	END DO ! i = 1, ncum
	1488	!
	1489	IF (cvflag_prec_eject) THEN
	1490	DO i = 1, ncum
	1491	IF (k>=(icbs(i)+1)) THEN ! convect3
	1492	! Specific precipitation (liquid and solid) and ice content
	1493	! before ejection of precipitation !!jygprl
	1494	elacrit = elcrit*min(max(1.-(tp(i,k)-T0)/Tlcrit, 0.), 1.) !!jygprl
	1495	!!!! qcld(i,k) = min(clw(i,k), elacrit) !!jygprl
	1496	qcld(i,k) = min(clw(i,k), elacrit*(1.-qta(i,k-1))/(1.-elacrit)) !!jygprl
[3617]	1497	phinu2p = max(qhsat(i,k-1) + qcld(i,k-1) - (qhsat(i,k) + qcld(i,k)),0.) !!jygprl
[3605]	1498	qpl(i,k) = qpl(i,k-1) + (1.-frac(i,k))*phinu2p !!jygprl
	1499	qps(i,k) = qps(i,k-1) + frac(i,k) *phinu2p !!jygprl
	1500	qi(i,k) = (1.-ejectliq)clw(i,k)frac(i,k) + & !!jygprl
	1501	ejectliq(qps(i,k-1) + frac(i,k)(phinu2p+qcld(i,k))) !!jygprl
	1502	!!
	1503	! =====================================================================================
	1504	! Ejection of precipitation from adiabatic ascents if requested (cvflag_prec_eject=True):
	1505	! Compute the steps of total water (qta), of moist static energy (ha), of specific
	1506	! precipitation (qpl and qps) and of specific cloud water (qcld) associated with precipitation
	1507	! ejection.
	1508	! =====================================================================================
	1509	!
	1510	! Verif
	1511	qpreca(i,k) = ejectliqqpl(i,k) + ejecticeqps(i,k) !!jygprl
	1512	frac_a(i,k) = ejectice*qps(i,k)/max(qpreca(i,k),smallestreal) !!jygprl
	1513	frac_s(i,k) = (1.-ejectliq)*frac(i,k) + & !!jygprl
	1514	ejectliq(1. - (qpl(i,k)+(1.-frac(i,k))qcld(i,k))/max(clw(i,k),smallestreal)) !!jygprl
	1515	!
	1516	denomm1 = 1./(1. - qpreca(i,k))
	1517	!
	1518	qta(i,k) = qta(i,k-1) - &
	1519	qpreca(i,k)(1.-qta(i,k-1))denomm1
	1520	ha(i,k) = ha(i,k-1) + &
	1521	( qpreca(i,k)(-(1.-qta(i,k-1))(cl-cpd)*tp(i,k) + &
	1522	lv(i,k)qhsat(i,k) - lf(i,k)(frac_s(i,k)*qcld(i,k)+qps(i,k))) + &
	1523	lf(i,k)ejecticeqps(i,k))*denomm1
	1524	hla(i,k) = hla(i,k-1) + &
	1525	( qpreca(i,k)(-(1.-qta(i,k-1))(cpv-cpd)*tp(i,k) - &
	1526	lv(i,k)((1.-frac_s(i,k))qcld(i,k)+qpl(i,k)) - &
	1527	(lv(i,k)+lf(i,k))(frac_s(i,k)qcld(i,k)+qps(i,k))) + &
	1528	lv(i,k)ejectliqqpl(i,k) + (lv(i,k)+lf(i,k))ejecticeqps(i,k))*denomm1
	1529	qpl(i,k) = qpl(i,k)(1.-ejectliq)denomm1
	1530	qps(i,k) = qps(i,k)(1.-ejectice)denomm1
	1531	qcld(i,k) = qcld(i,k)*denomm1
	1532	qhsat(i,k) = qhsat(i,k)*(1.-qta(i,k))/(1.-qta(i,k-1))
	1533	END IF ! (k>=(icbs(i)+1))
	1534	END DO ! i = 1, ncum
	1535	ENDIF ! (cvflag_prec_eject)
	1536	!
	1537	END DO ! k = minorig + 1, nl
	1538	!
	1539	!----------------------------------------------------------------------------
	1540	!
	1541	ELSE IF (icvflag_Tpa == 0) THEN! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
	1542	!
	1543	!----------------------------------------------------------------------------
	1544	!
[1992]	1545	DO k = minorig + 1, nl
	1546	DO i = 1, ncum
[2007]	1547	! ori if(k.ge.(icb(i)+1))then
	1548	IF (k>=(icbs(i)+1)) THEN ! convect3
[1992]	1549	tg = t(i, k)
	1550	qg = qs(i, k)
[2007]	1551	! debug alv=lv0-clmcpv*(t(i,k)-t0)
[1992]	1552	alv = lv0 - clmcpv*(t(i,k)-273.15)
[524]	1553
[2007]	1554	! First iteration.
[524]	1555
[2007]	1556	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	1557	s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
	1558	alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
[1992]	1559	s = 1./s
[2007]	1560	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
[1992]	1561	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
	1562	tg = tg + s*(ah0(i)-ahg)
[2007]	1563	! ori tg=max(tg,35.0)
	1564	! debug tc=tg-t0
[1992]	1565	tc = tg - 273.15
	1566	denom = 243.5 + tc
[2007]	1567	denom = max(denom, 1.0) ! convect3
	1568	! ori if(tc.ge.0.0)then
[1992]	1569	es = 6.112exp(17.67tc/denom)
[2007]	1570	! ori else
	1571	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	1572	! ori endif
[1992]	1573	qg = epses/(p(i,k)-es(1.-eps))
[524]	1574
[2007]	1575	! Second iteration.
[524]	1576
[2007]	1577	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	1578	! ori s=1./s
	1579	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
[1992]	1580	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
	1581	tg = tg + s*(ah0(i)-ahg)
[2007]	1582	! ori tg=max(tg,35.0)
	1583	! debug tc=tg-t0
[1992]	1584	tc = tg - 273.15
	1585	denom = 243.5 + tc
[2007]	1586	denom = max(denom, 1.0) ! convect3
	1587	! ori if(tc.ge.0.0)then
[1992]	1588	es = 6.112exp(17.67tc/denom)
[2007]	1589	! ori else
	1590	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	1591	! ori endif
[1992]	1592	qg = epses/(p(i,k)-es(1.-eps))
[524]	1593
[2007]	1594	! debug alv=lv0-clmcpv*(t(i,k)-t0)
[1992]	1595	alv = lv0 - clmcpv*(t(i,k)-273.15)
[2007]	1596	! print*,'cpd dans convect2 ',cpd
	1597	! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
	1598	! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
[524]	1599
[2007]	1600	! ori c approximation here:
	1601	! ori tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
[524]	1602
[2007]	1603	! convect3: no approximation:
[1992]	1604	IF (cvflag_ice) THEN
	1605	tp(i, k) = max(0., (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i)))
	1606	ELSE
	1607	tp(i, k) = (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
	1608	END IF
[524]	1609
[1992]	1610	clw(i, k) = qnk(i) - qg
	1611	clw(i, k) = max(0.0, clw(i,k))
	1612	rg = qg/(1.-qnk(i))
[2007]	1613	! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
	1614	! convect3: (qg utilise au lieu du vrai mixing ratio rg):
[1992]	1615	tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
	1616	IF (cvflag_ice) THEN
	1617	IF (clw(i,k)<1.E-11) THEN
	1618	tp(i, k) = tv(i, k)
	1619	tvp(i, k) = tv(i, k)
	1620	END IF
	1621	END IF
[2253]	1622	!jyg<
	1623	!! END IF ! Endif moved to the end of the loop
	1624	!>jyg
[524]	1625
[1992]	1626	IF (cvflag_ice) THEN
[2007]	1627	!CR:attention boucle en klon dans Icefrac
	1628	! Call Icefrac(t,clw,qi,nl,nloc)
[1992]	1629	IF (t(i,k)>263.15) THEN
	1630	qi(i, k) = 0.
	1631	ELSE
	1632	IF (t(i,k)<243.15) THEN
	1633	qi(i, k) = clw(i, k)
	1634	ELSE
	1635	fracg = (263.15-t(i,k))/20
	1636	qi(i, k) = clw(i, k)*fracg
	1637	END IF
	1638	END IF
[2007]	1639	!CR: fin test
[1992]	1640	IF (t(i,k)<263.15) THEN
[2007]	1641	!CR: on commente les calculs d'Arnaud car division par zero
	1642	! nouveau calcul propose par JYG
	1643	! alv=lv0-clmcpv*(t(i,k)-273.15)
	1644	! alf=lf0-clmci*(t(i,k)-273.15)
	1645	! tg=tp(i,k)
	1646	! tc=tp(i,k)-273.15
	1647	! denom=243.5+tc
	1648	! do j=1,3
	1649	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1650	! il faudra que esi vienne en argument de la convection
	1651	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1652	! tbis=t(i,k)+(tp(i,k)-tg)
	1653	! esi=exp(23.33086-(6111.72784/tbis) + &
	1654	! 0.15215*log(tbis))
	1655	! qsat_new=epsesi/(p(i,k)-esi(1.-eps))
	1656	! snew=cpd(1.-qnk(i))+clqnk(i)+alvalvqsat_new/ &
	1657	! (rrvtbistbis)
	1658	! snew=1./snew
	1659	! print*,esi,qsat_new,snew,'esi,qsat,snew'
	1660	! tp(i,k)=tg+(alfqi(i,k)+alvqg(1.-(esi/es)))snew
	1661	! print*,k,tp(i,k),qnk(i),'avec glace'
	1662	! print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
	1663	! enddo
[524]	1664
[1992]	1665	alv = lv0 - clmcpv*(t(i,k)-273.15)
	1666	alf = lf0 + clmci*(t(i,k)-273.15)
	1667	als = alf + alv
	1668	tg = tp(i, k)
	1669	tp(i, k) = t(i, k)
	1670	DO j = 1, 3
	1671	esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
	1672	qsat_new = epsesi/(p(i,k)-esi(1.-eps))
[2007]	1673	snew = cpd(1.-qnk(i)) + clqnk(i) + alvalsqsat_new/ &
	1674	(rrvtp(i,k)tp(i,k))
[1992]	1675	snew = 1./snew
[2007]	1676	! c print*,esi,qsat_new,snew,'esi,qsat,snew'
	1677	tp(i, k) = tp(i, k) + &
	1678	((cpd(1.-qnk(i))+clqnk(i))*(tg-tp(i,k)) + &
	1679	alv(qg-qsat_new)+alfqi(i,k))*snew
	1680	! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
	1681	! 'k,tp,q,qt,qi avec glace'
[1992]	1682	END DO
[524]	1683
[2007]	1684	!CR:reprise du code AJ
[1992]	1685	clw(i, k) = qnk(i) - qsat_new
	1686	clw(i, k) = max(0.0, clw(i,k))
	1687	tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
[2007]	1688	! print*,tvp(i,k),'tvp'
[1992]	1689	END IF
	1690	IF (clw(i,k)<1.E-11) THEN
	1691	tp(i, k) = tv(i, k)
	1692	tvp(i, k) = tv(i, k)
	1693	END IF
	1694	END IF ! (cvflag_ice)
[2253]	1695	!jyg<
	1696	END IF ! (k>=(icbs(i)+1))
	1697	!>jyg
[1992]	1698	END DO
	1699	END DO
[1849]	1700
[3605]	1701	!----------------------------------------------------------------------------
	1702	!
	1703	ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE (icvflag_Tpa == 0)
	1704	!
	1705	!----------------------------------------------------------------------------
	1706	!
[2007]	1707	! =====================================================================
[3605]	1708	! --- SET THE PRECIPITATION EFFICIENCIES
[2007]	1709	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
	1710	! =====================================================================
[2253]	1711	!
[1992]	1712	IF (flag_epkeorig/=1) THEN
	1713	DO k = 1, nl ! convect3
	1714	DO i = 1, ncum
[2253]	1715	!jyg<
	1716	IF(k>=icb(i)) THEN
	1717	!>jyg
	1718	pden = ptcrit - pbcrit
	1719	ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
	1720	ep(i, k) = max(ep(i,k), 0.0)
	1721	ep(i, k) = min(ep(i,k), epmax)
	1722	!! sigp(i, k) = spfac ! jyg
	1723	ENDIF ! (k>=icb(i))
[1992]	1724	END DO
	1725	END DO
	1726	ELSE
	1727	DO k = 1, nl
	1728	DO i = 1, ncum
[2253]	1729	IF(k>=icb(i)) THEN
	1730	!! IF (k>=(nk(i)+1)) THEN
	1731	!>jyg
[1992]	1732	tca = tp(i, k) - t0
	1733	IF (tca>=0.0) THEN
	1734	elacrit = elcrit
	1735	ELSE
	1736	elacrit = elcrit*(1.0-tca/tlcrit)
	1737	END IF
	1738	elacrit = max(elacrit, 0.0)
	1739	ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
	1740	ep(i, k) = max(ep(i,k), 0.0)
	1741	ep(i, k) = min(ep(i,k), epmax)
[2253]	1742	!! sigp(i, k) = spfac ! jyg
	1743	END IF ! (k>=icb(i))
[1992]	1744	END DO
	1745	END DO
	1746	END IF
[2253]	1747	!
[3605]	1748	! =========================================================================
	1749	IF (prt_level >= 10) THEN
	1750	print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
	1751	(k, tp(1,k), tvp(1,k), k = 1,nl)
	1752	ENDIF
	1753	!
[2007]	1754	! =====================================================================
	1755	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
	1756	! --- VIRTUAL TEMPERATURE
	1757	! =====================================================================
[1849]	1758
[2007]	1759	! dans convect3, tvp est calcule en une seule fois, et sans retirer
	1760	! l'eau condensee (~> reversible CAPE)
[1849]	1761
[2007]	1762	! ori do 340 k=minorig+1,nl
	1763	! ori do 330 i=1,ncum
	1764	! ori if(k.ge.(icb(i)+1))then
	1765	! ori tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
	1766	! oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
	1767	! oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
	1768	! ori endif
	1769	! ori 330 continue
	1770	! ori 340 continue
[524]	1771
[2007]	1772	! ori do 350 i=1,ncum
	1773	! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
	1774	! ori 350 continue
[524]	1775
[2007]	1776	DO i = 1, ncum ! convect3
	1777	tp(i, nlp) = tp(i, nl) ! convect3
	1778	END DO ! convect3
[524]	1779
[2007]	1780	! =====================================================================
	1781	! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
	1782	! =====================================================================
[524]	1783
[2007]	1784	! -- this is for convect3 only:
[524]	1785
[2007]	1786	! first estimate of buoyancy:
[879]	1787
[2257]	1788	!jyg : k-loop outside i-loop (07042015)
	1789	DO k = 1, nl
	1790	DO i = 1, ncum
[1992]	1791	buoy(i, k) = tvp(i, k) - tv(i, k)
	1792	END DO
	1793	END DO
[524]	1794
[2007]	1795	! set buoyancy=buoybase for all levels below base
	1796	! for safety, set buoy(icb)=buoybase
[524]	1797
[2257]	1798	!jyg : k-loop outside i-loop (07042015)
	1799	DO k = 1, nl
	1800	DO i = 1, ncum
[1992]	1801	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
	1802	buoy(i, k) = buoybase(i)
	1803	END IF
	1804	END DO
[2257]	1805	END DO
	1806	DO i = 1, ncum
[2007]	1807	! buoy(icb(i),k)=buoybase(i)
[1992]	1808	buoy(i, icb(i)) = buoybase(i)
	1809	END DO
[524]	1810
[2007]	1811	! -- end convect3
[524]	1812
[2007]	1813	! =====================================================================
	1814	! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
	1815	! --- LEVEL OF NEUTRAL BUOYANCY
	1816	! =====================================================================
[524]	1817
[2007]	1818	! -- this is for convect3 only:
[524]	1819
[1992]	1820	DO i = 1, ncum
	1821	inb(i) = nl - 1
	1822	iposit(i) = nl
	1823	END DO
[524]	1824
	1825
[2007]	1826	! -- iposit(i) = first level, above icb, with positive buoyancy
[1992]	1827	DO k = 1, nl - 1
	1828	DO i = 1, ncum
	1829	IF (k>=icb(i) .AND. buoy(i,k)>0.) THEN
	1830	iposit(i) = min(iposit(i), k)
	1831	END IF
	1832	END DO
	1833	END DO
[1849]	1834
[1992]	1835	DO i = 1, ncum
	1836	IF (iposit(i)==nl) THEN
	1837	iposit(i) = icb(i)
	1838	END IF
	1839	END DO
[1849]	1840
[1992]	1841	DO k = 1, nl - 1
	1842	DO i = 1, ncum
	1843	IF ((k>=iposit(i)) .AND. (buoy(i,k)<dtovsh)) THEN
	1844	inb(i) = min(inb(i), k)
	1845	END IF
	1846	END DO
	1847	END DO
[1849]	1848
[2420]	1849	!CR fix computation of inb
	1850	!keep flag or modify in all cases?
	1851	IF (iflag_mix_adiab.eq.1) THEN
	1852	DO i = 1, ncum
	1853	cape(i)=0.
	1854	inb(i)=icb(i)+1
	1855	ENDDO
	1856
	1857	DO k = 2, nl
	1858	DO i = 1, ncum
	1859	IF ((k>=iposit(i))) THEN
	1860	deltap = min(plcl(i), ph(i,k-1)) - min(plcl(i), ph(i,k))
	1861	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
	1862	IF (cape(i).gt.0.) THEN
	1863	inb(i) = max(inb(i), k)
	1864	END IF
	1865	ENDIF
	1866	ENDDO
	1867	ENDDO
	1868
	1869	! DO i = 1, ncum
	1870	! print*,"inb",inb(i)
	1871	! ENDDO
	1872
	1873	endif
	1874
[2007]	1875	! -- end convect3
[1849]	1876
[2007]	1877	! ori do 510 i=1,ncum
	1878	! ori cape(i)=0.0
	1879	! ori capem(i)=0.0
	1880	! ori inb(i)=icb(i)+1
	1881	! ori inb1(i)=inb(i)
	1882	! ori 510 continue
[524]	1883
[2007]	1884	! Originial Code
[524]	1885
[2007]	1886	! do 530 k=minorig+1,nl-1
	1887	! do 520 i=1,ncum
	1888	! if(k.ge.(icb(i)+1))then
	1889	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1890	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1891	! cape(i)=cape(i)+by
	1892	! if(by.ge.0.0)inb1(i)=k+1
	1893	! if(cape(i).gt.0.0)then
	1894	! inb(i)=k+1
	1895	! capem(i)=cape(i)
	1896	! endif
	1897	! endif
	1898	!520 continue
	1899	!530 continue
	1900	! do 540 i=1,ncum
	1901	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
	1902	! cape(i)=capem(i)+byp
	1903	! defrac=capem(i)-cape(i)
	1904	! defrac=max(defrac,0.001)
	1905	! frac(i)=-cape(i)/defrac
	1906	! frac(i)=min(frac(i),1.0)
	1907	! frac(i)=max(frac(i),0.0)
	1908	!540 continue
[524]	1909
[2007]	1910	! K Emanuel fix
[524]	1911
[2007]	1912	! call zilch(byp,ncum)
	1913	! do 530 k=minorig+1,nl-1
	1914	! do 520 i=1,ncum
	1915	! if(k.ge.(icb(i)+1))then
	1916	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1917	! cape(i)=cape(i)+by
	1918	! if(by.ge.0.0)inb1(i)=k+1
	1919	! if(cape(i).gt.0.0)then
	1920	! inb(i)=k+1
	1921	! capem(i)=cape(i)
	1922	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1923	! endif
	1924	! endif
	1925	!520 continue
	1926	!530 continue
	1927	! do 540 i=1,ncum
	1928	! inb(i)=max(inb(i),inb1(i))
	1929	! cape(i)=capem(i)+byp(i)
	1930	! defrac=capem(i)-cape(i)
	1931	! defrac=max(defrac,0.001)
	1932	! frac(i)=-cape(i)/defrac
	1933	! frac(i)=min(frac(i),1.0)
	1934	! frac(i)=max(frac(i),0.0)
	1935	!540 continue
[524]	1936
[2007]	1937	! J Teixeira fix
[524]	1938
[2007]	1939	! ori call zilch(byp,ncum)
	1940	! ori do 515 i=1,ncum
	1941	! ori lcape(i)=.true.
	1942	! ori 515 continue
	1943	! ori do 530 k=minorig+1,nl-1
	1944	! ori do 520 i=1,ncum
	1945	! ori if(cape(i).lt.0.0)lcape(i)=.false.
	1946	! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
	1947	! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1948	! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1949	! ori cape(i)=cape(i)+by
	1950	! ori if(by.ge.0.0)inb1(i)=k+1
	1951	! ori if(cape(i).gt.0.0)then
	1952	! ori inb(i)=k+1
	1953	! ori capem(i)=cape(i)
	1954	! ori endif
	1955	! ori endif
	1956	! ori 520 continue
	1957	! ori 530 continue
	1958	! ori do 540 i=1,ncum
	1959	! ori cape(i)=capem(i)+byp(i)
	1960	! ori defrac=capem(i)-cape(i)
	1961	! ori defrac=max(defrac,0.001)
	1962	! ori frac(i)=-cape(i)/defrac
	1963	! ori frac(i)=min(frac(i),1.0)
	1964	! ori frac(i)=max(frac(i),0.0)
	1965	! ori 540 continue
[524]	1966
[2761]	1967	! --------------------------------------------------------------------
	1968	! Prevent convection when top is too hot
	1969	! --------------------------------------------------------------------
	1970	DO i = 1,ncum
	1971	IF (t(i,inb(i)) > T_top_max) iflag(i) = 10
	1972	ENDDO
	1973
[2007]	1974	! =====================================================================
	1975	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
	1976	! =====================================================================
[524]	1977
[2393]	1978	DO k = 1, nl
[1992]	1979	DO i = 1, ncum
	1980	hp(i, k) = h(i, k)
	1981	END DO
	1982	END DO
[524]	1983
[2257]	1984	!jyg : cvflag_ice test outside the loops (07042015)
	1985	!
	1986	IF (cvflag_ice) THEN
	1987	!
[3605]	1988	IF (cvflag_prec_eject) THEN
	1989	!! DO k = minorig + 1, nl
	1990	!! DO i = 1, ncum
	1991	!! IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
	1992	!! frac_s(i,k) = qi(i,k)/max(clw(i,k),smallestreal)
	1993	!! frac_s(i,k) = 1. - (qpl(i,k)+(1.-frac_s(i,k))*qcld(i,k))/max(clw(i,k),smallestreal)
	1994	!! END IF
	1995	!! END DO
	1996	!! END DO
	1997	ELSE ! (cvflag_prec_eject)
[2257]	1998	DO k = minorig + 1, nl
	1999	DO i = 1, ncum
	2000	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
[3605]	2001	!jyg< frac computation moved to beginning of cv3_undilute2.
	2002	! kept here for compatibility test with CMip6 version
	2003	frac_s(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
	2004	frac_s(i, k) = min(max(frac_s(i,k),0.0), 1.0)
[2257]	2005	END IF
	2006	END DO
	2007	END DO
[3605]	2008	ENDIF ! (cvflag_prec_eject) ELSE
	2009	DO k = minorig + 1, nl
[2376]	2010	DO i = 1, ncum
[3605]	2011	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
	2012	!! hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac_s(i,k)lf(i,k))* & !!jygprl
	2013	!! ep(i, k)*clw(i, k) !!jygprl
	2014	hp(i, k) = hla(i,k-1) + (lv(i,k)+(cpd-cpv)t(i,k)+frac_s(i,k)lf(i,k))* & !!jygprl
	2015	ep(i, k)*clw(i, k) !!jygprl
[2376]	2016	END IF
	2017	END DO
	2018	END DO
[2257]	2019	!
[3605]	2020	ELSE ! (cvflag_ice)
[2257]	2021	!
	2022	DO k = minorig + 1, nl
	2023	DO i = 1, ncum
	2024	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
[2902]	2025	!jyg< (energy conservation tests)
	2026	!! hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)tp(i,k))ep(i, k)*clw(i, k)
	2027	!! hp(i, k) = ( hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k) ) / &
	2028	!! (1. - ep(i,k)*clw(i,k))
	2029	!! hp(i, k) = ( hnk(i) + (lv(i,k)+(cpd-cl)t(i,k))ep(i, k)*clw(i, k) ) / &
	2030	!! (1. - ep(i,k)*clw(i,k))
[1992]	2031	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k)
	2032	END IF
[2257]	2033	END DO
[1992]	2034	END DO
[2257]	2035	!
	2036	END IF ! (cvflag_ice)
[524]	2037
[1992]	2038	RETURN
	2039	END SUBROUTINE cv3_undilute2
[524]	2040
[2007]	2041	SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
	2042	pbase, p, ph, tv, buoy, &
	2043	sig, w0, cape, m, iflag)
[1992]	2044	IMPLICIT NONE
[524]	2045
[2007]	2046	! ===================================================================
	2047	! --- CLOSURE OF CONVECT3
	2048	!
	2049	! vectorization: S. Bony
	2050	! ===================================================================
[524]	2051
[1992]	2052	include "cvthermo.h"
	2053	include "cv3param.h"
[524]	2054
[2007]	2055	!input:
[1992]	2056	INTEGER ncum, nd, nloc
	2057	INTEGER icb(nloc), inb(nloc)
	2058	REAL pbase(nloc)
	2059	REAL p(nloc, nd), ph(nloc, nd+1)
	2060	REAL tv(nloc, nd), buoy(nloc, nd)
[524]	2061
[2007]	2062	!input/output:
[1992]	2063	REAL sig(nloc, nd), w0(nloc, nd)
	2064	INTEGER iflag(nloc)
[524]	2065
[2007]	2066	!output:
[1992]	2067	REAL cape(nloc)
	2068	REAL m(nloc, nd)
[524]	2069
[2007]	2070	!local variables:
[1992]	2071	INTEGER i, j, k, icbmax
	2072	REAL deltap, fac, w, amu
	2073	REAL dtmin(nloc, nd), sigold(nloc, nd)
	2074	REAL cbmflast(nloc)
[524]	2075
[776]	2076
[2007]	2077	! -------------------------------------------------------
	2078	! -- Initialization
	2079	! -------------------------------------------------------
[524]	2080
[1992]	2081	DO k = 1, nl
	2082	DO i = 1, ncum
	2083	m(i, k) = 0.0
	2084	END DO
	2085	END DO
[524]	2086
[2007]	2087	! -------------------------------------------------------
	2088	! -- Reset sig(i) and w0(i) for i>inb and i<icb
	2089	! -------------------------------------------------------
[524]	2090
[2007]	2091	! update sig and w0 above LNB:
[879]	2092
[1992]	2093	DO k = 1, nl - 1
	2094	DO i = 1, ncum
	2095	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
[2007]	2096	sig(i, k) = beta*sig(i, k) + &
	2097	2.alphabuoy(i, inb(i))*abs(buoy(i,inb(i)))
[1992]	2098	sig(i, k) = amax1(sig(i,k), 0.0)
	2099	w0(i, k) = beta*w0(i, k)
	2100	END IF
	2101	END DO
	2102	END DO
[1494]	2103
[2007]	2104	! compute icbmax:
[524]	2105
[1992]	2106	icbmax = 2
	2107	DO i = 1, ncum
	2108	icbmax = max(icbmax, icb(i))
	2109	END DO
[524]	2110
[2007]	2111	! update sig and w0 below cloud base:
[1494]	2112
[1992]	2113	DO k = 1, icbmax
	2114	DO i = 1, ncum
	2115	IF (k<=icb(i)) THEN
[2007]	2116	sig(i, k) = beta*sig(i, k) - &
	2117	2.alphabuoy(i, icb(i))*buoy(i, icb(i))
[1992]	2118	sig(i, k) = max(sig(i,k), 0.0)
	2119	w0(i, k) = beta*w0(i, k)
	2120	END IF
	2121	END DO
	2122	END DO
[524]	2123
[2007]	2124	!! if(inb.lt.(nl-1))then
	2125	!! do 85 i=inb+1,nl-1
	2126	!! sig(i)=betasig(i)+2.alphabuoy(inb)
	2127	!! 1 abs(buoy(inb))
	2128	!! sig(i)=max(sig(i),0.0)
	2129	!! w0(i)=beta*w0(i)
	2130	!! 85 continue
	2131	!! end if
[524]	2132
[2007]	2133	!! do 87 i=1,icb
	2134	!! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
	2135	!! sig(i)=max(sig(i),0.0)
	2136	!! w0(i)=beta*w0(i)
	2137	!! 87 continue
[524]	2138
[2007]	2139	! -------------------------------------------------------------
	2140	! -- Reset fractional areas of updrafts and w0 at initial time
	2141	! -- and after 10 time steps of no convection
	2142	! -------------------------------------------------------------
[524]	2143
[1992]	2144	DO k = 1, nl - 1
	2145	DO i = 1, ncum
	2146	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
	2147	sig(i, k) = 0.0
	2148	w0(i, k) = 0.0
	2149	END IF
	2150	END DO
	2151	END DO
[524]	2152
[2007]	2153	! -------------------------------------------------------------
	2154	! -- Calculate convective available potential energy (cape),
	2155	! -- vertical velocity (w), fractional area covered by
	2156	! -- undilute updraft (sig), and updraft mass flux (m)
	2157	! -------------------------------------------------------------
[1849]	2158
[1992]	2159	DO i = 1, ncum
	2160	cape(i) = 0.0
	2161	END DO
[1849]	2162
[2007]	2163	! compute dtmin (minimum buoyancy between ICB and given level k):
[1849]	2164
[1992]	2165	DO i = 1, ncum
	2166	DO k = 1, nl
	2167	dtmin(i, k) = 100.0
	2168	END DO
	2169	END DO
[524]	2170
[1992]	2171	DO i = 1, ncum
	2172	DO k = 1, nl
	2173	DO j = minorig, nl
[2007]	2174	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
[1992]	2175	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
	2176	END IF
	2177	END DO
	2178	END DO
	2179	END DO
[1849]	2180
[2007]	2181	! the interval on which cape is computed starts at pbase :
[1849]	2182
[1992]	2183	DO k = 1, nl
	2184	DO i = 1, ncum
[1849]	2185
[1992]	2186	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
[1849]	2187
[1992]	2188	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
	2189	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
	2190	cape(i) = amax1(0.0, cape(i))
	2191	sigold(i, k) = sig(i, k)
[1849]	2192
[2007]	2193	! dtmin(i,k)=100.0
	2194	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
	2195	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
	2196	! 97 continue
[1849]	2197
[1992]	2198	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
	2199	sig(i, k) = max(sig(i,k), 0.0)
	2200	sig(i, k) = amin1(sig(i,k), 0.01)
	2201	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
	2202	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
	2203	amu = 0.5(sig(i,k)+sigold(i,k))w
	2204	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
	2205	w0(i, k) = w
	2206	END IF
[1849]	2207
[1992]	2208	END DO
	2209	END DO
[1849]	2210
[1992]	2211	DO i = 1, ncum
	2212	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
[2007]	2213	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]	2214	sig(i, icb(i)) = sig(i, icb(i)+1)
	2215	sig(i, icb(i)-1) = sig(i, icb(i))
	2216	END DO
[1849]	2217
[2007]	2218	! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
	2219	! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
	2220	! ccc the final mass flux (cbmflast) is greater than the target mass flux
	2221	! ccc (cbmf) ??).
	2222	! cc
	2223	! c do i = 1,ncum
	2224	! c cbmflast(i) = 0.
	2225	! c enddo
	2226	! cc
	2227	! c do k= 1,nl
	2228	! c do i = 1,ncum
	2229	! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
	2230	! c cbmflast(i) = cbmflast(i)+M(i,k)
	2231	! c ENDIF
	2232	! c enddo
	2233	! c enddo
	2234	! cc
	2235	! c do i = 1,ncum
	2236	! c IF (cbmflast(i) .lt. 1.e-6) THEN
	2237	! c iflag(i) = 3
	2238	! c ENDIF
	2239	! c enddo
	2240	! cc
	2241	! c do k= 1,nl
	2242	! c do i = 1,ncum
	2243	! c IF (iflag(i) .ge. 3) THEN
	2244	! c M(i,k) = 0.
	2245	! c sig(i,k) = 0.
	2246	! c w0(i,k) = 0.
	2247	! c ENDIF
	2248	! c enddo
	2249	! c enddo
	2250	! cc
	2251	!! cape=0.0
	2252	!! do 98 i=icb+1,inb
	2253	!! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
	2254	!! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
	2255	!! dcape=rrdbuoy(i-1)deltap/p(i-1)
	2256	!! dlnp=deltap/p(i-1)
	2257	!! cape=max(0.0,cape)
	2258	!! sigold=sig(i)
[1849]	2259
[2007]	2260	!! dtmin=100.0
	2261	!! do 97 j=icb,i-1
	2262	!! dtmin=amin1(dtmin,buoy(j))
	2263	!! 97 continue
[1849]	2264
[2007]	2265	!! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
	2266	!! sig(i)=max(sig(i),0.0)
	2267	!! sig(i)=amin1(sig(i),0.01)
	2268	!! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
	2269	!! w=(1.-beta)facsqrt(cape)+beta*w0(i)
	2270	!! amu=0.5(sig(i)+sigold)w
	2271	!! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
	2272	!! w0(i)=w
	2273	!! 98 continue
	2274	!! w0(icb)=0.5*w0(icb+1)
	2275	!! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
	2276	!! sig(icb)=sig(icb+1)
	2277	!! sig(icb-1)=sig(icb)
[1849]	2278
[1992]	2279	RETURN
	2280	END SUBROUTINE cv3_closure
[1849]	2281
[2007]	2282	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
	2283	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
	2284	unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
	2285	ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
[1992]	2286	IMPLICIT NONE
[1849]	2287
[2007]	2288	! ---------------------------------------------------------------------
	2289	! a faire:
	2290	! - vectorisation de la partie normalisation des flux (do 789...)
	2291	! ---------------------------------------------------------------------
[524]	2292
[1992]	2293	include "cvthermo.h"
	2294	include "cv3param.h"
	2295	include "cvflag.h"
[524]	2296
[2007]	2297	!inputs:
[2253]	2298	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
	2299	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
	2300	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
	2301	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
	2302	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	2303	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
	2304	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
	2305	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
	2306	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, h, hp
	2307	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf, frac
	2308	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp, ep, clw
	2309	REAL, DIMENSION (nloc, na), INTENT (IN) :: m ! input of convect3
[524]	2310
[2007]	2311	!outputs:
[2253]	2312	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: ment, qent
	2313	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
	2314	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: sij, elij
	2315	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
	2316	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: ments, qents
	2317	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
[524]	2318
[2007]	2319	!local variables:
[1992]	2320	INTEGER i, j, k, il, im, jm
	2321	INTEGER num1, num2
	2322	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
	2323	REAL alt, smid, sjmin, sjmax, delp, delm
	2324	REAL asij(nloc), smax(nloc), scrit(nloc)
	2325	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
[2253]	2326	REAL sigij(nloc, nd, nd)
[1992]	2327	REAL wgh
	2328	REAL zm(nloc, na)
	2329	LOGICAL lwork(nloc)
[524]	2330
[2007]	2331	! =====================================================================
	2332	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	2333	! =====================================================================
[524]	2334
[2007]	2335	! ori do 360 i=1,ncum*nlp
[1992]	2336	DO j = 1, nl
	2337	DO i = 1, ncum
	2338	nent(i, j) = 0
[2007]	2339	! in convect3, m is computed in cv3_closure
	2340	! ori m(i,1)=0.0
[1992]	2341	END DO
	2342	END DO
[524]	2343
[2007]	2344	! ori do 400 k=1,nlp
	2345	! ori do 390 j=1,nlp
[1992]	2346	DO j = 1, nl
[2459]	2347	DO k = 1, nl
	2348	DO i = 1, ncum
	2349	qent(i, k, j) = rr(i, j)
	2350	uent(i, k, j) = u(i, j)
	2351	vent(i, k, j) = v(i, j)
	2352	elij(i, k, j) = 0.0
	2353	!ym ment(i,k,j)=0.0
	2354	!ym sij(i,k,j)=0.0
	2355	END DO
	2356	END DO
[1992]	2357	END DO
[524]	2358
[2007]	2359	!ym
[1992]	2360	ment(1:ncum, 1:nd, 1:nd) = 0.0
	2361	sij(1:ncum, 1:nd, 1:nd) = 0.0
[524]	2362
[2007]	2363	!AC! do k=1,ntra
	2364	!AC! do j=1,nd ! instead nlp
	2365	!AC! do i=1,nd ! instead nlp
	2366	!AC! do il=1,ncum
	2367	!AC! traent(il,i,j,k)=tra(il,j,k)
	2368	!AC! enddo
	2369	!AC! enddo
	2370	!AC! enddo
	2371	!AC! enddo
[1992]	2372	zm(:, :) = 0.
[524]	2373
[2007]	2374	! =====================================================================
	2375	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	2376	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	2377	! --- FRACTION (sij)
	2378	! =====================================================================
[524]	2379
[1992]	2380	DO i = minorig + 1, nl
[524]	2381
[1992]	2382	DO j = minorig, nl
	2383	DO il = 1, ncum
[2007]	2384	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
[524]	2385
[1992]	2386	rti = qnk(il) - ep(il, i)*clw(il, i)
	2387	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
[524]	2388
	2389
[1992]	2390	IF (cvflag_ice) THEN
[2007]	2391	! print*,cvflag_ice,'cvflag_ice dans do 700'
[1992]	2392	IF (t(il,j)<=263.15) THEN
[2007]	2393	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
	2394	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
[1992]	2395	END IF
	2396	END IF
[524]	2397
[1992]	2398	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
	2399	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
	2400	dei = denom
	2401	IF (abs(dei)<0.01) dei = 0.01
	2402	sij(il, i, j) = anum/dei
	2403	sij(il, i, i) = 1.0
	2404	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
	2405	altem = altem/bf2
	2406	cwat = clw(il, j)*(1.-ep(il,j))
	2407	stemp = sij(il, i, j)
	2408	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
[524]	2409
[1992]	2410	IF (cvflag_ice) THEN
[2007]	2411	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
[1992]	2412	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
	2413	ELSE
	2414	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
	2415	denom = denom + lv(il, j)*(rr(il,i)-rti)
	2416	END IF
[524]	2417
[1992]	2418	IF (abs(denom)<0.01) denom = 0.01
	2419	sij(il, i, j) = anum/denom
[2007]	2420	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
[1992]	2421	altem = altem - (bf2-1.)*cwat
	2422	END IF
	2423	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
	2424	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
[2007]	2425	uent(il, i, j) = sij(il, i, j)u(il, i) + (1.-sij(il,i,j))unk(il)
	2426	vent(il, i, j) = sij(il, i, j)v(il, i) + (1.-sij(il,i,j))vnk(il)
	2427	!!!! do k=1,ntra
	2428	!!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	2429	!!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	2430	!!!! end do
[1992]	2431	elij(il, i, j) = altem
	2432	elij(il, i, j) = max(0.0, elij(il,i,j))
	2433	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
	2434	nent(il, i) = nent(il, i) + 1
	2435	END IF
	2436	sij(il, i, j) = max(0.0, sij(il,i,j))
	2437	sij(il, i, j) = amin1(1.0, sij(il,i,j))
	2438	END IF ! new
	2439	END DO
	2440	END DO
[524]	2441
[2007]	2442	!AC! do k=1,ntra
	2443	!AC! do j=minorig,nl
	2444	!AC! do il=1,ncum
	2445	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
	2446	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
	2447	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	2448	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	2449	!AC! endif
	2450	!AC! enddo
	2451	!AC! enddo
	2452	!AC! enddo
[524]	2453
	2454
[2007]	2455	! * if no air can entrain at level i assume that updraft detrains *
	2456	! * at that level and calculate detrained air flux and properties *
[524]	2457
	2458
[2007]	2459	! @ do 170 i=icb(il),inb(il)
[524]	2460
[1992]	2461	DO il = 1, ncum
	2462	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
[2007]	2463	! @ if(nent(il,i).eq.0)then
[1992]	2464	ment(il, i, i) = m(il, i)
	2465	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	2466	uent(il, i, i) = unk(il)
	2467	vent(il, i, i) = vnk(il)
	2468	elij(il, i, i) = clw(il, i)
[2007]	2469	! MAF sij(il,i,i)=1.0
[1992]	2470	sij(il, i, i) = 0.0
	2471	END IF
	2472	END DO
	2473	END DO
[879]	2474
[2007]	2475	!AC! do j=1,ntra
	2476	!AC! do i=minorig+1,nl
	2477	!AC! do il=1,ncum
	2478	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
	2479	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
	2480	!AC! endif
	2481	!AC! enddo
	2482	!AC! enddo
	2483	!AC! enddo
[879]	2484
[1992]	2485	DO j = minorig, nl
	2486	DO i = minorig, nl
	2487	DO il = 1, ncum
[2007]	2488	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
[1992]	2489	sigij(il, i, j) = sij(il, i, j)
	2490	END IF
	2491	END DO
	2492	END DO
	2493	END DO
[2007]	2494	! @ enddo
[879]	2495
[2007]	2496	! @170 continue
[524]	2497
[2007]	2498	! =====================================================================
	2499	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	2500	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	2501	! =====================================================================
[970]	2502
[1992]	2503	CALL zilch(asum, nloc*nd)
	2504	CALL zilch(csum, nloc*nd)
	2505	CALL zilch(csum, nloc*nd)
[524]	2506
[1992]	2507	DO il = 1, ncum
	2508	lwork(il) = .FALSE.
	2509	END DO
[524]	2510
[1992]	2511	DO i = minorig + 1, nl
[524]	2512
[1992]	2513	num1 = 0
	2514	DO il = 1, ncum
	2515	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
	2516	END DO
	2517	IF (num1<=0) GO TO 789
[524]	2518
[879]	2519
[1992]	2520	DO il = 1, ncum
	2521	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	2522	lwork(il) = (nent(il,i)/=0)
	2523	qp = qnk(il) - ep(il, i)*clw(il, i)
[524]	2524
[1992]	2525	IF (cvflag_ice) THEN
[524]	2526
[2007]	2527	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
	2528	(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
	2529	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
	2530	(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
[1992]	2531	ELSE
[879]	2532
[1992]	2533	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
[2007]	2534	(cpv-cpd)t(il, i)(qp-rr(il,i))
[1992]	2535	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
[2007]	2536	(cpd-cpv)t(il, i)(rr(il,i)-qp)
[1992]	2537	END IF
[524]	2538
[1992]	2539	IF (abs(denom)<0.01) denom = 0.01
	2540	scrit(il) = anum/denom
	2541	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
	2542	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
	2543	smax(il) = 0.0
	2544	asij(il) = 0.0
	2545	END IF
	2546	END DO
[524]	2547
[1992]	2548	DO j = nl, minorig, -1
[524]	2549
[1992]	2550	num2 = 0
	2551	DO il = 1, ncum
[2007]	2552	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	2553	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	2554	lwork(il)) num2 = num2 + 1
[1992]	2555	END DO
	2556	IF (num2<=0) GO TO 175
[524]	2557
[1992]	2558	DO il = 1, ncum
[2007]	2559	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	2560	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	2561	lwork(il)) THEN
[524]	2562
[1992]	2563	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
	2564	wgh = 1.0
	2565	IF (j>i) THEN
	2566	sjmax = max(sij(il,i,j+1), smax(il))
	2567	sjmax = amin1(sjmax, scrit(il))
	2568	smax(il) = max(sij(il,i,j), smax(il))
	2569	sjmin = max(sij(il,i,j-1), smax(il))
	2570	sjmin = amin1(sjmin, scrit(il))
	2571	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
	2572	smid = amin1(sij(il,i,j), scrit(il))
	2573	ELSE
	2574	sjmax = max(sij(il,i,j+1), scrit(il))
	2575	smid = max(sij(il,i,j), scrit(il))
	2576	sjmin = 0.0
	2577	IF (j>1) sjmin = sij(il, i, j-1)
	2578	sjmin = max(sjmin, scrit(il))
	2579	END IF
	2580	delp = abs(sjmax-smid)
	2581	delm = abs(sjmin-smid)
	2582	asij(il) = asij(il) + wgh*(delp+delm)
	2583	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
	2584	END IF
	2585	END IF
	2586	END DO
[524]	2587
[1992]	2588	175 END DO
[524]	2589
[1992]	2590	DO il = 1, ncum
	2591	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	2592	asij(il) = max(1.0E-16, asij(il))
	2593	asij(il) = 1.0/asij(il)
	2594	asum(il, i) = 0.0
	2595	bsum(il, i) = 0.0
	2596	csum(il, i) = 0.0
	2597	END IF
	2598	END DO
[524]	2599
[1992]	2600	DO j = minorig, nl
	2601	DO il = 1, ncum
[2007]	2602	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2603	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[1992]	2604	ment(il, i, j) = ment(il, i, j)*asij(il)
	2605	END IF
	2606	END DO
	2607	END DO
[524]	2608
[1992]	2609	DO j = minorig, nl
	2610	DO il = 1, ncum
[2007]	2611	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2612	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[1992]	2613	asum(il, i) = asum(il, i) + ment(il, i, j)
	2614	ment(il, i, j) = ment(il, i, j)*sig(il, j)
	2615	bsum(il, i) = bsum(il, i) + ment(il, i, j)
	2616	END IF
	2617	END DO
	2618	END DO
[1849]	2619
[1992]	2620	DO il = 1, ncum
	2621	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	2622	bsum(il, i) = max(bsum(il,i), 1.0E-16)
	2623	bsum(il, i) = 1.0/bsum(il, i)
	2624	END IF
	2625	END DO
[1849]	2626
[1992]	2627	DO j = minorig, nl
	2628	DO il = 1, ncum
[2007]	2629	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2630	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[1992]	2631	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
	2632	END IF
	2633	END DO
	2634	END DO
[879]	2635
[1992]	2636	DO j = minorig, nl
	2637	DO il = 1, ncum
[2007]	2638	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2639	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[1992]	2640	csum(il, i) = csum(il, i) + ment(il, i, j)
	2641	END IF
	2642	END DO
	2643	END DO
[1849]	2644
[1992]	2645	DO il = 1, ncum
	2646	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2647	csum(il,i)<m(il,i)) THEN
	2648	nent(il, i) = 0
	2649	ment(il, i, i) = m(il, i)
	2650	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	2651	uent(il, i, i) = unk(il)
	2652	vent(il, i, i) = vnk(il)
	2653	elij(il, i, i) = clw(il, i)
[2007]	2654	! MAF sij(il,i,i)=1.0
[1992]	2655	sij(il, i, i) = 0.0
	2656	END IF
	2657	END DO ! il
[1849]	2658
[2007]	2659	!AC! do j=1,ntra
	2660	!AC! do il=1,ncum
	2661	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	2662	!AC! : .and. csum(il,i).lt.m(il,i) ) then
	2663	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
	2664	!AC! endif
	2665	!AC! enddo
	2666	!AC! enddo
[1992]	2667	789 END DO
[879]	2668
[2007]	2669	! MAF: renormalisation de MENT
[1992]	2670	CALL zilch(zm, nloc*na)
[2393]	2671	DO jm = 1, nl
[2459]	2672	DO im = 1, nl
	2673	DO il = 1, ncum
	2674	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
	2675	END DO
	2676	END DO
[1992]	2677	END DO
[524]	2678
[2393]	2679	DO jm = 1, nl
	2680	DO im = 1, nl
[1992]	2681	DO il = 1, ncum
	2682	IF (zm(il,im)/=0.) THEN
	2683	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
	2684	END IF
	2685	END DO
	2686	END DO
	2687	END DO
[524]	2688
[2459]	2689	DO jm = 1, nl
	2690	DO im = 1, nl
	2691	DO il = 1, ncum
	2692	qents(il, im, jm) = qent(il, im, jm)
	2693	ments(il, im, jm) = ment(il, im, jm)
	2694	END DO
	2695	END DO
	2696	END DO
[524]	2697
[1992]	2698	RETURN
	2699	END SUBROUTINE cv3_mixing
[879]	2700
[2007]	2701	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
	2702	t, rr, rs, gz, u, v, tra, p, ph, &
[3605]	2703	th, tv, lv, lf, cpn, ep, sigp, clw, frac_s, qpreca, frac_a, qta , & !!jygprl
[2007]	2704	m, ment, elij, delt, plcl, coef_clos, &
	2705	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
	2706	faci, b, sigd, &
[3605]	2707	wdtrainA, wdtrainS, wdtrainM) ! RomP
[2393]	2708	USE print_control_mod, ONLY: prt_level, lunout
[1992]	2709	IMPLICIT NONE
[879]	2710
	2711
[1992]	2712	include "cvthermo.h"
	2713	include "cv3param.h"
	2714	include "cvflag.h"
[2287]	2715	include "nuage.h"
[524]	2716
[2007]	2717	!inputs:
[2393]	2718	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
	2719	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
	2720	REAL, INTENT(IN) :: delt
	2721	REAL, DIMENSION (nloc), INTENT (IN) :: plcl
	2722	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
	2723	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz
	2724	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
[3605]	2725	REAL, DIMENSION (nloc, nd, ntra), INTENT(IN) :: tra
	2726	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
	2727	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	2728	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, sigp, clw !adiab ascent shedding
	2729	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_s !ice fraction in adiab ascent shedding !!jygprl
	2730	REAL, DIMENSION (nloc, na), INTENT (IN) :: qpreca !adiab ascent precip !!jygprl
	2731	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_a !ice fraction in adiab ascent precip !!jygprl
	2732	REAL, DIMENSION (nloc, na), INTENT (IN) :: qta !adiab ascent specific total water !!jygprl
[2393]	2733	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tv, lv, cpn
	2734	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
	2735	REAL, DIMENSION (nloc, na), INTENT (IN) :: m
	2736	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: ment, elij
	2737	REAL, DIMENSION (nloc), INTENT (IN) :: coef_clos
[524]	2738
[2007]	2739	!input/output
[2393]	2740	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag(nloc)
[524]	2741
[2007]	2742	!outputs:
[2393]	2743	REAL, DIMENSION (nloc, na), INTENT (OUT) :: mp, rp, up, vp
	2744	REAL, DIMENSION (nloc, na), INTENT (OUT) :: water, evap, wt
[3605]	2745	REAL, DIMENSION (nloc, na), INTENT (OUT) :: ice, fondue
	2746	REAL, DIMENSION (nloc, na), INTENT (OUT) :: faci ! ice fraction in precipitation
[2393]	2747	REAL, DIMENSION (nloc, na, ntra), INTENT (OUT) :: trap
	2748	REAL, DIMENSION (nloc, na), INTENT (OUT) :: b
	2749	REAL, DIMENSION (nloc), INTENT (OUT) :: sigd
[2007]	2750	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
	2751	! de l ascendance adiabatique et des flux melanges Pa et Pm.
	2752	! Distinction des wdtrain
	2753	! Pa = wdtrainA Pm = wdtrainM
[3605]	2754	REAL, DIMENSION (nloc, na), INTENT (OUT) :: wdtrainA, wdtrainS, wdtrainM
[879]	2755
[2007]	2756	!local variables
[1992]	2757	INTEGER i, j, k, il, num1, ndp1
[3605]	2758	REAL smallestreal
[1992]	2759	REAL tinv, delti, coef
	2760	REAL awat, afac, afac1, afac2, bfac
	2761	REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
	2762	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
	2763	REAL ampmax, thaw
	2764	REAL tevap(nloc)
[3605]	2765	REAL, DIMENSION (nloc, na) :: lvcp, lfcp
	2766	REAL, DIMENSION (nloc, na) :: h, hm
	2767	REAL, DIMENSION (nloc, na) :: ma
	2768	REAL, DIMENSION (nloc, na) :: frac ! ice fraction in precipitation source
	2769	REAL, DIMENSION (nloc, na) :: fraci ! provisionnal ice fraction in precipitation
	2770	REAL, DIMENSION (nloc, na) :: prec
[1992]	2771	REAL wdtrain(nloc)
	2772	LOGICAL lwork(nloc), mplus(nloc)
[524]	2773
	2774
[2007]	2775	! ------------------------------------------------------
[2671]	2776	IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1)
[524]	2777
[3605]	2778	smallestreal=tiny(smallestreal)
	2779
[2393]	2780	! =============================
	2781	! --- INITIALIZE OUTPUT ARRAYS
	2782	! =============================
	2783	! (loops up to nl+1)
[2671]	2784	mp(:,:) = 0.
	2785	rp(:,:) = 0.
	2786	up(:,:) = 0.
	2787	vp(:,:) = 0.
	2788	water(:,:) = 0.
	2789	evap(:,:) = 0.
	2790	wt(:,:) = 0.
	2791	ice(:,:) = 0.
	2792	fondue(:,:) = 0.
	2793	faci(:,:) = 0.
	2794	b(:,:) = 0.
	2795	sigd(:) = 0.
	2796	!! RomP >>>
	2797	wdtrainA(:,:) = 0.
[3605]	2798	wdtrainS(:,:) = 0.
[2671]	2799	wdtrainM(:,:) = 0.
	2800	!! RomP <<<
[524]	2801
[2393]	2802	DO i = 1, nlp
[1992]	2803	DO il = 1, ncum
	2804	rp(il, i) = rr(il, i)
	2805	up(il, i) = u(il, i)
	2806	vp(il, i) = v(il, i)
	2807	wt(il, i) = 0.001
[2393]	2808	END DO
	2809	END DO
	2810
	2811	! *** Set the fractionnal area sigd of precipitating downdraughts
	2812	DO il = 1, ncum
	2813	sigd(il) = sigdz*coef_clos(il)
	2814	END DO
	2815
	2816	! =====================================================================
	2817	! --- INITIALIZE VARIOUS ARRAYS AND PARAMETERS USED IN THE COMPUTATIONS
	2818	! =====================================================================
	2819	! (loops up to nl+1)
	2820
	2821	delti = 1./delt
	2822	tinv = 1./3.
	2823
	2824	DO i = 1, nlp
	2825	DO il = 1, ncum
	2826	frac(il, i) = 0.0
	2827	fraci(il, i) = 0.0
	2828	prec(il, i) = 0.0
[1992]	2829	lvcp(il, i) = lv(il, i)/cpn(il, i)
	2830	lfcp(il, i) = lf(il, i)/cpn(il, i)
	2831	END DO
	2832	END DO
[2393]	2833
[2007]	2834	!AC! do k=1,ntra
	2835	!AC! do i=1,nd
	2836	!AC! do il=1,ncum
	2837	!AC! trap(il,i,k)=tra(il,i,k)
	2838	!AC! enddo
	2839	!AC! enddo
	2840	!AC! enddo
[524]	2841
[2007]	2842	! * check whether ep(inb)=0, if so, skip precipitating *
	2843	! * downdraft calculation *
[524]	2844
	2845
[1992]	2846	DO il = 1, ncum
[2007]	2847	!! lwork(il)=.TRUE.
	2848	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
[2671]	2849	!jyg<
	2850	!! lwork(il) = ep(il, inb(il)) >= 0.0001
	2851	lwork(il) = ep(il, inb(il)) >= 0.0001 .AND. iflag(il) <= 2
[1992]	2852	END DO
[524]	2853
[3605]	2854	!
	2855	! Get adiabatic ascent mass flux
	2856	!
	2857	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2858	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
	2859	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2860	!!! Warning : this option leads to water conservation violation
	2861	!!! Expert only
	2862	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2863	DO il = 1, ncum
	2864	ma(il, nlp) = 0.
	2865	ma(il, 1) = 0.
	2866	END DO
[524]	2867
[3605]	2868	DO i = nl, 2, -1
	2869	DO il = 1, ncum
	2870	ma(il, i) = ma(il, i+1)*(1.-qta(il,i))/(1.-qta(il,i-1)) + m(il, i)
	2871	END DO
	2872	END DO
	2873	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2874	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
	2875	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2876	DO il = 1, ncum
	2877	ma(il, nlp) = 0.
	2878	ma(il, 1) = 0.
	2879	END DO
	2880
	2881	DO i = nl, 2, -1
	2882	DO il = 1, ncum
	2883	ma(il, i) = ma(il, i+1) + m(il, i)
	2884	END DO
	2885	END DO
	2886
	2887	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
	2888	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2889
[2007]	2890	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	2891	!
	2892	! * begin downdraft loop *
	2893	!
	2894	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[524]	2895
[1992]	2896	DO i = nl + 1, 1, -1
[524]	2897
[1992]	2898	num1 = 0
	2899	DO il = 1, ncum
	2900	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
	2901	END DO
	2902	IF (num1<=0) GO TO 400
[1146]	2903
[1992]	2904	CALL zilch(wdtrain, ncum)
[1146]	2905
	2906
[2007]	2907	! * integrate liquid water equation to find condensed water *
	2908	! * and condensed water flux *
	2909	!
	2910	!
	2911	! * calculate detrained precipitation *
[524]	2912
	2913
[3605]	2914	DO il = 1, ncum
	2915	IF (i<=inb(il) .AND. lwork(il)) THEN
	2916	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
	2917	wdtrainS(il, i) = wdtrain(il)/grav ! Ps jyg
	2918	!! wdtrainA(il, i) = wdtrain(il)/grav ! Ps RomP
	2919	END IF
	2920	END DO
	2921
[1992]	2922	IF (i>1) THEN
	2923	DO j = 1, i - 1
	2924	DO il = 1, ncum
	2925	IF (i<=inb(il) .AND. lwork(il)) THEN
	2926	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
	2927	awat = max(awat, 0.0)
[3605]	2928	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
	2929	wdtrainM(il, i) = wdtrain(il)/grav - wdtrainS(il, i) ! Pm jyg
	2930	!! wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP
[1992]	2931	END IF
	2932	END DO
	2933	END DO
	2934	END IF
[524]	2935
[3605]	2936	IF (cvflag_prec_eject) THEN
	2937	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2938	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
	2939	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2940	!!! Warning : this option leads to water conservation violation
	2941	!!! Expert only
	2942	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2943	IF ( i > 1) THEN
	2944	DO il = 1, ncum
	2945	IF (i<=inb(il) .AND. lwork(il)) THEN
	2946	wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))/(1. - qta(il, i-1)) ! Pa jygprl
	2947	wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
	2948	END IF
	2949	END DO
	2950	ENDIF ! ( i > 1)
	2951	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2952	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
	2953	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2954	IF ( i > 1) THEN
	2955	DO il = 1, ncum
	2956	IF (i<=inb(il) .AND. lwork(il)) THEN
	2957	wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i)) ! Pa jygprl
	2958	wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
	2959	END IF
	2960	END DO
	2961	ENDIF ! ( i > 1)
[1146]	2962
[3605]	2963	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
	2964	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2965	ENDIF ! (cvflag_prec_eject)
	2966
	2967
[2007]	2968	! * find rain water and evaporation using provisional *
	2969	! * estimates of rp(i)and rp(i-1) *
[1146]	2970
[1650]	2971
[3605]	2972	IF (cvflag_ice) THEN !!jygprl
	2973	IF (cvflag_prec_eject) THEN
	2974	DO il = 1, ncum !!jygprl
	2975	IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl
	2976	frac(il, i) = (frac_a(il,i)wdtrainA(il,i)+frac_s(il,i)(wdtrainS(il,i)+wdtrainM(il,i))) / & !!jygprl
	2977	max(wdtrainA(il,i)+wdtrainS(il,i)+wdtrainM(il,i),smallestreal) !!jygprl
	2978	fraci(il, i) = frac(il, i) !!jygprl
	2979	END IF !!jygprl
	2980	END DO !!jygprl
	2981	ELSE ! (cvflag_prec_eject)
	2982	DO il = 1, ncum !!jygprl
	2983	IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl
	2984	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2985	IF (keepbug_ice_frac) THEN
	2986	frac(il, i) = frac_s(il, i)
	2987	! Ice fraction computed again here as a function of the temperature seen by unsaturated downdraughts
	2988	! (i.e. the cold pool temperature) for compatibility with earlier versions.
	2989	fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
	2990	fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
	2991	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2992	ELSE ! (keepbug_ice_frac)
	2993	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2994	frac(il, i) = frac_s(il, i)
	2995	fraci(il, i) = frac(il, i) !!jygprl
	2996	ENDIF ! (keepbug_ice_frac)
	2997	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2998	END IF !!jygprl
	2999	END DO !!jygprl
	3000	ENDIF ! (cvflag_prec_eject)
	3001	END IF !!jygprl
	3002
	3003
[1992]	3004	DO il = 1, ncum
	3005	IF (i<=inb(il) .AND. lwork(il)) THEN
[524]	3006
[1992]	3007	wt(il, i) = 45.0
[524]	3008
[1992]	3009	IF (i<inb(il)) THEN
[2007]	3010	rp(il, i) = rp(il, i+1) + &
	3011	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
[1992]	3012	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
	3013	END IF
	3014	rp(il, i) = max(rp(il,i), 0.0)
	3015	rp(il, i) = amin1(rp(il,i), rs(il,i))
	3016	rp(il, inb(il)) = rr(il, inb(il))
[524]	3017
[1992]	3018	IF (i==1) THEN
	3019	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
	3020	IF (cvflag_ice) THEN
[2007]	3021	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
[1992]	3022	END IF
	3023	ELSE
[2007]	3024	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]	3025	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
	3026	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
	3027	rp(il, i-1) = max(rp(il,i-1), 0.0)
[2007]	3028	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
	3029	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]	3030	afac = 0.5*(afac1+afac2)
	3031	END IF
	3032	IF (i==inb(il)) afac = 0.0
	3033	afac = max(afac, 0.0)
	3034	bfac = 1./(sigd(il)*wt(il,i))
[524]	3035
[2393]	3036	!
	3037	IF (prt_level >= 20) THEN
	3038	Print*, 'cv3_unsat after provisional rp estimate: rp, afac, bfac ', &
	3039	i, rp(1, i), afac,bfac
	3040	ENDIF
	3041	!
[2007]	3042	!JYG1
	3043	! cc sigt=1.0
	3044	! cc if(i.ge.icb)sigt=sigp(i)
	3045	! prise en compte de la variation progressive de sigt dans
	3046	! les couches icb et icb-1:
	3047	! pour plcl<ph(i+1), pr1=0 & pr2=1
	3048	! pour plcl>ph(i), pr1=1 & pr2=0
	3049	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
	3050	! sur le nuage, et pr2 est la proportion sous la base du
	3051	! nuage.
[1992]	3052	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
	3053	pr1 = max(0., min(1.,pr1))
	3054	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
	3055	pr2 = max(0., min(1.,pr2))
	3056	sigt = sigp(il, i)*pr1 + pr2
[2007]	3057	!JYG2
[524]	3058
[2007]	3059	!JYG----
	3060	! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
	3061	! c6 = water(il,i+1) + wdtrain(il)*bfac
	3062	! c6 = prec(il,i+1) + wdtrain(il)*bfac
	3063	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
	3064	! evap(il,i)=sigtafacrevap
	3065	! water(il,i)=revap*revap
	3066	! prec(il,i)=revap*revap
	3067	!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
	3068	!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
	3069	!!---end jyg---
[879]	3070
[2007]	3071	! --------retour à la formulation originale d'Emanuel.
[1992]	3072	IF (cvflag_ice) THEN
[524]	3073
[2007]	3074	! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
	3075	! c6=prec(il,i+1)+bfac*wdtrain(il) &
	3076	! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
	3077	! if(c6.gt.0.0)then
	3078	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
[524]	3079
[2007]	3080	!JAM Attention: evap=sigt*E
	3081	! Modification: evap devient l'évaporation en milieu de couche
	3082	! car nécessaire dans cv3_yield
	3083	! Du coup, il faut modifier pas mal d'équations...
	3084	! et l'expression de afac qui devient afac1
	3085	! revap=sqrt((prec(i+1)+prec(i))/2)
[524]	3086
[1992]	3087	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
	3088	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
[2007]	3089	! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
	3090	! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
	3091	! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
[1992]	3092	IF (c6>b6*b6+1.E-20) THEN
	3093	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
	3094	ELSE
	3095	revap = (-b6+sqrt(b6b6+4.c6))/2.
	3096	END IF
	3097	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
[2007]	3098	! print*,prec(il,i),'neige'
[524]	3099
[2007]	3100	!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
	3101	! c evap(il,i)=sigtafacrevap
	3102	! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
	3103	! Ici,l'evaporation evap est simplement calculee par l'equation de
	3104	! conservation.
	3105	! prec(il,i)=revap*revap
	3106	! else
	3107	!JYG---- Correction : si c6 <= 0, water(il,i)=0.
	3108	! prec(il,i)=0.
	3109	! endif
[524]	3110
[2007]	3111	!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
	3112	! moins [tt ce qui sort de la couche i]
	3113	! print *, 'evap avec ice'
	3114	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
	3115	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
[2393]	3116	!
	3117	IF (prt_level >= 20) THEN
	3118	Print*, 'cv3_unsat after evap computation: wdtrain, sigd, wt, prec(i+1),prec(i) ', &
	3119	i, wdtrain(1), sigd(1), wt(1,i), prec(1,i+1),prec(1,i)
	3120	ENDIF
	3121	!
[524]	3122
[2490]	3123	!jyg<
	3124	d6 = prec(il,i)-prec(il,i+1)
[524]	3125
[2490]	3126	!! d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
	3127	!! e6 = bfac*wdtrain(il)
	3128	!! f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
	3129	!>jyg
[2287]	3130	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
	3131	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
[1992]	3132	thaw = min(max(thaw,0.0), 1.0)
[2490]	3133	!jyg<
[1992]	3134	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
[2490]	3135	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
	3136	water(il, i) = min(prec(il,i), max(water(il,i), 0.))
	3137	ice(il, i) = min(prec(il,i), max(ice(il,i), 0.))
	3138
	3139	!! water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
	3140	!! water(il, i) = max(water(il,i), 0.)
	3141	!! ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
	3142	!! ice(il, i) = max(ice(il,i), 0.)
	3143	!>jyg
[1992]	3144	fondue(il, i) = ice(il, i)*thaw
	3145	water(il, i) = water(il, i) + fondue(il, i)
	3146	ice(il, i) = ice(il, i) - fondue(il, i)
[524]	3147
[1992]	3148	IF (water(il,i)+ice(il,i)<1.E-30) THEN
	3149	faci(il, i) = 0.
	3150	ELSE
	3151	faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
	3152	END IF
[524]	3153
[2007]	3154	! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
	3155	! water(il,i)=max(water(il,i),0.)
	3156	! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
	3157	! ice(il,i)=max(ice(il,i),0.)
	3158	! fondue(il,i)=ice(il,i)*thaw
	3159	! water(il,i)=water(il,i)+fondue(il,i)
	3160	! ice(il,i)=ice(il,i)-fondue(il,i)
	3161
	3162	! if((water(il,i)+ice(il,i)).lt.1.e-30)then
	3163	! faci(il,i)=0.
	3164	! else
	3165	! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
	3166	! endif
[524]	3167
[1992]	3168	ELSE
	3169	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
[2007]	3170	c6 = water(il, i+1) + bfac*wdtrain(il) - &
	3171	50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
[1992]	3172	IF (c6>0.0) THEN
	3173	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
	3174	water(il, i) = revap*revap
	3175	ELSE
	3176	water(il, i) = 0.
	3177	END IF
[2007]	3178	! print *, 'evap sans ice'
	3179	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
	3180	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
[524]	3181
[1992]	3182	END IF
	3183	END IF !(i.le.inb(il) .and. lwork(il))
	3184	END DO
[2007]	3185	! ----------------------------------------------------------------
[524]	3186
[2007]	3187	! cc
	3188	! * calculate precipitating downdraft mass flux under *
	3189	! * hydrostatic approximation *
[524]	3190
[1992]	3191	DO il = 1, ncum
	3192	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
[524]	3193
[1992]	3194	tevap(il) = max(0.0, evap(il,i))
	3195	delth = max(0.001, (th(il,i)-th(il,i-1)))
	3196	IF (cvflag_ice) THEN
	3197	IF (cvflag_grav) THEN
[2007]	3198	mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
	3199	(p(il,i-1)-p(il,i))/delth + &
	3200	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
	3201	(p(il,i-1)-p(il,i))/delth + &
	3202	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
	3203	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
[1992]	3204	ELSE
[2007]	3205	mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
	3206	(p(il,i-1)-p(il,i))/delth + &
	3207	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
	3208	(p(il,i-1)-p(il,i))/delth + &
	3209	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
	3210	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
[524]	3211
[1992]	3212	END IF
	3213	ELSE
	3214	IF (cvflag_grav) THEN
	3215	mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
[2007]	3216	(p(il,i-1)-p(il,i))/delth
[1992]	3217	ELSE
	3218	mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
[2007]	3219	(p(il,i-1)-p(il,i))/delth
[1992]	3220	END IF
[524]	3221
[1992]	3222	END IF
[879]	3223
[1992]	3224	END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
[2671]	3225	IF (prt_level .GE. 20) THEN
	3226	PRINT *,'cv3_unsat, mp hydrostatic ', i, mp(il,i)
	3227	ENDIF
[1992]	3228	END DO
[2007]	3229	! ----------------------------------------------------------------
[524]	3230
[2007]	3231	! * if hydrostatic assumption fails, *
	3232	! * solve cubic difference equation for downdraft theta *
	3233	! * and mass flux from two simultaneous differential eqns *
[524]	3234
[1992]	3235	DO il = 1, ncum
	3236	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
[1742]	3237
[1992]	3238	amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
[2007]	3239	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
[1992]	3240	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
[1742]	3241
[1992]	3242	IF (amp2>(0.1*amfac)) THEN
	3243	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
[2007]	3244	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
	3245	(lvcp(il,i)sigd(il)th(il,i))
[1992]	3246	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
[1742]	3247
[1992]	3248	IF (cvflag_ice) THEN
	3249	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
[2007]	3250	50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
	3251	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
[1992]	3252	ELSE
[1774]	3253
[1992]	3254	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
[2007]	3255	50.(p(il,i-1)-p(il,i))xf*tevap(il)
[1992]	3256	END IF
[1742]	3257
[1992]	3258	fac2 = 1.0
	3259	IF (bf<0.0) fac2 = -1.0
	3260	bf = abs(bf)
	3261	ur = 0.25bfbf - afafaftinvtinv*tinv
	3262	IF (ur>=0.0) THEN
	3263	sru = sqrt(ur)
	3264	fac = 1.0
	3265	IF ((0.5*bf-sru)<0.0) fac = -1.0
	3266	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
[2007]	3267	fac(abs(0.5bf-sru))**tinv
[1992]	3268	ELSE
	3269	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
	3270	IF (fac2<0.0) d = 3.14159 - d
	3271	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
	3272	END IF
	3273	mp(il, i) = max(0.0, mp(il,i))
[2671]	3274	IF (prt_level .GE. 20) THEN
	3275	PRINT *,'cv3_unsat, mp cubic ', i, mp(il,i)
	3276	ENDIF
[524]	3277
[1992]	3278	IF (cvflag_ice) THEN
	3279	IF (cvflag_grav) THEN
[2007]	3280	!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
	3281	! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
	3282	! Et il faut bien revoir les facteurs 100.
	3283	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i)) &
	3284	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
	3285	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
	3286	(ph(il,i)-ph(il,i+1))) / &
	3287	(mp(il,i)+sigd(il)*0.1) - &
	3288	10.0(th(il,i)-th(il,i-1))t(il, i) / &
	3289	(lvcp(il,i)sigd(il)th(il,i))
[1992]	3290	ELSE
[2007]	3291	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
	3292	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
	3293	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
	3294	(ph(il,i)-ph(il,i+1))) / &
	3295	(mp(il,i)+sigd(il)*0.1) - &
	3296	10.0(th(il,i)-th(il,i-1))t(il, i) / &
	3297	(lvcp(il,i)sigd(il)th(il,i))
[1992]	3298	END IF
	3299	ELSE
	3300	IF (cvflag_grav) THEN
[2007]	3301	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
	3302	(mp(il,i)+sigd(il)*0.1) - &
	3303	10.0(th(il,i)-th(il,i-1))t(il, i) / &
	3304	(lvcp(il,i)sigd(il)th(il,i))
[1992]	3305	ELSE
[2007]	3306	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
	3307	(mp(il,i)+sigd(il)*0.1) - &
	3308	10.0(th(il,i)-th(il,i-1))t(il, i) / &
	3309	(lvcp(il,i)sigd(il)th(il,i))
[1992]	3310	END IF
	3311	END IF
	3312	b(il, i-1) = max(b(il,i-1), 0.0)
[524]	3313
[1992]	3314	END IF !(amp2.gt.(0.1*amfac))
[524]	3315
[2759]	3316	!jyg< This part shifted 10 lines farther
	3317	!!! * limit magnitude of mp(i) to meet cfl condition *
	3318	!!
	3319	!! ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
	3320	!! amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
	3321	!! ampmax = min(ampmax, amp2)
	3322	!! mp(il, i) = min(mp(il,i), ampmax)
	3323	!>jyg
[524]	3324
[2007]	3325	! * force mp to decrease linearly to zero *
	3326	! * between cloud base and the surface *
[524]	3327
	3328
[2007]	3329	! c if(p(il,i).gt.p(il,icb(il)))then
	3330	! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
	3331	! c endif
[1992]	3332	IF (ph(il,i)>0.9*plcl(il)) THEN
	3333	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
	3334	END IF
[524]	3335
[2759]	3336	!jyg< Shifted part
	3337	! * limit magnitude of mp(i) to meet cfl condition *
	3338
	3339	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
	3340	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
	3341	ampmax = min(ampmax, amp2)
	3342	mp(il, i) = min(mp(il,i), ampmax)
	3343	!>jyg
	3344
[1992]	3345	END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
	3346	END DO
[2007]	3347	! ----------------------------------------------------------------
[2393]	3348	!
	3349	IF (prt_level >= 20) THEN
	3350	Print*, 'cv3_unsat after mp computation: mp, b(i), b(i-1) ', &
	3351	i, mp(1, i), b(1,i), b(1,max(i-1,1))
	3352	ENDIF
	3353	!
[524]	3354
[2007]	3355	! * find mixing ratio of precipitating downdraft *
[524]	3356
[1992]	3357	DO il = 1, ncum
	3358	IF (i<inb(il) .AND. lwork(il)) THEN
	3359	mplus(il) = mp(il, i) > mp(il, i+1)
	3360	END IF ! (i.lt.inb(il) .and. lwork(il))
	3361	END DO
	3362
	3363	DO il = 1, ncum
	3364	IF (i<inb(il) .AND. lwork(il)) THEN
	3365
	3366	rp(il, i) = rr(il, i)
	3367
	3368	IF (mplus(il)) THEN
	3369
	3370	IF (cvflag_grav) THEN
[2007]	3371	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
	3372	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
[1992]	3373	ELSE
[2007]	3374	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
	3375	5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
[1992]	3376	END IF
	3377	rp(il, i) = rp(il, i)/mp(il, i)
[2007]	3378	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
[1992]	3379	up(il, i) = up(il, i)/mp(il, i)
[2007]	3380	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
[1992]	3381	vp(il, i) = vp(il, i)/mp(il, i)
	3382
	3383	ELSE ! if (mplus(il))
	3384
	3385	IF (mp(il,i+1)>1.0E-16) THEN
	3386	IF (cvflag_grav) THEN
[2007]	3387	rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
	3388	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
[1992]	3389	ELSE
[2007]	3390	rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
	3391	(evap(il,i+1)+evap(il,i))/mp(il, i+1)
[1992]	3392	END IF
	3393	up(il, i) = up(il, i+1)
	3394	vp(il, i) = vp(il, i+1)
	3395	END IF ! (mp(il,i+1).gt.1.0e-16)
	3396	END IF ! (mplus(il)) else if (.not.mplus(il))
	3397
	3398	rp(il, i) = amin1(rp(il,i), rs(il,i))
	3399	rp(il, i) = max(rp(il,i), 0.0)
	3400
	3401	END IF ! (i.lt.inb(il) .and. lwork(il))
	3402	END DO
[2007]	3403	! ----------------------------------------------------------------
[1992]	3404
[2007]	3405	! * find tracer concentrations in precipitating downdraft *
[1992]	3406
[2007]	3407	!AC! do j=1,ntra
	3408	!AC! do il = 1,ncum
	3409	!AC! if (i.lt.inb(il) .and. lwork(il)) then
	3410	!AC!c
	3411	!AC! if(mplus(il))then
	3412	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
	3413	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
	3414	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
	3415	!AC! else ! if (mplus(il))
	3416	!AC! if(mp(il,i+1).gt.1.0e-16)then
	3417	!AC! trap(il,i,j)=trap(il,i+1,j)
	3418	!AC! endif
	3419	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
	3420	!AC!c
	3421	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
	3422	!AC! enddo
	3423	!AC! end do
[1992]	3424
	3425	400 END DO
[2007]	3426	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1992]	3427
[2007]	3428	! * end of downdraft loop *
[1992]	3429
[2007]	3430	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1992]	3431
	3432
	3433	RETURN
[3605]	3434
[1992]	3435	END SUBROUTINE cv3_unsat
	3436
[2007]	3437	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
	3438	icb, inb, delt, &
	3439	t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
	3440	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
[3605]	3441	ep, clw, qpreca, m, tp, mp, rp, up, vp, trap, &
[2007]	3442	wt, water, ice, evap, fondue, faci, b, sigd, &
	3443	ment, qent, hent, iflag_mix, uent, vent, &
	3444	nent, elij, traent, sig, &
	3445	tv, tvp, wghti, &
[2306]	3446	iflag, precip, Vprecip, Vprecipi, & ! jyg: Vprecipi
	3447	ft, fr, fu, fv, ftra, & ! jyg
[2007]	3448	cbmf, upwd, dnwd, dnwd0, ma, mip, &
[2259]	3449	!! tls, tps, ! useless . jyg
	3450	qcondc, wd, &
[3605]	3451	ftd, fqd, qta, qtc, sigt, tau_cld_cv, coefw_cld_cv)
[1992]	3452
[2901]	3453	USE print_control_mod, ONLY: lunout, prt_level
[2908]	3454	USE add_phys_tend_mod, only : fl_cor_ebil
[2901]	3455
[1992]	3456	IMPLICIT NONE
	3457
	3458	include "cvthermo.h"
	3459	include "cv3param.h"
	3460	include "cvflag.h"
	3461	include "conema3.h"
	3462
[2007]	3463	!inputs:
[2327]	3464	INTEGER, INTENT (IN) :: iflag_mix
	3465	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
	3466	LOGICAL, INTENT (IN) :: ok_conserv_q
	3467	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
	3468	REAL, INTENT (IN) :: delt
	3469	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, u, v
	3470	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t_wake, rr_wake
	3471	REAL, DIMENSION (nloc), INTENT (IN) :: s_wake
	3472	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra
	3473	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
	3474	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	3475	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz, h, hp
	3476	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tp
	3477	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, cpn, ep, clw
	3478	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
	3479	REAL, DIMENSION (nloc, na), INTENT (IN) :: rp, up
	3480	REAL, DIMENSION (nloc, na), INTENT (IN) :: vp
	3481	REAL, DIMENSION (nloc, nd), INTENT (IN) :: wt
	3482	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: trap
	3483	REAL, DIMENSION (nloc, na), INTENT (IN) :: water, evap, b
	3484	REAL, DIMENSION (nloc, na), INTENT (IN) :: fondue, faci, ice
	3485	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: qent, uent
	3486	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: hent
	3487	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: vent, elij
	3488	INTEGER, DIMENSION (nloc, nd), INTENT (IN) :: nent
	3489	REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent
	3490	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti
[3605]	3491	REAL, DIMENSION (nloc, nd), INTENT (IN) :: qta
	3492	REAL, DIMENSION (nloc, na),INTENT(IN) :: qpreca
	3493	REAL, INTENT(IN) :: tau_cld_cv, coefw_cld_cv
[2007]	3494	!
	3495	!input/output:
[2327]	3496	REAL, DIMENSION (nloc, na), INTENT (INOUT) :: m, mp
	3497	REAL, DIMENSION (nloc, na, na), INTENT (INOUT) :: ment
	3498	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
	3499	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig
	3500	REAL, DIMENSION (nloc), INTENT (INOUT) :: sigd
[2007]	3501	!
	3502	!outputs:
[2327]	3503	REAL, DIMENSION (nloc), INTENT (OUT) :: precip
	3504	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ft, fr, fu, fv
	3505	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ftd, fqd
	3506	REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT) :: ftra
	3507	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: upwd, dnwd, ma
	3508	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: dnwd0, mip
	3509	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecip
	3510	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecipi
	3511	!! REAL tls(nloc, nd), tps(nloc, nd) ! useless . jyg
	3512	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qcondc ! cld
	3513	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qtc, sigt ! cld
	3514	REAL, DIMENSION (nloc), INTENT (OUT) :: wd ! gust
	3515	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf
[2007]	3516	!
	3517	!local variables:
[2508]	3518	INTEGER :: i, k, il, n, j, num1
	3519	REAL :: rat, delti
	3520	REAL :: ax, bx, cx, dx, ex
	3521	REAL :: cpinv, rdcp, dpinv
[3605]	3522	REAL :: sigaq
[2508]	3523	REAL, DIMENSION (nloc) :: awat
	3524	REAL, DIMENSION (nloc, nd) :: lvcp, lfcp ! , mke ! unused . jyg
	3525	REAL, DIMENSION (nloc) :: am, work, ad, amp1
[2007]	3526	!! real up1(nloc), dn1(nloc)
[2508]	3527	REAL, DIMENSION (nloc, nd, nd) :: up1, dn1
	3528	!jyg<
	3529	REAL, DIMENSION (nloc, nd) :: up_to, up_from
	3530	REAL, DIMENSION (nloc, nd) :: dn_to, dn_from
	3531	!>jyg
	3532	REAL, DIMENSION (nloc) :: asum, bsum, csum, dsum
	3533	REAL, DIMENSION (nloc) :: esum, fsum, gsum, hsum
	3534	REAL, DIMENSION (nloc, nd) :: th_wake
	3535	REAL, DIMENSION (nloc) :: alpha_qpos, alpha_qpos1
	3536	REAL, DIMENSION (nloc, nd) :: qcond, nqcond, wa ! cld
	3537	REAL, DIMENSION (nloc, nd) :: siga, sax, mac ! cld
	3538	REAL, DIMENSION (nloc) :: sument
	3539	REAL, DIMENSION (nloc, nd) :: sigment, qtment ! cld
[2007]	3540	REAL sumdq !jyg
	3541	!
	3542	! -------------------------------------------------------------
[1992]	3543
[2007]	3544	! initialization:
[1992]	3545
	3546	delti = 1.0/delt
[2007]	3547	! print*,'cv3_yield initialisation delt', delt
[2393]	3548
[1992]	3549	DO il = 1, ncum
	3550	precip(il) = 0.0
	3551	wd(il) = 0.0 ! gust
	3552	END DO
	3553
[2393]	3554	! Fluxes are on a staggered grid : loops extend up to nl+1
	3555	DO i = 1, nlp
[1992]	3556	DO il = 1, ncum
[2007]	3557	Vprecip(il, i) = 0.0
[2306]	3558	Vprecipi(il, i) = 0.0 ! jyg
[2393]	3559	upwd(il, i) = 0.0
	3560	dnwd(il, i) = 0.0
	3561	dnwd0(il, i) = 0.0
	3562	mip(il, i) = 0.0
	3563	END DO
	3564	END DO
	3565	DO i = 1, nl
	3566	DO il = 1, ncum
[1992]	3567	ft(il, i) = 0.0
	3568	fr(il, i) = 0.0
	3569	fu(il, i) = 0.0
	3570	fv(il, i) = 0.0
	3571	ftd(il, i) = 0.0
	3572	fqd(il, i) = 0.0
	3573	qcondc(il, i) = 0.0 ! cld
	3574	qcond(il, i) = 0.0 ! cld
[2205]	3575	qtc(il, i) = 0.0 ! cld
	3576	qtment(il, i) = 0.0 ! cld
	3577	sigment(il, i) = 0.0 ! cld
	3578	sigt(il, i) = 0.0 ! cld
[1992]	3579	nqcond(il, i) = 0.0 ! cld
	3580	END DO
	3581	END DO
[2007]	3582	! print*,'cv3_yield initialisation 2'
	3583	!AC! do j=1,ntra
	3584	!AC! do i=1,nd
	3585	!AC! do il=1,ncum
	3586	!AC! ftra(il,i,j)=0.0
	3587	!AC! enddo
	3588	!AC! enddo
	3589	!AC! enddo
	3590	! print*,'cv3_yield initialisation 3'
[1992]	3591	DO i = 1, nl
	3592	DO il = 1, ncum
	3593	lvcp(il, i) = lv(il, i)/cpn(il, i)
	3594	lfcp(il, i) = lf(il, i)/cpn(il, i)
	3595	END DO
	3596	END DO
	3597
	3598
	3599
[2007]	3600	! * calculate surface precipitation in mm/day *
[1992]	3601
	3602	DO il = 1, ncum
	3603	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
	3604	IF (cvflag_ice) THEN
[2007]	3605	precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
	3606	86400.1000./(rowl*grav)
[1992]	3607	ELSE
[2007]	3608	precip(il) = wt(il, 1)sigd(il)water(il, 1) &
	3609	86400.1000./(rowl*grav)
[1992]	3610	END IF
	3611	END IF
	3612	END DO
[2007]	3613	! print*,'cv3_yield apres calcul precip'
[1992]	3614
	3615
[2007]	3616	! === calculate vertical profile of precipitation in kg/m2/s ===
[1992]	3617
	3618	DO i = 1, nl
	3619	DO il = 1, ncum
	3620	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
	3621	IF (cvflag_ice) THEN
[2007]	3622	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
[2306]	3623	Vprecipi(il, i) = wt(il, i)sigd(il)ice(il,i)/grav ! jyg
[1992]	3624	ELSE
[2007]	3625	Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
[2306]	3626	Vprecipi(il, i) = 0. ! jyg
[1992]	3627	END IF
	3628	END IF
	3629	END DO
	3630	END DO
	3631
	3632
[2007]	3633	! * Calculate downdraft velocity scale *
	3634	! * NE PAS UTILISER POUR L'INSTANT *
[1992]	3635
[2007]	3636	!! do il=1,ncum
	3637	!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
	3638	!! /(sigd(il)*p(il,icb(il)))
	3639	!! enddo
[1992]	3640
	3641
[2007]	3642	! * calculate tendencies of lowest level potential temperature *
	3643	! * and mixing ratio *
[1992]	3644
	3645	DO il = 1, ncum
	3646	work(il) = 1.0/(ph(il,1)-ph(il,2))
	3647	cbmf(il) = 0.0
	3648	END DO
	3649
[3605]	3650	! - Adiabatic ascent mass flux "ma" and cloud base mass flux "cbmf"
	3651	!-----------------------------------------------------------------
	3652	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3653	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
	3654	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3655	!!! Warning : this option leads to water conservation violation
	3656	!!! Expert only
	3657	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3658	DO il = 1, ncum
	3659	ma(il, nlp) = 0.
	3660	ma(il, 1) = 0.
	3661	END DO
	3662	DO k = nl, 2, -1
	3663	DO il = 1, ncum
	3664	ma(il, k) = ma(il, k+1)*(1.-qta(il, k))/(1.-qta(il, k-1)) + m(il, k)
	3665	cbmf(il) = max(cbmf(il), ma(il,k))
	3666	END DO
	3667	END DO
	3668	DO k = 2,nl
	3669	DO il = 1, ncum
	3670	IF (k <icb(il)) THEN
	3671	ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
	3672	ENDIF
	3673	END DO
	3674	END DO
	3675	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3676	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
	3677	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3678	!! Line kept for compatibility with earlier versions
[1992]	3679	DO k = 2, nl
	3680	DO il = 1, ncum
	3681	IF (k>=icb(il)) THEN
	3682	cbmf(il) = cbmf(il) + m(il, k)
	3683	END IF
	3684	END DO
	3685	END DO
	3686
[3605]	3687	DO il = 1, ncum
	3688	ma(il, nlp) = 0.
	3689	ma(il, 1) = 0.
	3690	END DO
	3691	DO k = nl, 2, -1
	3692	DO il = 1, ncum
	3693	ma(il, k) = ma(il, k+1) + m(il, k)
	3694	END DO
	3695	END DO
	3696	DO k = 2,nl
	3697	DO il = 1, ncum
	3698	IF (k <icb(il)) THEN
	3699	ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
	3700	ENDIF
	3701	END DO
	3702	END DO
	3703
	3704	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
	3705	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3706	!
[2007]	3707	! print*,'cv3_yield avant ft'
	3708	! am is the part of cbmf taken from the first level
[1992]	3709	DO il = 1, ncum
	3710	am(il) = cbmf(il)*wghti(il, 1)
	3711	END DO
	3712
	3713	DO il = 1, ncum
	3714	IF (iflag(il)<=1) THEN
[2007]	3715	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
	3716	!JYG Correction pour conserver l'eau
	3717	! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
[1992]	3718	IF (cvflag_ice) THEN
	3719	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
[2007]	3720	lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
	3721	lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
	3722	(100.*(ph(il,1)-ph(il,2))) !precip
[1992]	3723	ELSE
	3724	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
	3725	END IF
	3726
[2007]	3727	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
[1992]	3728
	3729	IF (cvflag_ice) THEN
[2007]	3730	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
	3731	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
	3732	0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
	3733	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
[1992]	3734	ELSE
[2007]	3735	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
	3736	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
[1992]	3737	END IF
	3738
[2007]	3739	ftd(il, 1) = ft(il, 1) ! fin precip
[1992]	3740
[2007]	3741	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
[2908]	3742	!jyg<
	3743	IF (fl_cor_ebil >= 2) THEN
	3744	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
	3745	((t(il,2)-t(il,1))*cpn(il,2)+gz(il,2)-gz(il,1))/cpn(il,1)
	3746	ELSE
	3747	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
	3748	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
	3749	ENDIF
	3750	!>jyg
[1992]	3751	END IF ! iflag
	3752	END DO
	3753
	3754
	3755	DO j = 2, nl
	3756	IF (iflag_mix>0) THEN
	3757	DO il = 1, ncum
[2007]	3758	! FH WARNING a modifier :
[1992]	3759	cpinv = 0.
[2007]	3760	! cpinv=1.0/cpn(il,1)
[1992]	3761	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
[2007]	3762	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
	3763	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
[1992]	3764	END IF ! j
	3765	END DO
	3766	END IF
	3767	END DO
[2007]	3768	! fin sature
[1992]	3769
	3770
	3771	DO il = 1, ncum
	3772	IF (iflag(il)<=1) THEN
[2007]	3773	!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
	3774	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
	3775	sigd(il)*evap(il, 1)
	3776	!!! sigd(il)0.5(evap(il,1)+evap(il,2))
[1992]	3777
[2007]	3778	fqd(il, 1) = fr(il, 1) !precip
[1992]	3779
[2007]	3780	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
[1992]	3781
[2007]	3782	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
	3783	am(il)*(u(il,2)-u(il,1)))
	3784	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
	3785	am(il)*(v(il,2)-v(il,1)))
[1992]	3786	END IF ! iflag
	3787	END DO ! il
	3788
	3789
[2007]	3790	!AC! do j=1,ntra
	3791	!AC! do il=1,ncum
	3792	!AC! if (iflag(il) .le. 1) then
	3793	!AC! if (cvflag_grav) then
	3794	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
	3795	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	3796	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	3797	!AC! else
	3798	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
	3799	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	3800	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	3801	!AC! endif
	3802	!AC! endif ! iflag
	3803	!AC! enddo
	3804	!AC! enddo
[1992]	3805
	3806	DO j = 2, nl
	3807	DO il = 1, ncum
	3808	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
[2007]	3809	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
	3810	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
	3811	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
[1992]	3812	END IF ! j
	3813	END DO
	3814	END DO
	3815
[2007]	3816	!AC! do k=1,ntra
	3817	!AC! do j=2,nl
	3818	!AC! do il=1,ncum
	3819	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
	3820	!AC!
	3821	!AC! if (cvflag_grav) then
	3822	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
	3823	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
	3824	!AC! else
	3825	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
	3826	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
	3827	!AC! endif
	3828	!AC!
	3829	!AC! endif
	3830	!AC! enddo
	3831	!AC! enddo
	3832	!AC! enddo
	3833	! print*,'cv3_yield apres ft'
[1992]	3834
[2508]	3835	!jyg<
	3836	!-----------------------------------------------------------
	3837	IF (ok_optim_yield) THEN !\|
	3838	!-----------------------------------------------------------
	3839	!
	3840	!* *
	3841	!* Compute convective mass fluxes upwd and dnwd *
	3842
[3605]	3843	!
	3844	! =================================================
	3845	! upward fluxes \|
	3846	! ------------------------------------------------
	3847	!
[2508]	3848	upwd(:,:) = 0.
	3849	up_to(:,:) = 0.
	3850	up_from(:,:) = 0.
	3851	!
[3605]	3852	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3853	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
	3854	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3855	!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
	3856	!! is taken into account.
	3857	!! WARNING : in the present version, taking into account the mass-flux decrease due to
	3858	!! precipitation ejection leads to water conservation violation.
	3859	!
	3860	! - Upward mass flux of mixed draughts
	3861	!---------------------------------------
[2508]	3862	DO i = 2, nl
[3605]	3863	DO j = 1, i-1
	3864	DO il = 1, ncum
	3865	IF (i<=inb(il)) THEN
	3866	up_to(il,i) = up_to(il,i) + ment(il,j,i)
	3867	ENDIF
	3868	ENDDO
	3869	ENDDO
	3870	ENDDO
	3871	!
	3872	DO j = 3, nl
	3873	DO i = 2, j-1
	3874	DO il = 1, ncum
	3875	IF (j<=inb(il)) THEN
	3876	up_from(il,i) = up_from(il,i) + ment(il,i,j)
	3877	ENDIF
	3878	ENDDO
	3879	ENDDO
	3880	ENDDO
	3881	!
	3882	! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
	3883	!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
	3884	!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
	3885	!
	3886	DO i = 2, nlp
[2508]	3887	DO il = 1, ncum
[3605]	3888	IF (i<=inb(il)+1) THEN
	3889	upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
	3890	ENDIF
	3891	ENDDO
	3892	ENDDO
	3893	!
	3894	! - Total upward mass flux
	3895	!---------------------------
	3896	DO i = 2, nlp
	3897	DO il = 1, ncum
	3898	IF (i<=inb(il)+1) THEN
	3899	upwd(il,i) = upwd(il,i) + ma(il,i)
	3900	ENDIF
	3901	ENDDO
	3902	ENDDO
	3903	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3904	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
	3905	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3906	!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
	3907	!! is not taken into account.
	3908	!
	3909	! - Upward mass flux
	3910	!-------------------
	3911	DO i = 2, nl
	3912	DO il = 1, ncum
[2508]	3913	IF (i<=inb(il)) THEN
	3914	up_to(il,i) = m(il,i)
	3915	ENDIF
	3916	ENDDO
	3917	DO j = 1, i-1
	3918	DO il = 1, ncum
	3919	IF (i<=inb(il)) THEN
	3920	up_to(il,i) = up_to(il,i) + ment(il,j,i)
	3921	ENDIF
	3922	ENDDO
	3923	ENDDO
	3924	ENDDO
	3925	!
	3926	DO i = 1, nl
	3927	DO il = 1, ncum
	3928	IF (i<=inb(il)) THEN
	3929	up_from(il,i) = cbmf(il)*wghti(il,i)
	3930	ENDIF
	3931	ENDDO
	3932	ENDDO
[3605]	3933	!
[2508]	3934	DO j = 3, nl
	3935	DO i = 2, j-1
	3936	DO il = 1, ncum
	3937	IF (j<=inb(il)) THEN
	3938	up_from(il,i) = up_from(il,i) + ment(il,i,j)
	3939	ENDIF
	3940	ENDDO
	3941	ENDDO
	3942	ENDDO
	3943	!
	3944	! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
	3945	!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
	3946	!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
	3947	!
	3948	DO i = 2, nlp
	3949	DO il = 1, ncum
[3197]	3950	IF (i<=inb(il)+1) THEN
	3951	upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
	3952	ENDIF
[2508]	3953	ENDDO
	3954	ENDDO
[3605]	3955
	3956
	3957	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
	3958	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3959
[2508]	3960	!
	3961	! =================================================
	3962	! downward fluxes \|
	3963	! ------------------------------------------------
[3605]	3964	dnwd(:,:) = 0.
	3965	dn_to(:,:) = 0.
	3966	dn_from(:,:) = 0.
[2508]	3967	DO i = 1, nl
	3968	DO j = i+1, nl
	3969	DO il = 1, ncum
	3970	IF (j<=inb(il)) THEN
	3971	dn_to(il,i) = dn_to(il,i) + ment(il,j,i)
	3972	ENDIF
	3973	ENDDO
	3974	ENDDO
	3975	ENDDO
	3976	!
	3977	DO j = 1, nl
	3978	DO i = j+1, nl
	3979	DO il = 1, ncum
	3980	IF (i<=inb(il)) THEN
	3981	dn_from(il,i) = dn_from(il,i) + ment(il,i,j)
	3982	ENDIF
	3983	ENDDO
	3984	ENDDO
	3985	ENDDO
	3986	!
	3987	! The difference between dnwd(il,i) and dnwd(il,i+1) is due to downdrafts ending in layer
	3988	!(i) (theses drafts cross interface (i+1) but not interface(i)) and to downdrafts
	3989	!starting from layer (i) (theses drafts cross interface (i) but not interface(i+1)):
	3990	!
	3991	DO i = nl-1, 1, -1
	3992	DO il = 1, ncum
	3993	dnwd(il,i) = max(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i))
	3994	ENDDO
	3995	ENDDO
	3996	! =================================================
	3997	!
	3998	!-----------------------------------------------------------
	3999	ENDIF !(ok_optim_yield) !\|
	4000	!-----------------------------------------------------------
	4001	!>jyg
	4002
[2007]	4003	! * calculate tendencies of potential temperature and mixing ratio *
	4004	! * at levels above the lowest level *
[1992]	4005
[2007]	4006	! * first find the net saturated updraft and downdraft mass fluxes *
	4007	! * through each level *
[1992]	4008
	4009
[2508]	4010	!jyg<
	4011	!! DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
	4012	DO i = 2, nl
	4013	!>jyg
[1992]	4014
	4015	num1 = 0
	4016	DO il = 1, ncum
	4017	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
	4018	END DO
	4019	IF (num1<=0) GO TO 500
	4020
[2508]	4021	!
[2393]	4022	!jyg<
[2508]	4023	!-----------------------------------------------------------
	4024	IF (ok_optim_yield) THEN !\|
	4025	!-----------------------------------------------------------
	4026	DO il = 1, ncum
	4027	amp1(il) = upwd(il,i+1)
	4028	ad(il) = dnwd(il,i)
	4029	ENDDO
	4030	!-----------------------------------------------------------
	4031	ELSE !(ok_optim_yield) !\|
	4032	!-----------------------------------------------------------
	4033	!>jyg
[2393]	4034	DO il = 1,ncum
	4035	amp1(il) = 0.
	4036	ad(il) = 0.
	4037	ENDDO
[1992]	4038
	4039	DO k = 1, nl + 1
	4040	DO il = 1, ncum
	4041	IF (i>=icb(il)) THEN
	4042	IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
	4043	amp1(il) = amp1(il) + m(il, k)
	4044	END IF
	4045	ELSE
[2007]	4046	! AMP1 is the part of cbmf taken from layers I and lower
[1992]	4047	IF (k<=i) THEN
	4048	amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
	4049	END IF
	4050	END IF
	4051	END DO
	4052	END DO
	4053
[2508]	4054	DO j = i + 1, nl + 1
	4055	DO k = 1, i
	4056	!yor! reverted j and k loops
	4057	DO il = 1, ncum
	4058	!yor! IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN ! the second condition implies the first !
	4059	IF (j<=(inb(il)+1)) THEN
	4060	amp1(il) = amp1(il) + ment(il, k, j)
	4061	END IF
	4062	END DO
	4063	END DO
[1992]	4064	END DO
	4065
	4066	DO k = 1, i - 1
[2508]	4067	!jyg<
	4068	!! DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
	4069	DO j = i, nl
	4070	!>jyg
[1992]	4071	DO il = 1, ncum
[2508]	4072	!yor! IF (i<=inb(il) .AND. j<=inb(il)) THEN ! the second condition implies the 1st !
	4073	IF (j<=inb(il)) THEN
[1992]	4074	ad(il) = ad(il) + ment(il, j, k)
	4075	END IF
	4076	END DO
	4077	END DO
	4078	END DO
[2508]	4079	!
	4080	!-----------------------------------------------------------
	4081	ENDIF !(ok_optim_yield) !\|
	4082	!-----------------------------------------------------------
	4083	!
	4084	!! print *,'yield, i, amp1, ad', i, amp1(1), ad(1)
[1992]	4085
	4086	DO il = 1, ncum
	4087	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
	4088	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	4089	cpinv = 1.0/cpn(il, i)
	4090
[2007]	4091	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
	4092	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
[1992]	4093
[2007]	4094	! precip
	4095	! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
[1992]	4096	IF (cvflag_ice) THEN
	4097	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
[2007]	4098	sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
	4099	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
[1992]	4100	ELSE
	4101	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
	4102	END IF
	4103
	4104	rat = cpn(il, i-1)*cpinv
	4105
[2007]	4106	ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
	4107	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
	4108	IF (cvflag_ice) THEN
	4109	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
	4110	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
	4111	0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
	4112	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
	4113	ELSE
	4114	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
	4115	(t_wake(il,i+1)-t_wake(il,i))dpinv &
	4116	cpinv
	4117	END IF
[1992]	4118
[2007]	4119	ftd(il, i) = ft(il, i)
	4120	! fin precip
[1992]	4121
[2007]	4122	! sature
[2908]	4123	!jyg<
	4124	IF (fl_cor_ebil >= 2) THEN
	4125	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
	4126	( amp1(il)( (t(il,i+1)-t(il,i))cpn(il,i+1) + gz(il,i+1)-gz(il,i))*cpinv - &
	4127	ad(il)( (t(il,i)-t(il,i-1))cpn(il,i-1) + gz(il,i)-gz(il,i-1))*cpinv)
	4128	ELSE
	4129	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
[2007]	4130	(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
	4131	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
[2908]	4132	ENDIF
	4133	!>jyg
[1992]	4134
	4135
[2007]	4136	IF (iflag_mix==0) THEN
	4137	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
	4138	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
	4139	END IF
[2902]	4140	!
[2007]	4141	! sb: on ne fait pas encore la correction permettant de mieux
	4142	! conserver l'eau:
	4143	!JYG: correction permettant de mieux conserver l'eau:
	4144	! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
	4145	fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
	4146	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
	4147	fqd(il, i) = fr(il, i) ! precip
[1992]	4148
[2007]	4149	fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
	4150	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	4151	fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
	4152	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
[1992]	4153
	4154
[2007]	4155	fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
	4156	ad(il)*(rr(il,i)-rr(il,i-1)))
	4157	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
	4158	ad(il)*(u(il,i)-u(il,i-1)))
	4159	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
	4160	ad(il)*(v(il,i)-v(il,i-1)))
[1992]	4161
	4162	END IF ! i
	4163	END DO
	4164
[2007]	4165	!AC! do k=1,ntra
	4166	!AC! do il=1,ncum
	4167	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
	4168	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	4169	!AC! cpinv=1.0/cpn(il,i)
	4170	!AC! if (cvflag_grav) then
	4171	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
	4172	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	4173	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	4174	!AC! else
	4175	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
	4176	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	4177	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	4178	!AC! endif
	4179	!AC! endif
	4180	!AC! enddo
	4181	!AC! enddo
[1992]	4182
	4183	DO k = 1, i - 1
	4184
	4185	DO il = 1, ncum
	4186	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
	4187	awat(il) = max(awat(il), 0.0)
	4188	END DO
	4189
	4190	IF (iflag_mix/=0) THEN
	4191	DO il = 1, ncum
	4192	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
	4193	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	4194	cpinv = 1.0/cpn(il, i)
[2007]	4195	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
	4196	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
	4197	!
	4198	!
[1992]	4199	END IF ! i
	4200	END DO
	4201	END IF
	4202
	4203	DO il = 1, ncum
	4204	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
	4205	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	4206	cpinv = 1.0/cpn(il, i)
[2007]	4207	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
	4208	(qent(il,k,i)-awat(il)-rr(il,i))
	4209	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
	4210	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
[1992]	4211
[2007]	4212	! (saturated updrafts resulting from mixing) ! cld
	4213	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
[2205]	4214	qtment(il, i) = qtment(il, i) + qent(il,k,i) ! cld
	4215	nqcond(il, i) = nqcond(il, i) + 1. ! cld
[1992]	4216	END IF ! i
	4217	END DO
	4218	END DO
	4219
[2007]	4220	!AC! do j=1,ntra
	4221	!AC! do k=1,i-1
	4222	!AC! do il=1,ncum
	4223	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
	4224	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	4225	!AC! cpinv=1.0/cpn(il,i)
	4226	!AC! if (cvflag_grav) then
	4227	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	4228	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
	4229	!AC! else
	4230	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	4231	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
	4232	!AC! endif
	4233	!AC! endif
	4234	!AC! enddo
	4235	!AC! enddo
	4236	!AC! enddo
[1992]	4237
[2508]	4238	!jyg<
	4239	!! DO k = i, nl + 1
	4240	DO k = i, nl
	4241	!>jyg
[1992]	4242
	4243	IF (iflag_mix/=0) THEN
	4244	DO il = 1, ncum
	4245	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
	4246	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	4247	cpinv = 1.0/cpn(il, i)
[2007]	4248	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
	4249	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
[1992]	4250
	4251
	4252	END IF ! i
	4253	END DO
	4254	END IF
	4255
	4256	DO il = 1, ncum
	4257	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
	4258	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	4259	cpinv = 1.0/cpn(il, i)
	4260
[2007]	4261	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
	4262	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
	4263	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
[1992]	4264	END IF ! i and k
	4265	END DO
	4266	END DO
	4267
[2007]	4268	!AC! do j=1,ntra
	4269	!AC! do k=i,nl+1
	4270	!AC! do il=1,ncum
	4271	!AC! if (i.le.inb(il) .and. k.le.inb(il)
	4272	!AC! $ .and. iflag(il) .le. 1) then
	4273	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	4274	!AC! cpinv=1.0/cpn(il,i)
	4275	!AC! if (cvflag_grav) then
	4276	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	4277	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
	4278	!AC! else
	4279	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	4280	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
	4281	!AC! endif
	4282	!AC! endif ! i and k
	4283	!AC! enddo
	4284	!AC! enddo
	4285	!AC! enddo
[1992]	4286
[2007]	4287	! sb: interface with the cloud parameterization: ! cld
[1992]	4288
	4289	DO k = i + 1, nl
	4290	DO il = 1, ncum
[2007]	4291	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
	4292	! (saturated downdrafts resulting from mixing) ! cld
	4293	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
[2205]	4294	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
	4295	nqcond(il, i) = nqcond(il, i) + 1. ! cld
[1992]	4296	END IF ! cld
	4297	END DO ! cld
	4298	END DO ! cld
	4299
[3605]	4300	!ym BIG Warning : it seems that the k loop is missing !!!
	4301	!ym Strong advice to check this
	4302	!ym add a k loop temporary
	4303
[2007]	4304	! (particular case: no detraining level is found) ! cld
[3605]	4305	! Verif merge Dynamico<<<<<<< .working
[2007]	4306	DO il = 1, ncum ! cld
	4307	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
	4308	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
[3345]	4309	!jyg< Bug correction 20180620
	4310	! PROBLEM: Should not qent(il,i,i) be taken into account even if nent(il,i)/=0?
	4311	!! qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
	4312	qtment(il, i) = qent(il,i,i) + qtment(il,i) ! cld
	4313	!>jyg
[2007]	4314	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	4315	END IF ! cld
	4316	END DO ! cld
[3605]	4317	! Verif merge Dynamico =======
	4318	! Verif merge Dynamico DO k = i + 1, nl
	4319	! Verif merge Dynamico DO il = 1, ncum !ym k loop added ! cld
	4320	! Verif merge Dynamico IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
	4321	! Verif merge Dynamico qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
	4322	! Verif merge Dynamico qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
	4323	! Verif merge Dynamico nqcond(il, i) = nqcond(il, i) + 1. ! cld
	4324	! Verif merge Dynamico END IF ! cld
	4325	! Verif merge Dynamico END DO
	4326	! Verif merge Dynamico ENDDO ! cld
	4327	! Verif merge Dynamico >>>>>>> .merge-right.r3413
[1992]	4328
[2007]	4329	DO il = 1, ncum ! cld
	4330	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
	4331	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
[2205]	4332	qtment(il, i) = qtment(il,i)/nqcond(il, i) ! cld
[2007]	4333	END IF ! cld
[1992]	4334	END DO
	4335
[2007]	4336	!AC! do j=1,ntra
	4337	!AC! do il=1,ncum
	4338	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
	4339	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	4340	!AC! cpinv=1.0/cpn(il,i)
	4341	!AC!
	4342	!AC! if (cvflag_grav) then
	4343	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
	4344	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	4345	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
	4346	!AC! else
	4347	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
	4348	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	4349	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
	4350	!AC! endif
	4351	!AC! endif ! i
	4352	!AC! enddo
	4353	!AC! enddo
[1992]	4354
	4355
	4356	500 END DO
	4357
[2007]	4358	!JYG<
	4359	!Conservation de l'eau
	4360	! sumdq = 0.
	4361	! DO k = 1, nl
	4362	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
	4363	! END DO
	4364	! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
	4365	!JYG>
	4366	! * move the detrainment at level inb down to level inb-1 *
	4367	! * in such a way as to preserve the vertically *
	4368	! * integrated enthalpy and water tendencies *
[1992]	4369
[2007]	4370	! Correction bug le 18-03-09
[1992]	4371	DO il = 1, ncum
	4372	IF (iflag(il)<=1) THEN
[2007]	4373	ax = 0.01gravment(il, inb(il), inb(il))* &
	4374	(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
	4375	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
	4376	ft(il, inb(il)) = ft(il, inb(il)) - ax
	4377	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
	4378	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
[1992]	4379
[2007]	4380	bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
	4381	(ph(il,inb(il))-ph(il,inb(il)+1))
	4382	fr(il, inb(il)) = fr(il, inb(il)) - bx
	4383	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
	4384	(ph(il,inb(il)-1)-ph(il,inb(il)))
[1992]	4385
[2007]	4386	cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
	4387	(ph(il,inb(il))-ph(il,inb(il)+1))
	4388	fu(il, inb(il)) = fu(il, inb(il)) - cx
	4389	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
	4390	(ph(il,inb(il)-1)-ph(il,inb(il)))
[1992]	4391
[2007]	4392	dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
	4393	(ph(il,inb(il))-ph(il,inb(il)+1))
	4394	fv(il, inb(il)) = fv(il, inb(il)) - dx
	4395	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
	4396	(ph(il,inb(il)-1)-ph(il,inb(il)))
[1992]	4397	END IF !iflag
	4398	END DO
	4399
[2007]	4400	!JYG<
	4401	!Conservation de l'eau
	4402	! sumdq = 0.
	4403	! DO k = 1, nl
	4404	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
	4405	! END DO
	4406	! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
	4407	!JYG>
[1992]	4408
[2007]	4409	!AC! do j=1,ntra
	4410	!AC! do il=1,ncum
	4411	!AC! IF (iflag(il) .le. 1) THEN
	4412	!AC! IF (cvflag_grav) then
	4413	!AC! ex=0.01gravment(il,inb(il),inb(il))
	4414	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
	4415	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
	4416	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
	4417	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
	4418	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
	4419	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
	4420	!AC! else
	4421	!AC! ex=0.1*ment(il,inb(il),inb(il))
	4422	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
	4423	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
	4424	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
	4425	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
	4426	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
	4427	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
	4428	!AC! ENDIF !cvflag grav
	4429	!AC! ENDIF !iflag
	4430	!AC! enddo
	4431	!AC! enddo
[1992]	4432
	4433
[2007]	4434	! * homogenize tendencies below cloud base *
[1992]	4435
[2007]	4436
[1992]	4437	DO il = 1, ncum
	4438	asum(il) = 0.0
	4439	bsum(il) = 0.0
	4440	csum(il) = 0.0
	4441	dsum(il) = 0.0
	4442	esum(il) = 0.0
	4443	fsum(il) = 0.0
	4444	gsum(il) = 0.0
	4445	hsum(il) = 0.0
	4446	END DO
	4447
[2007]	4448	!do i=1,nl
	4449	!do il=1,ncum
	4450	!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
	4451	!enddo
	4452	!enddo
[1992]	4453
	4454	DO i = 1, nl
	4455	DO il = 1, ncum
	4456	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
[2007]	4457	!jyg Saturated part : use T profile
[1992]	4458	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
[2007]	4459	!jyg<20140311
	4460	!Correction pour conserver l eau
	4461	IF (ok_conserv_q) THEN
	4462	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
	4463	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
	4464
	4465	ELSE
	4466	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
	4467	(ph(il,i)-ph(il,i+1))
	4468	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
	4469	(ph(il,i)-ph(il,i+1))
	4470	ENDIF ! (ok_conserv_q)
	4471	!jyg>
[1992]	4472	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
[2007]	4473	!jyg Unsaturated part : use T_wake profile
[1992]	4474	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
[2007]	4475	!jyg<20140311
	4476	!Correction pour conserver l eau
	4477	IF (ok_conserv_q) THEN
	4478	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
	4479	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
	4480	ELSE
	4481	fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
	4482	(ph(il,i)-ph(il,i+1))
	4483	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
	4484	(ph(il,i)-ph(il,i+1))
	4485	ENDIF ! (ok_conserv_q)
	4486	!jyg>
	4487	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
[1992]	4488	END IF
	4489	END DO
	4490	END DO
	4491
[2007]	4492	!!!! do 700 i=1,icb(il)-1
[2901]	4493	IF (ok_homo_tend) THEN
	4494	DO i = 1, nl
	4495	DO il = 1, ncum
	4496	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
	4497	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
	4498	fqd(il, i) = fsum(il)/gsum(il)
	4499	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
	4500	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
	4501	END IF
	4502	END DO
[1992]	4503	END DO
[2901]	4504	ENDIF
[1992]	4505
[2007]	4506	!jyg<
	4507	!Conservation de l'eau
	4508	!! sumdq = 0.
	4509	!! DO k = 1, nl
	4510	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
	4511	!! END DO
	4512	!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
	4513	!jyg>
[1992]	4514
[2007]	4515
	4516	! *** Check that moisture stays positive. If not, scale tendencies
	4517	! in order to ensure moisture positivity
[1992]	4518	DO il = 1, ncum
	4519	alpha_qpos(il) = 1.
	4520	IF (iflag(il)<=1) THEN
	4521	IF (fr(il,1)<=0.) THEN
[2007]	4522	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]	4523	END IF
	4524	END IF
	4525	END DO
	4526	DO i = 2, nl
	4527	DO il = 1, ncum
	4528	IF (iflag(il)<=1) THEN
	4529	IF (fr(il,i)<=0.) THEN
[2007]	4530	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
	4531	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
[1992]	4532	END IF
	4533	END IF
	4534	END DO
	4535	END DO
	4536	DO il = 1, ncum
	4537	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
	4538	alpha_qpos(il) = alpha_qpos(il)*1.1
	4539	END IF
	4540	END DO
[2327]	4541	!
[2901]	4542	IF (prt_level .GE. 5) THEN
	4543	print *,' CV3_YIELD : alpha_qpos ',alpha_qpos(1)
	4544	ENDIF
[2327]	4545	!
[1992]	4546	DO il = 1, ncum
	4547	IF (iflag(il)<=1) THEN
	4548	sigd(il) = sigd(il)/alpha_qpos(il)
	4549	precip(il) = precip(il)/alpha_qpos(il)
[2579]	4550	cbmf(il) = cbmf(il)/alpha_qpos(il)
[1992]	4551	END IF
	4552	END DO
	4553	DO i = 1, nl
	4554	DO il = 1, ncum
	4555	IF (iflag(il)<=1) THEN
	4556	fr(il, i) = fr(il, i)/alpha_qpos(il)
	4557	ft(il, i) = ft(il, i)/alpha_qpos(il)
	4558	fqd(il, i) = fqd(il, i)/alpha_qpos(il)
	4559	ftd(il, i) = ftd(il, i)/alpha_qpos(il)
	4560	fu(il, i) = fu(il, i)/alpha_qpos(il)
	4561	fv(il, i) = fv(il, i)/alpha_qpos(il)
	4562	m(il, i) = m(il, i)/alpha_qpos(il)
	4563	mp(il, i) = mp(il, i)/alpha_qpos(il)
[2306]	4564	Vprecip(il, i) = Vprecip(il, i)/alpha_qpos(il)
	4565	Vprecipi(il, i) = Vprecipi(il, i)/alpha_qpos(il) ! jyg
[1992]	4566	END IF
	4567	END DO
	4568	END DO
[2508]	4569	!jyg<
	4570	!-----------------------------------------------------------
	4571	IF (ok_optim_yield) THEN !\|
	4572	!-----------------------------------------------------------
[2459]	4573	DO i = 1, nl
[2508]	4574	DO il = 1, ncum
	4575	IF (iflag(il)<=1) THEN
	4576	upwd(il, i) = upwd(il, i)/alpha_qpos(il)
	4577	dnwd(il, i) = dnwd(il, i)/alpha_qpos(il)
	4578	END IF
	4579	END DO
	4580	END DO
	4581	!-----------------------------------------------------------
	4582	ENDIF !(ok_optim_yield) !\|
	4583	!-----------------------------------------------------------
	4584	!>jyg
	4585	DO j = 1, nl !yor! inverted i and j loops
	4586	DO i = 1, nl
[1992]	4587	DO il = 1, ncum
	4588	IF (iflag(il)<=1) THEN
	4589	ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
	4590	END IF
	4591	END DO
	4592	END DO
	4593	END DO
	4594
[2007]	4595	!AC! DO j = 1,ntra
	4596	!AC! DO i = 1,nl
	4597	!AC! DO il = 1,ncum
	4598	!AC! IF (iflag(il) .le. 1) THEN
	4599	!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
	4600	!AC! ENDIF
	4601	!AC! ENDDO
	4602	!AC! ENDDO
	4603	!AC! ENDDO
[1992]	4604
	4605
[2007]	4606	! * reset counter and return *
[1992]	4607
[2253]	4608	! Reset counter only for points actually convective (jyg)
	4609	! In order take into account the possibility of changing the compression,
	4610	! reset m, sig and w0 to zero for non-convecting points.
[1992]	4611	DO il = 1, ncum
[2253]	4612	IF (iflag(il) < 3) THEN
	4613	sig(il, nd) = 2.0
	4614	ENDIF
[1992]	4615	END DO
	4616
	4617
[2393]	4618	DO i = 1, nl
[1992]	4619	DO il = 1, ncum
	4620	dnwd0(il, i) = -mp(il, i)
	4621	END DO
	4622	END DO
[2393]	4623	!jyg< (loops stop at nl)
	4624	!! DO i = nl + 1, nd
	4625	!! DO il = 1, ncum
	4626	!! dnwd0(il, i) = 0.
	4627	!! END DO
	4628	!! END DO
	4629	!>jyg
[1992]	4630
	4631
[2508]	4632	!jyg<
	4633	!-----------------------------------------------------------
	4634	IF (.NOT.ok_optim_yield) THEN !\|
	4635	!-----------------------------------------------------------
[1992]	4636	DO i = 1, nl
	4637	DO il = 1, ncum
[2508]	4638	upwd(il, i) = 0.0
	4639	dnwd(il, i) = 0.0
[1992]	4640	END DO
	4641	END DO
	4642
[2508]	4643	!! DO i = 1, nl ! useless; jyg
	4644	!! DO il = 1, ncum ! useless; jyg
	4645	!! IF (i>=icb(il) .AND. i<=inb(il)) THEN ! useless; jyg
	4646	!! upwd(il, i) = 0.0 ! useless; jyg
	4647	!! dnwd(il, i) = 0.0 ! useless; jyg
	4648	!! END IF ! useless; jyg
	4649	!! END DO ! useless; jyg
	4650	!! END DO ! useless; jyg
	4651
[1992]	4652	DO i = 1, nl
	4653	DO k = 1, nl
	4654	DO il = 1, ncum
	4655	up1(il, k, i) = 0.0
	4656	dn1(il, k, i) = 0.0
	4657	END DO
	4658	END DO
	4659	END DO
	4660
[2508]	4661	!yor! commented original
	4662	! DO i = 1, nl
	4663	! DO k = i, nl
	4664	! DO n = 1, i - 1
	4665	! DO il = 1, ncum
	4666	! IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
	4667	! up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
	4668	! dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
	4669	! END IF
	4670	! END DO
	4671	! END DO
	4672	! END DO
	4673	! END DO
	4674	!yor! replaced with
[1992]	4675	DO i = 1, nl
	4676	DO k = i, nl
	4677	DO n = 1, i - 1
	4678	DO il = 1, ncum
[2508]	4679	IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! as i always <= k
	4680	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
[1992]	4681	END IF
	4682	END DO
	4683	END DO
	4684	END DO
	4685	END DO
[2508]	4686	DO i = 1, nl
	4687	DO n = 1, i - 1
	4688	DO k = i, nl
	4689	DO il = 1, ncum
	4690	IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! i always <= k
	4691	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
	4692	END IF
	4693	END DO
	4694	END DO
	4695	END DO
	4696	END DO
	4697	!yor! end replace
[1992]	4698
	4699	DO i = 1, nl
	4700	DO k = 1, nl
	4701	DO il = 1, ncum
	4702	IF (i>=icb(il)) THEN
	4703	IF (k>=i .AND. k<=(inb(il))) THEN
	4704	upwd(il, i) = upwd(il, i) + m(il, k)
	4705	END IF
	4706	ELSE
	4707	IF (k<i) THEN
	4708	upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
	4709	END IF
	4710	END IF
[2007]	4711	! c print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
[1992]	4712	END DO
	4713	END DO
	4714	END DO
	4715
	4716	DO i = 2, nl
	4717	DO k = i, nl
	4718	DO il = 1, ncum
[2007]	4719	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
[1992]	4720	IF (i<=inb(il) .AND. k<=inb(il)) THEN
	4721	upwd(il, i) = upwd(il, i) + up1(il, k, i)
	4722	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
	4723	END IF
[2007]	4724	! c print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
[1992]	4725	END DO
	4726	END DO
	4727	END DO
	4728
	4729
[2007]	4730	!!!! DO il=1,ncum
	4731	!!!! do i=icb(il),inb(il)
	4732	!!!!
	4733	!!!! upwd(il,i)=0.0
	4734	!!!! dnwd(il,i)=0.0
	4735	!!!! do k=i,inb(il)
	4736	!!!! up1=0.0
	4737	!!!! dn1=0.0
	4738	!!!! do n=1,i-1
	4739	!!!! up1=up1+ment(il,n,k)
	4740	!!!! dn1=dn1-ment(il,k,n)
	4741	!!!! enddo
	4742	!!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
	4743	!!!! dnwd(il,i)=dnwd(il,i)+dn1
	4744	!!!! enddo
	4745	!!!! enddo
	4746	!!!!
	4747	!!!! ENDDO
[3605]	4748
	4749	!! DO i = 1, nlp
	4750	!! DO il = 1, ncum
	4751	!! ma(il, i) = 0
	4752	!! END DO
	4753	!! END DO
	4754	!!
	4755	!! DO i = 1, nl
	4756	!! DO j = i, nl
	4757	!! DO il = 1, ncum
	4758	!! ma(il, i) = ma(il, i) + m(il, j)
	4759	!! END DO
	4760	!! END DO
	4761	!! END DO
	4762
	4763	!jyg< (loops stop at nl)
	4764	!! DO i = nl + 1, nd
	4765	!! DO il = 1, ncum
	4766	!! ma(il, i) = 0.
	4767	!! END DO
	4768	!! END DO
	4769	!>jyg
	4770
	4771	!! DO i = 1, nl
	4772	!! DO il = 1, ncum
	4773	!! IF (i<=(icb(il)-1)) THEN
	4774	!! ma(il, i) = 0
	4775	!! END IF
	4776	!! END DO
	4777	!! END DO
	4778
[2508]	4779	!-----------------------------------------------------------
	4780	ENDIF !(.NOT.ok_optim_yield) !\|
	4781	!-----------------------------------------------------------
	4782	!>jyg
[1992]	4783
[2007]	4784	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	4785	! determination de la variation de flux ascendant entre
	4786	! deux niveau non dilue mip
	4787	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1992]	4788
	4789	DO i = 1, nl
	4790	DO il = 1, ncum
	4791	mip(il, i) = m(il, i)
	4792	END DO
	4793	END DO
	4794
[2393]	4795	!jyg< (loops stop at nl)
	4796	!! DO i = nl + 1, nd
	4797	!! DO il = 1, ncum
	4798	!! mip(il, i) = 0.
	4799	!! END DO
	4800	!! END DO
	4801	!>jyg
[1992]	4802
	4803
[2007]	4804	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	4805	! icb represente de niveau ou se trouve la
	4806	! base du nuage , et inb le top du nuage
	4807	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1992]	4808
[2259]	4809	!! DO i = 1, nd ! unused . jyg
	4810	!! DO il = 1, ncum ! unused . jyg
	4811	!! mke(il, i) = upwd(il, i) + dnwd(il, i) ! unused . jyg
	4812	!! END DO ! unused . jyg
	4813	!! END DO ! unused . jyg
[1992]	4814
[2259]	4815	!! DO i = 1, nd ! unused . jyg
	4816	!! DO il = 1, ncum ! unused . jyg
	4817	!! rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd(1.-rr(il,i))+rr(il,i)cpv) ! unused . jyg
	4818	!! tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp ! unused . jyg
	4819	!! tps(il, i) = tp(il, i) ! unused . jyg
	4820	!! END DO ! unused . jyg
	4821	!! END DO ! unused . jyg
[1992]	4822
	4823
[2007]	4824	! * diagnose the in-cloud mixing ratio * ! cld
	4825	! * of condensed water * ! cld
	4826	!! cld
	4827
[2393]	4828	DO i = 1, nl+1 ! cld
[2007]	4829	DO il = 1, ncum ! cld
	4830	mac(il, i) = 0.0 ! cld
	4831	wa(il, i) = 0.0 ! cld
	4832	siga(il, i) = 0.0 ! cld
	4833	sax(il, i) = 0.0 ! cld
	4834	END DO ! cld
	4835	END DO ! cld
	4836
	4837	DO i = minorig, nl ! cld
	4838	DO k = i + 1, nl + 1 ! cld
	4839	DO il = 1, ncum ! cld
[1992]	4840	IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
[2007]	4841	mac(il, i) = mac(il, i) + m(il, k) ! cld
	4842	END IF ! cld
	4843	END DO ! cld
	4844	END DO ! cld
	4845	END DO ! cld
[1992]	4846
[2007]	4847	DO i = 1, nl ! cld
	4848	DO j = 1, i ! cld
	4849	DO il = 1, ncum ! cld
	4850	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
	4851	.AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
	4852	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
	4853	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
	4854	END IF ! cld
	4855	END DO ! cld
	4856	END DO ! cld
	4857	END DO ! cld
[1992]	4858
[2007]	4859	DO i = 1, nl ! cld
	4860	DO il = 1, ncum ! cld
	4861	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
	4862	.AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
	4863	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
	4864	END IF ! cld
	4865	END DO ! cld
[2205]	4866	END DO
	4867	! cld
	4868	DO i = 1, nl
[1992]	4869
[2205]	4870	! 14/01/15 AJ je remets les parties manquantes cf JYG
	4871	! Initialize sument to 0
	4872
	4873	DO il = 1,ncum
	4874	sument(il) = 0.
	4875	ENDDO
	4876
	4877	! Sum mixed mass fluxes in sument
	4878
	4879	DO k = 1,nl
	4880	DO il = 1,ncum
	4881	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
	4882	sument(il) =sument(il) + abs(ment(il,k,i))
	4883	ENDIF
	4884	ENDDO ! il
	4885	ENDDO ! k
	4886
	4887	! 14/01/15 AJ delta n'a rien à faire là...
[2007]	4888	DO il = 1, ncum ! cld
[3605]	4889	!! IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
	4890	!! siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
	4891	!! rrdtvp(il, i)/p(il, i)/100. ! cld
	4892	!!
	4893	!! siga(il, i) = min(siga(il,i), 1.0) ! cld
	4894	sigaq = 0.
	4895	IF (wa(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
[2205]	4896	siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
[3605]	4897	rrdtvp(il, i)/p(il, i)/100. ! cld
	4898	siga(il, i) = min(siga(il,i), 1.0) ! cld
	4899	sigaq = siga(il,i)*qta(il,i-1) ! cld
	4900	ENDIF
[2205]	4901
	4902	! IM cf. FH
	4903	! 14/01/15 AJ ne correspond pas à ce qui a été codé par JYG et SB
	4904
[2007]	4905	IF (iflag_clw==0) THEN ! cld
	4906	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
	4907	+(1.-siga(il,i))*qcond(il, i) ! cld
[2205]	4908
	4909
	4910	sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1)) ! cld
	4911	sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i)) ! cld
[3605]	4912	!! qtc(il, i) = (siga(il,i)qta(il,i-1)+sigment(il,i)qtment(il,i)) & ! cld
	4913	qtc(il, i) = (sigaq+sigment(il,i)*qtment(il,i)) & ! cld
[2205]	4914	/(siga(il,i)+sigment(il,i)) ! cld
	4915	sigt(il,i) = sigment(il, i) + siga(il, i)
	4916
[3605]	4917	! qtc(il, i) = siga(il,i)qta(il,i-1)+(1.-siga(il,i))qtment(il,i) ! cld
[2208]	4918	! print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)
[2205]	4919
[2007]	4920	ELSE IF (iflag_clw==1) THEN ! cld
	4921	qcondc(il, i) = qcond(il, i) ! cld
[2205]	4922	qtc(il,i) = qtment(il,i) ! cld
[2007]	4923	END IF ! cld
[1992]	4924
[2007]	4925	END DO ! cld
[1992]	4926	END DO
[2007]	4927	! print*,'cv3_yield fin'
	4928
[1992]	4929	RETURN
	4930	END SUBROUTINE cv3_yield
	4931
[2007]	4932	!AC! et !RomP >>>
	4933	SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
	4934	ment, sigij, da, phi, phi2, d1a, dam, &
	4935	ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
	4936	icb, inb)
[1992]	4937	IMPLICIT NONE
	4938
	4939	include "cv3param.h"
	4940
[2007]	4941	!inputs:
[3126]	4942	INTEGER, INTENT (IN) :: ncum, nd, na, nloc, len
	4943	INTEGER, DIMENSION (len), INTENT (IN) :: icb, inb
	4944	REAL, DIMENSION (len, na, na), INTENT (IN) :: ment, sigij, elij
	4945	REAL, DIMENSION (len, nd), INTENT (IN) :: clw
	4946	REAL, DIMENSION (len, na), INTENT (IN) :: ep
	4947	REAL, DIMENSION (len, nd+1), INTENT (IN) :: Vprecip
[2007]	4948	!ouputs:
[3126]	4949	REAL, DIMENSION (len, na, na), INTENT (OUT) :: phi, phi2, epmlmMm
	4950	REAL, DIMENSION (len, na), INTENT (OUT) :: da, d1a, dam, eplaMm
	4951	!
[2007]	4952	! variables pour tracer dans precip de l'AA et des mel
	4953	!local variables:
[1992]	4954	INTEGER i, j, k
	4955	REAL epm(nloc, na, na)
	4956
[2007]	4957	! variables d'Emanuel : du second indice au troisieme
	4958	! ---> tab(i,k,j) -> de l origine k a l arrivee j
	4959	! ment, sigij, elij
	4960	! variables personnelles : du troisieme au second indice
	4961	! ---> tab(i,j,k) -> de k a j
	4962	! phi, phi2
[1992]	4963
[2007]	4964	! initialisations
[1992]	4965
	4966	da(:, :) = 0.
	4967	d1a(:, :) = 0.
	4968	dam(:, :) = 0.
	4969	epm(:, :, :) = 0.
[2007]	4970	eplaMm(:, :) = 0.
	4971	epmlmMm(:, :, :) = 0.
[1992]	4972	phi(:, :, :) = 0.
	4973	phi2(:, :, :) = 0.
	4974
[2007]	4975	! fraction deau condensee dans les melanges convertie en precip : epm
	4976	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
[2393]	4977	DO j = 1, nl
	4978	DO k = 1, nl
[1992]	4979	DO i = 1, ncum
[2007]	4980	IF (k>=icb(i) .AND. k<=inb(i) .AND. &
	4981	!!jyg j.ge.k.and.j.le.inb(i)) then
	4982	!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
[1992]	4983	j>k .AND. j<=inb(i)) THEN
	4984	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
[2007]	4985	!!
[1992]	4986	epm(i, j, k) = max(epm(i,j,k), 0.0)
	4987	END IF
	4988	END DO
	4989	END DO
	4990	END DO
	4991
	4992
[2393]	4993	DO j = 1, nl
	4994	DO k = 1, nl
[1992]	4995	DO i = 1, ncum
	4996	IF (k>=icb(i) .AND. k<=inb(i)) THEN
[2007]	4997	eplaMm(i, j) = eplamm(i, j) + &
	4998	ep(i, j)clw(i, j)ment(i, j, k)*(1.-sigij(i,j,k))
[1992]	4999	END IF
	5000	END DO
	5001	END DO
	5002	END DO
	5003
[2393]	5004	DO j = 1, nl
[1992]	5005	DO k = 1, j - 1
	5006	DO i = 1, ncum
	5007	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
[2007]	5008	epmlmMm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
[1992]	5009	END IF
	5010	END DO
	5011	END DO
	5012	END DO
	5013
[2007]	5014	! matrices pour calculer la tendance des concentrations dans cvltr.F90
[2393]	5015	DO j = 1, nl
	5016	DO k = 1, nl
[1992]	5017	DO i = 1, ncum
	5018	da(i, j) = da(i, j) + (1.-sigij(i,k,j))*ment(i, k, j)
	5019	phi(i, j, k) = sigij(i, k, j)*ment(i, k, j)
	5020	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sigij(i,k,j))
	5021	IF (k<=j) THEN
[2007]	5022	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, k, j)(1.-ep(i,k))*(1.-sigij(i,k,j))
[1992]	5023	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
	5024	END IF
	5025	END DO
	5026	END DO
	5027	END DO
	5028
	5029	RETURN
	5030	END SUBROUTINE cv3_tracer
[2007]	5031	!AC! et !RomP <<<
[1992]	5032
[2007]	5033	SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
	5034	iflag, &
	5035	precip, sig, w0, &
	5036	ft, fq, fu, fv, ftra, &
	5037	Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
[2481]	5038	epmax_diag, & ! epmax_cape
[2007]	5039	iflag1, &
	5040	precip1, sig1, w01, &
	5041	ft1, fq1, fu1, fv1, ftra1, &
[2481]	5042	Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1, &
	5043	epmax_diag1) ! epmax_cape
[1992]	5044	IMPLICIT NONE
	5045
	5046	include "cv3param.h"
	5047
[2007]	5048	!inputs:
[1992]	5049	INTEGER len, ncum, nd, ntra, nloc
	5050	INTEGER idcum(nloc)
	5051	INTEGER iflag(nloc)
	5052	REAL precip(nloc)
	5053	REAL sig(nloc, nd), w0(nloc, nd)
	5054	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	5055	REAL ftra(nloc, nd, ntra)
	5056	REAL ma(nloc, nd)
	5057	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
	5058	REAL qcondc(nloc, nd)
	5059	REAL wd(nloc), cape(nloc)
[2481]	5060	REAL epmax_diag(nloc)
[1992]	5061
[2007]	5062	!outputs:
[1992]	5063	INTEGER iflag1(len)
	5064	REAL precip1(len)
	5065	REAL sig1(len, nd), w01(len, nd)
	5066	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
	5067	REAL ftra1(len, nd, ntra)
	5068	REAL ma1(len, nd)
	5069	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
	5070	REAL qcondc1(nloc, nd)
	5071	REAL wd1(nloc), cape1(nloc)
[2481]	5072	REAL epmax_diag1(len) ! epmax_cape
[1992]	5073
[2007]	5074	!local variables:
[1992]	5075	INTEGER i, k, j
	5076
	5077	DO i = 1, ncum
	5078	precip1(idcum(i)) = precip(i)
	5079	iflag1(idcum(i)) = iflag(i)
	5080	wd1(idcum(i)) = wd(i)
	5081	cape1(idcum(i)) = cape(i)
[2481]	5082	epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
[1992]	5083	END DO
	5084
	5085	DO k = 1, nl
	5086	DO i = 1, ncum
	5087	sig1(idcum(i), k) = sig(i, k)
	5088	w01(idcum(i), k) = w0(i, k)
	5089	ft1(idcum(i), k) = ft(i, k)
	5090	fq1(idcum(i), k) = fq(i, k)
	5091	fu1(idcum(i), k) = fu(i, k)
	5092	fv1(idcum(i), k) = fv(i, k)
	5093	ma1(idcum(i), k) = ma(i, k)
	5094	upwd1(idcum(i), k) = upwd(i, k)
	5095	dnwd1(idcum(i), k) = dnwd(i, k)
	5096	dnwd01(idcum(i), k) = dnwd0(i, k)
	5097	qcondc1(idcum(i), k) = qcondc(i, k)
	5098	END DO
	5099	END DO
	5100
	5101	DO i = 1, ncum
	5102	sig1(idcum(i), nd) = sig(i, nd)
	5103	END DO
	5104
	5105
[2007]	5106	!AC! do 2100 j=1,ntra
	5107	!AC!c oct3 do 2110 k=1,nl
	5108	!AC! do 2110 k=1,nd ! oct3
	5109	!AC! do 2120 i=1,ncum
	5110	!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
	5111	!AC! 2120 continue
	5112	!AC! 2110 continue
	5113	!AC! 2100 continue
	5114	!
[1992]	5115	RETURN
	5116	END SUBROUTINE cv3_uncompress
[2481]	5117
	5118
	5119	subroutine cv3_epmax_fn_cape(nloc,ncum,nd &
	5120	, ep,hp,icb,inb,clw,nk,t,h,hnk,lv,lf,frac &
	5121	, pbase, p, ph, tv, buoy, sig, w0,iflag &
	5122	, epmax_diag)
	5123	implicit none
	5124
	5125	! On fait varier epmax en fn de la cape
	5126	! Il faut donc recalculer ep, et hp qui a déjà été calculé et
	5127	! qui en dépend
	5128	! Toutes les autres variables fn de ep sont calculées plus bas.
	5129
	5130	include "cvthermo.h"
	5131	include "cv3param.h"
	5132	include "conema3.h"
	5133	include "cvflag.h"
	5134
	5135	! inputs:
	5136	INTEGER, INTENT (IN) :: ncum, nd, nloc
	5137	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
	5138	REAL, DIMENSION (nloc), INTENT (IN) :: hnk,pbase
	5139	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, lv, lf, tv, h
	5140	REAL, DIMENSION (nloc, nd), INTENT (IN) :: clw, buoy,frac
	5141	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig,w0
	5142	INTEGER, DIMENSION (nloc), INTENT (IN) :: iflag(nloc)
	5143	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
	5144	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	5145	! inouts:
	5146	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: ep,hp
	5147	! outputs
	5148	REAL, DIMENSION (nloc), INTENT (OUT) :: epmax_diag
	5149
	5150	! local
	5151	integer i,k
	5152	! real hp_bak(nloc,nd)
	5153	! real ep_bak(nloc,nd)
	5154	real m_loc(nloc,nd)
	5155	real sig_loc(nloc,nd)
	5156	real w0_loc(nloc,nd)
	5157	integer iflag_loc(nloc)
	5158	real cape(nloc)
	5159
	5160	if (coef_epmax_cape.gt.1e-12) then
	5161
	5162	! il faut calculer la cape: on fait un calcule simple car tant qu'on ne
	5163	! connait pas ep, on ne connait pas les mélanges, ddfts etc... qui sont
	5164	! necessaires au calcul de la cape dans la nouvelle physique
	5165
	5166	! write(,) 'cv3_routines check 4303'
	5167	do i=1,ncum
	5168	do k=1,nd
	5169	sig_loc(i,k)=sig(i,k)
	5170	w0_loc(i,k)=w0(i,k)
	5171	iflag_loc(i)=iflag(i)
	5172	! ep_bak(i,k)=ep(i,k)
	5173	enddo ! do k=1,nd
	5174	enddo !do i=1,ncum
	5175
	5176	! write(,) 'cv3_routines check 4311'
	5177	! write(,) 'nl=',nl
	5178	CALL cv3_closure(nloc, ncum, nd, icb, inb, & ! na->nd
	5179	pbase, p, ph, tv, buoy, &
	5180	sig_loc, w0_loc, cape, m_loc,iflag_loc)
	5181
	5182	! write(,) 'cv3_routines check 4316'
	5183	! write(,) 'ep(1,:)=',ep(1,:)
	5184	do i=1,ncum
	5185	epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
	5186	epmax_diag(i)=amax1(epmax_diag(i),0.0)
	5187	! write(,) 'i,icb,inb,cape,epmax_diag=', &
	5188	! i,icb(i),inb(i),cape(i),epmax_diag(i)
	5189	do k=1,nl
	5190	ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
	5191	ep(i,k)=amax1(ep(i,k),0.0)
	5192	ep(i,k)=amin1(ep(i,k),epmax_diag(i))
	5193	enddo
	5194	enddo
	5195	! write(,) 'ep(1,:)=',ep(1,:)
	5196
	5197	!write(,) 'cv3_routines check 4326'
	5198	! On recalcule hp:
	5199	! do k=1,nl
	5200	! do i=1,ncum
	5201	! hp_bak(i,k)=hp(i,k)
	5202	! enddo
	5203	! enddo
	5204	do k=1,nl
	5205	do i=1,ncum
[3605]	5206	hp(i,k)=h(i,k)
	5207	enddo
[2481]	5208	enddo
	5209
	5210	IF (cvflag_ice) THEN
	5211
	5212	do k=minorig+1,nl
	5213	do i=1,ncum
	5214	if((k.ge.icb(i)).and.(k.le.inb(i)))then
	5215	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac(i,k)lf(i,k))* &
	5216	ep(i, k)*clw(i, k)
	5217	endif
	5218	enddo
	5219	enddo !do k=minorig+1,n
	5220	ELSE !IF (cvflag_ice) THEN
	5221
	5222	DO k = minorig + 1, nl
	5223	DO i = 1, ncum
	5224	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
	5225	hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
	5226	endif
	5227	enddo
	5228	enddo !do k=minorig+1,n
	5229
	5230	ENDIF !IF (cvflag_ice) THEN
	5231	!write(,) 'cv3_routines check 4345'
	5232	! do i=1,ncum
	5233	! do k=1,nl
	5234	! if ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-1).or. &
	5235	! ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-4).and. &
	5236	! (ep(i,k)-ep_bak(i,k).lt.1e-4))) then
	5237	! write(,) 'i,k=',i,k
	5238	! write(,) 'coef_epmax_cape=',coef_epmax_cape
	5239	! write(,) 'epmax_diag(i)=',epmax_diag(i)
	5240	! write(,) 'ep(i,k)=',ep(i,k)
	5241	! write(,) 'ep_bak(i,k)=',ep_bak(i,k)
	5242	! write(,) 'hp(i,k)=',hp(i,k)
	5243	! write(,) 'hp_bak(i,k)=',hp_bak(i,k)
	5244	! write(,) 'h(i,k)=',h(i,k)
	5245	! write(,) 'nk(i)=',nk(i)
	5246	! write(,) 'h(i,nk(i))=',h(i,nk(i))
	5247	! write(,) 'lv(i,k)=',lv(i,k)
	5248	! write(,) 't(i,k)=',t(i,k)
	5249	! write(,) 'clw(i,k)=',clw(i,k)
	5250	! write(,) 'cpd,cpv=',cpd,cpv
	5251	! stop
	5252	! endif
	5253	! enddo !do k=1,nl
	5254	! enddo !do i=1,ncum
	5255	endif !if (coef_epmax_cape.gt.1e-12) then
	5256	!write(,) 'cv3_routines check 4367'
	5257
	5258	return
	5259	end subroutine cv3_epmax_fn_cape
	5260
	5261
	5262

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

Download in other formats:

Original Format