Context Navigation

cv3_routines.F @ 5416

Last change on this file since 5416 was 1467, checked in by idelkadi, 14 years ago
Input file conv_param.data added to change convection parameters without changing them inside the code
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 107.1 KB

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

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

Download in other formats:

Original Format