Context Navigation

cv3_routines.F90 @ 5451

Last change on this file since 5451 was 4727, checked in by idelkadi, 15 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 180.1 KB

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

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

Download in other formats:

Original Format