Context Navigation

cv30_routines.F90 @ 4786

Last change on this file since 4786 was 2542, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2487:2541 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 98.3 KB

Rev	Line
[1992]	1
[1403]	2	! $Id: cv30_routines.F90 2542 2016-06-06 14:04:57Z idelkadi $
[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
[2542]	409	icb1(i) = min(max(icb(i), 2), nl)
[1992]	410	! if icb is below LCL, start loop at ICB+1:
	411	! (icbs est le premier niveau au-dessus du LCL)
	412	icbs(i) = icb1(i) !convect3
	413	IF (plcl(i)<p(i,icb1(i))) THEN
	414	icbs(i) = min(icbs(i)+1, nl) !convect3
	415	END IF
	416	END DO !convect3
[879]	417
[1992]	418	DO i = 1, len !convect3
	419	ticb(i) = t(i, icbs(i)) !convect3
	420	gzicb(i) = gz(i, icbs(i)) !convect3
	421	qsicb(i) = qs(i, icbs(i)) !convect3
	422	END DO !convect3
[879]	423
	424
[1992]	425	! Re-compute icbsmax (icbsmax2): !convect3
	426	! !convect3
	427	icbsmax2 = 2 !convect3
	428	DO i = 1, len !convect3
	429	icbsmax2 = max(icbsmax2, icbs(i)) !convect3
	430	END DO !convect3
[879]	431
[1992]	432	! initialization outputs:
[879]	433
[1992]	434	DO k = 1, icbsmax2 ! convect3
	435	DO i = 1, len ! convect3
	436	tp(i, k) = 0.0 ! convect3
	437	tvp(i, k) = 0.0 ! convect3
	438	clw(i, k) = 0.0 ! convect3
	439	END DO ! convect3
	440	END DO ! convect3
[879]	441
[1992]	442	! tp and tvp below cloud base:
[879]	443
[1992]	444	DO k = minorig, icbsmax2 - 1
	445	DO i = 1, len
	446	tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
	447	tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
	448	END DO
	449	END DO
[879]	450
[1992]	451	! * Find lifted parcel quantities above cloud base *
[879]	452
[1992]	453	DO i = 1, len
	454	tg = ticb(i)
	455	! ori qg=qs(i,icb(i))
	456	qg = qsicb(i) ! convect3
	457	! debug alv=lv0-clmcpv*(ticb(i)-t0)
	458	alv = lv0 - clmcpv*(ticb(i)-273.15)
[879]	459
[1992]	460	! First iteration.
[879]	461
[1992]	462	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	463	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
	464	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
	465	s = 1./s
	466	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	467	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
	468	tg = tg + s*(ah0(i)-ahg)
	469	! ori tg=max(tg,35.0)
	470	! debug tc=tg-t0
	471	tc = tg - 273.15
	472	denom = 243.5 + tc
	473	denom = max(denom, 1.0) ! convect3
	474	! ori if(tc.ge.0.0)then
	475	es = 6.112exp(17.67tc/denom)
	476	! ori else
	477	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	478	! ori endif
	479	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	480	qg = epses/(p(i,icbs(i))-es(1.-eps))
[879]	481
[1992]	482	! Second iteration.
[879]	483
	484
[1992]	485	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	486	! ori s=1./s
	487	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	488	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
	489	tg = tg + s*(ah0(i)-ahg)
	490	! ori tg=max(tg,35.0)
	491	! debug tc=tg-t0
	492	tc = tg - 273.15
	493	denom = 243.5 + tc
	494	denom = max(denom, 1.0) ! convect3
	495	! ori if(tc.ge.0.0)then
	496	es = 6.112exp(17.67tc/denom)
	497	! ori else
	498	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	499	! ori end if
	500	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	501	qg = epses/(p(i,icbs(i))-es(1.-eps))
[879]	502
[1992]	503	alv = lv0 - clmcpv*(ticb(i)-273.15)
[879]	504
[1992]	505	! ori c approximation here:
	506	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	507	! ori & -gz(i,icb(i))-alv*qg)/cpd
[879]	508
[1992]	509	! convect3: no approximation:
	510	tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alvqg)/(cpd+(cl-cpd)qnk(i))
[879]	511
[1992]	512	! ori clw(i,icb(i))=qnk(i)-qg
	513	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	514	clw(i, icbs(i)) = qnk(i) - qg
	515	clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
[879]	516
[1992]	517	rg = qg/(1.-qnk(i))
	518	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	519	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
	520	tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
[879]	521
[1992]	522	END DO
[879]	523
[1992]	524	! ori do 380 k=minorig,icbsmax2
	525	! ori do 370 i=1,len
	526	! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
	527	! ori 370 continue
	528	! ori 380 continue
[879]	529
	530
[1992]	531	! -- The following is only for convect3:
[879]	532
[1992]	533	! * icbs is the first level above the LCL:
	534	! if plcl<p(icb), then icbs=icb+1
	535	! if plcl>p(icb), then icbs=icb
[1403]	536
[1992]	537	! * the routine above computes tvp from minorig to icbs (included).
[879]	538
[1992]	539	! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
	540	! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
[879]	541
[1992]	542	! * therefore, in the case icbs=icb, we compute tvp at level icb+1
	543	! (tvp at other levels will be computed in cv3_undilute2.F)
[879]	544
	545
[1992]	546	DO i = 1, len
	547	ticb(i) = t(i, icb(i)+1)
	548	gzicb(i) = gz(i, icb(i)+1)
	549	qsicb(i) = qs(i, icb(i)+1)
	550	END DO
[879]	551
[1992]	552	DO i = 1, len
	553	tg = ticb(i)
	554	qg = qsicb(i) ! convect3
	555	! debug alv=lv0-clmcpv*(ticb(i)-t0)
	556	alv = lv0 - clmcpv*(ticb(i)-273.15)
[879]	557
[1992]	558	! First iteration.
[879]	559
[1992]	560	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	561	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
	562	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
	563	s = 1./s
	564	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	565	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
	566	tg = tg + s*(ah0(i)-ahg)
	567	! ori tg=max(tg,35.0)
	568	! debug tc=tg-t0
	569	tc = tg - 273.15
	570	denom = 243.5 + tc
	571	denom = max(denom, 1.0) ! convect3
	572	! ori if(tc.ge.0.0)then
	573	es = 6.112exp(17.67tc/denom)
	574	! ori else
	575	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	576	! ori endif
	577	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	578	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
[879]	579
[1992]	580	! Second iteration.
[879]	581
	582
[1992]	583	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	584	! ori s=1./s
	585	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	586	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
	587	tg = tg + s*(ah0(i)-ahg)
	588	! ori tg=max(tg,35.0)
	589	! debug tc=tg-t0
	590	tc = tg - 273.15
	591	denom = 243.5 + tc
	592	denom = max(denom, 1.0) ! convect3
	593	! ori if(tc.ge.0.0)then
	594	es = 6.112exp(17.67tc/denom)
	595	! ori else
	596	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	597	! ori end if
	598	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
	599	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
[879]	600
[1992]	601	alv = lv0 - clmcpv*(ticb(i)-273.15)
[879]	602
[1992]	603	! ori c approximation here:
	604	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	605	! ori & -gz(i,icb(i))-alv*qg)/cpd
[879]	606
[1992]	607	! convect3: no approximation:
	608	tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alvqg)/(cpd+(cl-cpd)qnk(i))
[879]	609
[1992]	610	! ori clw(i,icb(i))=qnk(i)-qg
	611	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	612	clw(i, icb(i)+1) = qnk(i) - qg
	613	clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
[879]	614
[1992]	615	rg = qg/(1.-qnk(i))
	616	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	617	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
	618	tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
[879]	619
[1992]	620	END DO
[879]	621
[1992]	622	RETURN
	623	END SUBROUTINE cv30_undilute1
[879]	624
[1992]	625	SUBROUTINE cv30_trigger(len, nd, icb, plcl, p, th, tv, tvp, pbase, buoybase, &
	626	iflag, sig, w0)
	627	IMPLICIT NONE
[879]	628
[1992]	629	! -------------------------------------------------------------------
	630	! --- TRIGGERING
[879]	631
[1992]	632	! - computes the cloud base
	633	! - triggering (crude in this version)
	634	! - relaxation of sig and w0 when no convection
[879]	635
[1992]	636	! Caution1: if no convection, we set iflag=4
	637	! (it used to be 0 in convect3)
[879]	638
[1992]	639	! Caution2: at this stage, tvp (and thus buoy) are know up
	640	! through icb only!
	641	! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
	642	! -------------------------------------------------------------------
[879]	643
[1992]	644	include "cv30param.h"
[879]	645
[1992]	646	! input:
	647	INTEGER len, nd
	648	INTEGER icb(len)
	649	REAL plcl(len), p(len, nd)
	650	REAL th(len, nd), tv(len, nd), tvp(len, nd)
[879]	651
[1992]	652	! output:
	653	REAL pbase(len), buoybase(len)
[879]	654
[1992]	655	! input AND output:
	656	INTEGER iflag(len)
	657	REAL sig(len, nd), w0(len, nd)
[879]	658
[1992]	659	! local variables:
	660	INTEGER i, k
	661	REAL tvpbase, tvbase, tdif, ath, ath1
[879]	662
	663
