Context Navigation

cv3_routines.F @ 433

Last change on this file since 433 was 433, checked in by lmdzadmin, 22 years ago

Convergence avec la version de Ionela dec 2002

YOMCST.? : suppression RI0 (IM)
albedo.F : facteur 1.2 sur le nouveau calcul (IM)
clesphys.h : rajout de différentes ctes (concentration des gaz) (IM)
clmain.F : separation des flux LW, SW (JLD)

remplace qsurf par yqsol (IM)

conf_phys.F90 : rajout de différentes ctes (gaz + orbite) (IM)
convect3.F : DPINV+SIGD*0.5*(EVAP(1)+EVAP(2)) (SBL)
cv3_routines.F:
cvparam3.h : compatibilite avec conema3 TEMPORAIRE (FH)
phyetat0.F : lecture de co2_ppm et solaire pour tests de coherence
phyredem.F : co2_ppm et solaire passé en common
physiq.F : separation flux LW, SW

rajout diagnostiques (slp, w500)
suppression iflag_con = 4
clwcon0=qcondc (FH)
position dU "ENDIF ! ok_cvl"

radlwsw.F : passage des concentrations gaz dans un common (IM)

PEMIS(i) = 1.0 (JLD pour cohérence ORCHIDEE)

stdlevvar.F90 :
suphec.F : suppression init. des ctes orbitales (IM)

nouvelles E/S (ini_hist..., write_hist...)

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format