Context Navigation

cv3_routines.F90 @ 3449

Last change on this file since 3449 was 3435, checked in by Laurent Fairhead, 5 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

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

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

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

Download in other formats:

Original Format