Context Navigation

cv3_routines.F @ 999

Last change on this file since 999 was 993, checked in by Laurent Fairhead, 16 years ago
Oubli de declaration de variables sur la derniere correction LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 101.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
[991]	1496	: ,ment,qent,uent,vent,nent,sij,elij,ments,qents,traent)
[524]	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)
[991]	1527	integer nent(nloc,nd)
[524]	1528
	1529	c local variables:
	1530	integer i, j, k, il, im, jm
	1531	integer num1, num2
	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)
[993]	2317	real esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
	2318	real th_wake(nloc,nd)
[524]	2319	real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
	2320	real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
	2321
[879]	2322	c print*,'cv3_yield declarations 3'
[524]	2323	c-------------------------------------------------------------
	2324
	2325	c initialization:
	2326
	2327	delti = 1.0/delt
[879]	2328	c print*,'cv3_yield initialisation delt', delt
	2329	cprecip,Vprecip,ft,fr,fu,fv,ftra
	2330	c : ,cbmf,upwd,dnwd,dnwd0,ma,mip
	2331	c : ,tls,tps,qcondc,wd
	2332	c : ,ftd,fqd )
[524]	2333	do il=1,ncum
	2334	precip(il)=0.0
[879]	2335	c Vprecip(il,nd+1)=0.0
[524]	2336	wd(il)=0.0 ! gust
	2337	enddo
	2338
	2339	do i=1,nd
	2340	do il=1,ncum
[879]	2341	Vprecip(il,i)=0.0
[524]	2342	ft(il,i)=0.0
	2343	fr(il,i)=0.0
	2344	fu(il,i)=0.0
	2345	fv(il,i)=0.0
[879]	2346	upwd(il,i)=0.0
	2347	dnwd(il,i)=0.0
	2348	dnwd0(il,i)=0.0
	2349	mip(il,i)=0.0
	2350	ftd(il,i)=0.0
	2351	fqd(il,i)=0.0
[524]	2352	qcondc(il,i)=0.0 ! cld
	2353	qcond(il,i)=0.0 ! cld
	2354	nqcond(il,i)=0.0 ! cld
	2355	enddo
	2356	enddo
[879]	2357	c print*,'cv3_yield initialisation 2'
	2358	do j=1,ntra
	2359	do i=1,nd
	2360	do il=1,ncum
	2361	ftra(il,i,j)=0.0
	2362	enddo
	2363	enddo
	2364	enddo
	2365	c print*,'cv3_yield initialisation 3'
[524]	2366	do i=1,nl
	2367	do il=1,ncum
	2368	lvcp(il,i)=lv(il,i)/cpn(il,i)
	2369	enddo
	2370	enddo
	2371
	2372
	2373	c
	2374	c * calculate surface precipitation in mm/day *
	2375	c
[879]	2376	do il=1,ncum
	2377	if(ep(il,inb(il)).ge.0.0001 .and. iflag(il) .le. 1)then
[524]	2378	if (cvflag_grav) then
[879]	2379	precip(il)=wt(il,1)sigd(il)water(il,1)86400.1000.
	2380	: /(rowl*grav)
[524]	2381	else
[879]	2382	precip(il)=wt(il,1)sigd(il)water(il,1)*8640.
[524]	2383	endif
[879]	2384	endif
	2385	enddo
	2386	c print*,'cv3_yield apres calcul precip'
[524]	2387
[619]	2388	C
[879]	2389	C === calculate vertical profile of precipitation in kg/m2/s ===
[619]	2390	C
[879]	2391	do i = 1,nl
	2392	do il=1,ncum
	2393	if(ep(il,inb(il)).ge.0.0001 .and. i.le.inb(il)
	2394	: .and. iflag(il) .le. 1)then
	2395	if (cvflag_grav) then
	2396	VPrecip(il,i) = wt(il,i)sigd(il)water(il,i)/grav
	2397	else
	2398	VPrecip(il,i) = wt(il,i)sigd(il)water(il,i)/10.
	2399	endif
	2400	endif
	2401	enddo
	2402	enddo
	2403	C
[524]	2404	c
	2405	c * Calculate downdraft velocity scale *
	2406	c * NE PAS UTILISER POUR L'INSTANT *
	2407	c
	2408	c! do il=1,ncum
	2409	c! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
[879]	2410	c! : /(sigd(il)*p(il,icb(il)))
[524]	2411	c! enddo
	2412
	2413	c
	2414	c * calculate tendencies of lowest level potential temperature *
	2415	c * and mixing ratio *
	2416	c
	2417	do il=1,ncum
	2418	work(il)=1.0/(ph(il,1)-ph(il,2))
[879]	2419	cbmf(il)=0.0
[524]	2420	enddo
	2421
	2422	do k=2,nl
	2423	do il=1,ncum
[879]	2424	if (k.ge.icb(il)) then
	2425	cbmf(il)=cbmf(il)+m(il,k)
[524]	2426	endif
	2427	enddo
	2428	enddo
	2429
[879]	2430	c print*,'cv3_yield avant ft'
	2431	c AM is the part of cbmf taken from the first level
[524]	2432	do il=1,ncum
[879]	2433	am(il)=cbmf(il)*wghti(il,1)
	2434	enddo
	2435	c
	2436	do il=1,ncum
	2437	if (iflag(il) .le. 1) then
	2438	c convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
	2439	cjyg Correction pour conserver l'eau
	2440	ccc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
	2441	ft(il,1)=-lvcp(il,1)sigd(il)evap(il,1) !precip
