Context Navigation

cv30_routines.F90 @ 2343

Last change on this file since 2343 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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

Rev	Line
[1992]	1
[1403]	2	! $Id: cv30_routines.F90 2311 2015-06-25 07:45:24Z emillour $
[879]	3
	4
	5
[1992]	6	SUBROUTINE cv30_param(nd, delt)
	7	IMPLICIT NONE
[879]	8
[1992]	9	! ------------------------------------------------------------
	10	! Set parameters for convectL for iflag_con = 3
	11	! ------------------------------------------------------------
[879]	12
	13
[1992]	14	! * PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *
	15	! * PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *
	16	! * PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *
	17	! * EFFICIENCY IS ASSUMED TO BE UNITY *
	18	! * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
	19	! * SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
	20	! * OF CLOUD *
[879]	21
[1992]	22	! [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
	23	! * ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *
	24	! * APPROACH TO QUASI-EQUILIBRIUM *
	25	! * (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *
	26	! * (BETA MUST BE LESS THAN OR EQUAL TO 1) *
[879]	27
[1992]	28	! * DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *
	29	! * APPROACH TO QUASI-EQUILIBRIUM *
	30	! * IT MUST BE LESS THAN 0 *
[879]	31
[1992]	32	include "cv30param.h"
	33	include "conema3.h"
[879]	34
[1992]	35	INTEGER nd
	36	REAL delt ! timestep (seconds)
[879]	37
[1992]	38	! noff: integer limit for convection (nd-noff)
	39	! minorig: First level of convection
[879]	40
[1992]	41	! -- limit levels for convection:
[879]	42
[1992]	43	noff = 1
	44	minorig = 1
	45	nl = nd - noff
	46	nlp = nl + 1
	47	nlm = nl - 1
[879]	48
[1992]	49	! -- "microphysical" parameters:
[879]	50
[1992]	51	sigd = 0.01
	52	spfac = 0.15
	53	pbcrit = 150.0
	54	ptcrit = 500.0
	55	! IM cf. FH epmax = 0.993
[879]	56
[1992]	57	omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
[879]	58
[1992]	59	! -- misc:
[879]	60
[1992]	61	dtovsh = -0.2 ! dT for overshoot
	62	dpbase = -40. ! definition cloud base (400m above LCL)
	63	dttrig = 5. ! (loose) condition for triggering
[879]	64
[1992]	65	! -- rate of approach to quasi-equilibrium:
[879]	66
[1992]	67	dtcrit = -2.0
	68	tau = 8000.
	69	beta = 1.0 - delt/tau
	70	alpha = 1.5E-3*delt/tau
	71	! increase alpha to compensate W decrease:
	72	alpha = alpha*1.5
[879]	73
[1992]	74	! -- interface cloud parameterization:
[879]	75
[1992]	76	delta = 0.01 ! cld
[879]	77
[1992]	78	! -- interface with boundary-layer (gust factor): (sb)
[879]	79
[1992]	80	betad = 10.0 ! original value (from convect 4.3)
[879]	81
[1992]	82	RETURN
	83	END SUBROUTINE cv30_param
[879]	84
[1992]	85	SUBROUTINE cv30_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm, &
	86	th)
	87	IMPLICIT NONE
[879]	88
[1992]	89	! =====================================================================
	90	! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
	91	! "ori": from convect4.3 (vectorized)
	92	! "convect3": to be exactly consistent with convect3
	93	! =====================================================================
[879]	94
[1992]	95	! inputs:
	96	INTEGER len, nd, ndp1
	97	REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
[879]	98
[1992]	99	! outputs:
	100	REAL lv(len, nd), cpn(len, nd), tv(len, nd)
	101	REAL gz(len, nd), h(len, nd), hm(len, nd)
	102	REAL th(len, nd)
[879]	103
[1992]	104	! local variables:
	105	INTEGER k, i
	106	REAL rdcp
	107	REAL tvx, tvy ! convect3
	108	REAL cpx(len, nd)
[879]	109
[1992]	110	include "cvthermo.h"
	111	include "cv30param.h"
[879]	112
	113
[1992]	114	! ori do 110 k=1,nlp
	115	DO k = 1, nl ! convect3
	116	DO i = 1, len
	117	! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
	118	lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
	119	cpn(i, k) = cpd(1.0-q(i,k)) + cpvq(i, k)
	120	cpx(i, k) = cpd(1.0-q(i,k)) + clq(i, k)
	121	! ori tv(i,k)=t(i,k)(1.0+q(i,k)epsim1)
	122	tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
	123	rdcp = (rrd(1.-q(i,k))+q(i,k)rrv)/cpn(i, k)
	124	th(i, k) = t(i, k)(1000.0/p(i,k))*rdcp
	125	END DO
	126	END DO
[879]	127
[1992]	128	! gz = phi at the full levels (same as p).
[879]	129
[1992]	130	DO i = 1, len
	131	gz(i, 1) = 0.0
	132	END DO
	133	! ori do 140 k=2,nlp
	134	DO k = 2, nl ! convect3
	135	DO i = 1, len
	136	tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3
	137	tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
	138	gz(i, k) = gz(i, k-1) + 0.5rrd(tvx+tvy) & !convect3
	139	*(p(i,k-1)-p(i,k))/ph(i, k) !convect3
[879]	140
[1992]	141	! ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
	142	! ori & *(p(i,k-1)-p(i,k))/ph(i,k)
	143	END DO
	144	END DO
[879]	145
[1992]	146	! h = phi + cpT (dry static energy).
	147	! hm = phi + cp(T-Tbase)+Lq
[879]	148
[1992]	149	! ori do 170 k=1,nlp
	150	DO k = 1, nl ! convect3
	151	DO i = 1, len
	152	h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
	153	hm(i, k) = gz(i, k) + cpx(i, k)(t(i,k)-t(i,1)) + lv(i, k)q(i, k)
	154	END DO
	155	END DO
[879]	156
[1992]	157	RETURN
	158	END SUBROUTINE cv30_prelim
[879]	159
[1992]	160	SUBROUTINE cv30_feed(len, nd, t, q, qs, p, ph, hm, gz, nk, icb, icbmax, &
	161	iflag, tnk, qnk, gznk, plcl)
	162	IMPLICIT NONE
[879]	163
[1992]	164	! ================================================================
	165	! Purpose: CONVECTIVE FEED
[879]	166
[1992]	167	! Main differences with cv_feed:
	168	! - ph added in input
	169	! - here, nk(i)=minorig
	170	! - icb defined differently (plcl compared with ph instead of p)
[879]	171
[1992]	172	! Main differences with convect3:
	173	! - we do not compute dplcldt and dplcldr of CLIFT anymore
	174	! - values iflag different (but tests identical)
	175	! - A,B explicitely defined (!...)
	176	! ================================================================
[879]	177
[1992]	178	include "cv30param.h"
[879]	179
[1992]	180	! inputs:
	181	INTEGER len, nd
	182	REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
	183	REAL hm(len, nd), gz(len, nd)
	184	REAL ph(len, nd+1)
[879]	185
[1992]	186	! outputs:
	187	INTEGER iflag(len), nk(len), icb(len), icbmax
	188	REAL tnk(len), qnk(len), gznk(len), plcl(len)
[879]	189
[1992]	190	! local variables:
	191	INTEGER i, k
	192	INTEGER ihmin(len)
	193	REAL work(len)
	194	REAL pnk(len), qsnk(len), rh(len), chi(len)
	195	REAL a, b ! convect3
	196	! ym
	197	plcl = 0.0
	198	! @ !-------------------------------------------------------------------
	199	! @ ! --- Find level of minimum moist static energy
	200	! @ ! --- If level of minimum moist static energy coincides with
	201	! @ ! --- or is lower than minimum allowable parcel origin level,
	202	! @ ! --- set iflag to 6.
	203	! @ !-------------------------------------------------------------------
	204	! @
	205	! @ do 180 i=1,len
	206	! @ work(i)=1.0e12
	207	! @ ihmin(i)=nl
	208	! @ 180 continue
	209	! @ do 200 k=2,nlp
	210	! @ do 190 i=1,len
	211	! @ if((hm(i,k).lt.work(i)).and.
	212	! @ & (hm(i,k).lt.hm(i,k-1)))then
	213	! @ work(i)=hm(i,k)
	214	! @ ihmin(i)=k
	215	! @ endif
	216	! @ 190 continue
	217	! @ 200 continue
	218	! @ do 210 i=1,len
	219	! @ ihmin(i)=min(ihmin(i),nlm)
	220	! @ if(ihmin(i).le.minorig)then
	221	! @ iflag(i)=6
	222	! @ endif
	223	! @ 210 continue
	224	! @ c
	225	! @ !-------------------------------------------------------------------
	226	! @ ! --- Find that model level below the level of minimum moist static
	227	! @ ! --- energy that has the maximum value of moist static energy
	228	! @ !-------------------------------------------------------------------
	229	! @
	230	! @ do 220 i=1,len
	231	! @ work(i)=hm(i,minorig)
	232	! @ nk(i)=minorig
	233	! @ 220 continue
	234	! @ do 240 k=minorig+1,nl
	235	! @ do 230 i=1,len
	236	! @ if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
	237	! @ work(i)=hm(i,k)
	238	! @ nk(i)=k
	239	! @ endif
	240	! @ 230 continue
	241	! @ 240 continue
[879]	242
[1992]	243	! -------------------------------------------------------------------
	244	! --- Origin level of ascending parcels for convect3:
	245	! -------------------------------------------------------------------
[879]	246
[1992]	247	DO i = 1, len
	248	nk(i) = minorig
	249	END DO
[879]	250
[1992]	251	! -------------------------------------------------------------------
	252	! --- Check whether parcel level temperature and specific humidity
	253	! --- are reasonable
	254	! -------------------------------------------------------------------
	255	DO i = 1, len
	256	IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0)) & ! @ & .or.(
	257	! p(i,ihmin(i)).lt.400.0
	258	! ) )
	259	.AND. (iflag(i)==0)) iflag(i) = 7
	260	END DO
	261	! -------------------------------------------------------------------
	262	! --- Calculate lifted condensation level of air at parcel origin level
	263	! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
	264	! -------------------------------------------------------------------
[879]	265
[1992]	266	a = 1669.0 ! convect3
	267	b = 122.0 ! convect3
[879]	268
[1992]	269	DO i = 1, len
[879]	270
[1992]	271	IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002
[879]	272
[1992]	273	tnk(i) = t(i, nk(i))
	274	qnk(i) = q(i, nk(i))
	275	gznk(i) = gz(i, nk(i))
	276	pnk(i) = p(i, nk(i))
	277	qsnk(i) = qs(i, nk(i))
[879]	278
[1992]	279	rh(i) = qnk(i)/qsnk(i)
	280	! ori rh(i)=min(1.0,rh(i)) ! removed for convect3
	281	! ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
	282	chi(i) = tnk(i)/(a-b*rh(i)-tnk(i)) ! convect3
	283	plcl(i) = pnk(i)(rh(i)*chi(i))
	284	IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag &
	285	(i) = 8
[879]	286
[1992]	287	END IF ! iflag=7
[879]	288
[1992]	289	END DO
[879]	290
[1992]	291	! -------------------------------------------------------------------
	292	! --- Calculate first level above lcl (=icb)
	293	! -------------------------------------------------------------------
[879]	294
[1992]	295	! @ do 270 i=1,len
	296	! @ icb(i)=nlm
	297	! @ 270 continue
	298	! @c
	299	! @ do 290 k=minorig,nl
	300	! @ do 280 i=1,len
	301	! @ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
	302	! @ & icb(i)=min(icb(i),k)
	303	! @ 280 continue
	304	! @ 290 continue
	305	! @c
	306	! @ do 300 i=1,len
	307	! @ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
	308	! @ 300 continue
[879]	309
[1992]	310	DO i = 1, len
	311	icb(i) = nlm
	312	END DO
[879]	313
[1992]	314	! la modification consiste a comparer plcl a ph et non a p:
	315	! icb est defini par : ph(icb)<plcl<ph(icb-1)
	316	! @ do 290 k=minorig,nl
	317	DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
	318	DO i = 1, len
	319	IF (ph(i,k)<plcl(i)) icb(i) = min(icb(i), k)
	320	END DO
	321	END DO
[879]	322
[1992]	323	DO i = 1, len
	324	! @ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
	325	IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
	326	END DO
[879]	327
[1992]	328	DO i = 1, len
	329	icb(i) = icb(i) - 1 ! icb sup ou egal a 2
	330	END DO
[879]	331
[1992]	332	! Compute icbmax.
[879]	333
[1992]	334	icbmax = 2
	335	DO i = 1, len
	336	! ! icbmax=max(icbmax,icb(i))
	337	IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
	338	END DO
[879]	339
[1992]	340	RETURN
	341	END SUBROUTINE cv30_feed
[879]	342
[1992]	343	SUBROUTINE cv30_undilute1(len, nd, t, q, qs, gz, plcl, p, nk, icb, tp, tvp, &
	344	clw, icbs)
	345	IMPLICIT NONE
[879]	346
[1992]	347	! ----------------------------------------------------------------
	348	! Equivalent de TLIFT entre NK et ICB+1 inclus
[879]	349
[1992]	350	! Differences with convect4:
	351	! - specify plcl in input
	352	! - icbs is the first level above LCL (may differ from icb)
	353	! - in the iterations, used x(icbs) instead x(icb)
	354	! - many minor differences in the iterations
	355	! - tvp is computed in only one time
	356	! - icbs: first level above Plcl (IMIN de TLIFT) in output
	357	! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
	358	! ----------------------------------------------------------------
[879]	359
[1992]	360	include "cvthermo.h"
	361	include "cv30param.h"
[879]	362
[1992]	363	! inputs:
	364	INTEGER len, nd
	365	INTEGER nk(len), icb(len)
	366	REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
	367	REAL p(len, nd)
	368	REAL plcl(len) ! convect3
[879]	369
[1992]	370	! outputs:
	371	REAL tp(len, nd), tvp(len, nd), clw(len, nd)
[879]	372
[1992]	373	! local variables:
	374	INTEGER i, k
	375	INTEGER icb1(len), icbs(len), icbsmax2 ! convect3
	376	REAL tg, qg, alv, s, ahg, tc, denom, es, rg
	377	REAL ah0(len), cpp(len)
	378	REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
	379	REAL qsicb(len) ! convect3
	380	REAL cpinv(len) ! convect3
[879]	381
[1992]	382	! -------------------------------------------------------------------
	383	! --- Calculates the lifted parcel virtual temperature at nk,
	384	! --- the actual temperature, and the adiabatic
	385	! --- liquid water content. The procedure is to solve the equation.
	386	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
	387	! -------------------------------------------------------------------
[879]	388
[1992]	389	DO i = 1, len
	390	tnk(i) = t(i, nk(i))
	391	qnk(i) = q(i, nk(i))
	392	gznk(i) = gz(i, nk(i))
	393	! ori ticb(i)=t(i,icb(i))
	394	! ori gzicb(i)=gz(i,icb(i))
	395	END DO
[879]	396
[1992]	397	! * Calculate certain parcel quantities, including static energy *
[879]	398
[1992]	399	DO i = 1, len
	400	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))tnk(i) + qnk(i)(lv0-clmcpv*(tnk(i)- &
	401	273.15)) + gznk(i)
	402	cpp(i) = cpd(1.-qnk(i)) + qnk(i)cpv
	403	cpinv(i) = 1./cpp(i)
	404	END DO