[1992]	664	! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
[879]	665
[1992]	666	DO i = 1, len
	667	pbase(i) = plcl(i) + dpbase
	668	tvpbase = tvp(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
	669	(p(i,icb(i))-p(i,icb(i)+1)) + tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/( &
	670	p(i,icb(i))-p(i,icb(i)+1))
	671	tvbase = tv(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
	672	(p(i,icb(i))-p(i,icb(i)+1)) + tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/(p &
	673	(i,icb(i))-p(i,icb(i)+1))
	674	buoybase(i) = tvpbase - tvbase
	675	END DO
[879]	676
	677
[1992]	678	! * make sure that column is dry adiabatic between the surface *
	679	! * and cloud base, and that lifted air is positively buoyant *
	680	! * at cloud base *
	681	! * if not, return to calling program after resetting *
	682	! * sig(i) and w0(i) *
[879]	683
	684
[1992]	685	! oct3 do 200 i=1,len
	686	! oct3
	687	! oct3 tdif = buoybase(i)
	688	! oct3 ath1 = th(i,1)
	689	! oct3 ath = th(i,icb(i)-1) - dttrig
	690	! oct3
	691	! oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
	692	! oct3 do 60 k=1,nl
	693	! oct3 sig(i,k) = betasig(i,k) - 2.alphatdiftdif
	694	! oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
	695	! oct3 w0(i,k) = beta*w0(i,k)
	696	! oct3 60 continue
	697	! oct3 iflag(i)=4 ! pour version vectorisee
	698	! oct3c convect3 iflag(i)=0
	699	! oct3cccc return
	700	! oct3 endif
	701	! oct3
	702	! oct3200 continue
[879]	703
[1992]	704	! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
[879]	705
[1992]	706	DO k = 1, nl
	707	DO i = 1, len
[879]	708
[1992]	709	tdif = buoybase(i)
	710	ath1 = th(i, 1)
	711	ath = th(i, icb(i)-1) - dttrig
[879]	712
[1992]	713	IF (tdif<dtcrit .OR. ath>ath1) THEN
	714	sig(i, k) = betasig(i, k) - 2.alphatdiftdif
	715	sig(i, k) = amax1(sig(i,k), 0.0)
	716	w0(i, k) = beta*w0(i, k)
	717	iflag(i) = 4 ! pour version vectorisee
	718	! convect3 iflag(i)=0
	719	END IF
[879]	720
[1992]	721	END DO
	722	END DO
[879]	723
[1992]	724	! fin oct3 --
[879]	725
[1992]	726	RETURN
	727	END SUBROUTINE cv30_trigger
[879]	728
[1992]	729	SUBROUTINE cv30_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, &
	730	plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, &
	731	th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, &
	732	iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, &
	733	v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0)
[2408]	734	USE print_control_mod, ONLY: lunout
[1992]	735	IMPLICIT NONE
[879]	736
[1992]	737	include "cv30param.h"
[879]	738
[1992]	739	! inputs:
	740	INTEGER len, ncum, nd, ntra, nloc
	741	INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
	742	REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
	743	REAL pbase1(len), buoybase1(len)
	744	REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
	745	REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
	746	REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
	747	REAL tvp1(len, nd), clw1(len, nd)
	748	REAL th1(len, nd)
	749	REAL sig1(len, nd), w01(len, nd)
	750	REAL tra1(len, nd, ntra)
[879]	751
[1992]	752	! outputs:
	753	! en fait, on a nloc=len pour l'instant (cf cv_driver)
	754	INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
	755	REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
	756	REAL pbase(nloc), buoybase(nloc)
	757	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
	758	REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
	759	REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
	760	REAL tvp(nloc, nd), clw(nloc, nd)
	761	REAL th(nloc, nd)
	762	REAL sig(nloc, nd), w0(nloc, nd)
	763	REAL tra(nloc, nd, ntra)
[879]	764
[1992]	765	! local variables:
	766	INTEGER i, k, nn, j
[879]	767
[1992]	768	CHARACTER (LEN=20) :: modname = 'cv30_compress'
	769	CHARACTER (LEN=80) :: abort_message
[879]	770
	771
[1992]	772	DO k = 1, nl + 1
	773	nn = 0
	774	DO i = 1, len
	775	IF (iflag1(i)==0) THEN
	776	nn = nn + 1
	777	sig(nn, k) = sig1(i, k)
	778	w0(nn, k) = w01(i, k)
	779	t(nn, k) = t1(i, k)
	780	q(nn, k) = q1(i, k)
	781	qs(nn, k) = qs1(i, k)
	782	u(nn, k) = u1(i, k)
	783	v(nn, k) = v1(i, k)
	784	gz(nn, k) = gz1(i, k)
	785	h(nn, k) = h1(i, k)
	786	lv(nn, k) = lv1(i, k)
	787	cpn(nn, k) = cpn1(i, k)
	788	p(nn, k) = p1(i, k)
	789	ph(nn, k) = ph1(i, k)
	790	tv(nn, k) = tv1(i, k)
	791	tp(nn, k) = tp1(i, k)
	792	tvp(nn, k) = tvp1(i, k)
	793	clw(nn, k) = clw1(i, k)
	794	th(nn, k) = th1(i, k)
	795	END IF
	796	END DO
	797	END DO
[879]	798
[1992]	799	! do 121 j=1,ntra
	800	! do 111 k=1,nd
	801	! nn=0
	802	! do 101 i=1,len
	803	! if(iflag1(i).eq.0)then
	804	! nn=nn+1
	805	! tra(nn,k,j)=tra1(i,k,j)
	806	! endif
	807	! 101 continue
	808	! 111 continue
	809	! 121 continue
[879]	810
[1992]	811	IF (nn/=ncum) THEN
	812	WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
	813	abort_message = ''
[2408]	814	CALL abort_physic(modname, abort_message, 1)
[1992]	815	END IF
[879]	816
[1992]	817	nn = 0
	818	DO i = 1, len
	819	IF (iflag1(i)==0) THEN
	820	nn = nn + 1
	821	pbase(nn) = pbase1(i)
	822	buoybase(nn) = buoybase1(i)
	823	plcl(nn) = plcl1(i)
	824	tnk(nn) = tnk1(i)
	825	qnk(nn) = qnk1(i)
	826	gznk(nn) = gznk1(i)
	827	nk(nn) = nk1(i)
	828	icb(nn) = icb1(i)
	829	icbs(nn) = icbs1(i)
	830	iflag(nn) = iflag1(i)
	831	END IF
	832	END DO
[879]	833
[1992]	834	RETURN
	835	END SUBROUTINE cv30_compress
[879]	836
[1992]	837	SUBROUTINE cv30_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
	838	q, qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
	839	ep, sigp, buoy)
[2488]	840	! epmax_cape: ajout arguments
[1992]	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
[1992]	1244	! -------------------------------------------------------
	1245	! -- Initialization
	1246	! -------------------------------------------------------
[879]	1247
[1992]	1248	DO k = 1, nl
	1249	DO i = 1, ncum
	1250	m(i, k) = 0.0
	1251	END DO
	1252	END DO
[879]	1253
[1992]	1254	! -------------------------------------------------------
	1255	! -- Reset sig(i) and w0(i) for i>inb and i<icb
	1256	! -------------------------------------------------------
[879]	1257
[1992]	1258	! update sig and w0 above LNB:
[879]	1259
[1992]	1260	DO k = 1, nl - 1
	1261	DO i = 1, ncum
	1262	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
	1263	sig(i, k) = betasig(i, k) + 2.alphabuoy(i, inb(i))abs(buoy(i,inb( &
	1264	i)))
	1265	sig(i, k) = amax1(sig(i,k), 0.0)
	1266	w0(i, k) = beta*w0(i, k)
	1267	END IF
	1268	END DO
	1269	END DO
[879]	1270
[1992]	1271	! compute icbmax:
[879]	1272
[1992]	1273	icbmax = 2
	1274	DO i = 1, ncum
	1275	icbmax = max(icbmax, icb(i))
	1276	END DO
[879]	1277
[1992]	1278	! update sig and w0 below cloud base:
[879]	1279
[1992]	1280	DO k = 1, icbmax
	1281	DO i = 1, ncum
	1282	IF (k<=icb(i)) THEN
	1283	sig(i, k) = betasig(i, k) - 2.alphabuoy(i, icb(i))buoy(i, icb(i))
	1284	sig(i, k) = amax1(sig(i,k), 0.0)
	1285	w0(i, k) = beta*w0(i, k)
	1286	END IF
	1287	END DO
	1288	END DO
[879]	1289
[1992]	1290	! ! if(inb.lt.(nl-1))then
	1291	! ! do 85 i=inb+1,nl-1
	1292	! ! sig(i)=betasig(i)+2.alphabuoy(inb)
	1293	! ! 1 abs(buoy(inb))
	1294	! ! sig(i)=amax1(sig(i),0.0)
	1295	! ! w0(i)=beta*w0(i)
	1296	! ! 85 continue
	1297	! ! end if
[879]	1298
[1992]	1299	! ! do 87 i=1,icb
	1300	! ! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
	1301	! ! sig(i)=amax1(sig(i),0.0)
	1302	! ! w0(i)=beta*w0(i)
	1303	! ! 87 continue
[1742]	1304
[1992]	1305	! -------------------------------------------------------------
	1306	! -- Reset fractional areas of updrafts and w0 at initial time
	1307	! -- and after 10 time steps of no convection
	1308	! -------------------------------------------------------------
[1742]	1309
[1992]	1310	DO k = 1, nl - 1
	1311	DO i = 1, ncum
	1312	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
	1313	sig(i, k) = 0.0
	1314	w0(i, k) = 0.0
	1315	END IF
	1316	END DO
	1317	END DO
[879]	1318
[1992]	1319	! -------------------------------------------------------------
	1320	! -- Calculate convective available potential energy (cape),
	1321	! -- vertical velocity (w), fractional area covered by
	1322	! -- undilute updraft (sig), and updraft mass flux (m)
	1323	! -------------------------------------------------------------
[879]	1324
[1992]	1325	DO i = 1, ncum
	1326	cape(i) = 0.0
	1327	END DO
[879]	1328
[1992]	1329	! compute dtmin (minimum buoyancy between ICB and given level k):
[879]	1330
[1992]	1331	DO i = 1, ncum
	1332	DO k = 1, nl
	1333	dtmin(i, k) = 100.0
	1334	END DO
	1335	END DO
[879]	1336
[1992]	1337	DO i = 1, ncum
	1338	DO k = 1, nl
	1339	DO j = minorig, nl
	1340	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- &
	1341	1))) THEN
	1342	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
	1343	END IF
	1344	END DO
	1345	END DO
	1346	END DO
[879]	1347
[1992]	1348	! the interval on which cape is computed starts at pbase :
	1349	DO k = 1, nl
	1350	DO i = 1, ncum
[879]	1351
[1992]	1352	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
[879]	1353
[1992]	1354	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
	1355	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
	1356	cape(i) = amax1(0.0, cape(i))
	1357	sigold(i, k) = sig(i, k)
[879]	1358
[1992]	1359	! dtmin(i,k)=100.0
	1360	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
	1361	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
	1362	! 97 continue
[879]	1363
[1992]	1364	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
	1365	sig(i, k) = amax1(sig(i,k), 0.0)
	1366	sig(i, k) = amin1(sig(i,k), 0.01)
	1367	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
	1368	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
	1369	amu = 0.5(sig(i,k)+sigold(i,k))w
	1370	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
	1371	w0(i, k) = w
	1372	END IF
[879]	1373
[1992]	1374	END DO
	1375	END DO
[879]	1376
[1992]	1377	DO i = 1, ncum
	1378	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
	1379	m(i, icb(i)) = 0.5m(i, icb(i)+1)(ph(i,icb(i))-ph(i,icb(i)+1))/ &
	1380	(ph(i,icb(i)+1)-ph(i,icb(i)+2))
	1381	sig(i, icb(i)) = sig(i, icb(i)+1)
	1382	sig(i, icb(i)-1) = sig(i, icb(i))
	1383	END DO
[879]	1384
	1385
[1992]	1386	! ! cape=0.0
	1387	! ! do 98 i=icb+1,inb
	1388	! ! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
	1389	! ! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
	1390	! ! dcape=rrdbuoy(i-1)deltap/p(i-1)
	1391	! ! dlnp=deltap/p(i-1)
	1392	! ! cape=amax1(0.0,cape)
	1393	! ! sigold=sig(i)
[879]	1394
[1992]	1395	! ! dtmin=100.0
	1396	! ! do 97 j=icb,i-1
	1397	! ! dtmin=amin1(dtmin,buoy(j))
	1398	! ! 97 continue
[879]	1399
[1992]	1400	! ! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
	1401	! ! sig(i)=amax1(sig(i),0.0)
	1402	! ! sig(i)=amin1(sig(i),0.01)
	1403	! ! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
	1404	! ! w=(1.-beta)facsqrt(cape)+beta*w0(i)
	1405	! ! amu=0.5(sig(i)+sigold)w
	1406	! ! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
	1407	! ! w0(i)=w
	1408	! ! 98 continue
	1409	! ! w0(icb)=0.5*w0(icb+1)
	1410	! ! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
	1411	! ! sig(icb)=sig(icb+1)
	1412	! ! sig(icb-1)=sig(icb)
[879]	1413
[1992]	1414	RETURN
	1415	END SUBROUTINE cv30_closure
[879]	1416
[1992]	1417	SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
	1418	u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, &
	1419	vent, sij, elij, ments, qents, traent)
	1420	IMPLICIT NONE
[879]	1421
[1992]	1422	! ---------------------------------------------------------------------
	1423	! a faire:
	1424	! - changer rr(il,1) -> qnk(il)
	1425	! - vectorisation de la partie normalisation des flux (do 789...)
	1426	! ---------------------------------------------------------------------
[879]	1427
[1992]	1428	include "cvthermo.h"
	1429	include "cv30param.h"
[879]	1430
[1992]	1431	! inputs:
	1432	INTEGER ncum, nd, na, ntra, nloc
	1433	INTEGER icb(nloc), inb(nloc), nk(nloc)
	1434	REAL sig(nloc, nd)
	1435	REAL qnk(nloc)
	1436	REAL ph(nloc, nd+1)
	1437	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
	1438	REAL u(nloc, nd), v(nloc, nd)
	1439	REAL tra(nloc, nd, ntra) ! input of convect3
	1440	REAL lv(nloc, na), h(nloc, na), hp(nloc, na)
	1441	REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
	1442	REAL m(nloc, na) ! input of convect3
[879]	1443
[1992]	1444	! outputs:
	1445	REAL ment(nloc, na, na), qent(nloc, na, na)
	1446	REAL uent(nloc, na, na), vent(nloc, na, na)
	1447	REAL sij(nloc, na, na), elij(nloc, na, na)
	1448	REAL traent(nloc, nd, nd, ntra)
	1449	REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
	1450	REAL sigij(nloc, nd, nd)