[524]	2442
	2443	if (cvflag_grav) then
[879]	2444	ft(il,1)=ft(il,1)-0.009gravsigd(il)*mp(il,2)
	2445	: t_wake(il,1)b(il,1)*work(il)
	2446	else
	2447	ft(il,1)=ft(il,1)-0.09sigd(il)mp(il,2)
	2448	: t_wake(il,1)b(il,1)*work(il)
	2449	endif
	2450
	2451	ft(il,1)=ft(il,1)+0.01sigd(il)wt(il,1)(cl-cpd)water(il,2)
	2452	: (t_wake(il,2)-t_wake(il,1))work(il)/cpn(il,1)
	2453
	2454	ftd(il,1) = ft(il,1) ! fin precip
	2455
	2456	if (cvflag_grav) then !sature
[524]	2457	if((0.01gravwork(il)*am(il)).ge.delti)iflag(il)=1!consist vect
[879]	2458	ft(il,1)=ft(il,1)+0.01gravwork(il)am(il)(t(il,2)-t(il,1)
[524]	2459	: +(gz(il,2)-gz(il,1))/cpn(il,1))
	2460	else
	2461	if((0.1work(il)am(il)).ge.delti)iflag(il)=1 !consistency vect
[879]	2462	ft(il,1)=ft(il,1)+0.1work(il)am(il)*(t(il,2)-t(il,1)
[524]	2463	: +(gz(il,2)-gz(il,1))/cpn(il,1))
	2464	endif
[879]	2465	endif ! iflag
	2466	enddo
[524]	2467
[970]	2468
[879]	2469	do j=2,nl
	2470	IF (iflag_mix .gt. 0) then
	2471	do il=1,ncum
[973]	2472	c FH WARNING a modifier :
	2473	cpinv=0.
	2474	c cpinv=1.0/cpn(il,1)
[879]	2475	if (j.le.inb(il) .and. iflag(il) .le. 1) then
	2476	if (cvflag_grav) then
	2477	ft(il,1)=ft(il,1)
	2478	: +0.01gravwork(il)ment(il,j,1)(hent(il,j,1)-h(il,1)
	2479	: +t(il,1)(cpv-cpd)(rr(il,1)-Qent(il,j,1)))*cpinv
	2480	else
	2481	ft(il,1)=ft(il,1)
	2482	: +0.1work(il)ment(il,j,1)*(hent(il,j,1)-h(il,1)
	2483	: +t(il,1)(cpv-cpd)(rr(il,1)-Qent(il,j,1)))*cpinv
	2484	endif ! cvflag_grav
	2485	endif ! j
	2486	enddo
	2487	ENDIF
	2488	enddo
	2489	! fin sature
[524]	2490
	2491
[879]	2492	do il=1,ncum
	2493	if (iflag(il) .le. 1) then
	2494	if (cvflag_grav) then
[524]	2495	Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
[879]	2496	fr(il,1)=0.01gravmp(il,2)(rp(il,2)-rr_wake(il,1))work(il)
	2497	: +sigd(il)*evap(il,1)
	2498	ccc : +sigd(il)0.5(evap(il,1)+evap(il,2))
[524]	2499
[879]	2500	fqd(il,1)=fr(il,1) !precip
[524]	2501
[879]	2502	fr(il,1)=fr(il,1)+0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
	2503
[524]	2504	fu(il,1)=fu(il,1)+0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1))
	2505	: +am(il)*(u(il,2)-u(il,1)))
	2506	fv(il,1)=fv(il,1)+0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1))
	2507	: +am(il)*(v(il,2)-v(il,1)))
	2508	else ! cvflag_grav
[879]	2509	fr(il,1)=0.1mp(il,2)(rp(il,2)-rr_wake(il,1))*work(il)
	2510	: +sigd(il)*evap(il,1)
	2511	ccc : +sigd(il)0.5(evap(il,1)+evap(il,2))
	2512	fqd(il,1)=fr(il,1) !precip
[524]	2513	fr(il,1)=fr(il,1)+0.1am(il)(rr(il,2)-rr(il,1))*work(il)
	2514	fu(il,1)=fu(il,1)+0.1work(il)(mp(il,2)*(up(il,2)-u(il,1))
	2515	: +am(il)*(u(il,2)-u(il,1)))
	2516	fv(il,1)=fv(il,1)+0.1work(il)(mp(il,2)*(vp(il,2)-v(il,1))
	2517	: +am(il)*(v(il,2)-v(il,1)))
[879]	2518	endif ! cvflag_grav
	2519	endif ! iflag
[524]	2520	enddo ! il
	2521
	2522