[879]	405
[1992]	406	! * Calculate lifted parcel quantities below cloud base *
[879]	407
[1992]	408	DO i = 1, len !convect3
	409	icb1(i) = max(icb(i), 2) !convect3
	410	icb1(i) = min(icb(i), nl) !convect3
	411	! if icb is below LCL, start loop at ICB+1:
	412	! (icbs est le premier niveau au-dessus du LCL)
	413	icbs(i) = icb1(i) !convect3
	414	IF (plcl(i)<p(i,icb1(i))) THEN
	415	icbs(i) = min(icbs(i)+1, nl) !convect3
	416	END IF
	417	END DO !convect3
[879]	418
[1992]	419	DO i = 1, len !convect3
	420	ticb(i) = t(i, icbs(i)) !convect3
	421	gzicb(i) = gz(i, icbs(i)) !convect3
	422	qsicb(i) = qs(i, icbs(i)) !convect3
	423	END DO !convect3
[879]	424
	425
[1992]	426	! Re-compute icbsmax (icbsmax2): !convect3
	427	! !convect3
	428	icbsmax2 = 2 !convect3
	429	DO i = 1, len !convect3
	430	icbsmax2 = max(icbsmax2, icbs(i)) !convect3
	431	END DO !convect3
[879]	432
[1992]	433	! initialization outputs:
[879]	434
[1992]	435	DO k = 1, icbsmax2 ! convect3
	436	DO i = 1, len ! convect3
	437	tp(i, k) = 0.0 ! convect3
	438	tvp(i, k) = 0.0 ! convect3
	439	clw(i, k) = 0.0 ! convect3
	440	END DO ! convect3
	441	END DO ! convect3
[879]	442
[1992]	443	! tp and tvp below cloud base:
[879]	444
[1992]	445	DO k = minorig, icbsmax2 - 1
	446	DO i = 1, len
	447	tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
	448	tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
	449	END DO
	450	END DO
[879]	451
[1992]	452	! * Find lifted parcel quantities above cloud base *
[879]	453
[1992]	454	DO i = 1, len
	455	tg = ticb(i)
	456	! ori qg=qs(i,icb(i))
	457	qg = qsicb(i) ! convect3
	458	! debug alv=lv0-clmcpv*(ticb(i)-t0)
	459	alv = lv0 - clmcpv*(ticb(i)-273.15)
[879]	460
[1992]	461	! First iteration.
[879]	462
[1992]	463	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	464	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
	465	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
	466	s = 1./s
	467	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	468	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
	469	tg = tg + s*(ah0(i)-ahg)
	470	! ori tg=max(tg,35.0)
	471	! debug tc=tg-t0
	472	tc = tg - 273.15
	473	denom = 243.5 + tc
	474	denom = max(denom, 1.0) ! convect3
	475	! ori if(tc.ge.0.0)then
	476	es = 6.112exp(17.67tc/denom)
	477	! ori else
	478	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	479	! ori endif
	480	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	481	qg = epses/(p(i,icbs(i))-es(1.-eps))
[879]	482
[1992]	483	! Second iteration.
[879]	484
	485
[1992]	486	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	487	! ori s=1./s
	488	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	489	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
	490	tg = tg + s*(ah0(i)-ahg)
	491	! ori tg=max(tg,35.0)
	492	! debug tc=tg-t0
	493	tc = tg - 273.15
	494	denom = 243.5 + tc
	495	denom = max(denom, 1.0) ! convect3
	496	! ori if(tc.ge.0.0)then
	497	es = 6.112exp(17.67tc/denom)
	498	! ori else
	499	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	500	! ori end if
	501	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	502	qg = epses/(p(i,icbs(i))-es(1.-eps))
[879]	503
[1992]	504	alv = lv0 - clmcpv*(ticb(i)-273.15)
[879]	505
[1992]	506	! ori c approximation here:
	507	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	508	! ori & -gz(i,icb(i))-alv*qg)/cpd
[879]	509
[1992]	510	! convect3: no approximation:
	511	tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alvqg)/(cpd+(cl-cpd)qnk(i))
[879]	512
[1992]	513	! ori clw(i,icb(i))=qnk(i)-qg
	514	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	515	clw(i, icbs(i)) = qnk(i) - qg
	516	clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
[879]	517
[1992]	518	rg = qg/(1.-qnk(i))
	519	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	520	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
	521	tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
[879]	522
[1992]	523	END DO
[879]	524
[1992]	525	! ori do 380 k=minorig,icbsmax2
	526	! ori do 370 i=1,len
	527	! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
	528	! ori 370 continue
	529	! ori 380 continue
[879]	530
	531
[1992]	532	! -- The following is only for convect3:
[879]	533
[1992]	534	! * icbs is the first level above the LCL:
	535	! if plcl<p(icb), then icbs=icb+1
	536	! if plcl>p(icb), then icbs=icb
[1403]	537
[1992]	538	! * the routine above computes tvp from minorig to icbs (included).
[879]	539
[1992]	540	! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
	541	! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
[879]	542
[1992]	543	! * therefore, in the case icbs=icb, we compute tvp at level icb+1
	544	! (tvp at other levels will be computed in cv3_undilute2.F)
[879]	545
	546
[1992]	547	DO i = 1, len
	548	ticb(i) = t(i, icb(i)+1)
	549	gzicb(i) = gz(i, icb(i)+1)
	550	qsicb(i) = qs(i, icb(i)+1)
	551	END DO
[879]	552
[1992]	553	DO i = 1, len
	554	tg = ticb(i)
	555	qg = qsicb(i) ! convect3
	556	! debug alv=lv0-clmcpv*(ticb(i)-t0)
	557	alv = lv0 - clmcpv*(ticb(i)-273.15)
[879]	558
[1992]	559	! First iteration.
[879]	560
[1992]	561	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	562	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
	563	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
	564	s = 1./s
	565	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	566	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
	567	tg = tg + s*(ah0(i)-ahg)
	568	! ori tg=max(tg,35.0)
	569	! debug tc=tg-t0
	570	tc = tg - 273.15
	571	denom = 243.5 + tc
	572	denom = max(denom, 1.0) ! convect3
	573	! ori if(tc.ge.0.0)then
	574	es = 6.112exp(17.67tc/denom)
	575	! ori else
	576	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	577	! ori endif
	578	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	579	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
[879]	580
[1992]	581	! Second iteration.
[879]	582
	583
[1992]	584	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	585	! ori s=1./s
	586	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	587	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
	588	tg = tg + s*(ah0(i)-ahg)
	589	! ori tg=max(tg,35.0)
	590	! debug tc=tg-t0
	591	tc = tg - 273.15
	592	denom = 243.5 + tc
	593	denom = max(denom, 1.0) ! convect3
	594	! ori if(tc.ge.0.0)then
	595	es = 6.112exp(17.67tc/denom)
	596	! ori else
	597	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	598	! ori end if
	599	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	600	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
[879]	601
[1992]	602	alv = lv0 - clmcpv*(ticb(i)-273.15)
[879]	603
[1992]	604	! ori c approximation here:
	605	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	606	! ori & -gz(i,icb(i))-alv*qg)/cpd
[879]	607
[1992]	608	! convect3: no approximation:
	609	tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alvqg)/(cpd+(cl-cpd)qnk(i))
[879]	610
[1992]	611	! ori clw(i,icb(i))=qnk(i)-qg
	612	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	613	clw(i, icb(i)+1) = qnk(i) - qg
	614	clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
[879]	615
[1992]	616	rg = qg/(1.-qnk(i))
	617	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	618	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
	619	tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
[879]	620
[1992]	621	END DO
[879]	622
[1992]	623	RETURN
	624	END SUBROUTINE cv30_undilute1
[879]	625
[1992]	626	SUBROUTINE cv30_trigger(len, nd, icb, plcl, p, th, tv, tvp, pbase, buoybase, &
	627	iflag, sig, w0)
	628	IMPLICIT NONE
[879]	629
[1992]	630	! -------------------------------------------------------------------
	631	! --- TRIGGERING
[879]	632
[1992]	633	! - computes the cloud base
	634	! - triggering (crude in this version)
	635	! - relaxation of sig and w0 when no convection
[879]	636
[1992]	637	! Caution1: if no convection, we set iflag=4
	638	! (it used to be 0 in convect3)
[879]	639
[1992]	640	! Caution2: at this stage, tvp (and thus buoy) are know up
	641	! through icb only!
	642	! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
	643	! -------------------------------------------------------------------
[879]	644
[1992]	645	include "cv30param.h"
[879]	646
[1992]	647	! input:
	648	INTEGER len, nd
	649	INTEGER icb(len)
	650	REAL plcl(len), p(len, nd)
	651	REAL th(len, nd), tv(len, nd), tvp(len, nd)
[879]	652
[1992]	653	! output:
	654	REAL pbase(len), buoybase(len)
[879]	655
[1992]	656	! input AND output:
	657	INTEGER iflag(len)
	658	REAL sig(len, nd), w0(len, nd)
[879]	659
[1992]	660	! local variables:
	661	INTEGER i, k
	662	REAL tvpbase, tvbase, tdif, ath, ath1
[879]	663
	664
[1992]	665	! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
[879]	666
[1992]	667	DO i = 1, len
	668	pbase(i) = plcl(i) + dpbase
	669	tvpbase = tvp(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
	670	(p(i,icb(i))-p(i,icb(i)+1)) + tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/( &
	671	p(i,icb(i))-p(i,icb(i)+1))
	672	tvbase = tv(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
	673	(p(i,icb(i))-p(i,icb(i)+1)) + tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/(p &
	674	(i,icb(i))-p(i,icb(i)+1))
	675	buoybase(i) = tvpbase - tvbase
	676	END DO
[879]	677
	678
[1992]	679	! * make sure that column is dry adiabatic between the surface *
	680	! * and cloud base, and that lifted air is positively buoyant *
	681	! * at cloud base *
	682	! * if not, return to calling program after resetting *
	683	! * sig(i) and w0(i) *
[879]	684
	685
[1992]	686	! oct3 do 200 i=1,len
	687	! oct3
	688	! oct3 tdif = buoybase(i)
	689	! oct3 ath1 = th(i,1)
	690	! oct3 ath = th(i,icb(i)-1) - dttrig
	691	! oct3
	692	! oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
	693	! oct3 do 60 k=1,nl
	694	! oct3 sig(i,k) = betasig(i,k) - 2.alphatdiftdif
	695	! oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
	696	! oct3 w0(i,k) = beta*w0(i,k)
	697	! oct3 60 continue
	698	! oct3 iflag(i)=4 ! pour version vectorisee
	699	! oct3c convect3 iflag(i)=0
	700	! oct3cccc return
	701	! oct3 endif
	702	! oct3
	703	! oct3200 continue
[879]	704
[1992]	705	! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
[879]	706
[1992]	707	DO k = 1, nl
	708	DO i = 1, len
[879]	709
[1992]	710	tdif = buoybase(i)
	711	ath1 = th(i, 1)
	712	ath = th(i, icb(i)-1) - dttrig
[879]	713
[1992]	714	IF (tdif<dtcrit .OR. ath>ath1) THEN
	715	sig(i, k) = betasig(i, k) - 2.alphatdiftdif
	716	sig(i, k) = amax1(sig(i,k), 0.0)
	717	w0(i, k) = beta*w0(i, k)
	718	iflag(i) = 4 ! pour version vectorisee
	719	! convect3 iflag(i)=0
	720	END IF
[879]	721
[1992]	722	END DO
	723	END DO
[879]	724
[1992]	725	! fin oct3 --
[879]	726
[1992]	727	RETURN
	728	END SUBROUTINE cv30_trigger
[879]	729
[1992]	730	SUBROUTINE cv30_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, &
	731	plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, &
	732	th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, &
	733	iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, &
	734	v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0)
[2311]	735	USE print_control_mod, ONLY: lunout
[1992]	736	IMPLICIT NONE
[879]	737
[1992]	738	include "cv30param.h"
[879]	739
[1992]	740	! inputs:
	741	INTEGER len, ncum, nd, ntra, nloc
	742	INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
	743	REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
	744	REAL pbase1(len), buoybase1(len)
	745	REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
	746	REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
	747	REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
	748	REAL tvp1(len, nd), clw1(len, nd)
	749	REAL th1(len, nd)
	750	REAL sig1(len, nd), w01(len, nd)
	751	REAL tra1(len, nd, ntra)
[879]	752
[1992]	753	! outputs:
	754	! en fait, on a nloc=len pour l'instant (cf cv_driver)
	755	INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
	756	REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
	757	REAL pbase(nloc), buoybase(nloc)
	758	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
	759	REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
	760	REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
	761	REAL tvp(nloc, nd), clw(nloc, nd)
	762	REAL th(nloc, nd)
	763	REAL sig(nloc, nd), w0(nloc, nd)
	764	REAL tra(nloc, nd, ntra)
[879]	765
[1992]	766	! local variables:
	767	INTEGER i, k, nn, j
[879]	768
[1992]	769	CHARACTER (LEN=20) :: modname = 'cv30_compress'
	770	CHARACTER (LEN=80) :: abort_message
[879]	771
	772
[1992]	773	DO k = 1, nl + 1
	774	nn = 0
	775	DO i = 1, len
	776	IF (iflag1(i)==0) THEN
	777	nn = nn + 1
	778	sig(nn, k) = sig1(i, k)
	779	w0(nn, k) = w01(i, k)
	780	t(nn, k) = t1(i, k)
	781	q(nn, k) = q1(i, k)
	782	qs(nn, k) = qs1(i, k)
	783	u(nn, k) = u1(i, k)
	784	v(nn, k) = v1(i, k)
	785	gz(nn, k) = gz1(i, k)
	786	h(nn, k) = h1(i, k)
	787	lv(nn, k) = lv1(i, k)
	788	cpn(nn, k) = cpn1(i, k)
	789	p(nn, k) = p1(i, k)
	790	ph(nn, k) = ph1(i, k)
	791	tv(nn, k) = tv1(i, k)
	792	tp(nn, k) = tp1(i, k)
	793	tvp(nn, k) = tvp1(i, k)
	794	clw(nn, k) = clw1(i, k)
	795	th(nn, k) = th1(i, k)
	796	END IF
	797	END DO
	798	END DO
[879]	799
[1992]	800	! do 121 j=1,ntra
	801	! do 111 k=1,nd
	802	! nn=0
	803	! do 101 i=1,len
	804	! if(iflag1(i).eq.0)then
	805	! nn=nn+1
	806	! tra(nn,k,j)=tra1(i,k,j)
	807	! endif
	808	! 101 continue
	809	! 111 continue
	810	! 121 continue
[879]	811
[1992]	812	IF (nn/=ncum) THEN
	813	WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
	814	abort_message = ''
[2311]	815	CALL abort_physic(modname, abort_message, 1)
[1992]	816	END IF
[879]	817
[1992]	818	nn = 0
	819	DO i = 1, len
	820	IF (iflag1(i)==0) THEN
	821	nn = nn + 1
	822	pbase(nn) = pbase1(i)
	823	buoybase(nn) = buoybase1(i)
	824	plcl(nn) = plcl1(i)
	825	tnk(nn) = tnk1(i)
	826	qnk(nn) = qnk1(i)
	827	gznk(nn) = gznk1(i)
	828	nk(nn) = nk1(i)
	829	icb(nn) = icb1(i)
	830	icbs(nn) = icbs1(i)
	831	iflag(nn) = iflag1(i)
	832	END IF
	833	END DO
[879]	834
[1992]	835	RETURN
	836	END SUBROUTINE cv30_compress
[879]	837
[1992]	838	SUBROUTINE cv30_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
	839	q, qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
	840	ep, sigp, buoy)
	841	IMPLICIT NONE
