Context Navigation

cv3_routines.F @ 1071

Last change on this file since 1071 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 93.6 KB

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

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

Download in other formats:

Original Format