[879]	2523	do j=1,ntra
[524]	2524	do il=1,ncum
[879]	2525	if (iflag(il) .le. 1) then
	2526	if (cvflag_grav) then
	2527	ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
	2528	: (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2529	: +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2530	else
	2531	ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
	2532	: (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2533	: +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2534	endif
	2535	endif ! iflag
	2536	enddo
	2537	enddo
	2538
	2539	do j=2,nl
	2540	do il=1,ncum
	2541	if (j.le.inb(il) .and. iflag(il) .le. 1) then
[524]	2542	if (cvflag_grav) then
	2543	fr(il,1)=fr(il,1)
	2544	: +0.01gravwork(il)ment(il,j,1)(qent(il,j,1)-rr(il,1))
	2545	fu(il,1)=fu(il,1)
	2546	: +0.01gravwork(il)ment(il,j,1)(uent(il,j,1)-u(il,1))
	2547	fv(il,1)=fv(il,1)
	2548	: +0.01gravwork(il)ment(il,j,1)(vent(il,j,1)-v(il,1))
	2549	else ! cvflag_grav
	2550	fr(il,1)=fr(il,1)
	2551	: +0.1work(il)ment(il,j,1)*(qent(il,j,1)-rr(il,1))
	2552	fu(il,1)=fu(il,1)
	2553	: +0.1work(il)ment(il,j,1)*(uent(il,j,1)-u(il,1))
	2554	fv(il,1)=fv(il,1)
[879]	2555	: +0.1work(il)ment(il,j,1)*(vent(il,j,1)-v(il,1)) ! fin sature
[524]	2556	endif ! cvflag_grav
	2557	endif ! j
	2558	enddo
	2559	enddo
	2560
[879]	2561	do k=1,ntra
	2562	do j=2,nl
	2563	do il=1,ncum
	2564	if (j.le.inb(il) .and. iflag(il) .le. 1) then
[524]	2565
[879]	2566	if (cvflag_grav) then
	2567	ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
	2568	: *(traent(il,j,1,k)-tra(il,1,k))
	2569	else
	2570	ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
	2571	: *(traent(il,j,1,k)-tra(il,1,k))
	2572	endif
[524]	2573
[879]	2574	endif
	2575	enddo
	2576	enddo
	2577	enddo
	2578	c print*,'cv3_yield apres ft'
[524]	2579	c
	2580	c * calculate tendencies of potential temperature and mixing ratio *
	2581	c * at levels above the lowest level *
	2582	c
	2583	c * first find the net saturated updraft and downdraft mass fluxes *
	2584	c * through each level *
	2585	c
	2586
	2587	do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
	2588
	2589	num1=0
	2590	do il=1,ncum
[879]	2591	if(i.le.inb(il) .and. iflag(il) .le. 1)num1=num1+1
[524]	2592	enddo
	2593	if(num1.le.0)go to 500
	2594
	2595	call zilch(amp1,ncum)
	2596	call zilch(ad,ncum)
	2597
[879]	2598	do 440 k=1,nl+1
[524]	2599	do 441 il=1,ncum
[879]	2600	if(i.ge.icb(il)) then
	2601	if(k.ge.i+1.and. k.le.(inb(il)+1)) then
	2602	amp1(il)=amp1(il)+m(il,k)
	2603	endif
	2604	else
	2605	c AMP1 is the part of cbmf taken from layers I and lower
	2606	if(k.le.i) then
	2607	amp1(il)=amp1(il)+cbmf(il)*wghti(il,k)
	2608	endif
[524]	2609	endif
	2610	441 continue
	2611	440 continue
	2612
	2613	do 450 k=1,i
	2614	do 451 j=i+1,nl+1
	2615	do 452 il=1,ncum
	2616	if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
	2617	amp1(il)=amp1(il)+ment(il,k,j)
	2618	endif
	2619	452 continue
	2620	451 continue
	2621	450 continue
	2622
	2623	do 470 k=1,i-1
	2624	do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
	2625	do 472 il=1,ncum
	2626	if (i.le.inb(il) .and. j.le.inb(il)) then
	2627	ad(il)=ad(il)+ment(il,j,k)
	2628	endif
	2629	472 continue
	2630	471 continue
	2631	470 continue
	2632
	2633	do 1350 il=1,ncum
[879]	2634	if (i.le.inb(il) .and. iflag(il) .le. 1) then
[524]	2635	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2636	cpinv=1.0/cpn(il,i)
	2637
	2638	c convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
	2639	if (cvflag_grav) then
	2640	if((0.01gravdpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
	2641	else
	2642	if((0.1dpinvamp1(il)).ge.delti)iflag(il)=1 ! vecto
	2643	endif
	2644
[879]	2645	! precip
	2646	ccc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
	2647	ft(il,i)= -sigd(il)lvcp(il,i)evap(il,i)
	2648	rat=cpn(il,i-1)*cpinv
	2649	c
[524]	2650	if (cvflag_grav) then
[879]	2651	ft(il,i)=ft(il,i)-0.009gravsigd(il)
	2652	: (mp(il,i+1)t_wake(il,i)*b(il,i)
	2653	: -mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
	2654	ft(il,i)=ft(il,i)+0.01sigd(il)wt(il,i)(cl-cpd)water(il,i+1)
	2655	: (t_wake(il,i+1)-t_wake(il,i))dpinv*cpinv
	2656	ftd(il,i)=ft(il,i)
	2657	! fin precip
	2658	c
	2659	! sature
	2660	ft(il,i)=ft(il,i)+0.01gravdpinv(amp1(il)(t(il,i+1)-t(il,i)
[524]	2661	: +(gz(il,i+1)-gz(il,i))*cpinv)
	2662	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
[879]	2663
	2664	c
	2665	IF (iflag_mix .eq. 0) then
[524]	2666	ft(il,i)=ft(il,i)+0.01gravdpinvment(il,i,i)(hp(il,i)-h(il,i)
	2667	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
[879]	2668	ENDIF
	2669	c
[524]	2670	else ! cvflag_grav
[879]	2671	ft(il,i)=ft(il,i)-0.09*sigd(il)
	2672	: (mp(il,i+1)t_wake(il,i)*b(il,i)
	2673	: -mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
	2674	ft(il,i)=ft(il,i)+0.01sigd(il)wt(il,i)(cl-cpd)water(il,i+1)
	2675	: (t_wake(il,i+1)-t_wake(il,i))dpinv*cpinv
	2676	ftd(il,i)=ft(il,i)
	2677	! fin precip
	2678	c
	2679	! sature
	2680	ft(il,i)=ft(il,i)+0.1dpinv(amp1(il)*(t(il,i+1)-t(il,i)
[524]	2681	: +(gz(il,i+1)-gz(il,i))*cpinv)
	2682	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
[879]	2683
	2684	c
	2685	IF (iflag_mix .eq. 0) then
[524]	2686	ft(il,i)=ft(il,i)+0.1dpinvment(il,i,i)*(hp(il,i)-h(il,i)
	2687	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
[879]	2688	ENDIF
[524]	2689	endif ! cvflag_grav
	2690
	2691
[879]	2692	if (cvflag_grav) then
	2693	c sb: on ne fait pas encore la correction permettant de mieux
	2694	c conserver l'eau:
	2695	c jyg: correction permettant de mieux conserver l'eau:
	2696	ccc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
	2697	fr(il,i)=sigd(il)*evap(il,i)
	2698	: +0.01grav(mp(il,i+1)*(rp(il,i+1)-rr_wake(il,i))
	2699	: -mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
	2700	fqd(il,i)=fr(il,i) ! precip
[524]	2701
[879]	2702	fu(il,i)=0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i))
	2703	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	2704	fv(il,i)=0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i))
	2705	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
	2706	else ! cvflag_grav
	2707	ccc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
	2708	fr(il,i)=sigd(il)*evap(il,i)
	2709	: +0.1(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i))
	2710	: -mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
	2711	fqd(il,i)=fr(il,i) ! precip
	2712
	2713	fu(il,i)=0.1(mp(il,i+1)(up(il,i+1)-u(il,i))
	2714	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	2715	fv(il,i)=0.1(mp(il,i+1)(vp(il,i+1)-v(il,i))
	2716	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
	2717	endif ! cvflag_grav
	2718
	2719
[524]	2720	if (cvflag_grav) then
[879]	2721	fr(il,i)=fr(il,i)+0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i))
[524]	2722	: -ad(il)*(rr(il,i)-rr(il,i-1)))
	2723	fu(il,i)=fu(il,i)+0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i))
	2724	: -ad(il)*(u(il,i)-u(il,i-1)))
	2725	fv(il,i)=fv(il,i)+0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i))
	2726	: -ad(il)*(v(il,i)-v(il,i-1)))
	2727	else ! cvflag_grav