[879]	842
[1992]	843	! ---------------------------------------------------------------------
	844	! Purpose:
	845	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	846	! &
	847	! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
	848	! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
	849	! &
	850	! FIND THE LEVEL OF NEUTRAL BUOYANCY
[879]	851
[1992]	852	! Main differences convect3/convect4:
	853	! - icbs (input) is the first level above LCL (may differ from icb)
	854	! - many minor differences in the iterations
	855	! - condensed water not removed from tvp in convect3
	856	! - vertical profile of buoyancy computed here (use of buoybase)
	857	! - the determination of inb is different
	858	! - no inb1, only inb in output
	859	! ---------------------------------------------------------------------
[879]	860
[1992]	861	include "cvthermo.h"
	862	include "cv30param.h"
	863	include "conema3.h"
[879]	864
[1992]	865	! inputs:
	866	INTEGER ncum, nd, nloc
	867	INTEGER icb(nloc), icbs(nloc), nk(nloc)
	868	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
	869	REAL p(nloc, nd)
	870	REAL tnk(nloc), qnk(nloc), gznk(nloc)
	871	REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
	872	REAL pbase(nloc), buoybase(nloc), plcl(nloc)
[879]	873
[1992]	874	! outputs:
	875	INTEGER inb(nloc)
	876	REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
	877	REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
	878	REAL buoy(nloc, nd)
[879]	879
[1992]	880	! local variables:
	881	INTEGER i, k
	882	REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
	883	REAL by, defrac, pden
	884	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
	885	LOGICAL lcape(nloc)
[879]	886
[1992]	887	! =====================================================================
	888	! --- SOME INITIALIZATIONS
	889	! =====================================================================
[879]	890
[1992]	891	DO k = 1, nl
	892	DO i = 1, ncum
	893	ep(i, k) = 0.0
	894	sigp(i, k) = spfac
	895	END DO
	896	END DO
[879]	897
[1992]	898	! =====================================================================
	899	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	900	! =====================================================================
[879]	901
[1992]	902	! --- The procedure is to solve the equation.
	903	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
[879]	904
[1992]	905	! * Calculate certain parcel quantities, including static energy *
[879]	906
	907
[1992]	908	DO i = 1, ncum
	909	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))*tnk(i) & ! debug &
	910	! +qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
	911	+qnk(i)(lv0-clmcpv(tnk(i)-273.15)) + gznk(i)
	912	END DO
[879]	913
	914
[1992]	915	! * Find lifted parcel quantities above cloud base *
[879]	916
	917
[1992]	918	DO k = minorig + 1, nl
	919	DO i = 1, ncum
	920	! ori if(k.ge.(icb(i)+1))then
	921	IF (k>=(icbs(i)+1)) THEN ! convect3
	922	tg = t(i, k)
	923	qg = qs(i, k)
	924	! debug alv=lv0-clmcpv*(t(i,k)-t0)
	925	alv = lv0 - clmcpv*(t(i,k)-273.15)
[879]	926
[1992]	927	! First iteration.
[879]	928
[1992]	929	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	930	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
	931	+alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
	932	s = 1./s
	933	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	934	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
	935	tg = tg + s*(ah0(i)-ahg)
	936	! ori tg=max(tg,35.0)
	937	! debug tc=tg-t0
	938	tc = tg - 273.15
	939	denom = 243.5 + tc
	940	denom = max(denom, 1.0) ! convect3
	941	! ori if(tc.ge.0.0)then
	942	es = 6.112exp(17.67tc/denom)
	943	! ori else
	944	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	945	! ori endif
	946	qg = epses/(p(i,k)-es(1.-eps))
[879]	947
[1992]	948	! Second iteration.
[879]	949
[1992]	950	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	951	! ori s=1./s
	952	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	953	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
	954	tg = tg + s*(ah0(i)-ahg)
	955	! ori tg=max(tg,35.0)
	956	! debug tc=tg-t0
	957	tc = tg - 273.15
	958	denom = 243.5 + tc
	959	denom = max(denom, 1.0) ! convect3
	960	! ori if(tc.ge.0.0)then
	961	es = 6.112exp(17.67tc/denom)
	962	! ori else
	963	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	964	! ori endif
	965	qg = epses/(p(i,k)-es(1.-eps))
[879]	966
[1992]	967	! debug alv=lv0-clmcpv*(t(i,k)-t0)
	968	alv = lv0 - clmcpv*(t(i,k)-273.15)
	969	! print*,'cpd dans convect2 ',cpd
	970	! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
	971	! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
[879]	972
[1992]	973	! ori c approximation here:
	974	! ori
	975	! tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
[879]	976
[1992]	977	! convect3: no approximation:
	978	tp(i, k) = (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
[879]	979
[1992]	980	clw(i, k) = qnk(i) - qg
	981	clw(i, k) = max(0.0, clw(i,k))
	982	rg = qg/(1.-qnk(i))
	983	! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
	984	! convect3: (qg utilise au lieu du vrai mixing ratio rg):
	985	tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
	986	END IF
	987	END DO
	988	END DO
[879]	989
[1992]	990	! =====================================================================
	991	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
	992	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
	993	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
	994	! =====================================================================
[879]	995
[1992]	996	! ori do 320 k=minorig+1,nl
	997	DO k = 1, nl ! convect3
	998	DO i = 1, ncum
	999	pden = ptcrit - pbcrit
	1000	ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
	1001	ep(i, k) = amax1(ep(i,k), 0.0)
	1002	ep(i, k) = amin1(ep(i,k), epmax)
	1003	sigp(i, k) = spfac
	1004	! ori if(k.ge.(nk(i)+1))then
	1005	! ori tca=tp(i,k)-t0
	1006	! ori if(tca.ge.0.0)then
	1007	! ori elacrit=elcrit
	1008	! ori else
	1009	! ori elacrit=elcrit*(1.0-tca/tlcrit)
	1010	! ori endif
	1011	! ori elacrit=max(elacrit,0.0)
	1012	! ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
	1013	! ori ep(i,k)=max(ep(i,k),0.0 )
	1014	! ori ep(i,k)=min(ep(i,k),1.0 )
	1015	! ori sigp(i,k)=sigs
	1016	! ori endif
	1017	END DO
	1018	END DO
[879]	1019
[1992]	1020	! =====================================================================
	1021	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
	1022	! --- VIRTUAL TEMPERATURE
	1023	! =====================================================================
[879]	1024
[1992]	1025	! dans convect3, tvp est calcule en une seule fois, et sans retirer
	1026	! l'eau condensee (~> reversible CAPE)
[879]	1027
[1992]	1028	! ori do 340 k=minorig+1,nl
	1029	! ori do 330 i=1,ncum
	1030	! ori if(k.ge.(icb(i)+1))then
	1031	! ori tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
	1032	! oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
	1033	! oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
	1034	! ori endif
	1035	! ori 330 continue
	1036	! ori 340 continue
[879]	1037
[1992]	1038	! ori do 350 i=1,ncum
	1039	! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
	1040	! ori 350 continue
[879]	1041
[1992]	1042	DO i = 1, ncum ! convect3
	1043	tp(i, nlp) = tp(i, nl) ! convect3
	1044	END DO ! convect3
[879]	1045
[1992]	1046	! =====================================================================
	1047	! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
	1048	! =====================================================================
[879]	1049
[1992]	1050	! -- this is for convect3 only:
[879]	1051
[1992]	1052	! first estimate of buoyancy:
[879]	1053
[1992]	1054	DO i = 1, ncum
	1055	DO k = 1, nl
	1056	buoy(i, k) = tvp(i, k) - tv(i, k)
	1057	END DO
	1058	END DO
[879]	1059
[1992]	1060	! set buoyancy=buoybase for all levels below base
	1061	! for safety, set buoy(icb)=buoybase
[879]	1062
[1992]	1063	DO i = 1, ncum
	1064	DO k = 1, nl
	1065	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
	1066	buoy(i, k) = buoybase(i)
	1067	END IF
	1068	END DO
	1069	! IM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i)
	1070	buoy(i, icb(i)) = buoybase(i)
	1071	END DO
[879]	1072
[1992]	1073	! -- end convect3
[879]	1074
[1992]	1075	! =====================================================================
	1076	! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
	1077	! --- LEVEL OF NEUTRAL BUOYANCY
	1078	! =====================================================================
[879]	1079
[1992]	1080	! -- this is for convect3 only:
[879]	1081
[1992]	1082	DO i = 1, ncum
	1083	inb(i) = nl - 1
	1084	END DO
[879]	1085
[1992]	1086	DO i = 1, ncum
	1087	DO k = 1, nl - 1
	1088	IF ((k>=icb(i)) .AND. (buoy(i,k)<dtovsh)) THEN
	1089	inb(i) = min(inb(i), k)
	1090	END IF
	1091	END DO
	1092	END DO
[879]	1093
[1992]	1094	! -- end convect3
[879]	1095
[1992]	1096	! ori do 510 i=1,ncum
	1097	! ori cape(i)=0.0
	1098	! ori capem(i)=0.0
	1099	! ori inb(i)=icb(i)+1
	1100	! ori inb1(i)=inb(i)
	1101	! ori 510 continue
[879]	1102
[1992]	1103	! Originial Code
[879]	1104
[1992]	1105	! do 530 k=minorig+1,nl-1
	1106	! do 520 i=1,ncum
	1107	! if(k.ge.(icb(i)+1))then
	1108	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1109	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1110	! cape(i)=cape(i)+by
	1111	! if(by.ge.0.0)inb1(i)=k+1
	1112	! if(cape(i).gt.0.0)then
	1113	! inb(i)=k+1
	1114	! capem(i)=cape(i)
	1115	! endif
	1116	! endif
	1117	! 520 continue
	1118	! 530 continue
	1119	! do 540 i=1,ncum
	1120	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
	1121	! cape(i)=capem(i)+byp
	1122	! defrac=capem(i)-cape(i)
	1123	! defrac=max(defrac,0.001)
	1124	! frac(i)=-cape(i)/defrac
	1125	! frac(i)=min(frac(i),1.0)
	1126	! frac(i)=max(frac(i),0.0)
	1127	! 540 continue
[879]	1128
[1992]	1129	! K Emanuel fix
[879]	1130
[1992]	1131	! call zilch(byp,ncum)
	1132	! do 530 k=minorig+1,nl-1
	1133	! do 520 i=1,ncum
	1134	! if(k.ge.(icb(i)+1))then
	1135	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1136	! cape(i)=cape(i)+by
	1137	! if(by.ge.0.0)inb1(i)=k+1
	1138	! if(cape(i).gt.0.0)then
	1139	! inb(i)=k+1
	1140	! capem(i)=cape(i)
	1141	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1142	! endif
	1143	! endif
	1144	! 520 continue
	1145	! 530 continue
	1146	! do 540 i=1,ncum
	1147	! inb(i)=max(inb(i),inb1(i))
	1148	! cape(i)=capem(i)+byp(i)
	1149	! defrac=capem(i)-cape(i)
	1150	! defrac=max(defrac,0.001)
	1151	! frac(i)=-cape(i)/defrac
	1152	! frac(i)=min(frac(i),1.0)
	1153	! frac(i)=max(frac(i),0.0)
	1154	! 540 continue
[879]	1155
[1992]	1156	! J Teixeira fix
[879]	1157
[1992]	1158	! ori call zilch(byp,ncum)
	1159	! ori do 515 i=1,ncum
	1160	! ori lcape(i)=.true.
	1161	! ori 515 continue
	1162	! ori do 530 k=minorig+1,nl-1
	1163	! ori do 520 i=1,ncum
	1164	! ori if(cape(i).lt.0.0)lcape(i)=.false.
	1165	! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
	1166	! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1167	! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1168	! ori cape(i)=cape(i)+by
	1169	! ori if(by.ge.0.0)inb1(i)=k+1
	1170	! ori if(cape(i).gt.0.0)then
	1171	! ori inb(i)=k+1
	1172	! ori capem(i)=cape(i)
	1173	! ori endif
	1174	! ori endif
	1175	! ori 520 continue
	1176	! ori 530 continue
	1177	! ori do 540 i=1,ncum
	1178	! ori cape(i)=capem(i)+byp(i)
	1179	! ori defrac=capem(i)-cape(i)
	1180	! ori defrac=max(defrac,0.001)
	1181	! ori frac(i)=-cape(i)/defrac
	1182	! ori frac(i)=min(frac(i),1.0)
	1183	! ori frac(i)=max(frac(i),0.0)
	1184	! ori 540 continue
[879]	1185
[1992]	1186	! =====================================================================
	1187	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
	1188	! =====================================================================
[879]	1189
[1992]	1190	! ym do i=1,ncum*nlp
	1191	! ym hp(i,1)=h(i,1)
	1192	! ym enddo
[879]	1193
[1992]	1194	DO k = 1, nlp
	1195	DO i = 1, ncum
	1196	hp(i, k) = h(i, k)
	1197	END DO
	1198	END DO
[879]	1199
[1992]	1200	DO k = minorig + 1, nl
	1201	DO i = 1, ncum
	1202	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
	1203	hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k &
	1204	)
	1205	END IF
	1206	END DO
	1207	END DO
[879]	1208
[1992]	1209	RETURN
	1210	END SUBROUTINE cv30_undilute2
[879]	1211
[1992]	1212	SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
	1213	sig, w0, cape, m)
	1214	IMPLICIT NONE
[879]	1215
[1992]	1216	! ===================================================================
	1217	! --- CLOSURE OF CONVECT3
[879]	1218
[1992]	1219	! vectorization: S. Bony
	1220	! ===================================================================
[879]	1221
[1992]	1222	include "cvthermo.h"
	1223	include "cv30param.h"
[879]	1224
[1992]	1225	! input:
	1226	INTEGER ncum, nd, nloc
	1227	INTEGER icb(nloc), inb(nloc)
	1228	REAL pbase(nloc)
	1229	REAL p(nloc, nd), ph(nloc, nd+1)
	1230	REAL tv(nloc, nd), buoy(nloc, nd)
[879]	1231
[1992]	1232	! input/output:
	1233	REAL sig(nloc, nd), w0(nloc, nd)
[879]	1234
[1992]	1235	! output:
	1236	REAL cape(nloc)
	1237	REAL m(nloc, nd)
[879]	1238
[1992]	1239	! local variables:
	1240	INTEGER i, j, k, icbmax
	1241	REAL deltap, fac, w, amu
	1242	REAL dtmin(nloc, nd), sigold(nloc, nd)
[879]	1243
	1244
[1992]	1245	! -------------------------------------------------------
	1246	! -- Initialization
	1247	! -------------------------------------------------------
[879]	1248
[1992]	1249	DO k = 1, nl
	1250	DO i = 1, ncum
	1251	m(i, k) = 0.0
	1252	END DO
	1253	END DO
