Context Navigation

cv3_routines.F @ 1227

Last change on this file since 1227 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 105.1 KB

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

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

Download in other formats:

Original Format