Context Navigation

cv3_routines.F90 @ 2094

Last change on this file since 2094 was 2078, checked in by idelkadi, 10 years ago
Rajout des directives OPENMP pour la lecture des cles et parametres de convection et wakes dans les fichier conv_param.data, ep_param.data et wake_param.data
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: 124.1 KB

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

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

Download in other formats:

Original Format