[879]	1254
[1992]	1255	! -------------------------------------------------------
	1256	! -- Reset sig(i) and w0(i) for i>inb and i<icb
	1257	! -------------------------------------------------------
[879]	1258
[1992]	1259	! update sig and w0 above LNB:
[879]	1260
[1992]	1261	DO k = 1, nl - 1
	1262	DO i = 1, ncum
	1263	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
	1264	sig(i, k) = betasig(i, k) + 2.alphabuoy(i, inb(i))abs(buoy(i,inb( &
	1265	i)))
	1266	sig(i, k) = amax1(sig(i,k), 0.0)
	1267	w0(i, k) = beta*w0(i, k)
	1268	END IF
	1269	END DO
	1270	END DO
[879]	1271
[1992]	1272	! compute icbmax:
[879]	1273
[1992]	1274	icbmax = 2
	1275	DO i = 1, ncum
	1276	icbmax = max(icbmax, icb(i))
	1277	END DO
[879]	1278
[1992]	1279	! update sig and w0 below cloud base:
[879]	1280
[1992]	1281	DO k = 1, icbmax
	1282	DO i = 1, ncum
	1283	IF (k<=icb(i)) THEN
	1284	sig(i, k) = betasig(i, k) - 2.alphabuoy(i, icb(i))buoy(i, icb(i))
	1285	sig(i, k) = amax1(sig(i,k), 0.0)
	1286	w0(i, k) = beta*w0(i, k)
	1287	END IF
	1288	END DO
	1289	END DO
[879]	1290
[1992]	1291	! ! if(inb.lt.(nl-1))then
	1292	! ! do 85 i=inb+1,nl-1
	1293	! ! sig(i)=betasig(i)+2.alphabuoy(inb)
	1294	! ! 1 abs(buoy(inb))
	1295	! ! sig(i)=amax1(sig(i),0.0)
	1296	! ! w0(i)=beta*w0(i)
	1297	! ! 85 continue
	1298	! ! end if
[879]	1299
[1992]	1300	! ! do 87 i=1,icb
	1301	! ! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
	1302	! ! sig(i)=amax1(sig(i),0.0)
	1303	! ! w0(i)=beta*w0(i)
	1304	! ! 87 continue
[1742]	1305
[1992]	1306	! -------------------------------------------------------------
	1307	! -- Reset fractional areas of updrafts and w0 at initial time
	1308	! -- and after 10 time steps of no convection
	1309	! -------------------------------------------------------------
[1742]	1310
[1992]	1311	DO k = 1, nl - 1
	1312	DO i = 1, ncum
	1313	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
	1314	sig(i, k) = 0.0
	1315	w0(i, k) = 0.0
	1316	END IF
	1317	END DO
	1318	END DO
[879]	1319
[1992]	1320	! -------------------------------------------------------------
	1321	! -- Calculate convective available potential energy (cape),
	1322	! -- vertical velocity (w), fractional area covered by
	1323	! -- undilute updraft (sig), and updraft mass flux (m)
	1324	! -------------------------------------------------------------
[879]	1325
[1992]	1326	DO i = 1, ncum
	1327	cape(i) = 0.0
	1328	END DO
[879]	1329
[1992]	1330	! compute dtmin (minimum buoyancy between ICB and given level k):
[879]	1331
[1992]	1332	DO i = 1, ncum
	1333	DO k = 1, nl
	1334	dtmin(i, k) = 100.0
	1335	END DO
	1336	END DO
[879]	1337
[1992]	1338	DO i = 1, ncum
	1339	DO k = 1, nl
	1340	DO j = minorig, nl
	1341	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- &
	1342	1))) THEN
	1343	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
	1344	END IF
	1345	END DO
	1346	END DO
	1347	END DO
[879]	1348
[1992]	1349	! the interval on which cape is computed starts at pbase :
[879]	1350
[1992]	1351	DO k = 1, nl
	1352	DO i = 1, ncum
[879]	1353
[1992]	1354	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
[879]	1355
[1992]	1356	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
	1357	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
	1358	cape(i) = amax1(0.0, cape(i))
	1359	sigold(i, k) = sig(i, k)
[879]	1360
[1992]	1361	! dtmin(i,k)=100.0
	1362	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
	1363	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
	1364	! 97 continue
[879]	1365
[1992]	1366	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
	1367	sig(i, k) = amax1(sig(i,k), 0.0)
	1368	sig(i, k) = amin1(sig(i,k), 0.01)
	1369	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
	1370	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
	1371	amu = 0.5(sig(i,k)+sigold(i,k))w
	1372	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
	1373	w0(i, k) = w
	1374	END IF
[879]	1375
[1992]	1376	END DO
	1377	END DO
[879]	1378
[1992]	1379	DO i = 1, ncum
	1380	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
	1381	m(i, icb(i)) = 0.5m(i, icb(i)+1)(ph(i,icb(i))-ph(i,icb(i)+1))/ &
	1382	(ph(i,icb(i)+1)-ph(i,icb(i)+2))
	1383	sig(i, icb(i)) = sig(i, icb(i)+1)
	1384	sig(i, icb(i)-1) = sig(i, icb(i))
	1385	END DO
[879]	1386
	1387
[1992]	1388	! ! cape=0.0
	1389	! ! do 98 i=icb+1,inb
	1390	! ! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
	1391	! ! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
	1392	! ! dcape=rrdbuoy(i-1)deltap/p(i-1)
	1393	! ! dlnp=deltap/p(i-1)
	1394	! ! cape=amax1(0.0,cape)
	1395	! ! sigold=sig(i)
[879]	1396
[1992]	1397	! ! dtmin=100.0
	1398	! ! do 97 j=icb,i-1
	1399	! ! dtmin=amin1(dtmin,buoy(j))
	1400	! ! 97 continue
[879]	1401
[1992]	1402	! ! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
	1403	! ! sig(i)=amax1(sig(i),0.0)
	1404	! ! sig(i)=amin1(sig(i),0.01)
	1405	! ! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
	1406	! ! w=(1.-beta)facsqrt(cape)+beta*w0(i)
	1407	! ! amu=0.5(sig(i)+sigold)w
	1408	! ! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
	1409	! ! w0(i)=w
	1410	! ! 98 continue
	1411	! ! w0(icb)=0.5*w0(icb+1)
	1412	! ! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
	1413	! ! sig(icb)=sig(icb+1)
	1414	! ! sig(icb-1)=sig(icb)
[879]	1415
[1992]	1416	RETURN
	1417	END SUBROUTINE cv30_closure
[879]	1418
[1992]	1419	SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
	1420	u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, &
	1421	vent, sij, elij, ments, qents, traent)
	1422	IMPLICIT NONE
[879]	1423
[1992]	1424	! ---------------------------------------------------------------------
	1425	! a faire:
	1426	! - changer rr(il,1) -> qnk(il)
	1427	! - vectorisation de la partie normalisation des flux (do 789...)
	1428	! ---------------------------------------------------------------------
[879]	1429
[1992]	1430	include "cvthermo.h"
	1431	include "cv30param.h"
[879]	1432
[1992]	1433	! inputs:
	1434	INTEGER ncum, nd, na, ntra, nloc
	1435	INTEGER icb(nloc), inb(nloc), nk(nloc)
	1436	REAL sig(nloc, nd)
	1437	REAL qnk(nloc)
	1438	REAL ph(nloc, nd+1)
	1439	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
	1440	REAL u(nloc, nd), v(nloc, nd)
	1441	REAL tra(nloc, nd, ntra) ! input of convect3
	1442	REAL lv(nloc, na), h(nloc, na), hp(nloc, na)
	1443	REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
	1444	REAL m(nloc, na) ! input of convect3
[879]	1445
[1992]	1446	! outputs:
	1447	REAL ment(nloc, na, na), qent(nloc, na, na)
	1448	REAL uent(nloc, na, na), vent(nloc, na, na)
	1449	REAL sij(nloc, na, na), elij(nloc, na, na)
	1450	REAL traent(nloc, nd, nd, ntra)
	1451	REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
	1452	REAL sigij(nloc, nd, nd)
[879]	1453
[1992]	1454	! local variables:
	1455	INTEGER i, j, k, il, im, jm
	1456	INTEGER num1, num2
	1457	INTEGER nent(nloc, na)
	1458	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
	1459	REAL alt, smid, sjmin, sjmax, delp, delm
	1460	REAL asij(nloc), smax(nloc), scrit(nloc)
	1461	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
	1462	REAL wgh
	1463	REAL zm(nloc, na)
	1464	LOGICAL lwork(nloc)
[879]	1465
[1992]	1466	! =====================================================================
	1467	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	1468	! =====================================================================
[879]	1469
[1992]	1470	! ori do 360 i=1,ncum*nlp
	1471	DO j = 1, nl
	1472	DO i = 1, ncum
	1473	nent(i, j) = 0
	1474	! in convect3, m is computed in cv3_closure
	1475	! ori m(i,1)=0.0
	1476	END DO
	1477	END DO
[879]	1478
[1992]	1479	! ori do 400 k=1,nlp
	1480	! ori do 390 j=1,nlp
	1481	DO j = 1, nl
	1482	DO k = 1, nl
	1483	DO i = 1, ncum
	1484	qent(i, k, j) = rr(i, j)
	1485	uent(i, k, j) = u(i, j)
	1486	vent(i, k, j) = v(i, j)
	1487	elij(i, k, j) = 0.0
	1488	! ym ment(i,k,j)=0.0
	1489	! ym sij(i,k,j)=0.0
	1490	END DO
	1491	END DO
	1492	END DO
[879]	1493
[1992]	1494	! ym
	1495	ment(1:ncum, 1:nd, 1:nd) = 0.0
	1496	sij(1:ncum, 1:nd, 1:nd) = 0.0
[879]	1497
[1992]	1498	! do k=1,ntra
	1499	! do j=1,nd ! instead nlp
	1500	! do i=1,nd ! instead nlp
	1501	! do il=1,ncum
	1502	! traent(il,i,j,k)=tra(il,j,k)
	1503	! enddo
	1504	! enddo
	1505	! enddo
	1506	! enddo
	1507	zm(:, :) = 0.
[879]	1508
[1992]	1509	! =====================================================================
	1510	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	1511	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	1512	! --- FRACTION (sij)
	1513	! =====================================================================
[879]	1514
[1992]	1515	DO i = minorig + 1, nl
[879]	1516
[1992]	1517	DO j = minorig, nl
	1518	DO il = 1, ncum
	1519	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
	1520	1)) .AND. (j<=inb(il))) THEN
[879]	1521
[1992]	1522	rti = rr(il, 1) - ep(il, i)*clw(il, i)
	1523	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
	1524	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
	1525	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
	1526	dei = denom
	1527	IF (abs(dei)<0.01) dei = 0.01
	1528	sij(il, i, j) = anum/dei
	1529	sij(il, i, i) = 1.0
	1530	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
	1531	altem = altem/bf2
	1532	cwat = clw(il, j)*(1.-ep(il,j))
	1533	stemp = sij(il, i, j)
	1534	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
	1535	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
	1536	denom = denom + lv(il, j)*(rr(il,i)-rti)
	1537	IF (abs(denom)<0.01) denom = 0.01
	1538	sij(il, i, j) = anum/denom
	1539	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - &
	1540	rs(il, j)
	1541	altem = altem - (bf2-1.)*cwat
	1542	END IF
	1543	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
	1544	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
	1545	uent(il, i, j) = sij(il, i, j)*u(il, i) + &
	1546	(1.-sij(il,i,j))*u(il, nk(il))
	1547	vent(il, i, j) = sij(il, i, j)*v(il, i) + &
	1548	(1.-sij(il,i,j))*v(il, nk(il))
	1549	! !!! do k=1,ntra
	1550	! !!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1551	! !!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1552	! !!! end do
	1553	elij(il, i, j) = altem
	1554	elij(il, i, j) = amax1(0.0, elij(il,i,j))
	1555	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
	1556	nent(il, i) = nent(il, i) + 1
	1557	END IF
	1558	sij(il, i, j) = amax1(0.0, sij(il,i,j))
	1559	sij(il, i, j) = amin1(1.0, sij(il,i,j))
	1560	END IF ! new
	1561	END DO
	1562	END DO
[879]	1563
[1992]	1564	! do k=1,ntra
	1565	! do j=minorig,nl
	1566	! do il=1,ncum
	1567	! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
	1568	! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
	1569	! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1570	! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1571	! endif
	1572	! enddo
	1573	! enddo
	1574	! enddo
[879]	1575
	1576
[1992]	1577	! *** if no air can entrain at level i assume that updraft detrains
	1578	! ***
	1579	! *** at that level and calculate detrained air flux and properties
	1580	! ***
[879]	1581
	1582
[1992]	1583	! @ do 170 i=icb(il),inb(il)
[879]	1584
[1992]	1585	DO il = 1, ncum
	1586	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
	1587	! @ if(nent(il,i).eq.0)then
	1588	ment(il, i, i) = m(il, i)
	1589	qent(il, i, i) = rr(il, nk(il)) - ep(il, i)*clw(il, i)
	1590	uent(il, i, i) = u(il, nk(il))
	1591	vent(il, i, i) = v(il, nk(il))
	1592	elij(il, i, i) = clw(il, i)
	1593	! MAF sij(il,i,i)=1.0
	1594	sij(il, i, i) = 0.0
	1595	END IF
	1596	END DO
	1597	END DO
[879]	1598
[1992]	1599	! do j=1,ntra
	1600	! do i=minorig+1,nl
	1601	! do il=1,ncum
	1602	! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
	1603	! traent(il,i,i,j)=tra(il,nk(il),j)
	1604	! endif
	1605	! enddo
	1606	! enddo
	1607	! enddo
[879]	1608
[1992]	1609	DO j = minorig, nl
	1610	DO i = minorig, nl
	1611	DO il = 1, ncum
	1612	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
	1613	inb(il))) THEN
	1614	sigij(il, i, j) = sij(il, i, j)
	1615	END IF
	1616	END DO
	1617	END DO
	1618	END DO
	1619	! @ enddo
[879]	1620
[1992]	1621	! @170 continue
[879]	1622
[1992]	1623	! =====================================================================
	1624	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	1625	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	1626	! =====================================================================
[879]	1627
[1992]	1628	! ym call zilch(asum,ncum*nd)
	1629	! ym call zilch(bsum,ncum*nd)
	1630	! ym call zilch(csum,ncum*nd)
	1631	CALL zilch(asum, nloc*nd)
	1632	CALL zilch(csum, nloc*nd)
	1633	CALL zilch(csum, nloc*nd)
[879]	1634
[1992]	1635	DO il = 1, ncum
	1636	lwork(il) = .FALSE.
	1637	END DO
[879]	1638
[1992]	1639	DO i = minorig + 1, nl
[879]	1640
[1992]	1641	num1 = 0
	1642	DO il = 1, ncum
	1643	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
	1644	END DO
	1645	IF (num1<=0) GO TO 789
[879]	1646
	1647
[1992]	1648	DO il = 1, ncum
	1649	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	1650	lwork(il) = (nent(il,i)/=0)
	1651	qp = rr(il, 1) - ep(il, i)*clw(il, i)
	1652	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
	1653	(cpv-cpd)t(il, i)(qp-rr(il,i))
	1654	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
	1655	(cpd-cpv)t(il, i)(rr(il,i)-qp)
	1656	IF (abs(denom)<0.01) denom = 0.01
	1657	scrit(il) = anum/denom
	1658	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
	1659	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
	1660	smax(il) = 0.0
	1661	asij(il) = 0.0
	1662	END IF
	1663	END DO
