Context Navigation

cv3_routines.F @ 2116

Last change on this file since 2116 was 704, checked in by Laurent Fairhead, 18 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 93.5 KB

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

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

Download in other formats:

Original Format