Context Navigation

cv3_routines.F90 @ 2343

Last change on this file since 2343 was 2327, checked in by jyg, 9 years ago
Changes in output levels of some variables and Fixes of units for Ale and Alp variables. Cleaning of variable declarations in cv3_yield.
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: 132.1 KB

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

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

Download in other formats:

Original Format