[879]	1451
[1992]	1452	! local variables:
	1453	INTEGER i, j, k, il, im, jm
	1454	INTEGER num1, num2
	1455	INTEGER nent(nloc, na)
	1456	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
	1457	REAL alt, smid, sjmin, sjmax, delp, delm
	1458	REAL asij(nloc), smax(nloc), scrit(nloc)
	1459	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
	1460	REAL wgh
	1461	REAL zm(nloc, na)
	1462	LOGICAL lwork(nloc)
[879]	1463
[1992]	1464	! =====================================================================
	1465	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	1466	! =====================================================================
[879]	1467
[1992]	1468	! ori do 360 i=1,ncum*nlp
	1469	DO j = 1, nl
	1470	DO i = 1, ncum
	1471	nent(i, j) = 0
	1472	! in convect3, m is computed in cv3_closure
	1473	! ori m(i,1)=0.0
	1474	END DO
	1475	END DO
[879]	1476
[1992]	1477	! ori do 400 k=1,nlp
	1478	! ori do 390 j=1,nlp
	1479	DO j = 1, nl
	1480	DO k = 1, nl
	1481	DO i = 1, ncum
	1482	qent(i, k, j) = rr(i, j)
	1483	uent(i, k, j) = u(i, j)
	1484	vent(i, k, j) = v(i, j)
	1485	elij(i, k, j) = 0.0
	1486	! ym ment(i,k,j)=0.0
	1487	! ym sij(i,k,j)=0.0
	1488	END DO
	1489	END DO
	1490	END DO
[879]	1491
[1992]	1492	! ym
	1493	ment(1:ncum, 1:nd, 1:nd) = 0.0
	1494	sij(1:ncum, 1:nd, 1:nd) = 0.0
[879]	1495
[1992]	1496	! do k=1,ntra
	1497	! do j=1,nd ! instead nlp
	1498	! do i=1,nd ! instead nlp
	1499	! do il=1,ncum
	1500	! traent(il,i,j,k)=tra(il,j,k)
	1501	! enddo
	1502	! enddo
	1503	! enddo
	1504	! enddo
	1505	zm(:, :) = 0.
[879]	1506
[1992]	1507	! =====================================================================
	1508	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	1509	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	1510	! --- FRACTION (sij)
	1511	! =====================================================================
[879]	1512
[1992]	1513	DO i = minorig + 1, nl
[879]	1514
[1992]	1515	DO j = minorig, nl
	1516	DO il = 1, ncum
	1517	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
	1518	1)) .AND. (j<=inb(il))) THEN
[879]	1519
[1992]	1520	rti = rr(il, 1) - ep(il, i)*clw(il, i)
	1521	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
	1522	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
	1523	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
	1524	dei = denom
	1525	IF (abs(dei)<0.01) dei = 0.01
	1526	sij(il, i, j) = anum/dei
	1527	sij(il, i, i) = 1.0
	1528	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
	1529	altem = altem/bf2
	1530	cwat = clw(il, j)*(1.-ep(il,j))
	1531	stemp = sij(il, i, j)
	1532	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
	1533	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
	1534	denom = denom + lv(il, j)*(rr(il,i)-rti)
	1535	IF (abs(denom)<0.01) denom = 0.01
	1536	sij(il, i, j) = anum/denom
	1537	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - &
	1538	rs(il, j)
	1539	altem = altem - (bf2-1.)*cwat
	1540	END IF
	1541	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
	1542	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
	1543	uent(il, i, j) = sij(il, i, j)*u(il, i) + &
	1544	(1.-sij(il,i,j))*u(il, nk(il))
	1545	vent(il, i, j) = sij(il, i, j)*v(il, i) + &
	1546	(1.-sij(il,i,j))*v(il, nk(il))
	1547	! !!! do k=1,ntra
	1548	! !!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1549	! !!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1550	! !!! end do
	1551	elij(il, i, j) = altem
	1552	elij(il, i, j) = amax1(0.0, elij(il,i,j))
	1553	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
	1554	nent(il, i) = nent(il, i) + 1
	1555	END IF
	1556	sij(il, i, j) = amax1(0.0, sij(il,i,j))
	1557	sij(il, i, j) = amin1(1.0, sij(il,i,j))
	1558	END IF ! new
	1559	END DO
	1560	END DO
[879]	1561
[1992]	1562	! do k=1,ntra
	1563	! do j=minorig,nl
	1564	! do il=1,ncum
	1565	! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
	1566	! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
	1567	! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1568	! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1569	! endif
	1570	! enddo
	1571	! enddo
	1572	! enddo
[879]	1573
	1574
[1992]	1575	! *** if no air can entrain at level i assume that updraft detrains
	1576	! ***
	1577	! *** at that level and calculate detrained air flux and properties
	1578	! ***
[879]	1579
	1580
[1992]	1581	! @ do 170 i=icb(il),inb(il)
[879]	1582
[1992]	1583	DO il = 1, ncum
	1584	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
	1585	! @ if(nent(il,i).eq.0)then
	1586	ment(il, i, i) = m(il, i)
	1587	qent(il, i, i) = rr(il, nk(il)) - ep(il, i)*clw(il, i)
	1588	uent(il, i, i) = u(il, nk(il))
	1589	vent(il, i, i) = v(il, nk(il))
	1590	elij(il, i, i) = clw(il, i)
	1591	! MAF sij(il,i,i)=1.0
	1592	sij(il, i, i) = 0.0
	1593	END IF
	1594	END DO
	1595	END DO
[879]	1596
[1992]	1597	! do j=1,ntra
	1598	! do i=minorig+1,nl
	1599	! do il=1,ncum
	1600	! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
	1601	! traent(il,i,i,j)=tra(il,nk(il),j)
	1602	! endif
	1603	! enddo
	1604	! enddo
	1605	! enddo
[879]	1606
[1992]	1607	DO j = minorig, nl
	1608	DO i = minorig, nl
	1609	DO il = 1, ncum
	1610	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
	1611	inb(il))) THEN
	1612	sigij(il, i, j) = sij(il, i, j)
	1613	END IF
	1614	END DO
	1615	END DO
	1616	END DO
	1617	! @ enddo
[879]	1618
[1992]	1619	! @170 continue
[879]	1620
[1992]	1621	! =====================================================================
	1622	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	1623	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	1624	! =====================================================================
[879]	1625
[1992]	1626	! ym call zilch(asum,ncum*nd)
	1627	! ym call zilch(bsum,ncum*nd)
	1628	! ym call zilch(csum,ncum*nd)
	1629	CALL zilch(asum, nloc*nd)
	1630	CALL zilch(csum, nloc*nd)
	1631	CALL zilch(csum, nloc*nd)
[879]	1632
[1992]	1633	DO il = 1, ncum
	1634	lwork(il) = .FALSE.
	1635	END DO
[879]	1636
[1992]	1637	DO i = minorig + 1, nl
[879]	1638
[1992]	1639	num1 = 0
	1640	DO il = 1, ncum
	1641	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
	1642	END DO
	1643	IF (num1<=0) GO TO 789
[879]	1644
	1645
[1992]	1646	DO il = 1, ncum
	1647	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	1648	lwork(il) = (nent(il,i)/=0)
	1649	qp = rr(il, 1) - ep(il, i)*clw(il, i)
	1650	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
	1651	(cpv-cpd)t(il, i)(qp-rr(il,i))
	1652	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
	1653	(cpd-cpv)t(il, i)(rr(il,i)-qp)
	1654	IF (abs(denom)<0.01) denom = 0.01
	1655	scrit(il) = anum/denom
	1656	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
	1657	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
	1658	smax(il) = 0.0
	1659	asij(il) = 0.0
	1660	END IF
	1661	END DO
[879]	1662
[1992]	1663	DO j = nl, minorig, -1
[879]	1664
[1992]	1665	num2 = 0
	1666	DO il = 1, ncum
	1667	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	1668	il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
	1669	END DO
	1670	IF (num2<=0) GO TO 175
[879]	1671
[1992]	1672	DO il = 1, ncum
	1673	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
	1674	il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
[879]	1675
[1992]	1676	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
	1677	wgh = 1.0
	1678	IF (j>i) THEN
	1679	sjmax = amax1(sij(il,i,j+1), smax(il))
	1680	sjmax = amin1(sjmax, scrit(il))
	1681	smax(il) = amax1(sij(il,i,j), smax(il))
	1682	sjmin = amax1(sij(il,i,j-1), smax(il))
	1683	sjmin = amin1(sjmin, scrit(il))
	1684	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
	1685	smid = amin1(sij(il,i,j), scrit(il))
	1686	ELSE
	1687	sjmax = amax1(sij(il,i,j+1), scrit(il))
	1688	smid = amax1(sij(il,i,j), scrit(il))
	1689	sjmin = 0.0
	1690	IF (j>1) sjmin = sij(il, i, j-1)
	1691	sjmin = amax1(sjmin, scrit(il))
	1692	END IF
	1693	delp = abs(sjmax-smid)
	1694	delm = abs(sjmin-smid)
	1695	asij(il) = asij(il) + wgh*(delp+delm)
	1696	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
	1697	END IF
	1698	END IF
	1699	END DO
[879]	1700
[1992]	1701	175 END DO
[879]	1702
[1992]	1703	DO il = 1, ncum
	1704	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	1705	asij(il) = amax1(1.0E-16, asij(il))
	1706	asij(il) = 1.0/asij(il)
	1707	asum(il, i) = 0.0
	1708	bsum(il, i) = 0.0
	1709	csum(il, i) = 0.0
	1710	END IF
	1711	END DO
[879]	1712
[1992]	1713	DO j = minorig, nl
	1714	DO il = 1, ncum
	1715	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	1716	il)-1) .AND. j<=inb(il)) THEN
	1717	ment(il, i, j) = ment(il, i, j)*asij(il)
	1718	END IF
	1719	END DO
	1720	END DO
[879]	1721
[1992]	1722	DO j = minorig, nl
	1723	DO il = 1, ncum
	1724	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	1725	il)-1) .AND. j<=inb(il)) THEN
	1726	asum(il, i) = asum(il, i) + ment(il, i, j)
	1727	ment(il, i, j) = ment(il, i, j)*sig(il, j)
	1728	bsum(il, i) = bsum(il, i) + ment(il, i, j)
	1729	END IF
	1730	END DO
	1731	END DO
[879]	1732
[1992]	1733	DO il = 1, ncum
	1734	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	1735	bsum(il, i) = amax1(bsum(il,i), 1.0E-16)
	1736	bsum(il, i) = 1.0/bsum(il, i)
	1737	END IF
	1738	END DO
[879]	1739
[1992]	1740	DO j = minorig, nl
	1741	DO il = 1, ncum
	1742	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	1743	il)-1) .AND. j<=inb(il)) THEN
	1744	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
	1745	END IF
	1746	END DO
	1747	END DO
[879]	1748
[1992]	1749	DO j = minorig, nl
	1750	DO il = 1, ncum
	1751	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
	1752	il)-1) .AND. j<=inb(il)) THEN
	1753	csum(il, i) = csum(il, i) + ment(il, i, j)
	1754	END IF
	1755	END DO
	1756	END DO
[879]	1757
[1992]	1758	DO il = 1, ncum
	1759	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	1760	csum(il,i)<m(il,i)) THEN
	1761	nent(il, i) = 0
	1762	ment(il, i, i) = m(il, i)
	1763	qent(il, i, i) = rr(il, 1) - ep(il, i)*clw(il, i)
	1764	uent(il, i, i) = u(il, nk(il))
	1765	vent(il, i, i) = v(il, nk(il))
	1766	elij(il, i, i) = clw(il, i)
	1767	! MAF sij(il,i,i)=1.0
	1768	sij(il, i, i) = 0.0
	1769	END IF
	1770	END DO ! il
[879]	1771
[1992]	1772	! do j=1,ntra
	1773	! do il=1,ncum
	1774	! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	1775	! : .and. csum(il,i).lt.m(il,i) ) then
	1776	! traent(il,i,i,j)=tra(il,nk(il),j)
	1777	! endif
	1778	! enddo
	1779	! enddo
	1780	789 END DO
