Context Navigation

cv3_routines.F @ 669

Last change on this file since 669 was 546, checked in by lmdzadmin, 20 years ago
Suite des adaptations au portage sur SGI et VPP pour Prism (Caubel & Demory) LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 92.3 KB

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

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

Download in other formats:

Original Format