[879]	2728	fr(il,i)=fr(il,i)+0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il,i))
[524]	2729	: -ad(il)*(rr(il,i)-rr(il,i-1)))
	2730	fu(il,i)=fu(il,i)+0.1dpinv(amp1(il)*(u(il,i+1)-u(il,i))
	2731	: -ad(il)*(u(il,i)-u(il,i-1)))
	2732	fv(il,i)=fv(il,i)+0.1dpinv(amp1(il)*(v(il,i+1)-v(il,i))
	2733	: -ad(il)*(v(il,i)-v(il,i-1)))
	2734	endif ! cvflag_grav
	2735
	2736	endif ! i
	2737	1350 continue
	2738
[879]	2739	do k=1,ntra
	2740	do il=1,ncum
	2741	if (i.le.inb(il) .and. iflag(il) .le. 1) then
	2742	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2743	cpinv=1.0/cpn(il,i)
	2744	if (cvflag_grav) then
	2745	ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
	2746	: (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2747	: -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2748	else
	2749	ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
	2750	: (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2751	: -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2752	endif
	2753	endif
	2754	enddo
	2755	enddo
[524]	2756
	2757	do 480 k=1,i-1
[879]	2758	c
	2759	do il = 1,ncum
	2760	awat(il)=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
	2761	awat(il)=amax1(awat(il),0.0)
	2762	enddo
	2763	c
	2764	IF (iflag_mix .ne. 0) then
	2765	do il=1,ncum
	2766	if (i.le.inb(il) .and. iflag(il) .le. 1) then
	2767	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2768	cpinv=1.0/cpn(il,i)
	2769	if (cvflag_grav) then
	2770	ft(il,i)=ft(il,i)
	2771	: +0.01gravdpinvment(il,k,i)(hent(il,k,i)-h(il,i)
	2772	: +t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-Qent(il,k,i)))*cpinv
	2773
	2774	c
	2775	c
	2776	else
	2777	ft(il,i)=ft(il,i)
	2778	: +0.1dpinvment(il,k,i)*(hent(il,k,i)-h(il,i)
	2779	: +t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-Qent(il,k,i)))*cpinv
	2780	endif !cvflag_grav
	2781	endif ! i
	2782	enddo
	2783	ENDIF
	2784	c
[524]	2785	do 1370 il=1,ncum
[879]	2786	if (i.le.inb(il) .and. iflag(il) .le. 1) then
[524]	2787	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2788	cpinv=1.0/cpn(il,i)
[879]	2789	if (cvflag_grav) then
[524]	2790	fr(il,i)=fr(il,i)
[879]	2791	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-awat(il)-rr(il,i))
[524]	2792	fu(il,i)=fu(il,i)
	2793	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
	2794	fv(il,i)=fv(il,i)
	2795	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
	2796	else ! cvflag_grav
	2797	fr(il,i)=fr(il,i)
