Context Navigation

cv3_routines.F90 @ 2393

Last change on this file since 2393 was 2393, checked in by jyg, 9 years ago
Add various intializations of arrays in lmdz1d.F90 and in the convection scheme. Add output variables for boundary layer splitting.
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: 134.2 KB

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

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

Download in other formats:

Original Format