[879]	1781
[1992]	1782	! MAF: renormalisation de MENT
	1783	DO jm = 1, nd
	1784	DO im = 1, nd
	1785	DO il = 1, ncum
	1786	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
	1787	END DO
	1788	END DO
	1789	END DO
[879]	1790
[1992]	1791	DO jm = 1, nd
	1792	DO im = 1, nd
	1793	DO il = 1, ncum
	1794	IF (zm(il,im)/=0.) THEN
	1795	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
	1796	END IF
	1797	END DO
	1798	END DO
	1799	END DO
[879]	1800
[1992]	1801	DO jm = 1, nd
	1802	DO im = 1, nd
	1803	DO il = 1, ncum
	1804	qents(il, im, jm) = qent(il, im, jm)
	1805	ments(il, im, jm) = ment(il, im, jm)
	1806	END DO
	1807	END DO
	1808	END DO
[879]	1809
[1992]	1810	RETURN
	1811	END SUBROUTINE cv30_mixing
[879]	1812
	1813
[1992]	1814	SUBROUTINE cv30_unsat(nloc, ncum, nd, na, ntra, icb, inb, t, rr, rs, gz, u, &
	1815	v, tra, p, ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, &
	1816	mp, rp, up, vp, trap, wt, water, evap, b & ! RomP-jyg
	1817	, wdtraina, wdtrainm) ! 26/08/10 RomP-jyg
	1818	IMPLICIT NONE
[879]	1819
	1820
[1992]	1821	include "cvthermo.h"
	1822	include "cv30param.h"
	1823	include "cvflag.h"
[879]	1824
[1992]	1825	! inputs:
	1826	INTEGER ncum, nd, na, ntra, nloc
	1827	INTEGER icb(nloc), inb(nloc)
	1828	REAL delt, plcl(nloc)
	1829	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
	1830	REAL u(nloc, nd), v(nloc, nd)
	1831	REAL tra(nloc, nd, ntra)
	1832	REAL p(nloc, nd), ph(nloc, nd+1)
	1833	REAL th(nloc, na), gz(nloc, na)
	1834	REAL lv(nloc, na), ep(nloc, na), sigp(nloc, na), clw(nloc, na)
	1835	REAL cpn(nloc, na), tv(nloc, na)
	1836	REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
[879]	1837
[1992]	1838	! outputs:
	1839	REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
	1840	REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
	1841	REAL trap(nloc, na, ntra)
	1842	REAL b(nloc, na)
	1843	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
	1844	! lascendance adiabatique et des flux melanges Pa et Pm.
	1845	! Distinction des wdtrain
	1846	! Pa = wdtrainA Pm = wdtrainM
	1847	REAL wdtraina(nloc, na), wdtrainm(nloc, na)
[879]	1848
[1992]	1849	! local variables
	1850	INTEGER i, j, k, il, num1
	1851	REAL tinv, delti
	1852	REAL awat, afac, afac1, afac2, bfac
	1853	REAL pr1, pr2, sigt, b6, c6, revap, tevap, delth
	1854	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
	1855	REAL ampmax
	1856	REAL lvcp(nloc, na)
	1857	REAL wdtrain(nloc)
	1858	LOGICAL lwork(nloc)
[879]	1859
	1860
[1992]	1861	! ------------------------------------------------------
[879]	1862
[1992]	1863	delti = 1./delt
	1864	tinv = 1./3.
[879]	1865
[1992]	1866	mp(:, :) = 0.
[879]	1867
[1992]	1868	DO i = 1, nl
	1869	DO il = 1, ncum
	1870	mp(il, i) = 0.0
	1871	rp(il, i) = rr(il, i)
	1872	up(il, i) = u(il, i)
	1873	vp(il, i) = v(il, i)
	1874	wt(il, i) = 0.001
	1875	water(il, i) = 0.0
	1876	evap(il, i) = 0.0
	1877	b(il, i) = 0.0
	1878	lvcp(il, i) = lv(il, i)/cpn(il, i)
	1879	END DO
	1880	END DO
[879]	1881
[1992]	1882	! do k=1,ntra
	1883	! do i=1,nd
	1884	! do il=1,ncum
	1885	! trap(il,i,k)=tra(il,i,k)
	1886	! enddo
	1887	! enddo
	1888	! enddo
	1889	! ! RomP >>>
	1890	DO i = 1, nd
	1891	DO il = 1, ncum
	1892	wdtraina(il, i) = 0.0
	1893	wdtrainm(il, i) = 0.0
	1894	END DO
	1895	END DO
	1896	! ! RomP <<<
[879]	1897
[1992]	1898	! * check whether ep(inb)=0, if so, skip precipitating *
	1899	! * downdraft calculation *
[879]	1900
	1901
[1992]	1902	DO il = 1, ncum
	1903	lwork(il) = .TRUE.
	1904	IF (ep(il,inb(il))<0.0001) lwork(il) = .FALSE.
	1905	END DO
[879]	1906
[1992]	1907	CALL zilch(wdtrain, ncum)
[879]	1908
[1992]	1909	DO i = nl + 1, 1, -1
[879]	1910
[1992]	1911	num1 = 0
	1912	DO il = 1, ncum
	1913	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
	1914	END DO
	1915	IF (num1<=0) GO TO 400
[879]	1916
	1917
[1992]	1918	! * integrate liquid water equation to find condensed water *
	1919	! * and condensed water flux *
[879]	1920
	1921
	1922
[1992]	1923	! * begin downdraft loop *
[879]	1924
	1925
	1926
[1992]	1927	! * calculate detrained precipitation *
[879]	1928
[1992]	1929	DO il = 1, ncum
	1930	IF (i<=inb(il) .AND. lwork(il)) THEN
	1931	IF (cvflag_grav) THEN
	1932	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
	1933	wdtraina(il, i) = wdtrain(il)/grav ! Pa 26/08/10 RomP
	1934	ELSE
	1935	wdtrain(il) = 10.0ep(il, i)m(il, i)*clw(il, i)
	1936	wdtraina(il, i) = wdtrain(il)/10. ! Pa 26/08/10 RomP
	1937	END IF
	1938	END IF
	1939	END DO
[879]	1940
[1992]	1941	IF (i>1) THEN
[879]	1942
[1992]	1943	DO j = 1, i - 1
	1944	DO il = 1, ncum
	1945	IF (i<=inb(il) .AND. lwork(il)) THEN
	1946	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
	1947	awat = amax1(awat, 0.0)
	1948	IF (cvflag_grav) THEN
	1949	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
	1950	ELSE
	1951	wdtrain(il) = wdtrain(il) + 10.0awatment(il, j, i)
	1952	END IF
	1953	END IF
	1954	END DO
	1955	END DO
	1956	DO il = 1, ncum
	1957	IF (cvflag_grav) THEN
	1958	wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) ! Pm 26/08/10 RomP
	1959	ELSE
	1960	wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) ! Pm 26/08/10 RomP
	1961	END IF
	1962	END DO
[879]	1963
[1992]	1964	END IF
[879]	1965
	1966
[1992]	1967	! * find rain water and evaporation using provisional *
	1968	! * estimates of rp(i)and rp(i-1) *
[879]	1969
	1970
[1992]	1971	DO il = 1, ncum
[879]	1972
[1992]	1973	IF (i<=inb(il) .AND. lwork(il)) THEN
[879]	1974
[1992]	1975	wt(il, i) = 45.0
[879]	1976
[1992]	1977	IF (i<inb(il)) THEN
	1978	rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, &
	1979	i))+gz(il,i+1)-gz(il,i))/lv(il, i)
	1980	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
	1981	END IF
	1982	rp(il, i) = amax1(rp(il,i), 0.0)
	1983	rp(il, i) = amin1(rp(il,i), rs(il,i))
	1984	rp(il, inb(il)) = rr(il, inb(il))
[879]	1985
[1992]	1986	IF (i==1) THEN
	1987	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
	1988	ELSE
	1989	rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, &
	1990	i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
	1991	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
	1992	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
	1993	rp(il, i-1) = amax1(rp(il,i-1), 0.0)
	1994	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i) &
	1995	)
	1996	afac2 = p(il, i-1)*(rs(il,i-1)-rp(il,i-1))/ &
	1997	(1.0E4+2000.0p(il,i-1)rs(il,i-1))
	1998	afac = 0.5*(afac1+afac2)
	1999	END IF
	2000	IF (i==inb(il)) afac = 0.0
	2001	afac = amax1(afac, 0.0)
	2002	bfac = 1./(sigd*wt(il,i))
[879]	2003
[1992]	2004	! jyg1
	2005	! cc sigt=1.0
	2006	! cc if(i.ge.icb)sigt=sigp(i)
	2007	! prise en compte de la variation progressive de sigt dans
	2008	! les couches icb et icb-1:
	2009	! pour plcl<ph(i+1), pr1=0 & pr2=1
	2010	! pour plcl>ph(i), pr1=1 & pr2=0
	2011	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
	2012	! sur le nuage, et pr2 est la proportion sous la base du
	2013	! nuage.
	2014	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
	2015	pr1 = max(0., min(1.,pr1))
	2016	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
	2017	pr2 = max(0., min(1.,pr2))
	2018	sigt = sigp(il, i)*pr1 + pr2
	2019	! jyg2
[879]	2020
[1992]	2021	b6 = bfac50.sigd(ph(il,i)-ph(il,i+1))sigt*afac
	2022	c6 = water(il, i+1) + bfacwdtrain(il) - 50.sigdbfac(ph(il,i)-ph( &
	2023	il,i+1))*evap(il, i+1)
	2024	IF (c6>0.0) THEN
	2025	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
	2026	evap(il, i) = sigtafacrevap
	2027	water(il, i) = revap*revap
	2028	ELSE
	2029	evap(il, i) = -evap(il, i+1) + 0.02(wdtrain(il)+sigdwt(il,i)* &
	2030	water(il,i+1))/(sigd*(ph(il,i)-ph(il,i+1)))
	2031	END IF
[879]	2032
[1992]	2033	! * calculate precipitating downdraft mass flux under *
	2034	! * hydrostatic approximation *
[879]	2035
[1992]	2036	IF (i/=1) THEN
[879]	2037
[1992]	2038	tevap = amax1(0.0, evap(il,i))
	2039	delth = amax1(0.001, (th(il,i)-th(il,i-1)))
	2040	IF (cvflag_grav) THEN
	2041	mp(il, i) = 100.ginvlvcp(il, i)sigdtevap*(p(il,i-1)-p(il,i))/ &
	2042	delth
	2043	ELSE
	2044	mp(il, i) = 10.lvcp(il, i)sigdtevap(p(il,i-1)-p(il,i))/delth
	2045	END IF
[879]	2046
[1992]	2047	! * if hydrostatic assumption fails, *
	2048	! * solve cubic difference equation for downdraft theta *
	2049	! * and mass flux from two simultaneous differential eqns *
[879]	2050
[1992]	2051	amfac = sigdsigd70.0ph(il, i)(p(il,i-1)-p(il,i))* &
	2052	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
	2053	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
	2054	IF (amp2>(0.1*amfac)) THEN
	2055	xf = 100.0sigdsigdsigd(ph(il,i)-ph(il,i+1))
	2056	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i)/(lvcp(il,i)* &
	2057	sigd*th(il,i))
	2058	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
	2059	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
	2060	50.(p(il,i-1)-p(il,i))xf*tevap
	2061	fac2 = 1.0
	2062	IF (bf<0.0) fac2 = -1.0
	2063	bf = abs(bf)
	2064	ur = 0.25bfbf - afafaftinvtinv*tinv
	2065	IF (ur>=0.0) THEN
	2066	sru = sqrt(ur)
	2067	fac = 1.0
	2068	IF ((0.5*bf-sru)<0.0) fac = -1.0
	2069	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
	2070	fac(abs(0.5bf-sru))**tinv
	2071	ELSE
	2072	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
	2073	IF (fac2<0.0) d = 3.14159 - d
	2074	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
	2075	END IF
	2076	mp(il, i) = amax1(0.0, mp(il,i))