[879]	2798	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-awat(il)-rr(il,i))
[524]	2799	fu(il,i)=fu(il,i)
	2800	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
	2801	fv(il,i)=fv(il,i)
	2802	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
	2803	endif ! cvflag_grav
	2804
	2805	c (saturated updrafts resulting from mixing) ! cld
[879]	2806	qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat(il)) ! cld
[524]	2807	nqcond(il,i)=nqcond(il,i)+1. ! cld
	2808	endif ! i
	2809	1370 continue
	2810	480 continue
	2811
[879]	2812	do j=1,ntra
	2813	do k=1,i-1
	2814	do il=1,ncum
	2815	if (i.le.inb(il) .and. iflag(il) .le. 1) then
	2816	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2817	cpinv=1.0/cpn(il,i)
	2818	if (cvflag_grav) then
	2819	ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2820	: *(traent(il,k,i,j)-tra(il,i,j))
	2821	else
	2822	ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2823	: *(traent(il,k,i,j)-tra(il,i,j))
	2824	endif
	2825	endif
	2826	enddo
	2827	enddo
	2828	enddo
[524]	2829
	2830	do 490 k=i,nl+1
[879]	2831	c
	2832	IF (iflag_mix .ne. 0) then
	2833	do il=1,ncum
	2834	if (i.le.inb(il) .and. k.le.inb(il)
	2835	$ .and. iflag(il) .le. 1) then
	2836	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2837	cpinv=1.0/cpn(il,i)
	2838	if (cvflag_grav) then
	2839	ft(il,i)=ft(il,i)
	2840	: +0.01gravdpinvment(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	c
	2843	c
	2844	else
	2845	ft(il,i)=ft(il,i)
	2846	: +0.1dpinvment(il,k,i)*(hent(il,k,i)-h(il,i)
	2847	: +t(il,i)(cpv-cpd)(rr(il,i)-Qent(il,k,i)))*cpinv
	2848	endif !cvflag_grav
	2849	endif ! i
	2850	enddo
	2851	ENDIF
	2852	c
[524]	2853	do 1380 il=1,ncum
[879]	2854	if (i.le.inb(il) .and. k.le.inb(il)
	2855	$ .and. iflag(il) .le. 1) then
[524]	2856	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2857	cpinv=1.0/cpn(il,i)
	2858
	2859	if (cvflag_grav) then
	2860	fr(il,i)=fr(il,i)
	2861	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-rr(il,i))
	2862	fu(il,i)=fu(il,i)
	2863	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
	2864	fv(il,i)=fv(il,i)
	2865	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
[879]	2866	else ! cvflag_grav
[524]	2867	fr(il,i)=fr(il,i)
	2868	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-rr(il,i))
	2869	fu(il,i)=fu(il,i)
	2870	: +0.1dpinvment(il,k,i)*(uent(il,k,i)-u(il,i))
	2871	fv(il,i)=fv(il,i)
	2872	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
[879]	2873	endif ! cvflag_grav
[524]	2874	endif ! i and k
	2875	1380 continue
	2876	490 continue
	2877
[879]	2878	do j=1,ntra
	2879	do k=i,nl+1
	2880	do il=1,ncum
	2881	if (i.le.inb(il) .and. k.le.inb(il)
	2882	$ .and. iflag(il) .le. 1) then
	2883	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2884	cpinv=1.0/cpn(il,i)
	2885	if (cvflag_grav) then
	2886	ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2887	: *(traent(il,k,i,j)-tra(il,i,j))
	2888	else
	2889	ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2890	: *(traent(il,k,i,j)-tra(il,i,j))
	2891	endif
	2892	endif ! i and k
	2893	enddo
	2894	enddo
	2895	enddo
[524]	2896
	2897	c sb: interface with the cloud parameterization: ! cld
	2898
	2899	do k=i+1,nl
[879]	2900	do il=1,ncum
	2901	if (k.le.inb(il) .and. i.le.inb(il)
	2902	$ .and. iflag(il) .le. 1) then ! cld
[524]	2903	C (saturated downdrafts resulting from mixing) ! cld
	2904	qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
	2905	nqcond(il,i)=nqcond(il,i)+1. ! cld
	2906	endif ! cld
	2907	enddo ! cld
	2908	enddo ! cld
	2909
	2910	C (particular case: no detraining level is found) ! cld
	2911	do il=1,ncum ! cld
[879]	2912	if (i.le.inb(il) .and. nent(il,i).eq.0
	2913	$ .and. iflag(il) .le. 1) then ! cld
[524]	2914	qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
	2915	nqcond(il,i)=nqcond(il,i)+1. ! cld
	2916	endif ! cld
	2917	enddo ! cld
	2918
	2919	do il=1,ncum ! cld
[879]	2920	if (i.le.inb(il) .and. nqcond(il,i).ne.0
	2921	$ .and. iflag(il) .le. 1) then ! cld
[524]	2922	qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
	2923	endif ! cld
	2924	enddo
	2925
[879]	2926	do j=1,ntra
	2927	do il=1,ncum
	2928	if (i.le.inb(il) .and. iflag(il) .le. 1) then
	2929	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2930	cpinv=1.0/cpn(il,i)