[879]	1664
[1992]	1665	DO j = nl, minorig, -1
[879]	1666
[1992]	1667	num2 = 0
	1668	DO il = 1, ncum
	1669	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	1670	il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
	1671	END DO
	1672	IF (num2<=0) GO TO 175
[879]	1673
[1992]	1674	DO il = 1, ncum
	1675	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	1676	il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
[879]	1677
[1992]	1678	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
	1679	wgh = 1.0
	1680	IF (j>i) THEN
	1681	sjmax = amax1(sij(il,i,j+1), smax(il))
	1682	sjmax = amin1(sjmax, scrit(il))
	1683	smax(il) = amax1(sij(il,i,j), smax(il))
	1684	sjmin = amax1(sij(il,i,j-1), smax(il))
	1685	sjmin = amin1(sjmin, scrit(il))
	1686	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
	1687	smid = amin1(sij(il,i,j), scrit(il))
	1688	ELSE
	1689	sjmax = amax1(sij(il,i,j+1), scrit(il))
	1690	smid = amax1(sij(il,i,j), scrit(il))
	1691	sjmin = 0.0
	1692	IF (j>1) sjmin = sij(il, i, j-1)
	1693	sjmin = amax1(sjmin, scrit(il))
	1694	END IF
	1695	delp = abs(sjmax-smid)
	1696	delm = abs(sjmin-smid)
	1697	asij(il) = asij(il) + wgh*(delp+delm)
	1698	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
	1699	END IF
	1700	END IF
	1701	END DO
[879]	1702
[1992]	1703	175 END DO
[879]	1704
[1992]	1705	DO il = 1, ncum
	1706	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	1707	asij(il) = amax1(1.0E-16, asij(il))
	1708	asij(il) = 1.0/asij(il)
	1709	asum(il, i) = 0.0
	1710	bsum(il, i) = 0.0
	1711	csum(il, i) = 0.0
	1712	END IF
	1713	END DO
[879]	1714
[1992]	1715	DO j = minorig, nl
	1716	DO il = 1, ncum
	1717	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	1718	il)-1) .AND. j<=inb(il)) THEN
	1719	ment(il, i, j) = ment(il, i, j)*asij(il)
	1720	END IF
	1721	END DO
	1722	END DO
[879]	1723
[1992]	1724	DO j = minorig, nl
	1725	DO il = 1, ncum
	1726	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	1727	il)-1) .AND. j<=inb(il)) THEN
	1728	asum(il, i) = asum(il, i) + ment(il, i, j)
	1729	ment(il, i, j) = ment(il, i, j)*sig(il, j)
	1730	bsum(il, i) = bsum(il, i) + ment(il, i, j)
	1731	END IF
	1732	END DO
	1733	END DO
[879]	1734
[1992]	1735	DO il = 1, ncum
	1736	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	1737	bsum(il, i) = amax1(bsum(il,i), 1.0E-16)
	1738	bsum(il, i) = 1.0/bsum(il, i)
	1739	END IF
	1740	END DO
[879]	1741
[1992]	1742	DO j = minorig, nl
	1743	DO il = 1, ncum
	1744	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	1745	il)-1) .AND. j<=inb(il)) THEN
	1746	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
	1747	END IF
	1748	END DO
	1749	END DO
[879]	1750
[1992]	1751	DO j = minorig, nl
	1752	DO il = 1, ncum
	1753	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	1754	il)-1) .AND. j<=inb(il)) THEN
	1755	csum(il, i) = csum(il, i) + ment(il, i, j)
	1756	END IF
	1757	END DO
	1758	END DO
[879]	1759
[1992]	1760	DO il = 1, ncum
	1761	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	1762	csum(il,i)<m(il,i)) THEN
	1763	nent(il, i) = 0
	1764	ment(il, i, i) = m(il, i)
	1765	qent(il, i, i) = rr(il, 1) - ep(il, i)*clw(il, i)
	1766	uent(il, i, i) = u(il, nk(il))
	1767	vent(il, i, i) = v(il, nk(il))
	1768	elij(il, i, i) = clw(il, i)
	1769	! MAF sij(il,i,i)=1.0
	1770	sij(il, i, i) = 0.0
	1771	END IF
	1772	END DO ! il
[879]	1773
[1992]	1774	! do j=1,ntra
	1775	! do il=1,ncum
	1776	! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	1777	! : .and. csum(il,i).lt.m(il,i) ) then
	1778	! traent(il,i,i,j)=tra(il,nk(il),j)
	1779	! endif
	1780	! enddo
	1781	! enddo
	1782	789 END DO
[879]	1783
[1992]	1784	! MAF: renormalisation de MENT
	1785	DO jm = 1, nd
	1786	DO im = 1, nd
	1787	DO il = 1, ncum
	1788	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
	1789	END DO
	1790	END DO
	1791	END DO
[879]	1792
[1992]	1793	DO jm = 1, nd
	1794	DO im = 1, nd
	1795	DO il = 1, ncum
	1796	IF (zm(il,im)/=0.) THEN
	1797	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
	1798	END IF
	1799	END DO
	1800	END DO
	1801	END DO
[879]	1802
[1992]	1803	DO jm = 1, nd
	1804	DO im = 1, nd
	1805	DO il = 1, ncum
	1806	qents(il, im, jm) = qent(il, im, jm)
	1807	ments(il, im, jm) = ment(il, im, jm)
	1808	END DO
	1809	END DO
	1810	END DO
[879]	1811
[1992]	1812	RETURN
	1813	END SUBROUTINE cv30_mixing
[879]	1814
	1815
[1992]	1816	SUBROUTINE cv30_unsat(nloc, ncum, nd, na, ntra, icb, inb, t, rr, rs, gz, u, &
	1817	v, tra, p, ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, &
	1818	mp, rp, up, vp, trap, wt, water, evap, b & ! RomP-jyg
	1819	, wdtraina, wdtrainm) ! 26/08/10 RomP-jyg
	1820	IMPLICIT NONE
[879]	1821
	1822
[1992]	1823	include "cvthermo.h"
	1824	include "cv30param.h"
	1825	include "cvflag.h"
[879]	1826
[1992]	1827	! inputs:
	1828	INTEGER ncum, nd, na, ntra, nloc
	1829	INTEGER icb(nloc), inb(nloc)
	1830	REAL delt, plcl(nloc)
	1831	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
	1832	REAL u(nloc, nd), v(nloc, nd)
	1833	REAL tra(nloc, nd, ntra)
	1834	REAL p(nloc, nd), ph(nloc, nd+1)
	1835	REAL th(nloc, na), gz(nloc, na)
	1836	REAL lv(nloc, na), ep(nloc, na), sigp(nloc, na), clw(nloc, na)
	1837	REAL cpn(nloc, na), tv(nloc, na)
	1838	REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
[879]	1839
[1992]	1840	! outputs:
	1841	REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
	1842	REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
	1843	REAL trap(nloc, na, ntra)
	1844	REAL b(nloc, na)
	1845	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
	1846	! lascendance adiabatique et des flux melanges Pa et Pm.
	1847	! Distinction des wdtrain
	1848	! Pa = wdtrainA Pm = wdtrainM
	1849	REAL wdtraina(nloc, na), wdtrainm(nloc, na)
[879]	1850
[1992]	1851	! local variables
	1852	INTEGER i, j, k, il, num1
	1853	REAL tinv, delti
	1854	REAL awat, afac, afac1, afac2, bfac
	1855	REAL pr1, pr2, sigt, b6, c6, revap, tevap, delth
	1856	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
	1857	REAL ampmax
	1858	REAL lvcp(nloc, na)
	1859	REAL wdtrain(nloc)
	1860	LOGICAL lwork(nloc)
[879]	1861
	1862
[1992]	1863	! ------------------------------------------------------
[879]	1864
[1992]	1865	delti = 1./delt
	1866	tinv = 1./3.
[879]	1867
[1992]	1868	mp(:, :) = 0.
[879]	1869
[1992]	1870	DO i = 1, nl
	1871	DO il = 1, ncum
	1872	mp(il, i) = 0.0
	1873	rp(il, i) = rr(il, i)
	1874	up(il, i) = u(il, i)
	1875	vp(il, i) = v(il, i)
	1876	wt(il, i) = 0.001
	1877	water(il, i) = 0.0
	1878	evap(il, i) = 0.0
	1879	b(il, i) = 0.0
	1880	lvcp(il, i) = lv(il, i)/cpn(il, i)
	1881	END DO
	1882	END DO
[879]	1883
[1992]	1884	! do k=1,ntra
	1885	! do i=1,nd
	1886	! do il=1,ncum
	1887	! trap(il,i,k)=tra(il,i,k)
	1888	! enddo
	1889	! enddo
	1890	! enddo
	1891	! ! RomP >>>
	1892	DO i = 1, nd
	1893	DO il = 1, ncum
	1894	wdtraina(il, i) = 0.0
	1895	wdtrainm(il, i) = 0.0
	1896	END DO
	1897	END DO
	1898	! ! RomP <<<
[879]	1899
[1992]	1900	! * check whether ep(inb)=0, if so, skip precipitating *
	1901	! * downdraft calculation *
[879]	1902
	1903
[1992]	1904	DO il = 1, ncum
	1905	lwork(il) = .TRUE.
	1906	IF (ep(il,inb(il))<0.0001) lwork(il) = .FALSE.
	1907	END DO
[879]	1908
[1992]	1909	CALL zilch(wdtrain, ncum)
[879]	1910
[1992]	1911	DO i = nl + 1, 1, -1
[879]	1912
[1992]	1913	num1 = 0
	1914	DO il = 1, ncum
	1915	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
	1916	END DO
	1917	IF (num1<=0) GO TO 400
[879]	1918
	1919
[1992]	1920	! * integrate liquid water equation to find condensed water *
	1921	! * and condensed water flux *
[879]	1922
	1923
	1924
[1992]	1925	! * begin downdraft loop *
[879]	1926
	1927
	1928
[1992]	1929	! * calculate detrained precipitation *
[879]	1930
[1992]	1931	DO il = 1, ncum
	1932	IF (i<=inb(il) .AND. lwork(il)) THEN
	1933	IF (cvflag_grav) THEN
	1934	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
	1935	wdtraina(il, i) = wdtrain(il)/grav ! Pa 26/08/10 RomP
	1936	ELSE
	1937	wdtrain(il) = 10.0ep(il, i)m(il, i)*clw(il, i)
	1938	wdtraina(il, i) = wdtrain(il)/10. ! Pa 26/08/10 RomP
	1939	END IF
	1940	END IF
	1941	END DO
[879]	1942
[1992]	1943	IF (i>1) THEN
[879]	1944
[1992]	1945	DO j = 1, i - 1
	1946	DO il = 1, ncum
	1947	IF (i<=inb(il) .AND. lwork(il)) THEN
	1948	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
	1949	awat = amax1(awat, 0.0)
	1950	IF (cvflag_grav) THEN
	1951	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
	1952	ELSE
	1953	wdtrain(il) = wdtrain(il) + 10.0awatment(il, j, i)
	1954	END IF
	1955	END IF
	1956	END DO
	1957	END DO
	1958	DO il = 1, ncum
	1959	IF (cvflag_grav) THEN
	1960	wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) ! Pm 26/08/10 RomP
	1961	ELSE
	1962	wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) ! Pm 26/08/10 RomP
	1963	END IF
	1964	END DO
[879]	1965
[1992]	1966	END IF
[879]	1967
	1968
[1992]	1969	! * find rain water and evaporation using provisional *
	1970	! * estimates of rp(i)and rp(i-1) *
[879]	1971
	1972
[1992]	1973	DO il = 1, ncum
[879]	1974
[1992]	1975	IF (i<=inb(il) .AND. lwork(il)) THEN
[879]	1976
[1992]	1977	wt(il, i) = 45.0
[879]	1978
[1992]	1979	IF (i<inb(il)) THEN
	1980	rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, &
	1981	i))+gz(il,i+1)-gz(il,i))/lv(il, i)
	1982	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
	1983	END IF
	1984	rp(il, i) = amax1(rp(il,i), 0.0)
	1985	rp(il, i) = amin1(rp(il,i), rs(il,i))
	1986	rp(il, inb(il)) = rr(il, inb(il))
[879]	1987
[1992]	1988	IF (i==1) THEN
	1989	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
	1990	ELSE
	1991	rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, &
	1992	i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
	1993	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
	1994	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
	1995	rp(il, i-1) = amax1(rp(il,i-1), 0.0)
	1996	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i) &
	1997	)
	1998	afac2 = p(il, i-1)*(rs(il,i-1)-rp(il,i-1))/ &
	1999	(1.0E4+2000.0p(il,i-1)rs(il,i-1))
	2000	afac = 0.5*(afac1+afac2)
	2001	END IF
	2002	IF (i==inb(il)) afac = 0.0
	2003	afac = amax1(afac, 0.0)
	2004	bfac = 1./(sigd*wt(il,i))
[879]	2005
[1992]	2006	! jyg1
	2007	! cc sigt=1.0
	2008	! cc if(i.ge.icb)sigt=sigp(i)
	2009	! prise en compte de la variation progressive de sigt dans
	2010	! les couches icb et icb-1:
	2011	! pour plcl<ph(i+1), pr1=0 & pr2=1
	2012	! pour plcl>ph(i), pr1=1 & pr2=0
	2013	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
	2014	! sur le nuage, et pr2 est la proportion sous la base du
	2015	! nuage.
	2016	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
	2017	pr1 = max(0., min(1.,pr1))
	2018	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
	2019	pr2 = max(0., min(1.,pr2))
	2020	sigt = sigp(il, i)*pr1 + pr2
	2021	! jyg2
[879]	2022
[1992]	2023	b6 = bfac50.sigd(ph(il,i)-ph(il,i+1))sigt*afac
	2024	c6 = water(il, i+1) + bfacwdtrain(il) - 50.sigdbfac(ph(il,i)-ph( &
	2025	il,i+1))*evap(il, i+1)
	2026	IF (c6>0.0) THEN
	2027	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
	2028	evap(il, i) = sigtafacrevap
	2029	water(il, i) = revap*revap
	2030	ELSE
	2031	evap(il, i) = -evap(il, i+1) + 0.02(wdtrain(il)+sigdwt(il,i)* &
	2032	water(il,i+1))/(sigd*(ph(il,i)-ph(il,i+1)))
	2033	END IF
[879]	2034
[1992]	2035	! * calculate precipitating downdraft mass flux under *
	2036	! * hydrostatic approximation *
[879]	2037
[1992]	2038	IF (i/=1) THEN
[879]	2039
[1992]	2040	tevap = amax1(0.0, evap(il,i))
	2041	delth = amax1(0.001, (th(il,i)-th(il,i-1)))
	2042	IF (cvflag_grav) THEN
	2043	mp(il, i) = 100.ginvlvcp(il, i)sigdtevap*(p(il,i-1)-p(il,i))/ &
	2044	delth
	2045	ELSE
	2046	mp(il, i) = 10.lvcp(il, i)sigdtevap(p(il,i-1)-p(il,i))/delth
	2047	END IF
[879]	2048
[1992]	2049	! * if hydrostatic assumption fails, *
	2050	! * solve cubic difference equation for downdraft theta *
	2051	! * and mass flux from two simultaneous differential eqns *