[879]	2077
[1992]	2078	IF (cvflag_grav) THEN
	2079	! jyg : il y a vraisemblablement une erreur dans la ligne 2
	2080	! suivante:
	2081	! il faut diviser par (mp(il,i)sigdgrav) et non par
	2082	! (mp(il,i)+sigd*0.1).
	2083	! Et il faut bien revoir les facteurs 100.
	2084	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap/(mp(il, &
	2085	i)+sigd0.1) - 10.0(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
	2086	)sigdth(il,i))
	2087	ELSE
	2088	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap/(mp(il, &
	2089	i)+sigd0.1) - 10.0(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
	2090	)sigdth(il,i))
	2091	END IF
	2092	b(il, i-1) = amax1(b(il,i-1), 0.0)
	2093	END IF
[879]	2094
[1992]	2095	! *** limit magnitude of mp(i) to meet cfl condition
	2096	! ***
[879]	2097
[1992]	2098	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
	2099	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
	2100	ampmax = amin1(ampmax, amp2)
	2101	mp(il, i) = amin1(mp(il,i), ampmax)
[879]	2102
[1992]	2103	! *** force mp to decrease linearly to zero
	2104	! ***
	2105	! *** between cloud base and the surface
	2106	! ***
[879]	2107
[1992]	2108	IF (p(il,i)>p(il,icb(il))) THEN
	2109	mp(il, i) = mp(il, icb(il))*(p(il,1)-p(il,i))/ &
	2110	(p(il,1)-p(il,icb(il)))
	2111	END IF
[879]	2112
[1992]	2113	END IF ! i.eq.1
[879]	2114
[1992]	2115	! * find mixing ratio of precipitating downdraft *
[879]	2116
	2117
[1992]	2118	IF (i/=inb(il)) THEN
[879]	2119
[1992]	2120	rp(il, i) = rr(il, i)
[879]	2121
[1992]	2122	IF (mp(il,i)>mp(il,i+1)) THEN
[879]	2123
[1992]	2124	IF (cvflag_grav) THEN
	2125	rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
	2126	rr(il, i)(mp(il,i)-mp(il,i+1)) + 100.ginv0.5sigd*(ph(il,i &
	2127	)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
	2128	ELSE
	2129	rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
	2130	rr(il, i)(mp(il,i)-mp(il,i+1)) + 5.sigd*(ph(il,i)-ph(il,i+1 &
	2131	))*(evap(il,i+1)+evap(il,i))
	2132	END IF
	2133	rp(il, i) = rp(il, i)/mp(il, i)
	2134	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+ &
	2135	1))
	2136	up(il, i) = up(il, i)/mp(il, i)
	2137	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+ &
	2138	1))
	2139	vp(il, i) = vp(il, i)/mp(il, i)
[879]	2140
[1992]	2141	! do j=1,ntra
	2142	! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
	2143	! testmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
	2144	! : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
	2145	! trap(il,i,j)=trap(il,i,j)/mp(il,i)
	2146	! end do
[879]	2147
[1992]	2148	ELSE
[1742]	2149
[1992]	2150	IF (mp(il,i+1)>1.0E-16) THEN
	2151	IF (cvflag_grav) THEN
	2152	rp(il, i) = rp(il, i+1) + 100.ginv0.5sigd(ph(il,i)-ph(il, &
	2153	i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1)
	2154	ELSE
	2155	rp(il, i) = rp(il, i+1) + 5.sigd(ph(il,i)-ph(il,i+1))*(evap &
	2156	(il,i+1)+evap(il,i))/mp(il, i+1)
	2157	END IF
	2158	up(il, i) = up(il, i+1)
	2159	vp(il, i) = vp(il, i+1)
[1742]	2160
[1992]	2161	! do j=1,ntra
	2162	! trap(il,i,j)=trap(il,i+1,j)
	2163	! end do
[1742]	2164
[1992]	2165	END IF
	2166	END IF
	2167	rp(il, i) = amin1(rp(il,i), rs(il,i))
	2168	rp(il, i) = amax1(rp(il,i), 0.0)
[1742]	2169
[1992]	2170	END IF
	2171	END IF
	2172	END DO
[1742]	2173
[1992]	2174	400 END DO
[879]	2175
[1992]	2176	RETURN
	2177	END SUBROUTINE cv30_unsat
[879]	2178
[1992]	2179	SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, &
	2180	tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, &
	2181	wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, &
	2182	tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, &
	2183	mike, tls, tps, qcondc, wd)
	2184	IMPLICIT NONE
[879]	2185
[1992]	2186	include "cvthermo.h"
	2187	include "cv30param.h"
	2188	include "cvflag.h"
	2189	include "conema3.h"
[879]	2190
[1992]	2191	! inputs:
	2192	INTEGER ncum, nd, na, ntra, nloc
	2193	INTEGER icb(nloc), inb(nloc)
	2194	REAL delt
	2195	REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
	2196	REAL tra(nloc, nd, ntra), sig(nloc, nd)
	2197	REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
	2198	REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
	2199	REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
	2200	REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
	2201	REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
	2202	REAL water(nloc, na), evap(nloc, na), b(nloc, na)
	2203	REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
	2204	! ym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
	2205	REAL vent(nloc, na, na), elij(nloc, na, na)
	2206	INTEGER nent(nloc, na)
	2207	REAL traent(nloc, na, na, ntra)
	2208	REAL tv(nloc, nd), tvp(nloc, nd)
[879]	2209
[1992]	2210	! input/output:
	2211	INTEGER iflag(nloc)
[879]	2212
[1992]	2213	! outputs:
	2214	REAL precip(nloc)
	2215	REAL vprecip(nloc, nd+1)
	2216	REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	2217	REAL ftra(nloc, nd, ntra)
	2218	REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
	2219	REAL dnwd0(nloc, nd), mike(nloc, nd)
	2220	REAL tls(nloc, nd), tps(nloc, nd)
	2221	REAL qcondc(nloc, nd) ! cld
	2222	REAL wd(nloc) ! gust
[879]	2223
[1992]	2224	! local variables:
	2225	INTEGER i, k, il, n, j, num1
	2226	REAL rat, awat, delti
	2227	REAL ax, bx, cx, dx, ex
	2228	REAL cpinv, rdcp, dpinv
	2229	REAL lvcp(nloc, na), mke(nloc, na)
	2230	REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
	2231	! !! real up1(nloc), dn1(nloc)
	2232	REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
	2233	REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
	2234	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
	2235	REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
[879]	2236
	2237