[524]	2931
[879]	2932	if (cvflag_grav) then
	2933	ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
	2934	: (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2935	: -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
	2936	else
	2937	ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
	2938	: (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2939	: -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
	2940	endif
	2941	endif ! i
	2942	enddo
	2943	enddo
[524]	2944
[879]	2945
[524]	2946	500 continue
	2947
	2948
	2949	c * move the detrainment at level inb down to level inb-1 *
	2950	c * in such a way as to preserve the vertically *
	2951	c * integrated enthalpy and water tendencies *
	2952	c
	2953	do 503 il=1,ncum
[879]	2954	IF (iflag(il) .le. 1) THEN
[524]	2955	ax=0.1ment(il,inb(il),inb(il))(hp(il,inb(il))-h(il,inb(il))
	2956	: +t(il,inb(il))*(cpv-cpd)
	2957	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
	2958	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
	2959	ft(il,inb(il))=ft(il,inb(il))-ax
	2960	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
	2961	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
	2962	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
	2963
	2964	bx=0.1ment(il,inb(il),inb(il))(qent(il,inb(il),inb(il))
	2965	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2966	fr(il,inb(il))=fr(il,inb(il))-bx
	2967	fr(il,inb(il)-1)=fr(il,inb(il)-1)
	2968	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
	2969	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2970
	2971	cx=0.1ment(il,inb(il),inb(il))(uent(il,inb(il),inb(il))
	2972	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2973	fu(il,inb(il))=fu(il,inb(il))-cx
	2974	fu(il,inb(il)-1)=fu(il,inb(il)-1)
	2975	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
	2976	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2977
	2978	dx=0.1ment(il,inb(il),inb(il))(vent(il,inb(il),inb(il))
	2979	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2980	fv(il,inb(il))=fv(il,inb(il))-dx
	2981	fv(il,inb(il)-1)=fv(il,inb(il)-1)
	2982	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
	2983	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
[879]	2984	ENDIF !iflag
[524]	2985	503 continue
	2986
[879]	2987	do j=1,ntra
	2988	do il=1,ncum
	2989	IF (iflag(il) .le. 1) THEN
	2990	ex=0.1*ment(il,inb(il),inb(il))
	2991	: *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
	2992	: /(ph(i l,inb(il))-ph(il,inb(il)+1))
	2993	ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
	2994	ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
	2995	: +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
	2996	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2997	ENDIF !iflag
	2998	enddo
	2999	enddo
[524]	3000
	3001	c
[879]	3002	c * homogenize tendencies below cloud base *
[524]	3003	c
	3004	c
	3005	do il=1,ncum
	3006	asum(il)=0.0
	3007	bsum(il)=0.0
	3008	csum(il)=0.0
	3009	dsum(il)=0.0
[992]	3010	esum(il)=0.0
	3011	fsum(il)=0.0
	3012	gsum(il)=0.0
	3013	hsum(il)=0.0
[524]	3014	enddo
[992]	3015	c
[524]	3016	do i=1,nl
	3017	do il=1,ncum
[992]	3018	th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
	3019	enddo
	3020	enddo
	3021	c
	3022	do i=1,nl
	3023	do il=1,ncum
[879]	3024	if (i.le.(icb(il)-1) .and. iflag(il) .le. 1) then
[992]	3025	cjyg Saturated part : use T profile
	3026	asum(il)=asum(il)+(ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
	3027	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))
	3028	: (lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))
[524]	3029	: *(ph(il,i)-ph(il,i+1))
	3030	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
	3031	: *(ph(il,i)-ph(il,i+1))
	3032	dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
[992]	3033	cjyg Unsaturated part : use T_wake profile
	3034	esum(il)=esum(il)+ftd(il,i)*(ph(il,i)-ph(il,i+1))
	3035	fsum(il)=fsum(il)+fqd(il,i)
	3036	: (lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))
	3037	: *(ph(il,i)-ph(il,i+1))
	3038	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))
	3039	: *(ph(il,i)-ph(il,i+1))
	3040	hsum(il)=hsum(il)+t_wake(il,i)
	3041	; *(ph(il,i)-ph(il,i+1))/th_wake(il,i)
[879]	3042	endif
[524]	3043	enddo
	3044	enddo
	3045
	3046	c!!! do 700 i=1,icb(il)-1
	3047	do i=1,nl
	3048	do il=1,ncum
[879]	3049	if (i.le.(icb(il)-1) .and. iflag(il) .le. 1) then
[992]	3050	ftd(il,i)=esum(il)t_wake(il,i)/(th_wake(il,i)hsum(il))
	3051	fqd(il,i)=fsum(il)/gsum(il)
	3052	ft(il,i)=ftd(il,i)+asum(il)t(il,i)/(th(il,i)dsum(il))
	3053	fr(il,i)=fqd(il,i)+bsum(il)/csum(il)
[524]	3054	endif
	3055	enddo
	3056	enddo
	3057	c
	3058	c * reset counter and return *
	3059	c
	3060	do il=1,ncum
	3061	sig(il,nd)=2.0
	3062	enddo
	3063
	3064
	3065	do i=1,nd
	3066	do il=1,ncum
	3067	upwd(il,i)=0.0
	3068	dnwd(il,i)=0.0
	3069	enddo
	3070	enddo
[879]	3071
[524]	3072	do i=1,nl
	3073	do il=1,ncum
	3074	dnwd0(il,i)=-mp(il,i)
	3075	enddo
	3076	enddo
	3077	do i=nl+1,nd
	3078	do il=1,ncum
	3079	dnwd0(il,i)=0.
	3080	enddo
	3081	enddo
	3082
	3083
	3084	do i=1,nl
	3085	do il=1,ncum
	3086	if (i.ge.icb(il) .and. i.le.inb(il)) then
	3087	upwd(il,i)=0.0
	3088	dnwd(il,i)=0.0
	3089	endif
	3090	enddo
	3091	enddo
	3092
	3093	do i=1,nl
	3094	do k=1,nl
	3095	do il=1,ncum
	3096	up1(il,k,i)=0.0
	3097	dn1(il,k,i)=0.0
	3098	enddo
	3099	enddo
	3100	enddo
	3101
	3102	do i=1,nl
	3103	do k=i,nl
	3104	do n=1,i-1
	3105	do il=1,ncum
	3106	if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
	3107	up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
	3108	dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
	3109	endif
	3110	enddo
	3111	enddo
	3112	enddo
	3113	enddo
	3114
[879]	3115	do i=1,nl
	3116	do k=1,nl
	3117	do il=1,ncum
	3118	if(i.ge.icb(il)) then
	3119	if(k.ge.i.and. k.le.(inb(il))) then
	3120	upwd(il,i)=upwd(il,i)+m(il,k)
	3121	endif
	3122	else
	3123	if(k.lt.i) then
	3124	upwd(il,i)=upwd(il,i)+cbmf(il)*wghti(il,k)
	3125	endif
	3126	endif
	3127	cc print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
	3128	end do
	3129	end do
	3130	end do
	3131
[524]	3132	do i=2,nl
	3133	do k=i,nl
	3134	do il=1,ncum
	3135	ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
	3136	if (i.le.inb(il).and.k.le.inb(il)) then
[879]	3137	upwd(il,i)=upwd(il,i)+up1(il,k,i)
[524]	3138	dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
	3139	endif
[879]	3140	cc print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
[524]	3141	enddo
	3142	enddo
	3143	enddo
	3144
	3145
	3146	c!!! DO il=1,ncum
	3147	c!!! do i=icb(il),inb(il)
[879]	3148	c!!!
[524]	3149	c!!! upwd(il,i)=0.0
	3150	c!!! dnwd(il,i)=0.0
	3151	c!!! do k=i,inb(il)
	3152	c!!! up1=0.0
	3153	c!!! dn1=0.0
	3154	c!!! do n=1,i-1
	3155	c!!! up1=up1+ment(il,n,k)
	3156	c!!! dn1=dn1-ment(il,k,n)
	3157	c!!! enddo
	3158	c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
	3159	c!!! dnwd(il,i)=dnwd(il,i)+dn1
	3160	c!!! enddo
	3161	c!!! enddo
	3162	c!!!
	3163	c!!! ENDDO
	3164
	3165	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	3166	c determination de la variation de flux ascendant entre
[879]	3167	c deux niveau non dilue mip
[524]	3168	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	3169
	3170	do i=1,nl
	3171	do il=1,ncum
[879]	3172	mip(il,i)=m(il,i)
[524]	3173	enddo
	3174	enddo
	3175
	3176	do i=nl+1,nd
	3177	do il=1,ncum
[879]	3178	mip(il,i)=0.
[524]	3179	enddo
	3180	enddo
	3181
	3182	do i=1,nd
	3183	do il=1,ncum
	3184	ma(il,i)=0
	3185	enddo
	3186	enddo
	3187
	3188	do i=1,nl
	3189	do j=i,nl
	3190	do il=1,ncum
	3191	ma(il,i)=ma(il,i)+m(il,j)
	3192	enddo
	3193	enddo
	3194	enddo
	3195
	3196	do i=nl+1,nd
	3197	do il=1,ncum
	3198	ma(il,i)=0.
	3199	enddo
	3200	enddo
	3201
	3202	do i=1,nl
	3203	do il=1,ncum
	3204	if (i.le.(icb(il)-1)) then
	3205	ma(il,i)=0
	3206	endif
	3207	enddo
	3208	enddo
	3209
	3210	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	3211	c icb represente de niveau ou se trouve la
	3212	c base du nuage , et inb le top du nuage
	3213	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	3214
	3215	do i=1,nd
	3216	do il=1,ncum
	3217	mke(il,i)=upwd(il,i)+dnwd(il,i)
	3218	enddo
	3219	enddo
	3220
	3221	do i=1,nd
	3222	DO 999 il=1,ncum
	3223	rdcp=(rrd(1.-rr(il,i))-rr(il,i)rrv)
	3224	: /(cpd(1.-rr(il,i))+rr(il,i)cpv)
	3225	tls(il,i)=t(il,i)(1000.0/p(il,i))*rdcp
	3226	tps(il,i)=tp(il,i)
	3227	999 CONTINUE
	3228	enddo
	3229
	3230	c
	3231	c * diagnose the in-cloud mixing ratio * ! cld
	3232	c * of condensed water * ! cld
	3233	c ! cld
	3234
	3235	do i=1,nd ! cld
	3236	do il=1,ncum ! cld
	3237	mac(il,i)=0.0 ! cld
	3238	wa(il,i)=0.0 ! cld
	3239	siga(il,i)=0.0 ! cld
	3240	sax(il,i)=0.0 ! cld
	3241	enddo ! cld
	3242	enddo ! cld
	3243
	3244	do i=minorig, nl ! cld
	3245	do k=i+1,nl+1 ! cld
	3246	do il=1,ncum ! cld
[879]	3247	if (i.le.inb(il) .and. k.le.(inb(il)+1)
	3248	$ .and. iflag(il) .le. 1) then ! cld
[524]	3249	mac(il,i)=mac(il,i)+m(il,k) ! cld
	3250	endif ! cld
	3251	enddo ! cld
	3252	enddo ! cld
	3253	enddo ! cld
	3254
	3255	do i=1,nl ! cld
	3256	do j=1,i ! cld
	3257	do il=1,ncum ! cld
	3258	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
[879]	3259	: .and. j.ge.icb(il)
	3260	: .and. iflag(il) .le. 1 ) then ! cld
