Context Navigation

cv3_routines.F @ 949

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

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

Download in other formats:

Original Format