[1992]	2238	! -------------------------------------------------------------
[879]	2239
[1992]	2240	! initialization:
[879]	2241
[1992]	2242	delti = 1.0/delt
[879]	2243
[1992]	2244	DO il = 1, ncum
	2245	precip(il) = 0.0
	2246	wd(il) = 0.0 ! gust
	2247	vprecip(il, nd+1) = 0.
	2248	END DO
	2249
	2250	DO i = 1, nd
	2251	DO il = 1, ncum
	2252	vprecip(il, i) = 0.0
	2253	ft(il, i) = 0.0
	2254	fr(il, i) = 0.0
	2255	fu(il, i) = 0.0
	2256	fv(il, i) = 0.0
	2257	qcondc(il, i) = 0.0 ! cld
	2258	qcond(il, i) = 0.0 ! cld
	2259	nqcond(il, i) = 0.0 ! cld
	2260	END DO
	2261	END DO
	2262
	2263	! do j=1,ntra
	2264	! do i=1,nd
	2265	! do il=1,ncum
	2266	! ftra(il,i,j)=0.0
	2267	! enddo
	2268	! enddo
	2269	! enddo
	2270
	2271	DO i = 1, nl
	2272	DO il = 1, ncum
	2273	lvcp(il, i) = lv(il, i)/cpn(il, i)
	2274	END DO
	2275	END DO
	2276
	2277
	2278
	2279	! * calculate surface precipitation in mm/day *
	2280
	2281	DO il = 1, ncum
	2282	IF (ep(il,inb(il))>=0.0001) THEN
	2283	IF (cvflag_grav) THEN
	2284	precip(il) = wt(il, 1)sigdwater(il, 1)86400.1000./(rowl*grav)
	2285	ELSE
	2286	precip(il) = wt(il, 1)sigdwater(il, 1)*8640.
	2287	END IF
	2288	END IF
	2289	END DO
	2290
	2291	! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
	2292
	2293	! MAF rajout pour lessivage
	2294	DO k = 1, nl
	2295	DO il = 1, ncum
	2296	IF (k<=inb(il)) THEN
	2297	IF (cvflag_grav) THEN
	2298	vprecip(il, k) = wt(il, k)sigdwater(il, k)/grav
	2299	ELSE
	2300	vprecip(il, k) = wt(il, k)sigdwater(il, k)/10.
	2301	END IF
	2302	END IF
	2303	END DO
	2304	END DO
	2305
	2306
	2307	! * Calculate downdraft velocity scale *
	2308	! * NE PAS UTILISER POUR L'INSTANT *
	2309
	2310	! ! do il=1,ncum
	2311	! ! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
	2312	! ! : /(sigd*p(il,icb(il)))
	2313	! ! enddo
	2314
	2315
	2316	! * calculate tendencies of lowest level potential temperature *
	2317	! * and mixing ratio *
	2318
	2319	DO il = 1, ncum
	2320	work(il) = 1.0/(ph(il,1)-ph(il,2))
	2321	am(il) = 0.0
	2322	END DO
	2323
	2324	DO k = 2, nl
	2325	DO il = 1, ncum
	2326	IF (k<=inb(il)) THEN
	2327	am(il) = am(il) + m(il, k)
	2328	END IF
	2329	END DO
	2330	END DO
	2331
	2332	DO il = 1, ncum
	2333
	2334	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
	2335	IF (cvflag_grav) THEN
	2336	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
	2337	ft(il, 1) = 0.01gravwork(il)am(il)(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
	2338	1))/cpn(il,1))
	2339	ELSE
	2340	IF ((0.1work(il)am(il))>=delti) iflag(il) = 1 !consistency vect
	2341	ft(il, 1) = 0.1work(il)am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
	2342	1))/cpn(il,1))
	2343	END IF
	2344
	2345	ft(il, 1) = ft(il, 1) - 0.5lvcp(il, 1)sigd*(evap(il,1)+evap(il,2))
	2346
	2347	IF (cvflag_grav) THEN
	2348	ft(il, 1) = ft(il, 1) - 0.009gravsigdmp(il, 2)t(il, 1)b(il, 1) &
	2349	work(il)
	2350	ELSE
	2351	ft(il, 1) = ft(il, 1) - 0.09sigdmp(il, 2)t(il, 1)b(il, 1)*work(il)
	2352	END IF
	2353
	2354	ft(il, 1) = ft(il, 1) + 0.01sigdwt(il, 1)(cl-cpd)water(il, 2)*(t(il,2 &
	2355	)-t(il,1))*work(il)/cpn(il, 1)
	2356
	2357	IF (cvflag_grav) THEN
	2358	! jyg1 Correction pour mieux conserver l'eau (conformite avec
	2359	! CONVECT4.3)
	2360	! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas
	2361	! evap)
	2362	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr(il,1))work(il) + &
	2363	sigd0.5(evap(il,1)+evap(il,2))
	2364	! +tard : +sigd*evap(il,1)
	2365
	2366	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il)
	2367
	2368	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il, &
	2369	1))+am(il)*(u(il,2)-u(il,1)))
	2370	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il, &
	2371	1))+am(il)*(v(il,2)-v(il,1)))
	2372	ELSE ! cvflag_grav
	2373	fr(il, 1) = 0.1mp(il, 2)(rp(il,2)-rr(il,1))*work(il) + &
	2374	sigd0.5(evap(il,1)+evap(il,2))
	2375	fr(il, 1) = fr(il, 1) + 0.1am(il)(rr(il,2)-rr(il,1))*work(il)
	2376	fu(il, 1) = fu(il, 1) + 0.1work(il)(mp(il,2)*(up(il,2)-u(il, &
	2377	1))+am(il)*(u(il,2)-u(il,1)))
	2378	fv(il, 1) = fv(il, 1) + 0.1work(il)(mp(il,2)*(vp(il,2)-v(il, &
	2379	1))+am(il)*(v(il,2)-v(il,1)))
	2380	END IF ! cvflag_grav
	2381
	2382	END DO ! il
	2383
	2384	! do j=1,ntra
	2385	! do il=1,ncum
	2386	! if (cvflag_grav) then
	2387	! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
	2388	! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2389	! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2390	! else
	2391	! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
	2392	! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2393	! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2394	! endif
	2395	! enddo
	2396	! enddo
	2397
	2398	DO j = 2, nl
	2399	DO il = 1, ncum
	2400	IF (j<=inb(il)) THEN
	2401	IF (cvflag_grav) THEN
	2402	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il, &
	2403	j,1)-rr(il,1))
	2404	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il, &
	2405	j,1)-u(il,1))
	2406	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il, &
	2407	j,1)-v(il,1))
	2408	ELSE ! cvflag_grav
	2409	fr(il, 1) = fr(il, 1) + 0.1work(il)ment(il, j, 1)*(qent(il,j,1)- &
	2410	rr(il,1))
	2411	fu(il, 1) = fu(il, 1) + 0.1work(il)ment(il, j, 1)*(uent(il,j,1)-u &
	2412	(il,1))
	2413	fv(il, 1) = fv(il, 1) + 0.1work(il)ment(il, j, 1)*(vent(il,j,1)-v &
	2414	(il,1))
	2415	END IF ! cvflag_grav
	2416	END IF ! j
	2417	END DO
	2418	END DO
	2419
	2420	! do k=1,ntra
	2421	! do j=2,nl
	2422	! do il=1,ncum
	2423	! if (j.le.inb(il)) then
	2424
	2425	! if (cvflag_grav) then
	2426	! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
	2427	! : *(traent(il,j,1,k)-tra(il,1,k))
	2428	! else
	2429	! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
	2430	! : *(traent(il,j,1,k)-tra(il,1,k))
	2431	! endif
	2432
	2433	! endif
	2434	! enddo
	2435	! enddo
	2436	! enddo
	2437
	2438
	2439	! * calculate tendencies of potential temperature and mixing ratio *
	2440	! * at levels above the lowest level *
	2441
	2442	! * first find the net saturated updraft and downdraft mass fluxes *
	2443	! * through each level *
	2444
	2445
	2446	DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
	2447
	2448	num1 = 0
	2449	DO il = 1, ncum
	2450	IF (i<=inb(il)) num1 = num1 + 1
	2451	END DO
	2452	IF (num1<=0) GO TO 500
	2453
	2454	CALL zilch(amp1, ncum)
	2455	CALL zilch(ad, ncum)
	2456
	2457	DO k = i + 1, nl + 1
	2458	DO il = 1, ncum
	2459	IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN
	2460	amp1(il) = amp1(il) + m(il, k)
	2461	END IF
	2462	END DO
	2463	END DO
	2464
	2465	DO k = 1, i
	2466	DO j = i + 1, nl + 1
	2467	DO il = 1, ncum
	2468	IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
	2469	amp1(il) = amp1(il) + ment(il, k, j)
	2470	END IF
	2471	END DO
	2472	END DO
	2473	END DO
	2474
	2475	DO k = 1, i - 1
	2476	DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
	2477	DO il = 1, ncum
	2478	IF (i<=inb(il) .AND. j<=inb(il)) THEN
	2479	ad(il) = ad(il) + ment(il, j, k)
	2480	END IF
	2481	END DO
	2482	END DO
	2483	END DO
	2484
	2485	DO il = 1, ncum
	2486	IF (i<=inb(il)) THEN
	2487	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	2488	cpinv = 1.0/cpn(il, i)
	2489
	2490	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
	2491	IF (cvflag_grav) THEN
	2492	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
	2493	ELSE
	2494	IF ((0.1dpinvamp1(il))>=delti) iflag(il) = 1 ! vecto
	2495	END IF
	2496
	2497	IF (cvflag_grav) THEN
	2498	ft(il, i) = 0.01gravdpinv(amp1(il)(t(il,i+1)-t(il, &
	2499	i)+(gz(il,i+1)-gz(il,i))cpinv)-ad(il)(t(il,i)-t(il, &
	2500	i-1)+(gz(il,i)-gz(il,i-1))cpinv)) - 0.5sigdlvcp(il, i)(evap( &
	2501	il,i)+evap(il,i+1))
	2502	rat = cpn(il, i-1)*cpinv
	2503	ft(il, i) = ft(il, i) - 0.009gravsigd(mp(il,i+1)t(il,i)*b(il,i) &
	2504	-mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
	2505	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h( &
	2506	il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
	2507	ELSE ! cvflag_grav
	2508	ft(il, i) = 0.1dpinv(amp1(il)*(t(il,i+1)-t(il, &
	2509	i)+(gz(il,i+1)-gz(il,i))cpinv)-ad(il)(t(il,i)-t(il, &
	2510	i-1)+(gz(il,i)-gz(il,i-1))cpinv)) - 0.5sigdlvcp(il, i)(evap( &
	2511	il,i)+evap(il,i+1))
	2512	rat = cpn(il, i-1)*cpinv
	2513	ft(il, i) = ft(il, i) - 0.09sigd(mp(il,i+1)t(il,i)b(il,i)-mp(il &
	2514	,i)t(il,i-1)ratb(il,i-1))dpinv
	2515	ft(il, i) = ft(il, i) + 0.1dpinvment(il, i, i)*(hp(il,i)-h(il,i)+ &
	2516	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
	2517	END IF ! cvflag_grav
	2518
	2519
	2520	ft(il, i) = ft(il, i) + 0.01sigdwt(il, i)(cl-cpd)water(il, i+1)*( &
	2521	t(il,i+1)-t(il,i))dpinvcpinv
	2522
	2523	IF (cvflag_grav) THEN
	2524	fr(il, i) = 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il, &
	2525	i))-ad(il)*(rr(il,i)-rr(il,i-1)))
	2526	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il, &
	2527	i))-ad(il)*(u(il,i)-u(il,i-1)))
	2528	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il, &
	2529	i))-ad(il)*(v(il,i)-v(il,i-1)))
	2530	ELSE ! cvflag_grav
	2531	fr(il, i) = 0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il, &
	2532	i))-ad(il)*(rr(il,i)-rr(il,i-1)))
	2533	fu(il, i) = fu(il, i) + 0.1dpinv(amp1(il)*(u(il,i+1)-u(il, &
	2534	i))-ad(il)*(u(il,i)-u(il,i-1)))
	2535	fv(il, i) = fv(il, i) + 0.1dpinv(amp1(il)*(v(il,i+1)-v(il, &
	2536	i))-ad(il)*(v(il,i)-v(il,i-1)))
	2537	END IF ! cvflag_grav
	2538
	2539	END IF ! i
	2540	END DO
	2541
	2542	! do k=1,ntra
	2543	! do il=1,ncum
	2544	! if (i.le.inb(il)) then
	2545	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2546	! cpinv=1.0/cpn(il,i)
	2547	! if (cvflag_grav) then
	2548	! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
	2549	! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2550	! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2551	! else
	2552	! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
	2553	! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2554	! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2555	! endif
	2556	! endif
	2557	! enddo
	2558	! enddo
	2559
	2560	DO k = 1, i - 1
	2561	DO il = 1, ncum
	2562	IF (i<=inb(il)) THEN
	2563	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	2564	cpinv = 1.0/cpn(il, i)
	2565
	2566	awat = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
	2567	awat = amax1(awat, 0.0)
	2568
	2569	IF (cvflag_grav) THEN
	2570	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k &
	2571	,i)-awat-rr(il,i))
	2572	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
	2573	,i)-u(il,i))
	2574	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k &
	2575	,i)-v(il,i))
	2576	ELSE ! cvflag_grav
	2577	fr(il, i) = fr(il, i) + 0.1dpinvment(il, k, i)*(qent(il,k,i)- &
	2578	awat-rr(il,i))
	2579	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
	2580	,i)-u(il,i))
	2581	fv(il, i) = fv(il, i) + 0.1dpinvment(il, k, i)*(vent(il,k,i)-v( &
	2582	il,i))
	2583	END IF ! cvflag_grav
	2584
	2585	! (saturated updrafts resulting from mixing) ! cld
	2586	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat) ! cld
	2587	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2588	END IF ! i
	2589	END DO
	2590	END DO
	2591
	2592	! do j=1,ntra
	2593	! do k=1,i-1
	2594	! do il=1,ncum
	2595	! if (i.le.inb(il)) then
	2596	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2597	! cpinv=1.0/cpn(il,i)
	2598	! if (cvflag_grav) then
	2599	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2600	! : *(traent(il,k,i,j)-tra(il,i,j))
	2601	! else
	2602	! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2603	! : *(traent(il,k,i,j)-tra(il,i,j))
	2604	! endif
	2605	! endif
	2606	! enddo
	2607	! enddo
	2608	! enddo
	2609
	2610	DO k = i, nl + 1
	2611	DO il = 1, ncum
	2612	IF (i<=inb(il) .AND. k<=inb(il)) THEN
	2613	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	2614	cpinv = 1.0/cpn(il, i)
	2615
	2616	IF (cvflag_grav) THEN
	2617	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k &
	2618	,i)-rr(il,i))
	2619	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
	2620	,i)-u(il,i))
	2621	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k &
	2622	,i)-v(il,i))
	2623	ELSE ! cvflag_grav
	2624	fr(il, i) = fr(il, i) + 0.1dpinvment(il, k, i)*(qent(il,k,i)-rr &
	2625	(il,i))
	2626	fu(il, i) = fu(il, i) + 0.1dpinvment(il, k, i)*(uent(il,k,i)-u( &
	2627	il,i))
	2628	fv(il, i) = fv(il, i) + 0.1dpinvment(il, k, i)*(vent(il,k,i)-v( &
	2629	il,i))
	2630	END IF ! cvflag_grav
	2631	END IF ! i and k
	2632	END DO
	2633	END DO
	2634
	2635	! do j=1,ntra
	2636	! do k=i,nl+1
	2637	! do il=1,ncum
	2638	! if (i.le.inb(il) .and. k.le.inb(il)) then
	2639	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2640	! cpinv=1.0/cpn(il,i)
	2641	! if (cvflag_grav) then
	2642	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2643	! : *(traent(il,k,i,j)-tra(il,i,j))
	2644	! else
	2645	! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2646	! : *(traent(il,k,i,j)-tra(il,i,j))
	2647	! endif
	2648	! endif ! i and k
	2649	! enddo
	2650	! enddo
	2651	! enddo
	2652
	2653	DO il = 1, ncum
	2654	IF (i<=inb(il)) THEN
	2655	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	2656	cpinv = 1.0/cpn(il, i)
	2657
	2658	IF (cvflag_grav) THEN
	2659	! sb: on ne fait pas encore la correction permettant de mieux
	2660	! conserver l'eau:
	2661	fr(il, i) = fr(il, i) + 0.5sigd(evap(il,i)+evap(il,i+1)) + &
	2662	0.01grav(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)(rp(il, &
	2663	i)-rr(il,i-1)))*dpinv
	2664
	2665	fu(il, i) = fu(il, i) + 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il, &
	2666	i))-mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	2667	fv(il, i) = fv(il, i) + 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il, &
	2668	i))-mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
	2669	ELSE ! cvflag_grav
	2670	fr(il, i) = fr(il, i) + 0.5sigd(evap(il,i)+evap(il,i+1)) + &
	2671	0.1(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il,i)-rr(il, &
	2672	i-1)))*dpinv
	2673	fu(il, i) = fu(il, i) + 0.1(mp(il,i+1)(up(il,i+1)-u(il, &
	2674	i))-mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	2675	fv(il, i) = fv(il, i) + 0.1(mp(il,i+1)(vp(il,i+1)-v(il, &
	2676	i))-mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
	2677	END IF ! cvflag_grav
	2678
	2679	END IF ! i
	2680	END DO
	2681
	2682	! sb: interface with the cloud parameterization: ! cld
	2683
	2684	DO k = i + 1, nl
	2685	DO il = 1, ncum
	2686	IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld
	2687	! (saturated downdrafts resulting from mixing) ! cld
	2688	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
	2689	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2690	END IF ! cld
	2691	END DO ! cld
	2692	END DO ! cld
	2693
	2694	! (particular case: no detraining level is found) ! cld
	2695	DO il = 1, ncum ! cld
	2696	IF (i<=inb(il) .AND. nent(il,i)==0) THEN ! cld
	2697	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
	2698	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	2699	END IF ! cld
	2700	END DO ! cld
	2701
	2702	DO il = 1, ncum ! cld
	2703	IF (i<=inb(il) .AND. nqcond(il,i)/=0.) THEN ! cld
	2704	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
	2705	END IF ! cld
	2706	END DO
	2707
	2708	! do j=1,ntra
	2709	! do il=1,ncum
	2710	! if (i.le.inb(il)) then
	2711	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2712	! cpinv=1.0/cpn(il,i)
	2713
	2714	! if (cvflag_grav) then
	2715	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
	2716	! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2717	! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
	2718	! else
	2719	! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
	2720	! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2721	! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
	2722	! endif
	2723	! endif ! i
	2724	! enddo
	2725	! enddo
	2726
	2727	500 END DO
	2728
	2729
	2730	! * move the detrainment at level inb down to level inb-1 *
	2731	! * in such a way as to preserve the vertically *
	2732	! * integrated enthalpy and water tendencies *
	2733
	2734	DO il = 1, ncum
	2735
	2736	ax = 0.1ment(il, inb(il), inb(il))(hp(il,inb(il))-h(il,inb(il))+t(il, &
	2737	inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il), &
	2738	inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
	2739	ft(il, inb(il)) = ft(il, inb(il)) - ax
	2740	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il &
	2741	))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il, &
	2742	inb(il))))
	2743
	2744	bx = 0.1ment(il, inb(il), inb(il))(qent(il,inb(il),inb(il))-rr(il,inb( &
	2745	il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2746	fr(il, inb(il)) = fr(il, inb(il)) - bx
	2747	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+ &
	2748	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
	2749
	2750	cx = 0.1ment(il, inb(il), inb(il))(uent(il,inb(il),inb(il))-u(il,inb(il &
	2751	)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2752	fu(il, inb(il)) = fu(il, inb(il)) - cx
	2753	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+ &
	2754	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
	2755
	2756	dx = 0.1ment(il, inb(il), inb(il))(vent(il,inb(il),inb(il))-v(il,inb(il &
	2757	)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2758	fv(il, inb(il)) = fv(il, inb(il)) - dx
	2759	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+ &
	2760	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
	2761
	2762	END DO
	2763
	2764	! do j=1,ntra
	2765	! do il=1,ncum
	2766	! ex=0.1*ment(il,inb(il),inb(il))
	2767	! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
	2768	! : /(ph(il,inb(il))-ph(il,inb(il)+1))
	2769	! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
	2770	! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
	2771	! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
	2772	! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2773	! enddo
	2774	! enddo
	2775
	2776
	2777	! * homoginize tendencies below cloud base *
	2778
	2779
	2780	DO il = 1, ncum
	2781	asum(il) = 0.0
	2782	bsum(il) = 0.0
	2783	csum(il) = 0.0
	2784	dsum(il) = 0.0
	2785	END DO
	2786
	2787	DO i = 1, nl
	2788	DO il = 1, ncum
	2789	IF (i<=(icb(il)-1)) THEN
	2790	asum(il) = asum(il) + ft(il, i)*(ph(il,i)-ph(il,i+1))
	2791	bsum(il) = bsum(il) + fr(il, i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il, &
	2792	1)))*(ph(il,i)-ph(il,i+1))
	2793	csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
	2794	1)))*(ph(il,i)-ph(il,i+1))
	2795	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
	2796	END IF
	2797	END DO
	2798	END DO
	2799
	2800	! !!! do 700 i=1,icb(il)-1
	2801	DO i = 1, nl
	2802	DO il = 1, ncum
	2803	IF (i<=(icb(il)-1)) THEN
	2804	ft(il, i) = asum(il)t(il, i)/(th(il,i)dsum(il))
	2805	fr(il, i) = bsum(il)/csum(il)
	2806	END IF
	2807	END DO
	2808	END DO
	2809
	2810
	2811	! * reset counter and return *
	2812
	2813	DO il = 1, ncum
	2814	sig(il, nd) = 2.0
	2815	END DO
	2816
	2817
	2818	DO i = 1, nd
	2819	DO il = 1, ncum
	2820	upwd(il, i) = 0.0
	2821	dnwd(il, i) = 0.0
	2822	END DO
	2823	END DO
	2824
	2825	DO i = 1, nl
	2826	DO il = 1, ncum
	2827	dnwd0(il, i) = -mp(il, i)
	2828	END DO
	2829	END DO
	2830	DO i = nl + 1, nd
	2831	DO il = 1, ncum
	2832	dnwd0(il, i) = 0.
	2833	END DO
	2834	END DO
	2835
	2836
	2837	DO i = 1, nl
	2838	DO il = 1, ncum
	2839	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	2840	upwd(il, i) = 0.0
	2841	dnwd(il, i) = 0.0
	2842	END IF
	2843	END DO
	2844	END DO
	2845
	2846	DO i = 1, nl
	2847	DO k = 1, nl
	2848	DO il = 1, ncum
	2849	up1(il, k, i) = 0.0
	2850	dn1(il, k, i) = 0.0
	2851	END DO
	2852	END DO
	2853	END DO
	2854
	2855	DO i = 1, nl
	2856	DO k = i, nl
	2857	DO n = 1, i - 1
	2858	DO il = 1, ncum
	2859	IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
	2860	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
	2861	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
	2862	END IF
	2863	END DO
	2864	END DO
	2865	END DO
	2866	END DO
	2867
	2868	DO i = 2, nl
	2869	DO k = i, nl
	2870	DO il = 1, ncum
	2871	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
	2872	! then
	2873	IF (i<=inb(il) .AND. k<=inb(il)) THEN
	2874	upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
	2875	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
	2876	END IF
	2877	END DO
	2878	END DO
	2879	END DO
	2880
	2881
	2882	! !!! DO il=1,ncum
	2883	! !!! do i=icb(il),inb(il)
	2884	! !!!
	2885	! !!! upwd(il,i)=0.0
	2886	! !!! dnwd(il,i)=0.0
	2887	! !!! do k=i,inb(il)
	2888	! !!! up1=0.0
	2889	! !!! dn1=0.0
	2890	! !!! do n=1,i-1
	2891	! !!! up1=up1+ment(il,n,k)
	2892	! !!! dn1=dn1-ment(il,k,n)
	2893	! !!! enddo
	2894	! !!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
	2895	! !!! dnwd(il,i)=dnwd(il,i)+dn1
	2896	! !!! enddo
	2897	! !!! enddo
	2898	! !!!
	2899	! !!! ENDDO
	2900
	2901	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2902	! determination de la variation de flux ascendant entre
	2903	! deux niveau non dilue mike
	2904	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2905
	2906	DO i = 1, nl
	2907	DO il = 1, ncum
	2908	mike(il, i) = m(il, i)
	2909	END DO
	2910	END DO
	2911
	2912	DO i = nl + 1, nd
	2913	DO il = 1, ncum
	2914	mike(il, i) = 0.
	2915	END DO
	2916	END DO
	2917
	2918	DO i = 1, nd
	2919	DO il = 1, ncum
	2920	ma(il, i) = 0
	2921	END DO
	2922	END DO
	2923
	2924	DO i = 1, nl
	2925	DO j = i, nl
	2926	DO il = 1, ncum
	2927	ma(il, i) = ma(il, i) + m(il, j)
	2928	END DO
	2929	END DO
	2930	END DO
	2931
	2932	DO i = nl + 1, nd
	2933	DO il = 1, ncum
	2934	ma(il, i) = 0.
	2935	END DO
	2936	END DO
	2937
	2938	DO i = 1, nl
	2939	DO il = 1, ncum
	2940	IF (i<=(icb(il)-1)) THEN
	2941	ma(il, i) = 0
	2942	END IF
	2943	END DO
	2944	END DO
	2945
	2946	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2947	! icb represente de niveau ou se trouve la
	2948	! base du nuage , et inb le top du nuage
	2949	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2950
	2951	DO i = 1, nd
	2952	DO il = 1, ncum
	2953	mke(il, i) = upwd(il, i) + dnwd(il, i)
	2954	END DO
	2955	END DO
	2956
	2957	DO i = 1, nd
	2958	DO il = 1, ncum
	2959	rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd*(1.-rr(il, &
	2960	i))+rr(il,i)*cpv)
	2961	tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp
	2962	tps(il, i) = tp(il, i)
	2963	END DO
	2964	END DO
	2965
	2966
	2967	! * diagnose the in-cloud mixing ratio * ! cld
	2968	! * of condensed water * ! cld
	2969	! ! cld
	2970
	2971	DO i = 1, nd ! cld
	2972	DO il = 1, ncum ! cld
	2973	mac(il, i) = 0.0 ! cld
	2974	wa(il, i) = 0.0 ! cld
	2975	siga(il, i) = 0.0 ! cld
	2976	sax(il, i) = 0.0 ! cld
	2977	END DO ! cld
	2978	END DO ! cld
	2979
	2980	DO i = minorig, nl ! cld
	2981	DO k = i + 1, nl + 1 ! cld
	2982	DO il = 1, ncum ! cld
	2983	IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN ! cld
	2984	mac(il, i) = mac(il, i) + m(il, k) ! cld
	2985	END IF ! cld
	2986	END DO ! cld
	2987	END DO ! cld
	2988	END DO ! cld
	2989
	2990	DO i = 1, nl ! cld
	2991	DO j = 1, i ! cld
	2992	DO il = 1, ncum ! cld
	2993	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
	2994	.AND. j>=icb(il)) THEN ! cld
	2995	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
	2996	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
	2997	END IF ! cld
	2998	END DO ! cld
	2999	END DO ! cld
	3000	END DO ! cld
	3001
	3002	DO i = 1, nl ! cld
	3003	DO il = 1, ncum ! cld
	3004	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
	3005	.AND. sax(il,i)>0.0) THEN ! cld
	3006	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
	3007	END IF ! cld
	3008	END DO ! cld
	3009	END DO ! cld
	3010
	3011	DO i = 1, nl ! cld
	3012	DO il = 1, ncum ! cld
	3013	IF (wa(il,i)>0.0) & ! cld
	3014	siga(il, i) = mac(il, i)/wa(il, i) & ! cld
	3015	rrdtvp(il, i)/p(il, i)/100./delta ! cld
	3016	siga(il, i) = min(siga(il,i), 1.0) ! cld
	3017	! IM cf. FH
	3018	IF (iflag_clw==0) THEN
	3019	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
	3020	+(1.-siga(il,i))*qcond(il, i) ! cld
	3021	ELSE IF (iflag_clw==1) THEN
	3022	qcondc(il, i) = qcond(il, i) ! cld
	3023	END IF
	3024
	3025	END DO ! cld
	3026	END DO ! cld
	3027
	3028	RETURN
	3029	END SUBROUTINE cv30_yield
	3030
	3031	! !RomP >>>
	3032	SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, &
	3033	d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
	3034	IMPLICIT NONE
	3035
	3036	include "cv30param.h"
	3037
	3038	! inputs:
	3039	INTEGER ncum, nd, na, nloc, len
	3040	REAL ment(nloc, na, na), sij(nloc, na, na)
	3041	REAL clw(nloc, nd), elij(nloc, na, na)
	3042	REAL ep(nloc, na)
	3043	INTEGER icb(nloc), inb(nloc)
	3044	REAL vprecip(nloc, nd+1)
	3045	! ouputs:
	3046	REAL da(nloc, na), phi(nloc, na, na)
	3047	REAL phi2(nloc, na, na)
	3048	REAL d1a(nloc, na), dam(nloc, na)
	3049	REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
	3050	! variables pour tracer dans precip de l'AA et des mel
	3051	! local variables:
	3052	INTEGER i, j, k
	3053	REAL epm(nloc, na, na)
	3054
	3055	! variables d'Emanuel : du second indice au troisieme
	3056	! ---> tab(i,k,j) -> de l origine k a l arrivee j
	3057	! ment, sij, elij
	3058	! variables personnelles : du troisieme au second indice
	3059	! ---> tab(i,j,k) -> de k a j
	3060	! phi, phi2
	3061
	3062	! initialisations
	3063	DO j = 1, na
	3064	DO i = 1, ncum
	3065	da(i, j) = 0.
	3066	d1a(i, j) = 0.
	3067	dam(i, j) = 0.
	3068	eplamm(i, j) = 0.
	3069	END DO
	3070	END DO
	3071	DO k = 1, na
	3072	DO j = 1, na
	3073	DO i = 1, ncum
	3074	epm(i, j, k) = 0.
	3075	epmlmmm(i, j, k) = 0.
	3076	phi(i, j, k) = 0.
	3077	phi2(i, j, k) = 0.
	3078	END DO
	3079	END DO
	3080	END DO
	3081
	3082	! fraction deau condensee dans les melanges convertie en precip : epm
	3083	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
	3084	DO j = 1, na
	3085	DO k = 1, j - 1
	3086	DO i = 1, ncum
	3087	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
	3088	! !jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
	3089	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
	3090	! !
	3091	epm(i, j, k) = max(epm(i,j,k), 0.0)
	3092	END IF
	3093	END DO
	3094	END DO
	3095	END DO
	3096
	3097	DO j = 1, na
	3098	DO k = 1, na
	3099	DO i = 1, ncum
	3100	IF (k>=icb(i) .AND. k<=inb(i)) THEN
	3101	eplamm(i, j) = eplamm(i, j) + ep(i, j)clw(i, j)ment(i, j, k)*(1.- &
	3102	sij(i,j,k))
	3103	END IF
	3104	END DO
	3105	END DO
	3106	END DO
	3107
	3108	DO j = 1, na
	3109	DO k = 1, j - 1
	3110	DO i = 1, ncum
	3111	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
	3112	epmlmmm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
	3113	END IF
	3114	END DO
	3115	END DO
	3116	END DO
	3117
	3118	! matrices pour calculer la tendance des concentrations dans cvltr.F90
	3119	DO j = 1, na
	3120	DO k = 1, na
	3121	DO i = 1, ncum
	3122	da(i, j) = da(i, j) + (1.-sij(i,k,j))*ment(i, k, j)
	3123	phi(i, j, k) = sij(i, k, j)*ment(i, k, j)
	3124	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sij(i,k,j))
	3125	END DO
	3126	END DO
	3127	END DO
	3128
	3129	DO j = 1, na
	3130	DO k = 1, j - 1
	3131	DO i = 1, ncum
	3132	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, j, k)(1.-ep(i,k))*(1.- &
	3133	sij(i,k,j))
	3134	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
	3135	END DO
	3136	END DO
	3137	END DO
	3138
	3139	RETURN
	3140	END SUBROUTINE cv30_tracer
	3141	! RomP <<<
	3142
	3143	SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
	3144	vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
	3145	dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