[524]	3261	sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
	3262	: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
	3263	endif ! cld
	3264	enddo ! cld
	3265	enddo ! cld
	3266	enddo ! cld
	3267
	3268	do i=1,nl ! cld
	3269	do il=1,ncum ! cld
	3270	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
[879]	3271	: .and. sax(il,i).gt.0.0
	3272	: .and. iflag(il) .le. 1 ) then ! cld
[524]	3273	wa(il,i)=sqrt(2.*sax(il,i)) ! cld
	3274	endif ! cld
	3275	enddo ! cld
	3276	enddo ! cld
[879]	3277
[524]	3278	do i=1,nl ! cld
	3279	do il=1,ncum ! cld
[879]	3280	if (wa(il,i).gt.0.0 .and. iflag(il) .le. 1) ! cld
[524]	3281	: siga(il,i)=mac(il,i)/wa(il,i) ! cld
	3282	: rrdtvp(il,i)/p(il,i)/100./delta ! cld
	3283	siga(il,i) = min(siga(il,i),1.0) ! cld
	3284	cIM cf. FH
	3285	if (iflag_clw.eq.0) then
	3286	qcondc(il,i)=siga(il,i)clw(il,i)(1.-ep(il,i)) ! cld
	3287	: + (1.-siga(il,i))*qcond(il,i) ! cld
	3288	else if (iflag_clw.eq.1) then
