Context Navigation

cv3_routines.F @ 5434

Last change on this file since 5434 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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

File size: 106.6 KB

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

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

Download in other formats:

Original Format