Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

cv3_routines.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 125.2 KB

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

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

Download in other formats:

Original Format