[879]	3289	qcondc(il,i)=qcond(il,i) ! cld
[524]	3290	endif
	3291
	3292	enddo ! cld
[879]	3293	enddo
	3294	c print*,'cv3_yield fin'
	3295	! cld
[524]	3296	return
	3297	end
	3298
	3299
	3300	SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
	3301	: ,iflag
[879]	3302	: ,precip,sig,w0
[524]	3303	: ,ft,fq,fu,fv,ftra
	3304	: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
	3305	: ,iflag1
[879]	3306	: ,precip1,sig1,w01
[524]	3307	: ,ft1,fq1,fu1,fv1,ftra1
	3308	: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
[879]	3309	: )
[524]	3310	implicit none
	3311
[879]	3312	#include "cv3param.h"
[524]	3313
	3314	c inputs:
	3315	integer len, ncum, nd, ntra, nloc
	3316	integer idcum(nloc)
	3317	integer iflag(nloc)
	3318	real precip(nloc)
	3319	real sig(nloc,nd), w0(nloc,nd)
	3320	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
	3321	real ftra(nloc,nd,ntra)
	3322	real Ma(nloc,nd)
	3323	real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
	3324	real qcondc(nloc,nd)
	3325	real wd(nloc),cape(nloc)
	3326
	3327	c outputs:
	3328	integer iflag1(len)
	3329	real precip1(len)
	3330	real sig1(len,nd), w01(len,nd)
	3331	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
	3332	real ftra1(len,nd,ntra)
	3333	real Ma1(len,nd)
	3334	real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
	3335	real qcondc1(nloc,nd)
	3336	real wd1(nloc),cape1(nloc)
	3337
	3338	c local variables:
	3339	integer i,k,j
	3340
	3341	do 2000 i=1,ncum
	3342	precip1(idcum(i))=precip(i)
	3343	iflag1(idcum(i))=iflag(i)
	3344	wd1(idcum(i))=wd(i)
	3345	cape1(idcum(i))=cape(i)
	3346	2000 continue
	3347
	3348	do 2020 k=1,nl
	3349	do 2010 i=1,ncum
	3350	sig1(idcum(i),k)=sig(i,k)
	3351	w01(idcum(i),k)=w0(i,k)
	3352	ft1(idcum(i),k)=ft(i,k)
	3353	fq1(idcum(i),k)=fq(i,k)
	3354	fu1(idcum(i),k)=fu(i,k)
	3355	fv1(idcum(i),k)=fv(i,k)
	3356	Ma1(idcum(i),k)=Ma(i,k)
	3357	upwd1(idcum(i),k)=upwd(i,k)
	3358	dnwd1(idcum(i),k)=dnwd(i,k)
	3359	dnwd01(idcum(i),k)=dnwd0(i,k)
	3360	qcondc1(idcum(i),k)=qcondc(i,k)
	3361	2010 continue
	3362	2020 continue
	3363
	3364	do 2200 i=1,ncum
	3365	sig1(idcum(i),nd)=sig(i,nd)
	3366	2200 continue
	3367
	3368
[879]	3369	do 2100 j=1,ntra
	3370	c oct3 do 2110 k=1,nl
	3371	do 2110 k=1,nd ! oct3
	3372	do 2120 i=1,ncum
	3373	ftra1(idcum(i),k,j)=ftra(i,k,j)
	3374	2120 continue
	3375	2110 continue
	3376	2100 continue
[524]	3377	return
	3378	end
	3379

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

Download in other formats:

Original Format