[2488]	3146	epmlmmm, eplamm, wdtraina, wdtrainm,epmax_diag, iflag1, precip1, vprecip1, evap1, &
[1992]	3147	ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
	3148	dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
[2488]	3149	elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1,epmax_diag1) ! epmax_cape
[1992]	3150	IMPLICIT NONE
	3151
	3152	include "cv30param.h"
	3153
	3154	! inputs:
	3155	INTEGER len, ncum, nd, ntra, nloc
	3156	INTEGER idcum(nloc)
	3157	INTEGER iflag(nloc)
	3158	INTEGER inb(nloc)
	3159	REAL precip(nloc)
	3160	REAL vprecip(nloc, nd+1), evap(nloc, nd)
	3161	REAL ep(nloc, nd)
	3162	REAL sig(nloc, nd), w0(nloc, nd)
	3163	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	3164	REAL ftra(nloc, nd, ntra)
	3165	REAL ma(nloc, nd)
	3166	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
	3167	REAL qcondc(nloc, nd)
	3168	REAL wd(nloc), cape(nloc)
	3169	REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
[2488]	3170	REAL epmax_diag(nloc) ! epmax_cape
[1992]	3171	! RomP >>>
	3172	REAL phi2(nloc, nd, nd)
	3173	REAL d1a(nloc, nd), dam(nloc, nd)
	3174	REAL wdtraina(nloc, nd), wdtrainm(nloc, nd)
	3175	REAL sij(nloc, nd, nd)
	3176	REAL elij(nloc, nd, nd), clw(nloc, nd)
	3177	REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd)
	3178	! RomP <<<
	3179
	3180	! outputs:
	3181	INTEGER iflag1(len)
	3182	INTEGER inb1(len)
	3183	REAL precip1(len)
	3184	REAL vprecip1(len, nd+1), evap1(len, nd) !<<< RomP
	3185	REAL ep1(len, nd) !<<< RomP
	3186	REAL sig1(len, nd), w01(len, nd)
	3187	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
	3188	REAL ftra1(len, nd, ntra)
	3189	REAL ma1(len, nd)
	3190	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
	3191	REAL qcondc1(nloc, nd)
	3192	REAL wd1(nloc), cape1(nloc)
	3193	REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd)
