Context Navigation

cv30_routines.F90 @ 2496

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

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

Download in other formats:

Original Format