Context Navigation

cv3_routines.F90 @ 2049

Last change on this file since 2049 was 2039, checked in by fhourdin, 10 years ago

Inclusion d'un appel supplementaire a geopot dans leapfrog dans

les versions dyn3d et dyn3dpar pour garantir la convergence
des 3 versions dyn3d/dyn3dpar/dyn3dmem.
La convergence fonctionne avec un compilateur gfortran 4.6.3 / openmpi
distribué sous ubuntu avec la nouvelle physique (NPv3.2) et le guidage.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Une seule exception : si on guide avec dyn3dpar et openMP.

Correction du guidage dans dyn3dmem.

un print supprime dans cv3_toutines.F90

Rm : On ne devrait perdre la convergence numérique avec les précédentes
versions que pour la nouvelle physique puisque le géopotentiel
n'est utilisé dans la physique que par les theermiques et sisvat.

=====================================================================

Added call to geopot in leapfrog in dyn3d and dyn3dpar in order to

insure numerical convergence with dyn3dmem.
The convergence dyn3d/dyn3dpar/dyn3dmem is satistied with
gfortran 4.6.3 / openmpi, new physics (NPv3.2) and nudging.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Only one exception : with nudging & dynd3dpar & openMP.

Bug fixing for nudging in dyn3dmem.

print supressed in cv3_routines.F90

Rm : Numerical convergence with previous releases should be lost for new
physics only since the geopotential is used only by thermals and
sisvat in the physics.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 123.9 KB

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

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

Download in other formats:

Original Format