[879]	2052
[1992]	2053	amfac = sigdsigd70.0ph(il, i)(p(il,i-1)-p(il,i))* &
	2054	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
	2055	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
	2056	IF (amp2>(0.1*amfac)) THEN
	2057	xf = 100.0sigdsigdsigd(ph(il,i)-ph(il,i+1))
	2058	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i)/(lvcp(il,i)* &
	2059	sigd*th(il,i))
	2060	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
	2061	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
	2062	50.(p(il,i-1)-p(il,i))xf*tevap
	2063	fac2 = 1.0
	2064	IF (bf<0.0) fac2 = -1.0
	2065	bf = abs(bf)
	2066	ur = 0.25bfbf - afafaftinvtinv*tinv
	2067	IF (ur>=0.0) THEN
	2068	sru = sqrt(ur)
	2069	fac = 1.0
	2070	IF ((0.5*bf-sru)<0.0) fac = -1.0
	2071	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
	2072	fac(abs(0.5bf-sru))**tinv
	2073	ELSE
	2074	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
	2075	IF (fac2<0.0) d = 3.14159 - d
	2076	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
	2077	END IF
	2078	mp(il, i) = amax1(0.0, mp(il,i))
[879]	2079
[1992]	2080	IF (cvflag_grav) THEN
	2081	! jyg : il y a vraisemblablement une erreur dans la ligne 2
	2082	! suivante:
	2083	! il faut diviser par (mp(il,i)sigdgrav) et non par
	2084	! (mp(il,i)+sigd*0.1).
	2085	! Et il faut bien revoir les facteurs 100.
	2086	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap/(mp(il, &
	2087	i)+sigd0.1) - 10.0(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
	2088	)sigdth(il,i))
	2089	ELSE
	2090	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap/(mp(il, &
	2091	i)+sigd0.1) - 10.0(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
	2092	)sigdth(il,i))
	2093	END IF
	2094	b(il, i-1) = amax1(b(il,i-1), 0.0)
	2095	END IF
[879]	2096
[1992]	2097	! *** limit magnitude of mp(i) to meet cfl condition
	2098	! ***
[879]	2099
[1992]	2100	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
	2101	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
	2102	ampmax = amin1(ampmax, amp2)
	2103	mp(il, i) = amin1(mp(il,i), ampmax)
[879]	2104
[1992]	2105	! *** force mp to decrease linearly to zero
	2106	! ***
	2107	! *** between cloud base and the surface
	2108	! ***
[879]	2109
[1992]	2110	IF (p(il,i)>p(il,icb(il))) THEN
	2111	mp(il, i) = mp(il, icb(il))*(p(il,1)-p(il,i))/ &
	2112	(p(il,1)-p(il,icb(il)))
	2113	END IF
[879]	2114
[1992]	2115	END IF ! i.eq.1
[879]	2116
[1992]	2117	! * find mixing ratio of precipitating downdraft *
[879]	2118
	2119
[1992]	2120	IF (i/=inb(il)) THEN
[879]	2121
[1992]	2122	rp(il, i) = rr(il, i)
[879]	2123
[1992]	2124	IF (mp(il,i)>mp(il,i+1)) THEN
[879]	2125
[1992]	2126	IF (cvflag_grav) THEN
	2127	rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
	2128	rr(il, i)(mp(il,i)-mp(il,i+1)) + 100.ginv0.5sigd*(ph(il,i &
	2129	)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
	2130	ELSE
	2131	rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
	2132	rr(il, i)(mp(il,i)-mp(il,i+1)) + 5.sigd*(ph(il,i)-ph(il,i+1 &
	2133	))*(evap(il,i+1)+evap(il,i))
	2134	END IF
	2135	rp(il, i) = rp(il, i)/mp(il, i)
	2136	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+ &
	2137	1))
	2138	up(il, i) = up(il, i)/mp(il, i)
	2139	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+ &
	2140	1))
	2141	vp(il, i) = vp(il, i)/mp(il, i)
[879]	2142
[1992]	2143	! do j=1,ntra
	2144	! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
	2145	! testmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
	2146	! : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
	2147	! trap(il,i,j)=trap(il,i,j)/mp(il,i)
	2148	! end do
[879]	2149
[1992]	2150	ELSE
[1742]	2151
[1992]	2152	IF (mp(il,i+1)>1.0E-16) THEN
	2153	IF (cvflag_grav) THEN
	2154	rp(il, i) = rp(il, i+1) + 100.ginv0.5sigd(ph(il,i)-ph(il, &
	2155	i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1)
	2156	ELSE
	2157	rp(il, i) = rp(il, i+1) + 5.sigd(ph(il,i)-ph(il,i+1))*(evap &
	2158	(il,i+1)+evap(il,i))/mp(il, i+1)
	2159	END IF
	2160	up(il, i) = up(il, i+1)
	2161	vp(il, i) = vp(il, i+1)
[1742]	2162
[1992]	2163	! do j=1,ntra
	2164	! trap(il,i,j)=trap(il,i+1,j)
	2165	! end do
[1742]	2166
[1992]	2167	END IF
	2168	END IF
	2169	rp(il, i) = amin1(rp(il,i), rs(il,i))
	2170	rp(il, i) = amax1(rp(il,i), 0.0)
[1742]	2171
[1992]	2172	END IF
	2173	END IF
	2174	END DO
[1742]	2175
[1992]	2176	400 END DO
[879]	2177
[1992]	2178	RETURN
	2179	END SUBROUTINE cv30_unsat
[879]	2180
[1992]	2181	SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, &
	2182	tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, &
	2183	wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, &
	2184	tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, &
	2185	mike, tls, tps, qcondc, wd)
	2186	IMPLICIT NONE
[879]	2187
[1992]	2188	include "cvthermo.h"
	2189	include "cv30param.h"
	2190	include "cvflag.h"
	2191	include "conema3.h"
[879]	2192
[1992]	2193	! inputs:
	2194	INTEGER ncum, nd, na, ntra, nloc
	2195	INTEGER icb(nloc), inb(nloc)
	2196	REAL delt
	2197	REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
	2198	REAL tra(nloc, nd, ntra), sig(nloc, nd)
	2199	REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
	2200	REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
	2201	REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
	2202	REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
	2203	REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
	2204	REAL water(nloc, na), evap(nloc, na), b(nloc, na)
	2205	REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
	2206	! ym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
	2207	REAL vent(nloc, na, na), elij(nloc, na, na)
	2208	INTEGER nent(nloc, na)
	2209	REAL traent(nloc, na, na, ntra)
	2210	REAL tv(nloc, nd), tvp(nloc, nd)
[879]	2211
[1992]	2212	! input/output:
	2213	INTEGER iflag(nloc)
[879]	2214
[1992]	2215	! outputs:
	2216	REAL precip(nloc)
	2217	REAL vprecip(nloc, nd+1)
	2218	REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	2219	REAL ftra(nloc, nd, ntra)
	2220	REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
	2221	REAL dnwd0(nloc, nd), mike(nloc, nd)
	2222	REAL tls(nloc, nd), tps(nloc, nd)
	2223	REAL qcondc(nloc, nd) ! cld
	2224	REAL wd(nloc) ! gust
[879]	2225
[1992]	2226	! local variables:
	2227	INTEGER i, k, il, n, j, num1
	2228	REAL rat, awat, delti
	2229	REAL ax, bx, cx, dx, ex
	2230	REAL cpinv, rdcp, dpinv
	2231	REAL lvcp(nloc, na), mke(nloc, na)
	2232	REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
	2233	! !! real up1(nloc), dn1(nloc)
	2234	REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
	2235	REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
	2236	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
	2237	REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
[879]	2238
	2239