[2488]	3194	REAL epmax_diag1(len) ! epmax_cape
[1992]	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)
[2488]	3213	epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
[1992]	3214	END DO
	3215
	3216	DO k = 1, nl
	3217	DO i = 1, ncum
	3218	vprecip1(idcum(i), k) = vprecip(i, k)
	3219	evap1(idcum(i), k) = evap(i, k) !<<< RomP
	3220	sig1(idcum(i), k) = sig(i, k)
	3221	w01(idcum(i), k) = w0(i, k)
	3222	ft1(idcum(i), k) = ft(i, k)
	3223	fq1(idcum(i), k) = fq(i, k)
	3224	fu1(idcum(i), k) = fu(i, k)
	3225	fv1(idcum(i), k) = fv(i, k)
	3226	ma1(idcum(i), k) = ma(i, k)
	3227	upwd1(idcum(i), k) = upwd(i, k)
	3228	dnwd1(idcum(i), k) = dnwd(i, k)
	3229	dnwd01(idcum(i), k) = dnwd0(i, k)
	3230	qcondc1(idcum(i), k) = qcondc(i, k)
	3231	da1(idcum(i), k) = da(i, k)
	3232	mp1(idcum(i), k) = mp(i, k)
	3233	! RomP >>>
	3234	ep1(idcum(i), k) = ep(i, k)
	3235	d1a1(idcum(i), k) = d1a(i, k)
	3236	dam1(idcum(i), k) = dam(i, k)
	3237	clw1(idcum(i), k) = clw(i, k)
	3238	eplamm1(idcum(i), k) = eplamm(i, k)
	3239	wdtraina1(idcum(i), k) = wdtraina(i, k)
	3240	wdtrainm1(idcum(i), k) = wdtrainm(i, k)
	3241	! RomP <<<
	3242	END DO
	3243	END DO
	3244
	3245	DO i = 1, ncum
	3246	sig1(idcum(i), nd) = sig(i, nd)
	3247	END DO
	3248
	3249
	3250	! do 2100 j=1,ntra
	3251	! do 2110 k=1,nd ! oct3
	3252	! do 2120 i=1,ncum
	3253	! ftra1(idcum(i),k,j)=ftra(i,k,j)
	3254	! 2120 continue
	3255	! 2110 continue
	3256	! 2100 continue
	3257	DO j = 1, nd
	3258	DO k = 1, nd
	3259	DO i = 1, ncum
	3260	sij1(idcum(i), k, j) = sij(i, k, j)
	3261	phi1(idcum(i), k, j) = phi(i, k, j)
	3262	phi21(idcum(i), k, j) = phi2(i, k, j)
	3263	elij1(idcum(i), k, j) = elij(i, k, j)
	3264	epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
	3265	END DO
	3266	END DO
	3267	END DO
	3268
	3269	RETURN
	3270	END SUBROUTINE cv30_uncompress
	3271
[2488]	3272	subroutine cv30_epmax_fn_cape(nloc,ncum,nd &
	3273	,cape,ep,hp,icb,inb,clw,nk,t,h,lv &
	3274	,epmax_diag)
	3275	implicit none
	3276
	3277	! On fait varier epmax en fn de la cape
	3278	! Il faut donc recalculer ep, et hp qui a déjà été calculé et
	3279	! qui en dépend
	3280	! Toutes les autres variables fn de ep sont calculées plus bas.
	3281
	3282	#include "cvthermo.h"
	3283	#include "cv30param.h"
	3284	#include "conema3.h"
	3285
	3286	! inputs:
	3287	integer ncum, nd, nloc
	3288	integer icb(nloc), inb(nloc)
	3289	real cape(nloc)
	3290	real clw(nloc,nd),lv(nloc,nd),t(nloc,nd),h(nloc,nd)
	3291	integer nk(nloc)
	3292	! inouts:
	3293	real ep(nloc,nd)
	3294	real hp(nloc,nd)
	3295	! outputs ou local
	3296	real epmax_diag(nloc)
	3297	! locals
	3298	integer i,k
	3299	real hp_bak(nloc,nd)
	3300
	3301	! on recalcule ep et hp
	3302
	3303	if (coef_epmax_cape.gt.1e-12) then
	3304	do i=1,ncum
	3305	epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
	3306	do k=1,nl
	3307	ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
	3308	ep(i,k)=amax1(ep(i,k),0.0)
	3309	ep(i,k)=amin1(ep(i,k),epmax_diag(i))
	3310	enddo
	3311	enddo
	3312
	3313	! On recalcule hp:
	3314	do k=1,nl
	3315	do i=1,ncum
	3316	hp_bak(i,k)=hp(i,k)
	3317	enddo
	3318	enddo
	3319	do k=1,nlp
	3320	do i=1,ncum
	3321	hp(i,k)=h(i,k)
	3322	enddo
	3323	enddo
	3324	do k=minorig+1,nl
	3325	do i=1,ncum
	3326	if((k.ge.icb(i)).and.(k.le.inb(i)))then
	3327	hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
	3328	endif
	3329	enddo
	3330	enddo !do k=minorig+1,n
	3331	! write(,) 'cv30_routines 6218: hp(1,20)=',hp(1,20)
	3332	do i=1,ncum
	3333	do k=1,nl
	3334	if (abs(hp_bak(i,k)-hp(i,k)).gt.0.01) then
	3335	write(,) 'i,k=',i,k
	3336	write(,) 'coef_epmax_cape=',coef_epmax_cape
	3337	write(,) 'epmax_diag(i)=',epmax_diag(i)
	3338	write(,) 'ep(i,k)=',ep(i,k)
	3339	write(,) 'hp(i,k)=',hp(i,k)
	3340	write(,) 'hp_bak(i,k)=',hp_bak(i,k)
	3341	write(,) 'h(i,k)=',h(i,k)
	3342	write(,) 'nk(i)=',nk(i)
	3343	write(,) 'h(i,nk(i))=',h(i,nk(i))
	3344	write(,) 'lv(i,k)=',lv(i,k)
	3345	write(,) 't(i,k)=',t(i,k)
	3346	write(,) 'clw(i,k)=',clw(i,k)
	3347	write(,) 'cpd,cpv=',cpd,cpv
	3348	stop
	3349	endif
	3350	enddo !do k=1,nl
	3351	enddo !do i=1,ncum
	3352	endif !if (coef_epmax_cape.gt.1e-12) then
	3353
	3354	return
	3355	end subroutine cv30_epmax_fn_cape
	3356
	3357

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

Download in other formats:

Original Format