Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

cv3_routines.F @ 5392

Last change on this file since 5392 was 416, checked in by lmdzadmin, 22 years ago
Inclusion initiale
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 90.5 KB

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

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

Download in other formats:

Original Format