[1992]	2240	! -------------------------------------------------------------
[879]	2241
[1992]	2242	! initialization:
[879]	2243
[1992]	2244	delti = 1.0/delt
[879]	2245
[1992]	2246	DO il = 1, ncum
	2247	precip(il) = 0.0
	2248	wd(il) = 0.0 ! gust
	2249	vprecip(il, nd+1) = 0.
	2250	END DO
	2251
	2252	DO i = 1, nd
	2253	DO il = 1, ncum
	2254	vprecip(il, i) = 0.0
	2255	ft(il, i) = 0.0
	2256	fr(il, i) = 0.0
	2257	fu(il, i) = 0.0
	2258	fv(il, i) = 0.0
	2259	qcondc(il, i) = 0.0 ! cld
	2260	qcond(il, i) = 0.0 ! cld
	2261	nqcond(il, i) = 0.0 ! cld
	2262	END DO
	2263	END DO
	2264
	2265	! do j=1,ntra
	2266	! do i=1,nd
	2267	! do il=1,ncum
	2268	! ftra(il,i,j)=0.0
	2269	! enddo
	2270	! enddo
	2271	! enddo
	2272
	2273	DO i = 1, nl
	2274	DO il = 1, ncum
	2275	lvcp(il, i) = lv(il, i)/cpn(il, i)
	2276	END DO
	2277	END DO
	2278
	2279
	2280
	2281	! * calculate surface precipitation in mm/day *
	2282
	2283	DO il = 1, ncum
	2284	IF (ep(il,inb(il))>=0.0001) THEN
	2285	IF (cvflag_grav) THEN
	2286	precip(il) = wt(il, 1)sigdwater(il, 1)86400.1000./(rowl*grav)
	2287	ELSE
	2288	precip(il) = wt(il, 1)sigdwater(il, 1)*8640.
	2289	END IF
	2290	END IF
	2291	END DO
	2292
	2293	! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
	2294
	2295	! MAF rajout pour lessivage
	2296	DO k = 1, nl
	2297	DO il = 1, ncum
	2298	IF (k<=inb(il)) THEN
	2299	IF (cvflag_grav) THEN
	2300	vprecip(il, k) = wt(il, k)sigdwater(il, k)/grav
	2301	ELSE
	2302	vprecip(il, k) = wt(il, k)sigdwater(il, k)/10.
	2303	END IF
	2304	END IF
	2305	END DO
	2306	END DO
	2307
	2308
	2309	! * Calculate downdraft velocity scale *
	2310	! * NE PAS UTILISER POUR L'INSTANT *
	2311
	2312	! ! do il=1,ncum
	2313	! ! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
	2314	! ! : /(sigd*p(il,icb(il)))
	2315	! ! enddo
	2316
	2317
	2318	! * calculate tendencies of lowest level potential temperature *
	2319	! * and mixing ratio *
	2320
	2321	DO il = 1, ncum
	2322	work(il) = 1.0/(ph(il,1)-ph(il,2))
	2323	am(il) = 0.0
	2324	END DO
	2325
	2326	DO k = 2, nl
	2327	DO il = 1, ncum
	2328	IF (k<=inb(il)) THEN
	2329	am(il) = am(il) + m(il, k)
	2330	END IF
	2331	END DO
	2332	END DO
	2333
	2334	DO il = 1, ncum
	2335
	2336	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
	2337	IF (cvflag_grav) THEN
	2338	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
	2339	ft(il, 1) = 0.01gravwork(il)am(il)(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
	2340	1))/cpn(il,1))
	2341	ELSE
	2342	IF ((0.1work(il)am(il))>=delti) iflag(il) = 1 !consistency vect
	2343	ft(il, 1) = 0.1work(il)am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
	2344	1))/cpn(il,1))
	2345	END IF
	2346
	2347	ft(il, 1) = ft(il, 1) - 0.5lvcp(il, 1)sigd*(evap(il,1)+evap(il,2))
	2348
	2349	IF (cvflag_grav) THEN
	2350	ft(il, 1) = ft(il, 1) - 0.009gravsigdmp(il, 2)t(il, 1)b(il, 1) &
	2351	work(il)
	2352	ELSE
	2353	ft(il, 1) = ft(il, 1) - 0.09sigdmp(il, 2)t(il, 1)b(il, 1)*work(il)
	2354	END IF
	2355
	2356	ft(il, 1) = ft(il, 1) + 0.01sigdwt(il, 1)(cl-cpd)water(il, 2)*(t(il,2 &
	2357	)-t(il,1))*work(il)/cpn(il, 1)
	2358
	2359	IF (cvflag_grav) THEN
	2360	! jyg1 Correction pour mieux conserver l'eau (conformite avec
	2361	! CONVECT4.3)
	2362	! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas
	2363	! evap)
	2364	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr(il,1))work(il) + &
	2365	sigd0.5(evap(il,1)+evap(il,2))
	2366	! +tard : +sigd*evap(il,1)
	2367
	2368	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il)
	2369
	2370	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il, &
	2371	1))+am(il)*(u(il,2)-u(il,1)))
	2372	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il, &
	2373	1))+am(il)*(v(il,2)-v(il,1)))
	2374	ELSE ! cvflag_grav
	2375	fr(il, 1) = 0.1mp(il, 2)(rp(il,2)-rr(il,1))*work(il) + &
	2376	sigd0.5(evap(il,1)+evap(il,2))
	2377	fr(il, 1) = fr(il, 1) + 0.1am(il)(rr(il,2)-rr(il,1))*work(il)
	2378	fu(il, 1) = fu(il, 1) + 0.1work(il)(mp(il,2)*(up(il,2)-u(il, &
	2379	1))+am(il)*(u(il,2)-u(il,1)))
	2380	fv(il, 1) = fv(il, 1) + 0.1work(il)(mp(il,2)*(vp(il,2)-v(il, &
	2381	1))+am(il)*(v(il,2)-v(il,1)))
	2382	END IF ! cvflag_grav
	2383
	2384	END DO ! il
	2385
	2386	! do j=1,ntra
	2387	! do il=1,ncum
	2388	! if (cvflag_grav) then
	2389	! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
	2390	! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2391	! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2392	! else
	2393	! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
	2394	! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2395	! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2396	! endif
	2397	! enddo
	2398	! enddo
	2399
	2400	DO j = 2, nl
	2401	DO il = 1, ncum
	2402	IF (j<=inb(il)) THEN
	2403	IF (cvflag_grav) THEN
	2404	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il, &
	2405	j,1)-rr(il,1))
	2406	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il, &
	2407	j,1)-u(il,1))
	2408	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il, &
	2409	j,1)-v(il,1))
	2410	ELSE ! cvflag_grav
	2411	fr(il, 1) = fr(il, 1) + 0.1work(il)ment(il, j, 1)*(qent(il,j,1)- &
	2412	rr(il,1))
	2413	fu(il, 1) = fu(il, 1) + 0.1work(il)ment(il, j, 1)*(uent(il,j,1)-u &
	2414	(il,1))
	2415	fv(il, 1) = fv(il, 1) + 0.1work(il)ment(il, j, 1)*(vent(il,j,1)-v &
	2416	(il,1))
	2417	END IF ! cvflag_grav
	2418	END IF ! j
	2419	END DO
	2420	END DO
	2421
	2422	! do k=1,ntra
	2423	! do j=2,nl
	2424	! do il=1,ncum
	2425	! if (j.le.inb(il)) then
	2426
	2427	! if (cvflag_grav) then
	2428	! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
	2429	! : *(traent(il,j,1,k)-tra(il,1,k))
	2430	! else
	2431	! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
	2432	! : *(traent(il,j,1,k)-tra(il,1,k))
	2433	! endif
	2434
	2435	! endif
	2436	! enddo
	2437	! enddo
	2438	! enddo
	2439
	2440
	2441	! * calculate tendencies of potential temperature and mixing ratio *
	2442	! * at levels above the lowest level *
	2443
	2444	! * first find the net saturated updraft and downdraft mass fluxes *
	2445	! * through each level *
	2446
	2447
	2448	DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
	2449
	2450	num1 = 0
	2451	DO il = 1, ncum
	2452	IF (i<=inb(il)) num1 = num1 + 1
	2453	END DO
	2454	IF (num1<=0) GO TO 500
	2455
	2456	CALL zilch(amp1, ncum)
	2457	CALL zilch(ad, ncum)
	2458
	2459	DO k = i + 1, nl + 1
	2460	DO il = 1, ncum
	2461	IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN
	2462	amp1(il) = amp1(il) + m(il, k)
	2463	END IF
	2464	END DO
	2465	END DO
	2466
	2467	DO k = 1, i
	2468	DO j = i + 1, nl + 1
	2469	DO il = 1, ncum
	2470	IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
	2471	amp1(il) = amp1(il) + ment(il, k, j)
	2472	END IF
	2473	END DO
	2474	END DO
	2475	END DO
	2476
	2477	DO k = 1, i - 1
	2478	DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
	2479	DO il = 1, ncum
	2480	IF (i<=inb(il) .AND. j<=inb(il)) THEN
	2481	ad(il) = ad(il) + ment(il, j, k)
	2482	END IF
	2483	END DO
	2484	END DO
	2485	END DO
	2486
	2487	DO il = 1, ncum
	2488	IF (i<=inb(il)) THEN
	2489	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	2490	cpinv = 1.0/cpn(il, i)
	2491
	2492	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
	2493	IF (cvflag_grav) THEN
	2494	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
	2495	ELSE
	2496	IF ((0.1dpinvamp1(il))>=delti) iflag(il) = 1 ! vecto
	2497	END IF
	2498
	2499	IF (cvflag_grav) THEN
	2500	ft(il, i) = 0.01gravdpinv(amp1(il)(t(il,i+1)-t(il, &
	2501	i)+(gz(il,i+1)-gz(il,i))cpinv)-ad(il)(t(il,i)-t(il, &
	2502	i-1)+(gz(il,i)-gz(il,i-1))cpinv)) - 0.5sigdlvcp(il, i)(evap( &
	2503	il,i)+evap(il,i+1))
	2504	rat = cpn(il, i-1)*cpinv
	2505	ft(il, i) = ft(il, i) - 0.009gravsigd(mp(il,i+1)t(il,i)*b(il,i) &
	2506	-mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
	2507	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h( &
	2508	il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
	2509	ELSE ! cvflag_grav
	2510	ft(il, i) = 0.1dpinv(amp1(il)*(t(il,i+1)-t(il, &
	2511	i)+(gz(il,i+1)-gz(il,i))cpinv)-ad(il)(t(il,i)-t(il, &
	2512	i-1)+(gz(il,i)-gz(il,i-1))cpinv)) - 0.5sigdlvcp(il, i)(evap( &
	2513	il,i)+evap(il,i+1))
	2514	rat = cpn(il, i-1)*cpinv
	2515	ft(il, i) = ft(il, i) - 0.09sigd(mp(il,i+1)t(il,i)b(il,i)-mp(il &
	2516	,i)t(il,i-1)ratb(il,i-1))dpinv
	2517	ft(il, i) = ft(il, i) + 0.1dpinvment(il, i, i)*(hp(il,i)-h(il,i)+ &
	2518	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
	2519	END IF ! cvflag_grav
	2520
	2521
	2522	ft(il, i) = ft(il, i) + 0.01sigdwt(il, i)(cl-cpd)water(il, i+1)*( &
	2523	t(il,i+1)-t(il,i))dpinvcpinv
	2524
	2525	IF (cvflag_grav) THEN
	2526	fr(il, i) = 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il, &
	2527	i))-ad(il)*(rr(il,i)-rr(il,i-1)))
	2528	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il, &
	2529	i))-ad(il)*(u(il,i)-u(il,i-1)))
	2530	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il, &
	2531	i))-ad(il)*(v(il,i)-v(il,i-1)))
	2532	ELSE ! cvflag_grav
	2533	fr(il, i) = 0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il, &
	2534	i))-ad(il)*(rr(il,i)-rr(il,i-1)))
	2535	fu(il, i) = fu(il, i) + 0.1dpinv(amp1(il)*(u(il,i+1)-u(il, &
	2536	i))-ad(il)*(u(il,i)-u(il,i-1)))
	2537	fv(il, i) = fv(il, i) + 0.1dpinv(amp1(il)*(v(il,i+1)-v(il, &
	2538	i))-ad(il)*(v(il,i)-v(il,i-1)))
	2539	END IF ! cvflag_grav
	2540
	2541	END IF ! i
	2542	END DO
	2543
	2544	! do k=1,ntra
	2545	! do il=1,ncum
	2546	! if (i.le.inb(il)) then
	2547	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2548	! cpinv=1.0/cpn(il,i)
	2549	! if (cvflag_grav) then
	2550	! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
	2551	! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2552	! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2553	! else
	2554	! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
	2555	! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2556	! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2557	! endif
	2558	! endif
	2559	! enddo
	2560	! enddo
	2561
	2562	DO k = 1, i - 1
	2563	DO il = 1, ncum
	2564	IF (i<=inb(il)) THEN
	2565	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	2566	cpinv = 1.0/cpn(il, i)
	2567
	2568	awat = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
	2569	awat = amax1(awat, 0.0)
	2570
	2571	IF (cvflag_grav) THEN
	2572	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k &
	2573	,i)-awat-rr(il,i))
	2574	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
	2575	,i)-u(il,i))
	2576	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k &
	2577	,i)-v(il,i))
	2578	ELSE ! cvflag_grav
	2579	fr(il, i) = fr(il, i) + 0.1dpinvment(il, k, i)*(qent(il,k,i)- &
	2580	awat-rr(il,i))
	2581	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
	2582	,i)-u(il,i))
	2583	fv(il, i) = fv(il, i) + 0.1dpinvment(il, k, i)*(vent(il,k,i)-v( &
	2584	il,i))
	2585	END IF ! cvflag_grav
	2586
	2587	! (saturated updrafts resulting from mixing) ! cld
	2588	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat) ! cld
	2589	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2590	END IF ! i
	2591	END DO
	2592	END DO
	2593
	2594	! do j=1,ntra
	2595	! do k=1,i-1
	2596	! do il=1,ncum
	2597	! if (i.le.inb(il)) then
	2598	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2599	! cpinv=1.0/cpn(il,i)
	2600	! if (cvflag_grav) then
	2601	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2602	! : *(traent(il,k,i,j)-tra(il,i,j))
	2603	! else
	2604	! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2605	! : *(traent(il,k,i,j)-tra(il,i,j))
	2606	! endif
	2607	! endif
	2608	! enddo
	2609	! enddo
	2610	! enddo
	2611
	2612	DO k = i, nl + 1
	2613	DO il = 1, ncum
	2614	IF (i<=inb(il) .AND. k<=inb(il)) THEN
	2615	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	2616	cpinv = 1.0/cpn(il, i)
	2617
	2618	IF (cvflag_grav) THEN
	2619	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k &
	2620	,i)-rr(il,i))
	2621	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
	2622	,i)-u(il,i))
	2623	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k &
	2624	,i)-v(il,i))
	2625	ELSE ! cvflag_grav
	2626	fr(il, i) = fr(il, i) + 0.1dpinvment(il, k, i)*(qent(il,k,i)-rr &
	2627	(il,i))
	2628	fu(il, i) = fu(il, i) + 0.1dpinvment(il, k, i)*(uent(il,k,i)-u( &
	2629	il,i))
	2630	fv(il, i) = fv(il, i) + 0.1dpinvment(il, k, i)*(vent(il,k,i)-v( &
	2631	il,i))
	2632	END IF ! cvflag_grav
	2633	END IF ! i and k
	2634	END DO
	2635	END DO
	2636
	2637	! do j=1,ntra
	2638	! do k=i,nl+1
	2639	! do il=1,ncum
	2640	! if (i.le.inb(il) .and. k.le.inb(il)) then
	2641	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2642	! cpinv=1.0/cpn(il,i)
	2643	! if (cvflag_grav) then
	2644	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2645	! : *(traent(il,k,i,j)-tra(il,i,j))
	2646	! else
	2647	! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2648	! : *(traent(il,k,i,j)-tra(il,i,j))
	2649	! endif
	2650	! endif ! i and k
	2651	! enddo
	2652	! enddo
	2653	! enddo
	2654
	2655	DO il = 1, ncum
	2656	IF (i<=inb(il)) THEN
	2657	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	2658	cpinv = 1.0/cpn(il, i)
	2659
	2660	IF (cvflag_grav) THEN
	2661	! sb: on ne fait pas encore la correction permettant de mieux
	2662	! conserver l'eau:
	2663	fr(il, i) = fr(il, i) + 0.5sigd(evap(il,i)+evap(il,i+1)) + &
	2664	0.01grav(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)(rp(il, &
	2665	i)-rr(il,i-1)))*dpinv
	2666
	2667	fu(il, i) = fu(il, i) + 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il, &
	2668	i))-mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	2669	fv(il, i) = fv(il, i) + 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il, &
	2670	i))-mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
	2671	ELSE ! cvflag_grav
	2672	fr(il, i) = fr(il, i) + 0.5sigd(evap(il,i)+evap(il,i+1)) + &
	2673	0.1(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il,i)-rr(il, &
	2674	i-1)))*dpinv
	2675	fu(il, i) = fu(il, i) + 0.1(mp(il,i+1)(up(il,i+1)-u(il, &
	2676	i))-mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	2677	fv(il, i) = fv(il, i) + 0.1(mp(il,i+1)(vp(il,i+1)-v(il, &
	2678	i))-mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
	2679	END IF ! cvflag_grav
	2680
	2681	END IF ! i
	2682	END DO
	2683
	2684	! sb: interface with the cloud parameterization: ! cld
	2685
	2686	DO k = i + 1, nl
	2687	DO il = 1, ncum
	2688	IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld
	2689	! (saturated downdrafts resulting from mixing) ! cld
	2690	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
	2691	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2692	END IF ! cld
	2693	END DO ! cld
	2694	END DO ! cld
	2695
	2696	! (particular case: no detraining level is found) ! cld
	2697	DO il = 1, ncum ! cld
	2698	IF (i<=inb(il) .AND. nent(il,i)==0) THEN ! cld
	2699	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
	2700	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2701	END IF ! cld
	2702	END DO ! cld
	2703
	2704	DO il = 1, ncum ! cld
	2705	IF (i<=inb(il) .AND. nqcond(il,i)/=0.) THEN ! cld
	2706	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
	2707	END IF ! cld
	2708	END DO
	2709
	2710	! do j=1,ntra
	2711	! do il=1,ncum
	2712	! if (i.le.inb(il)) then
	2713	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2714	! cpinv=1.0/cpn(il,i)
	2715
	2716	! if (cvflag_grav) then
	2717	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
	2718	! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2719	! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
	2720	! else
	2721	! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
	2722	! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2723	! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
	2724	! endif
	2725	! endif ! i
	2726	! enddo
	2727	! enddo
	2728
	2729	500 END DO
	2730
	2731
	2732	! * move the detrainment at level inb down to level inb-1 *
	2733	! * in such a way as to preserve the vertically *
	2734	! * integrated enthalpy and water tendencies *
	2735
	2736	DO il = 1, ncum
	2737
	2738	ax = 0.1ment(il, inb(il), inb(il))(hp(il,inb(il))-h(il,inb(il))+t(il, &
	2739	inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il), &
	2740	inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
	2741	ft(il, inb(il)) = ft(il, inb(il)) - ax
	2742	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il &
	2743	))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il, &
	2744	inb(il))))
	2745
	2746	bx = 0.1ment(il, inb(il), inb(il))(qent(il,inb(il),inb(il))-rr(il,inb( &
	2747	il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2748	fr(il, inb(il)) = fr(il, inb(il)) - bx
	2749	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+ &
	2750	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
	2751
	2752	cx = 0.1ment(il, inb(il), inb(il))(uent(il,inb(il),inb(il))-u(il,inb(il &
	2753	)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2754	fu(il, inb(il)) = fu(il, inb(il)) - cx
	2755	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+ &
	2756	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
	2757
	2758	dx = 0.1ment(il, inb(il), inb(il))(vent(il,inb(il),inb(il))-v(il,inb(il &
	2759	)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2760	fv(il, inb(il)) = fv(il, inb(il)) - dx
	2761	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+ &
	2762	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
	2763
	2764	END DO
	2765
	2766	! do j=1,ntra
	2767	! do il=1,ncum
	2768	! ex=0.1*ment(il,inb(il),inb(il))
	2769	! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
	2770	! : /(ph(il,inb(il))-ph(il,inb(il)+1))
	2771	! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
	2772	! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
	2773	! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
	2774	! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2775	! enddo
	2776	! enddo
	2777
	2778
	2779	! * homoginize tendencies below cloud base *
	2780
	2781
	2782	DO il = 1, ncum
	2783	asum(il) = 0.0
	2784	bsum(il) = 0.0
	2785	csum(il) = 0.0
	2786	dsum(il) = 0.0
	2787	END DO
	2788
	2789	DO i = 1, nl
	2790	DO il = 1, ncum
	2791	IF (i<=(icb(il)-1)) THEN
	2792	asum(il) = asum(il) + ft(il, i)*(ph(il,i)-ph(il,i+1))
	2793	bsum(il) = bsum(il) + fr(il, i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il, &
	2794	1)))*(ph(il,i)-ph(il,i+1))
	2795	csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
	2796	1)))*(ph(il,i)-ph(il,i+1))
	2797	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
	2798	END IF
	2799	END DO
	2800	END DO
	2801
	2802	! !!! do 700 i=1,icb(il)-1
	2803	DO i = 1, nl
	2804	DO il = 1, ncum
	2805	IF (i<=(icb(il)-1)) THEN
	2806	ft(il, i) = asum(il)t(il, i)/(th(il,i)dsum(il))
	2807	fr(il, i) = bsum(il)/csum(il)
	2808	END IF
	2809	END DO
	2810	END DO
	2811
	2812
	2813	! * reset counter and return *
	2814
	2815	DO il = 1, ncum
	2816	sig(il, nd) = 2.0
	2817	END DO
	2818
	2819
	2820	DO i = 1, nd
	2821	DO il = 1, ncum
	2822	upwd(il, i) = 0.0
	2823	dnwd(il, i) = 0.0
	2824	END DO
	2825	END DO
	2826
	2827	DO i = 1, nl
	2828	DO il = 1, ncum
	2829	dnwd0(il, i) = -mp(il, i)
	2830	END DO
	2831	END DO
	2832	DO i = nl + 1, nd
	2833	DO il = 1, ncum
	2834	dnwd0(il, i) = 0.
	2835	END DO
	2836	END DO
	2837
	2838
	2839	DO i = 1, nl
	2840	DO il = 1, ncum
	2841	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	2842	upwd(il, i) = 0.0
	2843	dnwd(il, i) = 0.0
	2844	END IF
	2845	END DO
	2846	END DO
	2847
	2848	DO i = 1, nl
	2849	DO k = 1, nl
	2850	DO il = 1, ncum
	2851	up1(il, k, i) = 0.0
	2852	dn1(il, k, i) = 0.0
	2853	END DO
	2854	END DO
	2855	END DO
	2856
	2857	DO i = 1, nl
	2858	DO k = i, nl
	2859	DO n = 1, i - 1
	2860	DO il = 1, ncum
	2861	IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
	2862	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
	2863	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
	2864	END IF
	2865	END DO
	2866	END DO
	2867	END DO
	2868	END DO
	2869
	2870	DO i = 2, nl
	2871	DO k = i, nl
	2872	DO il = 1, ncum
	2873	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
	2874	! then
	2875	IF (i<=inb(il) .AND. k<=inb(il)) THEN
	2876	upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
	2877	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
	2878	END IF
	2879	END DO
	2880	END DO
	2881	END DO
	2882
	2883
	2884	! !!! DO il=1,ncum
	2885	! !!! do i=icb(il),inb(il)
	2886	! !!!
	2887	! !!! upwd(il,i)=0.0
	2888	! !!! dnwd(il,i)=0.0
	2889	! !!! do k=i,inb(il)
	2890	! !!! up1=0.0
	2891	! !!! dn1=0.0
	2892	! !!! do n=1,i-1
	2893	! !!! up1=up1+ment(il,n,k)
	2894	! !!! dn1=dn1-ment(il,k,n)
	2895	! !!! enddo
	2896	! !!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
	2897	! !!! dnwd(il,i)=dnwd(il,i)+dn1
	2898	! !!! enddo
	2899	! !!! enddo
	2900	! !!!
	2901	! !!! ENDDO
	2902
	2903	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2904	! determination de la variation de flux ascendant entre
	2905	! deux niveau non dilue mike
	2906	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2907
	2908	DO i = 1, nl
	2909	DO il = 1, ncum
	2910	mike(il, i) = m(il, i)
	2911	END DO
	2912	END DO
	2913
	2914	DO i = nl + 1, nd
	2915	DO il = 1, ncum
	2916	mike(il, i) = 0.
	2917	END DO
	2918	END DO
	2919
	2920	DO i = 1, nd
	2921	DO il = 1, ncum
	2922	ma(il, i) = 0
	2923	END DO
	2924	END DO
	2925
	2926	DO i = 1, nl
	2927	DO j = i, nl
	2928	DO il = 1, ncum
	2929	ma(il, i) = ma(il, i) + m(il, j)
	2930	END DO
	2931	END DO
	2932	END DO
	2933
	2934	DO i = nl + 1, nd
	2935	DO il = 1, ncum
	2936	ma(il, i) = 0.
	2937	END DO
	2938	END DO
	2939
	2940	DO i = 1, nl
	2941	DO il = 1, ncum
	2942	IF (i<=(icb(il)-1)) THEN
	2943	ma(il, i) = 0
	2944	END IF
	2945	END DO
	2946	END DO
	2947
	2948	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2949	! icb represente de niveau ou se trouve la
	2950	! base du nuage , et inb le top du nuage
	2951	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2952
	2953	DO i = 1, nd
	2954	DO il = 1, ncum
	2955	mke(il, i) = upwd(il, i) + dnwd(il, i)
	2956	END DO
	2957	END DO
	2958
	2959	DO i = 1, nd
	2960	DO il = 1, ncum
	2961	rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd*(1.-rr(il, &
	2962	i))+rr(il,i)*cpv)
	2963	tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp
	2964	tps(il, i) = tp(il, i)
	2965	END DO
	2966	END DO
	2967
	2968
	2969	! * diagnose the in-cloud mixing ratio * ! cld
	2970	! * of condensed water * ! cld
	2971	! ! cld
	2972
	2973	DO i = 1, nd ! cld
	2974	DO il = 1, ncum ! cld
	2975	mac(il, i) = 0.0 ! cld
	2976	wa(il, i) = 0.0 ! cld
	2977	siga(il, i) = 0.0 ! cld
	2978	sax(il, i) = 0.0 ! cld
	2979	END DO ! cld
	2980	END DO ! cld
	2981
	2982	DO i = minorig, nl ! cld
	2983	DO k = i + 1, nl + 1 ! cld
	2984	DO il = 1, ncum ! cld
	2985	IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN ! cld
	2986	mac(il, i) = mac(il, i) + m(il, k) ! cld
	2987	END IF ! cld
	2988	END DO ! cld
	2989	END DO ! cld
	2990	END DO ! cld
	2991
	2992	DO i = 1, nl ! cld
	2993	DO j = 1, i ! cld
	2994	DO il = 1, ncum ! cld
	2995	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
	2996	.AND. j>=icb(il)) THEN ! cld
	2997	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
	2998	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
	2999	END IF ! cld
	3000	END DO ! cld
	3001	END DO ! cld
	3002	END DO ! cld
	3003
	3004	DO i = 1, nl ! cld
	3005	DO il = 1, ncum ! cld
	3006	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
	3007	.AND. sax(il,i)>0.0) THEN ! cld
	3008	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
	3009	END IF ! cld
	3010	END DO ! cld
	3011	END DO ! cld
	3012
	3013	DO i = 1, nl ! cld
	3014	DO il = 1, ncum ! cld
	3015	IF (wa(il,i)>0.0) & ! cld
	3016	siga(il, i) = mac(il, i)/wa(il, i) & ! cld
	3017	rrdtvp(il, i)/p(il, i)/100./delta ! cld
	3018	siga(il, i) = min(siga(il,i), 1.0) ! cld
	3019	! IM cf. FH
	3020	IF (iflag_clw==0) THEN
	3021	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
	3022	+(1.-siga(il,i))*qcond(il, i) ! cld
	3023	ELSE IF (iflag_clw==1) THEN
	3024	qcondc(il, i) = qcond(il, i) ! cld
	3025	END IF
	3026
	3027	END DO ! cld
	3028	END DO ! cld
	3029
	3030	RETURN
	3031	END SUBROUTINE cv30_yield
	3032
	3033	! !RomP >>>
	3034	SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, &
	3035	d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
	3036	IMPLICIT NONE
	3037
	3038	include "cv30param.h"
	3039
	3040	! inputs:
	3041	INTEGER ncum, nd, na, nloc, len
	3042	REAL ment(nloc, na, na), sij(nloc, na, na)
	3043	REAL clw(nloc, nd), elij(nloc, na, na)
	3044	REAL ep(nloc, na)
	3045	INTEGER icb(nloc), inb(nloc)
	3046	REAL vprecip(nloc, nd+1)
	3047	! ouputs:
	3048	REAL da(nloc, na), phi(nloc, na, na)
	3049	REAL phi2(nloc, na, na)
	3050	REAL d1a(nloc, na), dam(nloc, na)
	3051	REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
	3052	! variables pour tracer dans precip de l'AA et des mel
	3053	! local variables:
	3054	INTEGER i, j, k
	3055	REAL epm(nloc, na, na)
	3056
	3057	! variables d'Emanuel : du second indice au troisieme
	3058	! ---> tab(i,k,j) -> de l origine k a l arrivee j
	3059	! ment, sij, elij
	3060	! variables personnelles : du troisieme au second indice
	3061	! ---> tab(i,j,k) -> de k a j
	3062	! phi, phi2
	3063
	3064	! initialisations
	3065	DO j = 1, na
	3066	DO i = 1, ncum
	3067	da(i, j) = 0.
	3068	d1a(i, j) = 0.
	3069	dam(i, j) = 0.
	3070	eplamm(i, j) = 0.
	3071	END DO
	3072	END DO
	3073	DO k = 1, na
	3074	DO j = 1, na
	3075	DO i = 1, ncum
	3076	epm(i, j, k) = 0.
	3077	epmlmmm(i, j, k) = 0.
	3078	phi(i, j, k) = 0.
	3079	phi2(i, j, k) = 0.
	3080	END DO
	3081	END DO
	3082	END DO
	3083
	3084	! fraction deau condensee dans les melanges convertie en precip : epm
	3085	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
	3086	DO j = 1, na
	3087	DO k = 1, j - 1
	3088	DO i = 1, ncum
	3089	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
	3090	! !jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
	3091	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
	3092	! !
	3093	epm(i, j, k) = max(epm(i,j,k), 0.0)
	3094	END IF
	3095	END DO
	3096	END DO
	3097	END DO
	3098
	3099	DO j = 1, na
	3100	DO k = 1, na
	3101	DO i = 1, ncum
	3102	IF (k>=icb(i) .AND. k<=inb(i)) THEN
	3103	eplamm(i, j) = eplamm(i, j) + ep(i, j)clw(i, j)ment(i, j, k)*(1.- &
	3104	sij(i,j,k))
	3105	END IF
	3106	END DO
	3107	END DO
	3108	END DO
	3109
	3110	DO j = 1, na
	3111	DO k = 1, j - 1
	3112	DO i = 1, ncum
	3113	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
	3114	epmlmmm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
	3115	END IF
	3116	END DO
	3117	END DO
	3118	END DO
	3119
	3120	! matrices pour calculer la tendance des concentrations dans cvltr.F90
	3121	DO j = 1, na
	3122	DO k = 1, na
	3123	DO i = 1, ncum
	3124	da(i, j) = da(i, j) + (1.-sij(i,k,j))*ment(i, k, j)
	3125	phi(i, j, k) = sij(i, k, j)*ment(i, k, j)
	3126	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sij(i,k,j))
	3127	END DO
	3128	END DO
	3129	END DO
	3130
	3131	DO j = 1, na
	3132	DO k = 1, j - 1
	3133	DO i = 1, ncum
	3134	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, j, k)(1.-ep(i,k))*(1.- &
	3135	sij(i,k,j))
	3136	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
	3137	END DO
	3138	END DO
	3139	END DO
	3140
	3141	RETURN
	3142	END SUBROUTINE cv30_tracer
	3143	! RomP <<<
	3144
	3145	SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
	3146	vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
	3147	dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
	3148	epmlmmm, eplamm, wdtraina, wdtrainm, iflag1, precip1, vprecip1, evap1, &
	3149	ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
	3150	dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
	3151	elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1)
	3152	IMPLICIT NONE
	3153
	3154	include "cv30param.h"
	3155
	3156	! inputs:
	3157	INTEGER len, ncum, nd, ntra, nloc
	3158	INTEGER idcum(nloc)
	3159	INTEGER iflag(nloc)
	3160	INTEGER inb(nloc)
	3161	REAL precip(nloc)
	3162	REAL vprecip(nloc, nd+1), evap(nloc, nd)
	3163	REAL ep(nloc, nd)
	3164	REAL sig(nloc, nd), w0(nloc, nd)
	3165	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	3166	REAL ftra(nloc, nd, ntra)
	3167	REAL ma(nloc, nd)
	3168	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
	3169	REAL qcondc(nloc, nd)
	3170	REAL wd(nloc), cape(nloc)
	3171	REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
	3172	! RomP >>>
	3173	REAL phi2(nloc, nd, nd)
	3174	REAL d1a(nloc, nd), dam(nloc, nd)
	3175	REAL wdtraina(nloc, nd), wdtrainm(nloc, nd)
	3176	REAL sij(nloc, nd, nd)
	3177	REAL elij(nloc, nd, nd), clw(nloc, nd)
	3178	REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd)
	3179	! RomP <<<
	3180
	3181	! outputs:
	3182	INTEGER iflag1(len)
	3183	INTEGER inb1(len)
	3184	REAL precip1(len)
	3185	REAL vprecip1(len, nd+1), evap1(len, nd) !<<< RomP
	3186	REAL ep1(len, nd) !<<< RomP
	3187	REAL sig1(len, nd), w01(len, nd)
	3188	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
	3189	REAL ftra1(len, nd, ntra)
	3190	REAL ma1(len, nd)
	3191	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
	3192	REAL qcondc1(nloc, nd)
	3193	REAL wd1(nloc), cape1(nloc)
	3194	REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd)
	3195	! RomP >>>
	3196	REAL phi21(len, nd, nd)
	3197	REAL d1a1(len, nd), dam1(len, nd)
	3198	REAL wdtraina1(len, nd), wdtrainm1(len, nd)
	3199	REAL sij1(len, nd, nd)
	3200	REAL elij1(len, nd, nd), clw1(len, nd)
	3201	REAL epmlmmm1(len, nd, nd), eplamm1(len, nd)
	3202	! RomP <<<
	3203
	3204	! local variables:
	3205	INTEGER i, k, j
	3206
	3207	DO i = 1, ncum
	3208	precip1(idcum(i)) = precip(i)
	3209	iflag1(idcum(i)) = iflag(i)
	3210	wd1(idcum(i)) = wd(i)
	3211	inb1(idcum(i)) = inb(i)
	3212	cape1(idcum(i)) = cape(i)
	3213	END DO
	3214
	3215	DO k = 1, nl
	3216	DO i = 1, ncum
	3217	vprecip1(idcum(i), k) = vprecip(i, k)
	3218	evap1(idcum(i), k) = evap(i, k) !<<< RomP
	3219	sig1(idcum(i), k) = sig(i, k)
	3220	w01(idcum(i), k) = w0(i, k)
	3221	ft1(idcum(i), k) = ft(i, k)
	3222	fq1(idcum(i), k) = fq(i, k)
	3223	fu1(idcum(i), k) = fu(i, k)
	3224	fv1(idcum(i), k) = fv(i, k)
	3225	ma1(idcum(i), k) = ma(i, k)
	3226	upwd1(idcum(i), k) = upwd(i, k)
	3227	dnwd1(idcum(i), k) = dnwd(i, k)
	3228	dnwd01(idcum(i), k) = dnwd0(i, k)
	3229	qcondc1(idcum(i), k) = qcondc(i, k)
	3230	da1(idcum(i), k) = da(i, k)
	3231	mp1(idcum(i), k) = mp(i, k)
	3232	! RomP >>>
	3233	ep1(idcum(i), k) = ep(i, k)
	3234	d1a1(idcum(i), k) = d1a(i, k)
	3235	dam1(idcum(i), k) = dam(i, k)
	3236	clw1(idcum(i), k) = clw(i, k)
	3237	eplamm1(idcum(i), k) = eplamm(i, k)
	3238	wdtraina1(idcum(i), k) = wdtraina(i, k)
	3239	wdtrainm1(idcum(i), k) = wdtrainm(i, k)
	3240	! RomP <<<
	3241	END DO
	3242	END DO
	3243
	3244	DO i = 1, ncum
	3245	sig1(idcum(i), nd) = sig(i, nd)
	3246	END DO
	3247
	3248
	3249	! do 2100 j=1,ntra
	3250	! do 2110 k=1,nd ! oct3
	3251	! do 2120 i=1,ncum
	3252	! ftra1(idcum(i),k,j)=ftra(i,k,j)
	3253	! 2120 continue
	3254	! 2110 continue
	3255	! 2100 continue
	3256	DO j = 1, nd
	3257	DO k = 1, nd
	3258	DO i = 1, ncum
	3259	sij1(idcum(i), k, j) = sij(i, k, j)
	3260	phi1(idcum(i), k, j) = phi(i, k, j)
	3261	phi21(idcum(i), k, j) = phi2(i, k, j)
	3262	elij1(idcum(i), k, j) = elij(i, k, j)
	3263	epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
	3264	END DO
	3265	END DO
	3266	END DO
	3267
	3268	RETURN
	3269	END SUBROUTINE cv30_uncompress
	3270

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

Download in other formats:

Original Format