Context Navigation

cv3_routines.F90 @ 2390

Last change on this file since 2390 was 2376, checked in by jyg, 10 years ago
Add initialization of ice fraction below cloud base in cv3_undilute2 in cv3_routines.F90.
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.3 KB

Rev	Line
[1992]	1
[1403]	2	! $Id: cv3_routines.F90 2376 2015-10-20 09:57:39Z musat $
[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
[2376]	1076	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: frac
[524]	1077
[2007]	1078	!local variables:
[2253]	1079	INTEGER i, j, k
[1992]	1080	REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
	1081	REAL als
	1082	REAL qsat_new, snew, qi(nloc, nd)
	1083	REAL by, defrac, pden, tbis
	1084	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
	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
[2376]	1537	! Below cloud base, set ice fraction to cloud base value
	1538	DO k = 1, nl
	1539	DO i = 1, ncum
	1540	IF (k<icb(i)) THEN
	1541	frac(i,k) = frac(i,icb(i))
	1542	END IF
	1543	END DO
	1544	END DO
[2257]	1545	!
	1546	ELSE
	1547	!
	1548	DO k = minorig + 1, nl
	1549	DO i = 1, ncum
	1550	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
[1992]	1551	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k)
	1552	END IF
[2257]	1553	END DO
[1992]	1554	END DO
[2257]	1555	!
	1556	END IF ! (cvflag_ice)
[524]	1557
[1992]	1558	RETURN
	1559	END SUBROUTINE cv3_undilute2
[524]	1560
[2007]	1561	SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
	1562	pbase, p, ph, tv, buoy, &
	1563	sig, w0, cape, m, iflag)
[1992]	1564	IMPLICIT NONE
[524]	1565
[2007]	1566	! ===================================================================
	1567	! --- CLOSURE OF CONVECT3
	1568	!
	1569	! vectorization: S. Bony
	1570	! ===================================================================
[524]	1571
[1992]	1572	include "cvthermo.h"
	1573	include "cv3param.h"
[524]	1574
[2007]	1575	!input:
[1992]	1576	INTEGER ncum, nd, nloc
	1577	INTEGER icb(nloc), inb(nloc)
	1578	REAL pbase(nloc)
	1579	REAL p(nloc, nd), ph(nloc, nd+1)
	1580	REAL tv(nloc, nd), buoy(nloc, nd)
[524]	1581
[2007]	1582	!input/output:
[1992]	1583	REAL sig(nloc, nd), w0(nloc, nd)
	1584	INTEGER iflag(nloc)
[524]	1585
[2007]	1586	!output:
[1992]	1587	REAL cape(nloc)
	1588	REAL m(nloc, nd)
[524]	1589
[2007]	1590	!local variables:
[1992]	1591	INTEGER i, j, k, icbmax
	1592	REAL deltap, fac, w, amu
	1593	REAL dtmin(nloc, nd), sigold(nloc, nd)
	1594	REAL cbmflast(nloc)
[524]	1595
[776]	1596
[2007]	1597	! -------------------------------------------------------
	1598	! -- Initialization
	1599	! -------------------------------------------------------
[524]	1600
[1992]	1601	DO k = 1, nl
	1602	DO i = 1, ncum
	1603	m(i, k) = 0.0
	1604	END DO
	1605	END DO
[524]	1606
[2007]	1607	! -------------------------------------------------------
	1608	! -- Reset sig(i) and w0(i) for i>inb and i<icb
	1609	! -------------------------------------------------------
[524]	1610
[2007]	1611	! update sig and w0 above LNB:
[879]	1612
[1992]	1613	DO k = 1, nl - 1
	1614	DO i = 1, ncum
	1615	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
[2007]	1616	sig(i, k) = beta*sig(i, k) + &
	1617	2.alphabuoy(i, inb(i))*abs(buoy(i,inb(i)))
[1992]	1618	sig(i, k) = amax1(sig(i,k), 0.0)
	1619	w0(i, k) = beta*w0(i, k)
	1620	END IF
	1621	END DO
	1622	END DO
[1494]	1623
[2007]	1624	! compute icbmax:
[524]	1625
[1992]	1626	icbmax = 2
	1627	DO i = 1, ncum
	1628	icbmax = max(icbmax, icb(i))
	1629	END DO
[524]	1630
[2007]	1631	! update sig and w0 below cloud base:
[1494]	1632
[1992]	1633	DO k = 1, icbmax
	1634	DO i = 1, ncum
	1635	IF (k<=icb(i)) THEN
[2007]	1636	sig(i, k) = beta*sig(i, k) - &
	1637	2.alphabuoy(i, icb(i))*buoy(i, icb(i))
[1992]	1638	sig(i, k) = max(sig(i,k), 0.0)
	1639	w0(i, k) = beta*w0(i, k)
	1640	END IF
	1641	END DO
	1642	END DO
[524]	1643
[2007]	1644	!! if(inb.lt.(nl-1))then
	1645	!! do 85 i=inb+1,nl-1
	1646	!! sig(i)=betasig(i)+2.alphabuoy(inb)
	1647	!! 1 abs(buoy(inb))
	1648	!! sig(i)=max(sig(i),0.0)
	1649	!! w0(i)=beta*w0(i)
	1650	!! 85 continue
	1651	!! end if
[524]	1652
[2007]	1653	!! do 87 i=1,icb
	1654	!! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
	1655	!! sig(i)=max(sig(i),0.0)
	1656	!! w0(i)=beta*w0(i)
	1657	!! 87 continue
[524]	1658
[2007]	1659	! -------------------------------------------------------------
	1660	! -- Reset fractional areas of updrafts and w0 at initial time
	1661	! -- and after 10 time steps of no convection
	1662	! -------------------------------------------------------------
[524]	1663
[1992]	1664	DO k = 1, nl - 1
	1665	DO i = 1, ncum
	1666	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
	1667	sig(i, k) = 0.0
	1668	w0(i, k) = 0.0
	1669	END IF
	1670	END DO
	1671	END DO
[524]	1672
[2007]	1673	! -------------------------------------------------------------
	1674	! -- Calculate convective available potential energy (cape),
	1675	! -- vertical velocity (w), fractional area covered by
	1676	! -- undilute updraft (sig), and updraft mass flux (m)
	1677	! -------------------------------------------------------------
[1849]	1678
[1992]	1679	DO i = 1, ncum
	1680	cape(i) = 0.0
	1681	END DO
[1849]	1682
[2007]	1683	! compute dtmin (minimum buoyancy between ICB and given level k):
[1849]	1684
[1992]	1685	DO i = 1, ncum
	1686	DO k = 1, nl
	1687	dtmin(i, k) = 100.0
	1688	END DO
	1689	END DO
[524]	1690
[1992]	1691	DO i = 1, ncum
	1692	DO k = 1, nl
	1693	DO j = minorig, nl
[2007]	1694	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
[1992]	1695	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
	1696	END IF
	1697	END DO
	1698	END DO
	1699	END DO
[1849]	1700
[2007]	1701	! the interval on which cape is computed starts at pbase :
[1849]	1702
[1992]	1703	DO k = 1, nl
	1704	DO i = 1, ncum
[1849]	1705
[1992]	1706	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
[1849]	1707
[1992]	1708	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
	1709	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
	1710	cape(i) = amax1(0.0, cape(i))
	1711	sigold(i, k) = sig(i, k)
[1849]	1712
[2007]	1713	! dtmin(i,k)=100.0
	1714	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
	1715	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
	1716	! 97 continue
[1849]	1717
[1992]	1718	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
	1719	sig(i, k) = max(sig(i,k), 0.0)
	1720	sig(i, k) = amin1(sig(i,k), 0.01)
	1721	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
	1722	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
	1723	amu = 0.5(sig(i,k)+sigold(i,k))w
	1724	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
	1725	w0(i, k) = w
	1726	END IF
[1849]	1727
[1992]	1728	END DO
	1729	END DO
[1849]	1730
[1992]	1731	DO i = 1, ncum
	1732	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
[2007]	1733	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]	1734	sig(i, icb(i)) = sig(i, icb(i)+1)
	1735	sig(i, icb(i)-1) = sig(i, icb(i))
	1736	END DO
[1849]	1737
[2007]	1738	! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
	1739	! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
	1740	! ccc the final mass flux (cbmflast) is greater than the target mass flux
	1741	! ccc (cbmf) ??).
	1742	! cc
	1743	! c do i = 1,ncum
	1744	! c cbmflast(i) = 0.
	1745	! c enddo
	1746	! cc
	1747	! c do k= 1,nl
	1748	! c do i = 1,ncum
	1749	! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
	1750	! c cbmflast(i) = cbmflast(i)+M(i,k)
	1751	! c ENDIF
	1752	! c enddo
	1753	! c enddo
	1754	! cc
	1755	! c do i = 1,ncum
	1756	! c IF (cbmflast(i) .lt. 1.e-6) THEN
	1757	! c iflag(i) = 3
	1758	! c ENDIF
	1759	! c enddo
	1760	! cc
	1761	! c do k= 1,nl
	1762	! c do i = 1,ncum
	1763	! c IF (iflag(i) .ge. 3) THEN
	1764	! c M(i,k) = 0.
	1765	! c sig(i,k) = 0.
	1766	! c w0(i,k) = 0.
	1767	! c ENDIF
	1768	! c enddo
	1769	! c enddo
	1770	! cc
	1771	!! cape=0.0
	1772	!! do 98 i=icb+1,inb
	1773	!! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
	1774	!! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
	1775	!! dcape=rrdbuoy(i-1)deltap/p(i-1)
	1776	!! dlnp=deltap/p(i-1)
	1777	!! cape=max(0.0,cape)
	1778	!! sigold=sig(i)
[1849]	1779
[2007]	1780	!! dtmin=100.0
	1781	!! do 97 j=icb,i-1
	1782	!! dtmin=amin1(dtmin,buoy(j))
	1783	!! 97 continue
[1849]	1784
[2007]	1785	!! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
	1786	!! sig(i)=max(sig(i),0.0)
	1787	!! sig(i)=amin1(sig(i),0.01)
	1788	!! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
	1789	!! w=(1.-beta)facsqrt(cape)+beta*w0(i)
	1790	!! amu=0.5(sig(i)+sigold)w
	1791	!! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
	1792	!! w0(i)=w
	1793	!! 98 continue
	1794	!! w0(icb)=0.5*w0(icb+1)
	1795	!! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
	1796	!! sig(icb)=sig(icb+1)
	1797	!! sig(icb-1)=sig(icb)
[1849]	1798
[1992]	1799	RETURN
	1800	END SUBROUTINE cv3_closure
[1849]	1801
[2007]	1802	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
	1803	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
	1804	unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
	1805	ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
[1992]	1806	IMPLICIT NONE
[1849]	1807
[2007]	1808	! ---------------------------------------------------------------------
	1809	! a faire:
	1810	! - vectorisation de la partie normalisation des flux (do 789...)
	1811	! ---------------------------------------------------------------------
[524]	1812
[1992]	1813	include "cvthermo.h"
	1814	include "cv3param.h"
	1815	include "cvflag.h"
[524]	1816
[2007]	1817	!inputs:
[2253]	1818	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
	1819	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
	1820	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
	1821	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
	1822	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	1823	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
	1824	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
	1825	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
	1826	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, h, hp
	1827	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf, frac
	1828	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp, ep, clw
	1829	REAL, DIMENSION (nloc, na), INTENT (IN) :: m ! input of convect3
[524]	1830
[2007]	1831	!outputs:
[2253]	1832	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: ment, qent
	1833	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
	1834	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: sij, elij
	1835	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
	1836	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: ments, qents
	1837	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
[524]	1838
[2007]	1839	!local variables:
[1992]	1840	INTEGER i, j, k, il, im, jm
	1841	INTEGER num1, num2
	1842	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
	1843	REAL alt, smid, sjmin, sjmax, delp, delm
	1844	REAL asij(nloc), smax(nloc), scrit(nloc)
	1845	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
[2253]	1846	REAL sigij(nloc, nd, nd)
[1992]	1847	REAL wgh
	1848	REAL zm(nloc, na)
	1849	LOGICAL lwork(nloc)
[524]	1850
[2007]	1851	! =====================================================================
	1852	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	1853	! =====================================================================
[524]	1854
[2007]	1855	! ori do 360 i=1,ncum*nlp
[1992]	1856	DO j = 1, nl
	1857	DO i = 1, ncum
	1858	nent(i, j) = 0
[2007]	1859	! in convect3, m is computed in cv3_closure
	1860	! ori m(i,1)=0.0
[1992]	1861	END DO
	1862	END DO
[524]	1863
[2007]	1864	! ori do 400 k=1,nlp
	1865	! ori do 390 j=1,nlp
[1992]	1866	DO j = 1, nl
	1867	DO k = 1, nl
	1868	DO i = 1, ncum
	1869	qent(i, k, j) = rr(i, j)
	1870	uent(i, k, j) = u(i, j)
	1871	vent(i, k, j) = v(i, j)
	1872	elij(i, k, j) = 0.0
[2007]	1873	!ym ment(i,k,j)=0.0
	1874	!ym sij(i,k,j)=0.0
[1992]	1875	END DO
	1876	END DO
	1877	END DO
[524]	1878
[2007]	1879	!ym
[1992]	1880	ment(1:ncum, 1:nd, 1:nd) = 0.0
	1881	sij(1:ncum, 1:nd, 1:nd) = 0.0
[524]	1882
[2007]	1883	!AC! do k=1,ntra
	1884	!AC! do j=1,nd ! instead nlp
	1885	!AC! do i=1,nd ! instead nlp
	1886	!AC! do il=1,ncum
	1887	!AC! traent(il,i,j,k)=tra(il,j,k)
	1888	!AC! enddo
	1889	!AC! enddo
	1890	!AC! enddo
	1891	!AC! enddo
[1992]	1892	zm(:, :) = 0.
[524]	1893
[2007]	1894	! =====================================================================
	1895	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	1896	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	1897	! --- FRACTION (sij)
	1898	! =====================================================================
[524]	1899
[1992]	1900	DO i = minorig + 1, nl
[524]	1901
[1992]	1902	DO j = minorig, nl
	1903	DO il = 1, ncum
[2007]	1904	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
[524]	1905
[1992]	1906	rti = qnk(il) - ep(il, i)*clw(il, i)
	1907	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
[524]	1908
	1909
[1992]	1910	IF (cvflag_ice) THEN
[2007]	1911	! print*,cvflag_ice,'cvflag_ice dans do 700'
[1992]	1912	IF (t(il,j)<=263.15) THEN
[2007]	1913	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
	1914	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
[1992]	1915	END IF
	1916	END IF
[524]	1917
[1992]	1918	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
	1919	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
	1920	dei = denom
	1921	IF (abs(dei)<0.01) dei = 0.01
	1922	sij(il, i, j) = anum/dei
	1923	sij(il, i, i) = 1.0
	1924	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
	1925	altem = altem/bf2
	1926	cwat = clw(il, j)*(1.-ep(il,j))
	1927	stemp = sij(il, i, j)
	1928	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
[524]	1929
[1992]	1930	IF (cvflag_ice) THEN
[2007]	1931	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
[1992]	1932	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
	1933	ELSE
	1934	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
	1935	denom = denom + lv(il, j)*(rr(il,i)-rti)
	1936	END IF
[524]	1937
[1992]	1938	IF (abs(denom)<0.01) denom = 0.01
	1939	sij(il, i, j) = anum/denom
[2007]	1940	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
[1992]	1941	altem = altem - (bf2-1.)*cwat
	1942	END IF
	1943	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
	1944	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
[2007]	1945	uent(il, i, j) = sij(il, i, j)u(il, i) + (1.-sij(il,i,j))unk(il)
	1946	vent(il, i, j) = sij(il, i, j)v(il, i) + (1.-sij(il,i,j))vnk(il)
	1947	!!!! do k=1,ntra
	1948	!!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1949	!!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1950	!!!! end do
[1992]	1951	elij(il, i, j) = altem
	1952	elij(il, i, j) = max(0.0, elij(il,i,j))
	1953	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
	1954	nent(il, i) = nent(il, i) + 1
	1955	END IF
	1956	sij(il, i, j) = max(0.0, sij(il,i,j))
	1957	sij(il, i, j) = amin1(1.0, sij(il,i,j))
	1958	END IF ! new
	1959	END DO
	1960	END DO
[524]	1961
[2007]	1962	!AC! do k=1,ntra
	1963	!AC! do j=minorig,nl
	1964	!AC! do il=1,ncum
	1965	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
	1966	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
	1967	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1968	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1969	!AC! endif
	1970	!AC! enddo
	1971	!AC! enddo
	1972	!AC! enddo
[524]	1973
	1974
[2007]	1975	! * if no air can entrain at level i assume that updraft detrains *
	1976	! * at that level and calculate detrained air flux and properties *
[524]	1977
	1978
[2007]	1979	! @ do 170 i=icb(il),inb(il)
[524]	1980
[1992]	1981	DO il = 1, ncum
	1982	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
[2007]	1983	! @ if(nent(il,i).eq.0)then
[1992]	1984	ment(il, i, i) = m(il, i)
	1985	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	1986	uent(il, i, i) = unk(il)
	1987	vent(il, i, i) = vnk(il)
	1988	elij(il, i, i) = clw(il, i)
[2007]	1989	! MAF sij(il,i,i)=1.0
[1992]	1990	sij(il, i, i) = 0.0
	1991	END IF
	1992	END DO
	1993	END DO
[879]	1994
[2007]	1995	!AC! do j=1,ntra
	1996	!AC! do i=minorig+1,nl
	1997	!AC! do il=1,ncum
	1998	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
	1999	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
	2000	!AC! endif
	2001	!AC! enddo
	2002	!AC! enddo
	2003	!AC! enddo
[879]	2004
[1992]	2005	DO j = minorig, nl
	2006	DO i = minorig, nl
	2007	DO il = 1, ncum
[2007]	2008	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
[1992]	2009	sigij(il, i, j) = sij(il, i, j)
	2010	END IF
	2011	END DO
	2012	END DO
	2013	END DO
[2007]	2014	! @ enddo
[879]	2015
[2007]	2016	! @170 continue
[524]	2017
[2007]	2018	! =====================================================================
	2019	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	2020	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	2021	! =====================================================================
[970]	2022
[1992]	2023	CALL zilch(asum, nloc*nd)
	2024	CALL zilch(csum, nloc*nd)
	2025	CALL zilch(csum, nloc*nd)
[524]	2026
[1992]	2027	DO il = 1, ncum
	2028	lwork(il) = .FALSE.
	2029	END DO
[524]	2030
[1992]	2031	DO i = minorig + 1, nl
[524]	2032
[1992]	2033	num1 = 0
	2034	DO il = 1, ncum
	2035	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
	2036	END DO
	2037	IF (num1<=0) GO TO 789
[524]	2038
[879]	2039
[1992]	2040	DO il = 1, ncum
	2041	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	2042	lwork(il) = (nent(il,i)/=0)
	2043	qp = qnk(il) - ep(il, i)*clw(il, i)
[524]	2044
[1992]	2045	IF (cvflag_ice) THEN
[524]	2046
[2007]	2047	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
	2048	(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
	2049	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
	2050	(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
[1992]	2051	ELSE
[879]	2052
[1992]	2053	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
[2007]	2054	(cpv-cpd)t(il, i)(qp-rr(il,i))
[1992]	2055	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
[2007]	2056	(cpd-cpv)t(il, i)(rr(il,i)-qp)
[1992]	2057	END IF
[524]	2058
[1992]	2059	IF (abs(denom)<0.01) denom = 0.01
	2060	scrit(il) = anum/denom
	2061	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
	2062	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
	2063	smax(il) = 0.0
	2064	asij(il) = 0.0
	2065	END IF
	2066	END DO
[524]	2067
[1992]	2068	DO j = nl, minorig, -1
[524]	2069
[1992]	2070	num2 = 0
	2071	DO il = 1, ncum
[2007]	2072	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	2073	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	2074	lwork(il)) num2 = num2 + 1
[1992]	2075	END DO
	2076	IF (num2<=0) GO TO 175
[524]	2077
[1992]	2078	DO il = 1, ncum
[2007]	2079	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
	2080	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
	2081	lwork(il)) THEN
[524]	2082
[1992]	2083	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
	2084	wgh = 1.0
	2085	IF (j>i) THEN
	2086	sjmax = max(sij(il,i,j+1), smax(il))
	2087	sjmax = amin1(sjmax, scrit(il))
	2088	smax(il) = max(sij(il,i,j), smax(il))
	2089	sjmin = max(sij(il,i,j-1), smax(il))
	2090	sjmin = amin1(sjmin, scrit(il))
	2091	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
	2092	smid = amin1(sij(il,i,j), scrit(il))
	2093	ELSE
	2094	sjmax = max(sij(il,i,j+1), scrit(il))
	2095	smid = max(sij(il,i,j), scrit(il))
	2096	sjmin = 0.0
	2097	IF (j>1) sjmin = sij(il, i, j-1)
	2098	sjmin = max(sjmin, scrit(il))
	2099	END IF
	2100	delp = abs(sjmax-smid)
	2101	delm = abs(sjmin-smid)
	2102	asij(il) = asij(il) + wgh*(delp+delm)
	2103	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
	2104	END IF
	2105	END IF
	2106	END DO
[524]	2107
[1992]	2108	175 END DO
[524]	2109
[1992]	2110	DO il = 1, ncum
	2111	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	2112	asij(il) = max(1.0E-16, asij(il))
	2113	asij(il) = 1.0/asij(il)
	2114	asum(il, i) = 0.0
	2115	bsum(il, i) = 0.0
	2116	csum(il, i) = 0.0
	2117	END IF
	2118	END DO
[524]	2119
[1992]	2120	DO j = minorig, nl
	2121	DO il = 1, ncum
[2007]	2122	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2123	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[1992]	2124	ment(il, i, j) = ment(il, i, j)*asij(il)
	2125	END IF
	2126	END DO
	2127	END DO
[524]	2128
[1992]	2129	DO j = minorig, nl
	2130	DO il = 1, ncum
[2007]	2131	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2132	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[1992]	2133	asum(il, i) = asum(il, i) + ment(il, i, j)
	2134	ment(il, i, j) = ment(il, i, j)*sig(il, j)
	2135	bsum(il, i) = bsum(il, i) + ment(il, i, j)
	2136	END IF
	2137	END DO
	2138	END DO
[1849]	2139
[1992]	2140	DO il = 1, ncum
	2141	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
	2142	bsum(il, i) = max(bsum(il,i), 1.0E-16)
	2143	bsum(il, i) = 1.0/bsum(il, i)
	2144	END IF
	2145	END DO
[1849]	2146
[1992]	2147	DO j = minorig, nl
	2148	DO il = 1, ncum
[2007]	2149	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2150	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[1992]	2151	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
	2152	END IF
	2153	END DO
	2154	END DO
[879]	2155
[1992]	2156	DO j = minorig, nl
	2157	DO il = 1, ncum
[2007]	2158	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2159	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[1992]	2160	csum(il, i) = csum(il, i) + ment(il, i, j)
	2161	END IF
	2162	END DO
	2163	END DO
[1849]	2164
[1992]	2165	DO il = 1, ncum
	2166	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
	2167	csum(il,i)<m(il,i)) THEN
	2168	nent(il, i) = 0
	2169	ment(il, i, i) = m(il, i)
	2170	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	2171	uent(il, i, i) = unk(il)
	2172	vent(il, i, i) = vnk(il)
	2173	elij(il, i, i) = clw(il, i)
[2007]	2174	! MAF sij(il,i,i)=1.0
[1992]	2175	sij(il, i, i) = 0.0
	2176	END IF
	2177	END DO ! il
[1849]	2178
[2007]	2179	!AC! do j=1,ntra
	2180	!AC! do il=1,ncum
	2181	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	2182	!AC! : .and. csum(il,i).lt.m(il,i) ) then
	2183	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
	2184	!AC! endif
	2185	!AC! enddo
	2186	!AC! enddo
[1992]	2187	789 END DO
[879]	2188
[2007]	2189	! MAF: renormalisation de MENT
[1992]	2190	CALL zilch(zm, nloc*na)
	2191	DO jm = 1, nd
	2192	DO im = 1, nd
	2193	DO il = 1, ncum
	2194	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
	2195	END DO
	2196	END DO
	2197	END DO
[524]	2198
[1992]	2199	DO jm = 1, nd
	2200	DO im = 1, nd
	2201	DO il = 1, ncum
	2202	IF (zm(il,im)/=0.) THEN
	2203	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
	2204	END IF
	2205	END DO
	2206	END DO
	2207	END DO
[524]	2208
[1992]	2209	DO jm = 1, nd
	2210	DO im = 1, nd
	2211	DO il = 1, ncum
	2212	qents(il, im, jm) = qent(il, im, jm)
	2213	ments(il, im, jm) = ment(il, im, jm)
	2214	END DO
	2215	END DO
	2216	END DO
[524]	2217
[1992]	2218	RETURN
	2219	END SUBROUTINE cv3_mixing
[879]	2220
[2007]	2221	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
	2222	t, rr, rs, gz, u, v, tra, p, ph, &
	2223	th, tv, lv, lf, cpn, ep, sigp, clw, &
	2224	m, ment, elij, delt, plcl, coef_clos, &
	2225	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
	2226	faci, b, sigd, &
	2227	wdtrainA, wdtrainM) ! RomP
[1992]	2228	IMPLICIT NONE
[879]	2229
	2230
[1992]	2231	include "cvthermo.h"
	2232	include "cv3param.h"
	2233	include "cvflag.h"
[2287]	2234	include "nuage.h"
[524]	2235
[2007]	2236	!inputs:
[1992]	2237	INTEGER ncum, nd, na, ntra, nloc
	2238	INTEGER icb(nloc), inb(nloc)
	2239	REAL delt, plcl(nloc)
	2240	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd), gz(nloc, na)
	2241	REAL u(nloc, nd), v(nloc, nd)
	2242	REAL tra(nloc, nd, ntra)
	2243	REAL p(nloc, nd), ph(nloc, nd+1)
	2244	REAL ep(nloc, na), sigp(nloc, na), clw(nloc, na)
	2245	REAL th(nloc, na), tv(nloc, na), lv(nloc, na), cpn(nloc, na)
	2246	REAL lf(nloc, na)
	2247	REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
	2248	REAL coef_clos(nloc)
[524]	2249
[2007]	2250	!input/output
[1992]	2251	INTEGER iflag(nloc)
[524]	2252
[2007]	2253	!outputs:
[1992]	2254	REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
	2255	REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
	2256	REAL ice(nloc, na), fondue(nloc, na), faci(nloc, na)
	2257	REAL trap(nloc, na, ntra)
	2258	REAL b(nloc, na), sigd(nloc)
[2007]	2259	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
	2260	! de l ascendance adiabatique et des flux melanges Pa et Pm.
	2261	! Distinction des wdtrain
	2262	! Pa = wdtrainA Pm = wdtrainM
	2263	REAL wdtrainA(nloc, na), wdtrainM(nloc, na)
[879]	2264
[2007]	2265	!local variables
[1992]	2266	INTEGER i, j, k, il, num1, ndp1
	2267	REAL tinv, delti, coef
	2268	REAL awat, afac, afac1, afac2, bfac
	2269	REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
	2270	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
	2271	REAL ampmax, thaw
	2272	REAL tevap(nloc)
	2273	REAL lvcp(nloc, na), lfcp(nloc, na)
	2274	REAL h(nloc, na), hm(nloc, na)
	2275	REAL frac(nloc, na)
	2276	REAL fraci(nloc, na), prec(nloc, na)
	2277	REAL wdtrain(nloc)
	2278	LOGICAL lwork(nloc), mplus(nloc)
[524]	2279
	2280
[2007]	2281	! ------------------------------------------------------
[524]	2282
[1992]	2283	delti = 1./delt
	2284	tinv = 1./3.
[524]	2285
[1992]	2286	mp(:, :) = 0.
[524]	2287
[1992]	2288	DO i = 1, nl
	2289	DO il = 1, ncum
	2290	mp(il, i) = 0.0
	2291	rp(il, i) = rr(il, i)
	2292	up(il, i) = u(il, i)
	2293	vp(il, i) = v(il, i)
	2294	wt(il, i) = 0.001
	2295	water(il, i) = 0.0
	2296	frac(il, i) = 0.0
	2297	faci(il, i) = 0.0
	2298	fraci(il, i) = 0.0
	2299	ice(il, i) = 0.0
	2300	prec(il, i) = 0.0
	2301	fondue(il, i) = 0.0
	2302	evap(il, i) = 0.0
	2303	b(il, i) = 0.0
	2304	lvcp(il, i) = lv(il, i)/cpn(il, i)
	2305	lfcp(il, i) = lf(il, i)/cpn(il, i)
	2306	END DO
	2307	END DO
[2007]	2308	!AC! do k=1,ntra
	2309	!AC! do i=1,nd
	2310	!AC! do il=1,ncum
	2311	!AC! trap(il,i,k)=tra(il,i,k)
	2312	!AC! enddo
	2313	!AC! enddo
	2314	!AC! enddo
	2315	!! RomP >>>
[1992]	2316	DO i = 1, nd
	2317	DO il = 1, ncum
[2007]	2318	wdtrainA(il, i) = 0.0
	2319	wdtrainM(il, i) = 0.0
[1992]	2320	END DO
	2321	END DO
[2007]	2322	!! RomP <<<
[524]	2323
[2007]	2324	! * check whether ep(inb)=0, if so, skip precipitating *
	2325	! * downdraft calculation *
[524]	2326
	2327
[1992]	2328	DO il = 1, ncum
[2007]	2329	!! lwork(il)=.TRUE.
	2330	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
[1992]	2331	lwork(il) = ep(il, inb(il)) >= 0.0001
	2332	END DO
[524]	2333
[2007]	2334	! *** Set the fractionnal area sigd of precipitating downdraughts
[1992]	2335	DO il = 1, ncum
	2336	sigd(il) = sigdz*coef_clos(il)
	2337	END DO
[524]	2338
	2339
[2007]	2340	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	2341	!
	2342	! * begin downdraft loop *
	2343	!
	2344	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[524]	2345
[1992]	2346	DO i = nl + 1, 1, -1
[524]	2347
[1992]	2348	num1 = 0
	2349	DO il = 1, ncum
	2350	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
	2351	END DO
	2352	IF (num1<=0) GO TO 400
[1146]	2353
[1992]	2354	CALL zilch(wdtrain, ncum)
[1146]	2355
	2356
[2007]	2357	! * integrate liquid water equation to find condensed water *
	2358	! * and condensed water flux *
	2359	!
	2360	!
	2361	! * calculate detrained precipitation *
[524]	2362
[1992]	2363	DO il = 1, ncum
	2364	IF (i<=inb(il) .AND. lwork(il)) THEN
	2365	IF (cvflag_grav) THEN
	2366	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
[2007]	2367	wdtrainA(il, i) = wdtrain(il)/grav ! Pa RomP
[1992]	2368	ELSE
	2369	wdtrain(il) = 10.0ep(il, i)m(il, i)*clw(il, i)
[2007]	2370	wdtrainA(il, i) = wdtrain(il)/10. ! Pa RomP
[1992]	2371	END IF
	2372	END IF
	2373	END DO
[524]	2374
[1992]	2375	IF (i>1) THEN
	2376	DO j = 1, i - 1
	2377	DO il = 1, ncum
	2378	IF (i<=inb(il) .AND. lwork(il)) THEN
	2379	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
	2380	awat = max(awat, 0.0)
	2381	IF (cvflag_grav) THEN
	2382	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
[2007]	2383	wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP
[1992]	2384	ELSE
	2385	wdtrain(il) = wdtrain(il) + 10.0awatment(il, j, i)
[2007]	2386	wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i) ! Pm RomP
[1992]	2387	END IF
	2388	END IF
	2389	END DO
	2390	END DO
	2391	END IF
[524]	2392
[1146]	2393
[2007]	2394	! * find rain water and evaporation using provisional *
	2395	! * estimates of rp(i)and rp(i-1) *
[1146]	2396
[1650]	2397
[1992]	2398	DO il = 1, ncum
	2399	IF (i<=inb(il) .AND. lwork(il)) THEN
[524]	2400
[1992]	2401	wt(il, i) = 45.0
[524]	2402
[1992]	2403	IF (cvflag_ice) THEN
	2404	frac(il, inb(il)) = 1. - (t(il,inb(il))-243.15)/(263.15-243.15)
	2405	frac(il, inb(il)) = min(max(frac(il,inb(il)),0.), 1.)
	2406	fraci(il, inb(il)) = frac(il, inb(il))
	2407	ELSE
	2408	CONTINUE
	2409	END IF
[879]	2410
[1992]	2411	IF (i<inb(il)) THEN
[524]	2412
[1992]	2413	IF (cvflag_ice) THEN
[2287]	2414	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
	2415	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
[1992]	2416	thaw = min(max(thaw,0.0), 1.0)
	2417	frac(il, i) = frac(il, i)*(1.-thaw)
	2418	ELSE
	2419	CONTINUE
	2420	END IF
[524]	2421
[2007]	2422	rp(il, i) = rp(il, i+1) + &
	2423	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
[1992]	2424	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
	2425	END IF
	2426	fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
	2427	fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
	2428	rp(il, i) = max(rp(il,i), 0.0)
	2429	rp(il, i) = amin1(rp(il,i), rs(il,i))
	2430	rp(il, inb(il)) = rr(il, inb(il))
[524]	2431
[1992]	2432	IF (i==1) THEN
	2433	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
	2434	IF (cvflag_ice) THEN
[2007]	2435	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
[1992]	2436	END IF
	2437	ELSE
[2007]	2438	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]	2439	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
	2440	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
	2441	rp(il, i-1) = max(rp(il,i-1), 0.0)
[2007]	2442	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
	2443	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]	2444	afac = 0.5*(afac1+afac2)
	2445	END IF
	2446	IF (i==inb(il)) afac = 0.0
	2447	afac = max(afac, 0.0)
	2448	bfac = 1./(sigd(il)*wt(il,i))
[524]	2449
[2007]	2450	!JYG1
	2451	! cc sigt=1.0
	2452	! cc if(i.ge.icb)sigt=sigp(i)
	2453	! prise en compte de la variation progressive de sigt dans
	2454	! les couches icb et icb-1:
	2455	! pour plcl<ph(i+1), pr1=0 & pr2=1
	2456	! pour plcl>ph(i), pr1=1 & pr2=0
	2457	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
	2458	! sur le nuage, et pr2 est la proportion sous la base du
	2459	! nuage.
[1992]	2460	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
	2461	pr1 = max(0., min(1.,pr1))
	2462	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
	2463	pr2 = max(0., min(1.,pr2))
	2464	sigt = sigp(il, i)*pr1 + pr2
[2007]	2465	!JYG2
[524]	2466
[2007]	2467	!JYG----
	2468	! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
	2469	! c6 = water(il,i+1) + wdtrain(il)*bfac
	2470	! c6 = prec(il,i+1) + wdtrain(il)*bfac
	2471	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
	2472	! evap(il,i)=sigtafacrevap
	2473	! water(il,i)=revap*revap
	2474	! prec(il,i)=revap*revap
	2475	!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
	2476	!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
	2477	!!---end jyg---
[879]	2478
[2007]	2479	! --------retour à la formulation originale d'Emanuel.
[1992]	2480	IF (cvflag_ice) THEN
[524]	2481
[2007]	2482	! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
	2483	! c6=prec(il,i+1)+bfac*wdtrain(il) &
	2484	! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
	2485	! if(c6.gt.0.0)then
	2486	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
[524]	2487
[2007]	2488	!JAM Attention: evap=sigt*E
	2489	! Modification: evap devient l'évaporation en milieu de couche
	2490	! car nécessaire dans cv3_yield
	2491	! Du coup, il faut modifier pas mal d'équations...
	2492	! et l'expression de afac qui devient afac1
	2493	! revap=sqrt((prec(i+1)+prec(i))/2)
[524]	2494
[1992]	2495	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
	2496	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
[2007]	2497	! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
	2498	! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
	2499	! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
[1992]	2500	IF (c6>b6*b6+1.E-20) THEN
	2501	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
	2502	ELSE
	2503	revap = (-b6+sqrt(b6b6+4.c6))/2.
	2504	END IF
	2505	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
[2007]	2506	! print*,prec(il,i),'neige'
[524]	2507
[2007]	2508	!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
	2509	! c evap(il,i)=sigtafacrevap
	2510	! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
	2511	! Ici,l'evaporation evap est simplement calculee par l'equation de
	2512	! conservation.
	2513	! prec(il,i)=revap*revap
	2514	! else
	2515	!JYG---- Correction : si c6 <= 0, water(il,i)=0.
	2516	! prec(il,i)=0.
	2517	! endif
[524]	2518
[2007]	2519	!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
	2520	! moins [tt ce qui sort de la couche i]
	2521	! print *, 'evap avec ice'
	2522	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
	2523	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
[524]	2524
[2007]	2525	d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
[1992]	2526	e6 = bfac*wdtrain(il)
	2527	f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
[524]	2528
[2287]	2529	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
	2530	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
[1992]	2531	thaw = min(max(thaw,0.0), 1.0)
	2532	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
	2533	water(il, i) = max(water(il,i), 0.)
	2534	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
	2535	ice(il, i) = max(ice(il,i), 0.)
	2536	fondue(il, i) = ice(il, i)*thaw
	2537	water(il, i) = water(il, i) + fondue(il, i)
	2538	ice(il, i) = ice(il, i) - fondue(il, i)
[524]	2539
[1992]	2540	IF (water(il,i)+ice(il,i)<1.E-30) THEN
	2541	faci(il, i) = 0.
	2542	ELSE
	2543	faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
	2544	END IF
[524]	2545
[2007]	2546	! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
	2547	! water(il,i)=max(water(il,i),0.)
	2548	! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
	2549	! ice(il,i)=max(ice(il,i),0.)
	2550	! fondue(il,i)=ice(il,i)*thaw
	2551	! water(il,i)=water(il,i)+fondue(il,i)
	2552	! ice(il,i)=ice(il,i)-fondue(il,i)
	2553
	2554	! if((water(il,i)+ice(il,i)).lt.1.e-30)then
	2555	! faci(il,i)=0.
	2556	! else
	2557	! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
	2558	! endif
[524]	2559
[1992]	2560	ELSE
	2561	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
[2007]	2562	c6 = water(il, i+1) + bfac*wdtrain(il) - &
	2563	50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
[1992]	2564	IF (c6>0.0) THEN
	2565	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
	2566	water(il, i) = revap*revap
	2567	ELSE
	2568	water(il, i) = 0.
	2569	END IF
[2007]	2570	! print *, 'evap sans ice'
	2571	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
	2572	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
[524]	2573
[1992]	2574	END IF
	2575	END IF !(i.le.inb(il) .and. lwork(il))
	2576	END DO
[2007]	2577	! ----------------------------------------------------------------
[524]	2578
[2007]	2579	! cc
	2580	! * calculate precipitating downdraft mass flux under *
	2581	! * hydrostatic approximation *
[524]	2582
[1992]	2583	DO il = 1, ncum
	2584	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
[524]	2585
[1992]	2586	tevap(il) = max(0.0, evap(il,i))
	2587	delth = max(0.001, (th(il,i)-th(il,i-1)))
	2588	IF (cvflag_ice) THEN
	2589	IF (cvflag_grav) THEN
[2007]	2590	mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
	2591	(p(il,i-1)-p(il,i))/delth + &
	2592	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
	2593	(p(il,i-1)-p(il,i))/delth + &
	2594	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
	2595	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
[1992]	2596	ELSE
[2007]	2597	mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
	2598	(p(il,i-1)-p(il,i))/delth + &
	2599	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
	2600	(p(il,i-1)-p(il,i))/delth + &
	2601	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
	2602	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
[524]	2603
[1992]	2604	END IF
	2605	ELSE
	2606	IF (cvflag_grav) THEN
	2607	mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
[2007]	2608	(p(il,i-1)-p(il,i))/delth
[1992]	2609	ELSE
	2610	mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
[2007]	2611	(p(il,i-1)-p(il,i))/delth
[1992]	2612	END IF
[524]	2613
[1992]	2614	END IF
[879]	2615
[1992]	2616	END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
	2617	END DO
[2007]	2618	! ----------------------------------------------------------------
[524]	2619
[2007]	2620	! * if hydrostatic assumption fails, *
	2621	! * solve cubic difference equation for downdraft theta *
	2622	! * and mass flux from two simultaneous differential eqns *
[524]	2623
[1992]	2624	DO il = 1, ncum
	2625	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
[1742]	2626
[1992]	2627	amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
[2007]	2628	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
[1992]	2629	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
[1742]	2630
[1992]	2631	IF (amp2>(0.1*amfac)) THEN
	2632	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
[2007]	2633	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
	2634	(lvcp(il,i)sigd(il)th(il,i))
[1992]	2635	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
[1742]	2636
[1992]	2637	IF (cvflag_ice) THEN
	2638	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
[2007]	2639	50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
	2640	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
[1992]	2641	ELSE
[1774]	2642
[1992]	2643	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
[2007]	2644	50.(p(il,i-1)-p(il,i))xf*tevap(il)
[1992]	2645	END IF
[1742]	2646
[1992]	2647	fac2 = 1.0
	2648	IF (bf<0.0) fac2 = -1.0
	2649	bf = abs(bf)
	2650	ur = 0.25bfbf - afafaftinvtinv*tinv
	2651	IF (ur>=0.0) THEN
	2652	sru = sqrt(ur)
	2653	fac = 1.0
	2654	IF ((0.5*bf-sru)<0.0) fac = -1.0
	2655	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
[2007]	2656	fac(abs(0.5bf-sru))**tinv
[1992]	2657	ELSE
	2658	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
	2659	IF (fac2<0.0) d = 3.14159 - d
	2660	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
	2661	END IF
	2662	mp(il, i) = max(0.0, mp(il,i))
[524]	2663
[1992]	2664	IF (cvflag_ice) THEN
	2665	IF (cvflag_grav) THEN
[2007]	2666	!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
	2667	! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
	2668	! Et il faut bien revoir les facteurs 100.
	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	ELSE
[2007]	2677	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
	2678	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
	2679	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
	2680	(ph(il,i)-ph(il,i+1))) / &
	2681	(mp(il,i)+sigd(il)*0.1) - &
	2682	10.0(th(il,i)-th(il,i-1))t(il, i) / &
	2683	(lvcp(il,i)sigd(il)th(il,i))
[1992]	2684	END IF
	2685	ELSE
	2686	IF (cvflag_grav) THEN
[2007]	2687	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
	2688	(mp(il,i)+sigd(il)*0.1) - &
	2689	10.0(th(il,i)-th(il,i-1))t(il, i) / &
	2690	(lvcp(il,i)sigd(il)th(il,i))
[1992]	2691	ELSE
[2007]	2692	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
	2693	(mp(il,i)+sigd(il)*0.1) - &
	2694	10.0(th(il,i)-th(il,i-1))t(il, i) / &
	2695	(lvcp(il,i)sigd(il)th(il,i))
[1992]	2696	END IF
	2697	END IF
	2698	b(il, i-1) = max(b(il,i-1), 0.0)
[524]	2699
[1992]	2700	END IF !(amp2.gt.(0.1*amfac))
[524]	2701
[2007]	2702	! * limit magnitude of mp(i) to meet cfl condition *
[524]	2703
[1992]	2704	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
	2705	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
	2706	ampmax = min(ampmax, amp2)
	2707	mp(il, i) = min(mp(il,i), ampmax)
[524]	2708
[2007]	2709	! * force mp to decrease linearly to zero *
	2710	! * between cloud base and the surface *
[524]	2711
	2712
[2007]	2713	! c if(p(il,i).gt.p(il,icb(il)))then
	2714	! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
	2715	! c endif
[1992]	2716	IF (ph(il,i)>0.9*plcl(il)) THEN
	2717	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
	2718	END IF
[524]	2719
[1992]	2720	END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
	2721	END DO
[2007]	2722	! ----------------------------------------------------------------
[524]	2723
[2007]	2724	! * find mixing ratio of precipitating downdraft *
[524]	2725
[1992]	2726	DO il = 1, ncum
	2727	IF (i<inb(il) .AND. lwork(il)) THEN
	2728	mplus(il) = mp(il, i) > mp(il, i+1)
	2729	END IF ! (i.lt.inb(il) .and. lwork(il))
	2730	END DO
	2731
	2732	DO il = 1, ncum
	2733	IF (i<inb(il) .AND. lwork(il)) THEN
	2734
	2735	rp(il, i) = rr(il, i)
	2736
	2737	IF (mplus(il)) THEN
	2738
	2739	IF (cvflag_grav) THEN
[2007]	2740	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
	2741	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
[1992]	2742	ELSE
[2007]	2743	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
	2744	5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
[1992]	2745	END IF
	2746	rp(il, i) = rp(il, i)/mp(il, i)
[2007]	2747	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
[1992]	2748	up(il, i) = up(il, i)/mp(il, i)
[2007]	2749	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
[1992]	2750	vp(il, i) = vp(il, i)/mp(il, i)
	2751
	2752	ELSE ! if (mplus(il))
	2753
	2754	IF (mp(il,i+1)>1.0E-16) THEN
	2755	IF (cvflag_grav) THEN
[2007]	2756	rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
	2757	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
[1992]	2758	ELSE
[2007]	2759	rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
	2760	(evap(il,i+1)+evap(il,i))/mp(il, i+1)
[1992]	2761	END IF
	2762	up(il, i) = up(il, i+1)
	2763	vp(il, i) = vp(il, i+1)
	2764	END IF ! (mp(il,i+1).gt.1.0e-16)
	2765	END IF ! (mplus(il)) else if (.not.mplus(il))
	2766
	2767	rp(il, i) = amin1(rp(il,i), rs(il,i))
	2768	rp(il, i) = max(rp(il,i), 0.0)
	2769
	2770	END IF ! (i.lt.inb(il) .and. lwork(il))
	2771	END DO
[2007]	2772	! ----------------------------------------------------------------
[1992]	2773
[2007]	2774	! * find tracer concentrations in precipitating downdraft *
[1992]	2775
[2007]	2776	!AC! do j=1,ntra
	2777	!AC! do il = 1,ncum
	2778	!AC! if (i.lt.inb(il) .and. lwork(il)) then
	2779	!AC!c
	2780	!AC! if(mplus(il))then
	2781	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
	2782	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
	2783	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
	2784	!AC! else ! if (mplus(il))
	2785	!AC! if(mp(il,i+1).gt.1.0e-16)then
	2786	!AC! trap(il,i,j)=trap(il,i+1,j)
	2787	!AC! endif
	2788	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
	2789	!AC!c
	2790	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
	2791	!AC! enddo
	2792	!AC! end do
[1992]	2793
	2794	400 END DO
[2007]	2795	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1992]	2796
[2007]	2797	! * end of downdraft loop *
[1992]	2798
[2007]	2799	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1992]	2800
	2801
	2802	RETURN
	2803	END SUBROUTINE cv3_unsat
	2804
[2007]	2805	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
	2806	icb, inb, delt, &
	2807	t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
	2808	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
	2809	ep, clw, m, tp, mp, rp, up, vp, trap, &
	2810	wt, water, ice, evap, fondue, faci, b, sigd, &
	2811	ment, qent, hent, iflag_mix, uent, vent, &
	2812	nent, elij, traent, sig, &
	2813	tv, tvp, wghti, &
[2306]	2814	iflag, precip, Vprecip, Vprecipi, & ! jyg: Vprecipi
	2815	ft, fr, fu, fv, ftra, & ! jyg
[2007]	2816	cbmf, upwd, dnwd, dnwd0, ma, mip, &
[2259]	2817	!! tls, tps, ! useless . jyg
	2818	qcondc, wd, &
[2205]	2819	ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
[1992]	2820
	2821	IMPLICIT NONE
	2822
	2823	include "cvthermo.h"
	2824	include "cv3param.h"
	2825	include "cvflag.h"
	2826	include "conema3.h"
	2827
[2007]	2828	!inputs:
[2327]	2829	INTEGER, INTENT (IN) :: iflag_mix
	2830	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
	2831	LOGICAL, INTENT (IN) :: ok_conserv_q
	2832	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
	2833	REAL, INTENT (IN) :: delt
	2834	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, u, v
	2835	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t_wake, rr_wake
	2836	REAL, DIMENSION (nloc), INTENT (IN) :: s_wake
	2837	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra
	2838	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
	2839	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
	2840	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz, h, hp
	2841	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tp
	2842	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, cpn, ep, clw
	2843	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
	2844	REAL, DIMENSION (nloc, na), INTENT (IN) :: rp, up
	2845	REAL, DIMENSION (nloc, na), INTENT (IN) :: vp
	2846	REAL, DIMENSION (nloc, nd), INTENT (IN) :: wt
	2847	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: trap
	2848	REAL, DIMENSION (nloc, na), INTENT (IN) :: water, evap, b
	2849	REAL, DIMENSION (nloc, na), INTENT (IN) :: fondue, faci, ice
	2850	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: qent, uent
	2851	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: hent
	2852	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: vent, elij
	2853	INTEGER, DIMENSION (nloc, nd), INTENT (IN) :: nent
	2854	REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent
	2855	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti
	2856	REAL,INTENT(IN) :: tau_cld_cv, coefw_cld_cv
[2007]	2857	!
	2858	!input/output:
[2327]	2859	REAL, DIMENSION (nloc, na), INTENT (INOUT) :: m, mp
	2860	REAL, DIMENSION (nloc, na, na), INTENT (INOUT) :: ment
	2861	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
	2862	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig
	2863	REAL, DIMENSION (nloc), INTENT (INOUT) :: sigd
[2007]	2864	!
	2865	!outputs:
[2327]	2866	REAL, DIMENSION (nloc), INTENT (OUT) :: precip
	2867	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ft, fr, fu, fv
	2868	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ftd, fqd
	2869	REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT) :: ftra
	2870	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: upwd, dnwd, ma
	2871	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: dnwd0, mip
	2872	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecip
	2873	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecipi
	2874	!! REAL tls(nloc, nd), tps(nloc, nd) ! useless . jyg
	2875	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qcondc ! cld
	2876	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qtc, sigt ! cld
	2877	REAL, DIMENSION (nloc), INTENT (OUT) :: wd ! gust
	2878	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf
[2007]	2879	!
	2880	!local variables:
	2881	INTEGER i, k, il, n, j, num1
	2882	REAL rat, delti
	2883	REAL ax, bx, cx, dx, ex
	2884	REAL cpinv, rdcp, dpinv
	2885	REAL awat(nloc)
[2259]	2886	REAL lvcp(nloc, na), lfcp(nloc, na) ! , mke(nloc, na) ! unused . jyg
[2007]	2887	REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
	2888	!! real up1(nloc), dn1(nloc)
	2889	REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
	2890	REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
	2891	REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
	2892	REAL th_wake(nloc, nd)
	2893	REAL alpha_qpos(nloc), alpha_qpos1(nloc)
	2894	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
	2895	REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
[2205]	2896	REAL sument(nloc), sigment(nloc,nd), qtment(nloc,nd) ! cld
	2897	REAL qnk(nloc)
[2007]	2898	REAL sumdq !jyg
	2899	!
	2900	! -------------------------------------------------------------
[1992]	2901
[2007]	2902	! initialization:
[1992]	2903
	2904	delti = 1.0/delt
[2007]	2905	! print*,'cv3_yield initialisation delt', delt
	2906	!
[1992]	2907	DO il = 1, ncum
	2908	precip(il) = 0.0
[2007]	2909	Vprecip(il, nd+1) = 0.0
[2306]	2910	Vprecipi(il, nd+1) = 0.0 ! jyg: Vprecipi
[1992]	2911	wd(il) = 0.0 ! gust
	2912	END DO
	2913
	2914	DO i = 1, nd
	2915	DO il = 1, ncum
[2007]	2916	Vprecip(il, i) = 0.0
[2306]	2917	Vprecipi(il, i) = 0.0 ! jyg
[1992]	2918	ft(il, i) = 0.0
	2919	fr(il, i) = 0.0
	2920	fu(il, i) = 0.0
	2921	fv(il, i) = 0.0
	2922	upwd(il, i) = 0.0
	2923	dnwd(il, i) = 0.0
	2924	dnwd0(il, i) = 0.0
	2925	mip(il, i) = 0.0
	2926	ftd(il, i) = 0.0
	2927	fqd(il, i) = 0.0
	2928	qcondc(il, i) = 0.0 ! cld
	2929	qcond(il, i) = 0.0 ! cld
[2205]	2930	qtc(il, i) = 0.0 ! cld
	2931	qtment(il, i) = 0.0 ! cld
	2932	sigment(il, i) = 0.0 ! cld
	2933	sigt(il, i) = 0.0 ! cld
[1992]	2934	nqcond(il, i) = 0.0 ! cld
	2935	END DO
	2936	END DO
[2007]	2937	! print*,'cv3_yield initialisation 2'
	2938	!AC! do j=1,ntra
	2939	!AC! do i=1,nd
	2940	!AC! do il=1,ncum
	2941	!AC! ftra(il,i,j)=0.0
	2942	!AC! enddo
	2943	!AC! enddo
	2944	!AC! enddo
	2945	! print*,'cv3_yield initialisation 3'
[1992]	2946	DO i = 1, nl
	2947	DO il = 1, ncum
	2948	lvcp(il, i) = lv(il, i)/cpn(il, i)
	2949	lfcp(il, i) = lf(il, i)/cpn(il, i)
	2950	END DO
	2951	END DO
	2952
	2953
	2954
[2007]	2955	! * calculate surface precipitation in mm/day *
[1992]	2956
	2957	DO il = 1, ncum
	2958	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
	2959	IF (cvflag_ice) THEN
[2007]	2960	precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
	2961	86400.1000./(rowl*grav)
[1992]	2962	ELSE
[2007]	2963	precip(il) = wt(il, 1)sigd(il)water(il, 1) &
	2964	86400.1000./(rowl*grav)
[1992]	2965	END IF
	2966	END IF
	2967	END DO
[2007]	2968	! print*,'cv3_yield apres calcul precip'
[1992]	2969
	2970
[2007]	2971	! === calculate vertical profile of precipitation in kg/m2/s ===
[1992]	2972
	2973	DO i = 1, nl
	2974	DO il = 1, ncum
	2975	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
	2976	IF (cvflag_ice) THEN
[2007]	2977	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
[2306]	2978	Vprecipi(il, i) = wt(il, i)sigd(il)ice(il,i)/grav ! jyg
[1992]	2979	ELSE
[2007]	2980	Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
[2306]	2981	Vprecipi(il, i) = 0. ! jyg
[1992]	2982	END IF
	2983	END IF
	2984	END DO
	2985	END DO
	2986
	2987
[2007]	2988	! * Calculate downdraft velocity scale *
	2989	! * NE PAS UTILISER POUR L'INSTANT *
[1992]	2990
[2007]	2991	!! do il=1,ncum
	2992	!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
	2993	!! /(sigd(il)*p(il,icb(il)))
	2994	!! enddo
[1992]	2995
	2996
[2007]	2997	! * calculate tendencies of lowest level potential temperature *
	2998	! * and mixing ratio *
[1992]	2999
	3000	DO il = 1, ncum
	3001	work(il) = 1.0/(ph(il,1)-ph(il,2))
	3002	cbmf(il) = 0.0
	3003	END DO
	3004
	3005	DO k = 2, nl
	3006	DO il = 1, ncum
	3007	IF (k>=icb(il)) THEN
	3008	cbmf(il) = cbmf(il) + m(il, k)
	3009	END IF
	3010	END DO
	3011	END DO
	3012
[2007]	3013	! print*,'cv3_yield avant ft'
	3014	! am is the part of cbmf taken from the first level
[1992]	3015	DO il = 1, ncum
	3016	am(il) = cbmf(il)*wghti(il, 1)
	3017	END DO
	3018
	3019	DO il = 1, ncum
	3020	IF (iflag(il)<=1) THEN
[2007]	3021	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
	3022	!JYG Correction pour conserver l'eau
	3023	! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
[1992]	3024	IF (cvflag_ice) THEN
	3025	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
[2007]	3026	lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
	3027	lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
	3028	(100.*(ph(il,1)-ph(il,2))) !precip
[1992]	3029	ELSE
	3030	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
	3031	END IF
	3032
[2007]	3033	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
[1992]	3034
	3035	IF (cvflag_ice) THEN
[2007]	3036	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
	3037	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
	3038	0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
	3039	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
[1992]	3040	ELSE
[2007]	3041	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
	3042	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
[1992]	3043	END IF
	3044
[2007]	3045	ftd(il, 1) = ft(il, 1) ! fin precip
[1992]	3046
[2007]	3047	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
	3048	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
	3049	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
[1992]	3050	END IF ! iflag
	3051	END DO
	3052
	3053
	3054	DO j = 2, nl
	3055	IF (iflag_mix>0) THEN
	3056	DO il = 1, ncum
[2007]	3057	! FH WARNING a modifier :
[1992]	3058	cpinv = 0.
[2007]	3059	! cpinv=1.0/cpn(il,1)
[1992]	3060	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
[2007]	3061	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
	3062	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
[1992]	3063	END IF ! j
	3064	END DO
	3065	END IF
	3066	END DO
[2007]	3067	! fin sature
[1992]	3068
	3069
	3070	DO il = 1, ncum
	3071	IF (iflag(il)<=1) THEN
[2007]	3072	!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
	3073	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
	3074	sigd(il)*evap(il, 1)
	3075	!!! sigd(il)0.5(evap(il,1)+evap(il,2))
[1992]	3076
[2007]	3077	fqd(il, 1) = fr(il, 1) !precip
[1992]	3078
[2007]	3079	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
[1992]	3080
[2007]	3081	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
	3082	am(il)*(u(il,2)-u(il,1)))
	3083	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
	3084	am(il)*(v(il,2)-v(il,1)))
[1992]	3085	END IF ! iflag
	3086	END DO ! il
	3087
	3088
[2007]	3089	!AC! do j=1,ntra
	3090	!AC! do il=1,ncum
	3091	!AC! if (iflag(il) .le. 1) then
	3092	!AC! if (cvflag_grav) then
	3093	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
	3094	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	3095	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	3096	!AC! else
	3097	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
	3098	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	3099	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	3100	!AC! endif
	3101	!AC! endif ! iflag
	3102	!AC! enddo
	3103	!AC! enddo
[1992]	3104
	3105	DO j = 2, nl
	3106	DO il = 1, ncum
	3107	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
[2007]	3108	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
	3109	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
	3110	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
[1992]	3111	END IF ! j
	3112	END DO
	3113	END DO
	3114
[2007]	3115	!AC! do k=1,ntra
	3116	!AC! do j=2,nl
	3117	!AC! do il=1,ncum
	3118	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
	3119	!AC!
	3120	!AC! if (cvflag_grav) then
	3121	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
	3122	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
	3123	!AC! else
	3124	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
	3125	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
	3126	!AC! endif
	3127	!AC!
	3128	!AC! endif
	3129	!AC! enddo
	3130	!AC! enddo
	3131	!AC! enddo
	3132	! print*,'cv3_yield apres ft'
[1992]	3133
[2007]	3134	! * calculate tendencies of potential temperature and mixing ratio *
	3135	! * at levels above the lowest level *
[1992]	3136
[2007]	3137	! * first find the net saturated updraft and downdraft mass fluxes *
	3138	! * through each level *
[1992]	3139
	3140
	3141	DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
	3142
	3143	num1 = 0
	3144	DO il = 1, ncum
	3145	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
	3146	END DO
	3147	IF (num1<=0) GO TO 500
	3148
	3149	CALL zilch(amp1, ncum)
	3150	CALL zilch(ad, ncum)
	3151
	3152	DO k = 1, nl + 1
	3153	DO il = 1, ncum
	3154	IF (i>=icb(il)) THEN
	3155	IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
	3156	amp1(il) = amp1(il) + m(il, k)
	3157	END IF
	3158	ELSE
[2007]	3159	! AMP1 is the part of cbmf taken from layers I and lower
[1992]	3160	IF (k<=i) THEN
	3161	amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
	3162	END IF
	3163	END IF
	3164	END DO
	3165	END DO
	3166
	3167	DO k = 1, i
	3168	DO j = i + 1, nl + 1
	3169	DO il = 1, ncum
	3170	IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
	3171	amp1(il) = amp1(il) + ment(il, k, j)
	3172	END IF
	3173	END DO
	3174	END DO
	3175	END DO
	3176
	3177	DO k = 1, i - 1
	3178	DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
	3179	DO il = 1, ncum
	3180	IF (i<=inb(il) .AND. j<=inb(il)) THEN
	3181	ad(il) = ad(il) + ment(il, j, k)
	3182	END IF
	3183	END DO
	3184	END DO
	3185	END DO
	3186
	3187	DO il = 1, ncum
	3188	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
	3189	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	3190	cpinv = 1.0/cpn(il, i)
	3191
[2007]	3192	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
	3193	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
[1992]	3194
[2007]	3195	! precip
	3196	! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
[1992]	3197	IF (cvflag_ice) THEN
	3198	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
[2007]	3199	sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
	3200	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
[1992]	3201	ELSE
	3202	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
	3203	END IF
	3204
	3205	rat = cpn(il, i-1)*cpinv
	3206
[2007]	3207	ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
	3208	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
	3209	IF (cvflag_ice) THEN
	3210	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
	3211	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
	3212	0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
	3213	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
	3214	ELSE
	3215	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
	3216	(t_wake(il,i+1)-t_wake(il,i))dpinv &
	3217	cpinv
	3218	END IF
[1992]	3219
[2007]	3220	ftd(il, i) = ft(il, i)
	3221	! fin precip
[1992]	3222
[2007]	3223	! sature
	3224	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
	3225	(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
	3226	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
[1992]	3227
	3228
[2007]	3229	IF (iflag_mix==0) THEN
	3230	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
	3231	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
	3232	END IF
[1992]	3233
	3234
	3235
[2007]	3236	! sb: on ne fait pas encore la correction permettant de mieux
	3237	! conserver l'eau:
	3238	!JYG: correction permettant de mieux conserver l'eau:
	3239	! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
	3240	fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
	3241	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
	3242	fqd(il, i) = fr(il, i) ! precip
[1992]	3243
[2007]	3244	fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
	3245	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	3246	fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
	3247	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
[1992]	3248
	3249
[2007]	3250	fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
	3251	ad(il)*(rr(il,i)-rr(il,i-1)))
	3252	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
	3253	ad(il)*(u(il,i)-u(il,i-1)))
	3254	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
	3255	ad(il)*(v(il,i)-v(il,i-1)))
[1992]	3256
	3257	END IF ! i
	3258	END DO
	3259
[2007]	3260	!AC! do k=1,ntra
	3261	!AC! do il=1,ncum
	3262	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
	3263	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	3264	!AC! cpinv=1.0/cpn(il,i)
	3265	!AC! if (cvflag_grav) then
	3266	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
	3267	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	3268	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	3269	!AC! else
	3270	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
	3271	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	3272	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	3273	!AC! endif
	3274	!AC! endif
	3275	!AC! enddo
	3276	!AC! enddo
[1992]	3277
	3278	DO k = 1, i - 1
	3279
	3280	DO il = 1, ncum
	3281	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
	3282	awat(il) = max(awat(il), 0.0)
	3283	END DO
	3284
	3285	IF (iflag_mix/=0) THEN
	3286	DO il = 1, ncum
	3287	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
	3288	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	3289	cpinv = 1.0/cpn(il, i)
[2007]	3290	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
	3291	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
	3292	!
	3293	!
[1992]	3294	END IF ! i
	3295	END DO
	3296	END IF
	3297
	3298	DO il = 1, ncum
	3299	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
	3300	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	3301	cpinv = 1.0/cpn(il, i)
[2007]	3302	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
	3303	(qent(il,k,i)-awat(il)-rr(il,i))
	3304	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
	3305	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
[1992]	3306
[2007]	3307	! (saturated updrafts resulting from mixing) ! cld
	3308	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
[2205]	3309	qtment(il, i) = qtment(il, i) + qent(il,k,i) ! cld
	3310	nqcond(il, i) = nqcond(il, i) + 1. ! cld
[1992]	3311	END IF ! i
	3312	END DO
	3313	END DO
	3314
[2007]	3315	!AC! do j=1,ntra
	3316	!AC! do k=1,i-1
	3317	!AC! do il=1,ncum
	3318	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
	3319	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	3320	!AC! cpinv=1.0/cpn(il,i)
	3321	!AC! if (cvflag_grav) then
	3322	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	3323	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
	3324	!AC! else
	3325	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	3326	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
	3327	!AC! endif
	3328	!AC! endif
	3329	!AC! enddo
	3330	!AC! enddo
	3331	!AC! enddo
[1992]	3332
	3333	DO k = i, nl + 1
	3334
	3335	IF (iflag_mix/=0) THEN
	3336	DO il = 1, ncum
	3337	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
	3338	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	3339	cpinv = 1.0/cpn(il, i)
[2007]	3340	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
	3341	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
[1992]	3342
	3343
	3344	END IF ! i
	3345	END DO
	3346	END IF
	3347
	3348	DO il = 1, ncum
	3349	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
	3350	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
	3351	cpinv = 1.0/cpn(il, i)
	3352
[2007]	3353	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
	3354	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
	3355	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
[1992]	3356	END IF ! i and k
	3357	END DO
	3358	END DO
	3359
[2007]	3360	!AC! do j=1,ntra
	3361	!AC! do k=i,nl+1
	3362	!AC! do il=1,ncum
	3363	!AC! if (i.le.inb(il) .and. k.le.inb(il)
	3364	!AC! $ .and. iflag(il) .le. 1) then
	3365	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	3366	!AC! cpinv=1.0/cpn(il,i)
	3367	!AC! if (cvflag_grav) then
	3368	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	3369	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
	3370	!AC! else
	3371	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	3372	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
	3373	!AC! endif
	3374	!AC! endif ! i and k
	3375	!AC! enddo
	3376	!AC! enddo
	3377	!AC! enddo
[1992]	3378
[2007]	3379	! sb: interface with the cloud parameterization: ! cld
[1992]	3380
	3381	DO k = i + 1, nl
	3382	DO il = 1, ncum
[2007]	3383	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
	3384	! (saturated downdrafts resulting from mixing) ! cld
	3385	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
[2205]	3386	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
	3387	nqcond(il, i) = nqcond(il, i) + 1. ! cld
[1992]	3388	END IF ! cld
	3389	END DO ! cld
	3390	END DO ! cld
	3391
[2007]	3392	! (particular case: no detraining level is found) ! cld
	3393	DO il = 1, ncum ! cld
	3394	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
	3395	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
[2205]	3396	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
[2007]	3397	nqcond(il, i) = nqcond(il, i) + 1. ! cld
	3398	END IF ! cld
	3399	END DO ! cld
[1992]	3400
[2007]	3401	DO il = 1, ncum ! cld
	3402	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
	3403	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
[2205]	3404	qtment(il, i) = qtment(il,i)/nqcond(il, i) ! cld
[2007]	3405	END IF ! cld
[1992]	3406	END DO
	3407
[2007]	3408	!AC! do j=1,ntra
	3409	!AC! do il=1,ncum
	3410	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
	3411	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
	3412	!AC! cpinv=1.0/cpn(il,i)
	3413	!AC!
	3414	!AC! if (cvflag_grav) then
	3415	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
	3416	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	3417	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
	3418	!AC! else
	3419	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
	3420	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	3421	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
	3422	!AC! endif
	3423	!AC! endif ! i
	3424	!AC! enddo
	3425	!AC! enddo
[1992]	3426
	3427
	3428	500 END DO
	3429
[2007]	3430	!JYG<
	3431	!Conservation de l'eau
	3432	! sumdq = 0.
	3433	! DO k = 1, nl
	3434	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
	3435	! END DO
	3436	! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
	3437	!JYG>
	3438	! * move the detrainment at level inb down to level inb-1 *
	3439	! * in such a way as to preserve the vertically *
	3440	! * integrated enthalpy and water tendencies *
[1992]	3441
[2007]	3442	! Correction bug le 18-03-09
[1992]	3443	DO il = 1, ncum
	3444	IF (iflag(il)<=1) THEN
[2007]	3445	ax = 0.01gravment(il, inb(il), inb(il))* &
	3446	(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
	3447	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
	3448	ft(il, inb(il)) = ft(il, inb(il)) - ax
	3449	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
	3450	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
[1992]	3451
[2007]	3452	bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
	3453	(ph(il,inb(il))-ph(il,inb(il)+1))
	3454	fr(il, inb(il)) = fr(il, inb(il)) - bx
	3455	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
	3456	(ph(il,inb(il)-1)-ph(il,inb(il)))
[1992]	3457
[2007]	3458	cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
	3459	(ph(il,inb(il))-ph(il,inb(il)+1))
	3460	fu(il, inb(il)) = fu(il, inb(il)) - cx
	3461	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
	3462	(ph(il,inb(il)-1)-ph(il,inb(il)))
[1992]	3463
[2007]	3464	dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
	3465	(ph(il,inb(il))-ph(il,inb(il)+1))
	3466	fv(il, inb(il)) = fv(il, inb(il)) - dx
	3467	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
	3468	(ph(il,inb(il)-1)-ph(il,inb(il)))
[1992]	3469	END IF !iflag
	3470	END DO
	3471
[2007]	3472	!JYG<
	3473	!Conservation de l'eau
	3474	! sumdq = 0.
	3475	! DO k = 1, nl
	3476	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
	3477	! END DO
	3478	! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
	3479	!JYG>
[1992]	3480
[2007]	3481	!AC! do j=1,ntra
	3482	!AC! do il=1,ncum
	3483	!AC! IF (iflag(il) .le. 1) THEN
	3484	!AC! IF (cvflag_grav) then
	3485	!AC! ex=0.01gravment(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! else
	3493	!AC! ex=0.1*ment(il,inb(il),inb(il))
	3494	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
	3495	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
	3496	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
	3497	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
	3498	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
	3499	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
	3500	!AC! ENDIF !cvflag grav
	3501	!AC! ENDIF !iflag
	3502	!AC! enddo
	3503	!AC! enddo
[1992]	3504
	3505
[2007]	3506	! * homogenize tendencies below cloud base *
[1992]	3507
[2007]	3508
[1992]	3509	DO il = 1, ncum
	3510	asum(il) = 0.0
	3511	bsum(il) = 0.0
	3512	csum(il) = 0.0
	3513	dsum(il) = 0.0
	3514	esum(il) = 0.0
	3515	fsum(il) = 0.0
	3516	gsum(il) = 0.0
	3517	hsum(il) = 0.0
	3518	END DO
	3519
[2007]	3520	!do i=1,nl
	3521	!do il=1,ncum
	3522	!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
	3523	!enddo
	3524	!enddo
[1992]	3525
	3526	DO i = 1, nl
	3527	DO il = 1, ncum
	3528	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
[2007]	3529	!jyg Saturated part : use T profile
[1992]	3530	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
[2007]	3531	!jyg<20140311
	3532	!Correction pour conserver l eau
	3533	IF (ok_conserv_q) THEN
	3534	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
	3535	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
	3536
	3537	ELSE
	3538	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
	3539	(ph(il,i)-ph(il,i+1))
	3540	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
	3541	(ph(il,i)-ph(il,i+1))
	3542	ENDIF ! (ok_conserv_q)
	3543	!jyg>
[1992]	3544	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
[2007]	3545	!jyg Unsaturated part : use T_wake profile
[1992]	3546	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
[2007]	3547	!jyg<20140311
	3548	!Correction pour conserver l eau
	3549	IF (ok_conserv_q) THEN
	3550	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
	3551	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
	3552	ELSE
	3553	fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
	3554	(ph(il,i)-ph(il,i+1))
	3555	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
	3556	(ph(il,i)-ph(il,i+1))
	3557	ENDIF ! (ok_conserv_q)
	3558	!jyg>
	3559	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
[1992]	3560	END IF
	3561	END DO
	3562	END DO
	3563
[2007]	3564	!!!! do 700 i=1,icb(il)-1
[1992]	3565	DO i = 1, nl
	3566	DO il = 1, ncum
	3567	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
	3568	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
	3569	fqd(il, i) = fsum(il)/gsum(il)
	3570	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
	3571	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
	3572	END IF
	3573	END DO
	3574	END DO
	3575
[2007]	3576	!jyg<
	3577	!Conservation de l'eau
	3578	!! sumdq = 0.
	3579	!! DO k = 1, nl
	3580	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
	3581	!! END DO
	3582	!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
	3583	!jyg>
[1992]	3584
[2007]	3585
	3586	! *** Check that moisture stays positive. If not, scale tendencies
	3587	! in order to ensure moisture positivity
[1992]	3588	DO il = 1, ncum
	3589	alpha_qpos(il) = 1.
	3590	IF (iflag(il)<=1) THEN
	3591	IF (fr(il,1)<=0.) THEN
[2007]	3592	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]	3593	END IF
	3594	END IF
	3595	END DO
	3596	DO i = 2, nl
	3597	DO il = 1, ncum
	3598	IF (iflag(il)<=1) THEN
	3599	IF (fr(il,i)<=0.) THEN
[2007]	3600	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
	3601	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
[1992]	3602	END IF
	3603	END IF
	3604	END DO
	3605	END DO
	3606	DO il = 1, ncum
	3607	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
	3608	alpha_qpos(il) = alpha_qpos(il)*1.1
	3609	END IF
	3610	END DO
[2327]	3611	!
[2364]	3612	! print *,' YIELD : alpha_qpos ',alpha_qpos(1)
[2327]	3613	!
[1992]	3614	DO il = 1, ncum
	3615	IF (iflag(il)<=1) THEN
	3616	sigd(il) = sigd(il)/alpha_qpos(il)
	3617	precip(il) = precip(il)/alpha_qpos(il)
	3618	END IF
	3619	END DO
	3620	DO i = 1, nl
	3621	DO il = 1, ncum
	3622	IF (iflag(il)<=1) THEN
	3623	fr(il, i) = fr(il, i)/alpha_qpos(il)
	3624	ft(il, i) = ft(il, i)/alpha_qpos(il)
	3625	fqd(il, i) = fqd(il, i)/alpha_qpos(il)
	3626	ftd(il, i) = ftd(il, i)/alpha_qpos(il)
	3627	fu(il, i) = fu(il, i)/alpha_qpos(il)
	3628	fv(il, i) = fv(il, i)/alpha_qpos(il)
	3629	m(il, i) = m(il, i)/alpha_qpos(il)
	3630	mp(il, i) = mp(il, i)/alpha_qpos(il)
[2306]	3631	Vprecip(il, i) = Vprecip(il, i)/alpha_qpos(il)
	3632	Vprecipi(il, i) = Vprecipi(il, i)/alpha_qpos(il) ! jyg
[1992]	3633	END IF
	3634	END DO
	3635	END DO
	3636	DO i = 1, nl
	3637	DO j = 1, nl
	3638	DO il = 1, ncum
	3639	IF (iflag(il)<=1) THEN
	3640	ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
	3641	END IF
	3642	END DO
	3643	END DO
	3644	END DO
	3645
[2007]	3646	!AC! DO j = 1,ntra
	3647	!AC! DO i = 1,nl
	3648	!AC! DO il = 1,ncum
	3649	!AC! IF (iflag(il) .le. 1) THEN
	3650	!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
	3651	!AC! ENDIF
	3652	!AC! ENDDO
	3653	!AC! ENDDO
	3654	!AC! ENDDO
[1992]	3655
	3656
[2007]	3657	! * reset counter and return *
[1992]	3658
[2253]	3659	! Reset counter only for points actually convective (jyg)
	3660	! In order take into account the possibility of changing the compression,
	3661	! reset m, sig and w0 to zero for non-convecting points.
[1992]	3662	DO il = 1, ncum
[2253]	3663	IF (iflag(il) < 3) THEN
	3664	sig(il, nd) = 2.0
	3665	ENDIF
[1992]	3666	END DO
	3667
	3668
	3669	DO i = 1, nd
	3670	DO il = 1, ncum
	3671	upwd(il, i) = 0.0
	3672	dnwd(il, i) = 0.0
	3673	END DO
	3674	END DO
	3675
	3676	DO i = 1, nl
	3677	DO il = 1, ncum
	3678	dnwd0(il, i) = -mp(il, i)
	3679	END DO
	3680	END DO
	3681	DO i = nl + 1, nd
	3682	DO il = 1, ncum
	3683	dnwd0(il, i) = 0.
	3684	END DO
	3685	END DO
	3686
	3687
	3688	DO i = 1, nl
	3689	DO il = 1, ncum
	3690	IF (i>=icb(il) .AND. i<=inb(il)) THEN
	3691	upwd(il, i) = 0.0
	3692	dnwd(il, i) = 0.0
	3693	END IF
	3694	END DO
	3695	END DO
	3696
	3697	DO i = 1, nl
	3698	DO k = 1, nl
	3699	DO il = 1, ncum
	3700	up1(il, k, i) = 0.0
	3701	dn1(il, k, i) = 0.0
	3702	END DO
	3703	END DO
	3704	END DO
	3705
	3706	DO i = 1, nl
	3707	DO k = i, nl
	3708	DO n = 1, i - 1
	3709	DO il = 1, ncum
	3710	IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
	3711	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
	3712	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
	3713	END IF
	3714	END DO
	3715	END DO
	3716	END DO
	3717	END DO
	3718
	3719	DO i = 1, nl
	3720	DO k = 1, nl
	3721	DO il = 1, ncum
	3722	IF (i>=icb(il)) THEN
	3723	IF (k>=i .AND. k<=(inb(il))) THEN
	3724	upwd(il, i) = upwd(il, i) + m(il, k)
	3725	END IF
	3726	ELSE
	3727	IF (k<i) THEN
	3728	upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
	3729	END IF
	3730	END IF
[2007]	3731	! c print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
[1992]	3732	END DO
	3733	END DO
	3734	END DO
	3735
	3736	DO i = 2, nl
	3737	DO k = i, nl
	3738	DO il = 1, ncum
[2007]	3739	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
[1992]	3740	IF (i<=inb(il) .AND. k<=inb(il)) THEN
	3741	upwd(il, i) = upwd(il, i) + up1(il, k, i)
	3742	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
	3743	END IF
[2007]	3744	! c print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
[1992]	3745	END DO
	3746	END DO
	3747	END DO
	3748
	3749
[2007]	3750	!!!! DO il=1,ncum
	3751	!!!! do i=icb(il),inb(il)
	3752	!!!!
	3753	!!!! upwd(il,i)=0.0
	3754	!!!! dnwd(il,i)=0.0
	3755	!!!! do k=i,inb(il)
	3756	!!!! up1=0.0
	3757	!!!! dn1=0.0
	3758	!!!! do n=1,i-1
	3759	!!!! up1=up1+ment(il,n,k)
	3760	!!!! dn1=dn1-ment(il,k,n)
	3761	!!!! enddo
	3762	!!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
	3763	!!!! dnwd(il,i)=dnwd(il,i)+dn1
	3764	!!!! enddo
	3765	!!!! enddo
	3766	!!!!
	3767	!!!! ENDDO
[1992]	3768
[2007]	3769	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	3770	! determination de la variation de flux ascendant entre
	3771	! deux niveau non dilue mip
	3772	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1992]	3773
	3774	DO i = 1, nl
	3775	DO il = 1, ncum
	3776	mip(il, i) = m(il, i)
	3777	END DO
	3778	END DO
	3779
	3780	DO i = nl + 1, nd
	3781	DO il = 1, ncum
	3782	mip(il, i) = 0.
	3783	END DO
	3784	END DO
	3785
	3786	DO i = 1, nd
	3787	DO il = 1, ncum
	3788	ma(il, i) = 0
	3789	END DO
	3790	END DO
	3791
	3792	DO i = 1, nl
	3793	DO j = i, nl
	3794	DO il = 1, ncum
	3795	ma(il, i) = ma(il, i) + m(il, j)
	3796	END DO
	3797	END DO
	3798	END DO
	3799
	3800	DO i = nl + 1, nd
	3801	DO il = 1, ncum
	3802	ma(il, i) = 0.
	3803	END DO
	3804	END DO
	3805
	3806	DO i = 1, nl
	3807	DO il = 1, ncum
	3808	IF (i<=(icb(il)-1)) THEN
	3809	ma(il, i) = 0
	3810	END IF
	3811	END DO
	3812	END DO
	3813
[2007]	3814	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	3815	! icb represente de niveau ou se trouve la
	3816	! base du nuage , et inb le top du nuage
	3817	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1992]	3818
[2259]	3819	!! DO i = 1, nd ! unused . jyg
	3820	!! DO il = 1, ncum ! unused . jyg
	3821	!! mke(il, i) = upwd(il, i) + dnwd(il, i) ! unused . jyg
	3822	!! END DO ! unused . jyg
	3823	!! END DO ! unused . jyg
[1992]	3824
[2259]	3825	!! DO i = 1, nd ! unused . jyg
	3826	!! DO il = 1, ncum ! unused . jyg
	3827	!! rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd(1.-rr(il,i))+rr(il,i)cpv) ! unused . jyg
	3828	!! tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp ! unused . jyg
	3829	!! tps(il, i) = tp(il, i) ! unused . jyg
	3830	!! END DO ! unused . jyg
	3831	!! END DO ! unused . jyg
[1992]	3832
	3833
[2007]	3834	! * diagnose the in-cloud mixing ratio * ! cld
	3835	! * of condensed water * ! cld
	3836	!! cld
	3837
	3838	DO i = 1, nd ! cld
	3839	DO il = 1, ncum ! cld
	3840	mac(il, i) = 0.0 ! cld
	3841	wa(il, i) = 0.0 ! cld
	3842	siga(il, i) = 0.0 ! cld
	3843	sax(il, i) = 0.0 ! cld
	3844	END DO ! cld
	3845	END DO ! cld
	3846
	3847	DO i = minorig, nl ! cld
	3848	DO k = i + 1, nl + 1 ! cld
	3849	DO il = 1, ncum ! cld
[1992]	3850	IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
[2007]	3851	mac(il, i) = mac(il, i) + m(il, k) ! cld
	3852	END IF ! cld
	3853	END DO ! cld
	3854	END DO ! cld
	3855	END DO ! cld
[1992]	3856
[2007]	3857	DO i = 1, nl ! cld
	3858	DO j = 1, i ! cld
	3859	DO il = 1, ncum ! cld
	3860	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
	3861	.AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
	3862	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
	3863	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
	3864	END IF ! cld
	3865	END DO ! cld
	3866	END DO ! cld
	3867	END DO ! cld
[1992]	3868
[2007]	3869	DO i = 1, nl ! cld
	3870	DO il = 1, ncum ! cld
	3871	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
	3872	.AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
	3873	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
	3874	END IF ! cld
	3875	END DO ! cld
[2205]	3876	END DO
	3877	! cld
	3878	DO i = 1, nl
[1992]	3879
[2205]	3880	! 14/01/15 AJ je remets les parties manquantes cf JYG
	3881	! Initialize sument to 0
	3882
	3883	DO il = 1,ncum
	3884	sument(il) = 0.
	3885	ENDDO
	3886
	3887	! Sum mixed mass fluxes in sument
	3888
	3889	DO k = 1,nl
	3890	DO il = 1,ncum
	3891	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
	3892	sument(il) =sument(il) + abs(ment(il,k,i))
	3893	ENDIF
	3894	ENDDO ! il
	3895	ENDDO ! k
	3896
	3897	! 14/01/15 AJ delta n'a rien à faire là...
[2007]	3898	DO il = 1, ncum ! cld
	3899	IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
[2205]	3900	siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
	3901	rrdtvp(il, i)/p(il, i)/100. ! cld
	3902
[2007]	3903	siga(il, i) = min(siga(il,i), 1.0) ! cld
[2205]	3904
	3905	! IM cf. FH
	3906	! 14/01/15 AJ ne correspond pas à ce qui a été codé par JYG et SB
	3907
[2007]	3908	IF (iflag_clw==0) THEN ! cld
	3909	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
	3910	+(1.-siga(il,i))*qcond(il, i) ! cld
[2205]	3911
	3912
	3913	sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1)) ! cld
	3914	sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i)) ! cld
	3915	qtc(il, i) = (siga(il,i)qnk(il)+sigment(il,i)qtment(il,i)) & ! cld
	3916	/(siga(il,i)+sigment(il,i)) ! cld
	3917	sigt(il,i) = sigment(il, i) + siga(il, i)
	3918
	3919	! qtc(il, i) = siga(il,i)qnk(il)+(1.-siga(il,i))qtment(il,i) ! cld
[2208]	3920	! print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)
[2205]	3921
[2007]	3922	ELSE IF (iflag_clw==1) THEN ! cld
	3923	qcondc(il, i) = qcond(il, i) ! cld
[2205]	3924	qtc(il,i) = qtment(il,i) ! cld
[2007]	3925	END IF ! cld
[1992]	3926
[2007]	3927	END DO ! cld
[1992]	3928	END DO
[2007]	3929	! print*,'cv3_yield fin'
	3930
[1992]	3931	RETURN
	3932	END SUBROUTINE cv3_yield
	3933
[2007]	3934	!AC! et !RomP >>>
	3935	SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
	3936	ment, sigij, da, phi, phi2, d1a, dam, &
	3937	ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
	3938	icb, inb)
[1992]	3939	IMPLICIT NONE
	3940
	3941	include "cv3param.h"
	3942
[2007]	3943	!inputs:
[1992]	3944	INTEGER ncum, nd, na, nloc, len
	3945	REAL ment(nloc, na, na), sigij(nloc, na, na)
	3946	REAL clw(nloc, nd), elij(nloc, na, na)
	3947	REAL ep(nloc, na)
	3948	INTEGER icb(nloc), inb(nloc)
[2007]	3949	REAL Vprecip(nloc, nd+1)
	3950	!ouputs:
[1992]	3951	REAL da(nloc, na), phi(nloc, na, na)
	3952	REAL phi2(nloc, na, na)
	3953	REAL d1a(nloc, na), dam(nloc, na)
[2007]	3954	REAL epmlmMm(nloc, na, na), eplaMm(nloc, na)
	3955	! variables pour tracer dans precip de l'AA et des mel
	3956	!local variables:
[1992]	3957	INTEGER i, j, k
	3958	REAL epm(nloc, na, na)
	3959
[2007]	3960	! variables d'Emanuel : du second indice au troisieme
	3961	! ---> tab(i,k,j) -> de l origine k a l arrivee j
	3962	! ment, sigij, elij
	3963	! variables personnelles : du troisieme au second indice
	3964	! ---> tab(i,j,k) -> de k a j
	3965	! phi, phi2
[1992]	3966
[2007]	3967	! initialisations
[1992]	3968
	3969	da(:, :) = 0.
	3970	d1a(:, :) = 0.
	3971	dam(:, :) = 0.
	3972	epm(:, :, :) = 0.
[2007]	3973	eplaMm(:, :) = 0.
	3974	epmlmMm(:, :, :) = 0.
[1992]	3975	phi(:, :, :) = 0.
	3976	phi2(:, :, :) = 0.
	3977
[2007]	3978	! fraction deau condensee dans les melanges convertie en precip : epm
	3979	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
[1992]	3980	DO j = 1, na
	3981	DO k = 1, na
	3982	DO i = 1, ncum
[2007]	3983	IF (k>=icb(i) .AND. k<=inb(i) .AND. &
	3984	!!jyg j.ge.k.and.j.le.inb(i)) then
	3985	!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
[1992]	3986	j>k .AND. j<=inb(i)) THEN
	3987	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
[2007]	3988	!!
[1992]	3989	epm(i, j, k) = max(epm(i,j,k), 0.0)
	3990	END IF
	3991	END DO
	3992	END DO
	3993	END DO
	3994
	3995
	3996	DO j = 1, na
	3997	DO k = 1, na
	3998	DO i = 1, ncum
	3999	IF (k>=icb(i) .AND. k<=inb(i)) THEN
[2007]	4000	eplaMm(i, j) = eplamm(i, j) + &
	4001	ep(i, j)clw(i, j)ment(i, j, k)*(1.-sigij(i,j,k))
[1992]	4002	END IF
	4003	END DO
	4004	END DO
	4005	END DO
	4006
	4007	DO j = 1, na
	4008	DO k = 1, j - 1
	4009	DO i = 1, ncum
	4010	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
[2007]	4011	epmlmMm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
[1992]	4012	END IF
	4013	END DO
	4014	END DO
	4015	END DO
	4016
[2007]	4017	! matrices pour calculer la tendance des concentrations dans cvltr.F90
[1992]	4018	DO j = 1, na
	4019	DO k = 1, na
	4020	DO i = 1, ncum
	4021	da(i, j) = da(i, j) + (1.-sigij(i,k,j))*ment(i, k, j)
	4022	phi(i, j, k) = sigij(i, k, j)*ment(i, k, j)
	4023	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sigij(i,k,j))
	4024	IF (k<=j) THEN
[2007]	4025	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, k, j)(1.-ep(i,k))*(1.-sigij(i,k,j))
[1992]	4026	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
	4027	END IF
	4028	END DO
	4029	END DO
	4030	END DO
	4031
	4032	RETURN
	4033	END SUBROUTINE cv3_tracer
[2007]	4034	!AC! et !RomP <<<
[1992]	4035
[2007]	4036	SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
	4037	iflag, &
	4038	precip, sig, w0, &
	4039	ft, fq, fu, fv, ftra, &
	4040	Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
	4041	iflag1, &
	4042	precip1, sig1, w01, &
	4043	ft1, fq1, fu1, fv1, ftra1, &
	4044	Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1)
[1992]	4045	IMPLICIT NONE
	4046
	4047	include "cv3param.h"
	4048
[2007]	4049	!inputs:
[1992]	4050	INTEGER len, ncum, nd, ntra, nloc
	4051	INTEGER idcum(nloc)
	4052	INTEGER iflag(nloc)
	4053	REAL precip(nloc)
	4054	REAL sig(nloc, nd), w0(nloc, nd)
	4055	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	4056	REAL ftra(nloc, nd, ntra)
	4057	REAL ma(nloc, nd)
	4058	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
	4059	REAL qcondc(nloc, nd)
	4060	REAL wd(nloc), cape(nloc)
	4061
[2007]	4062	!outputs:
[1992]	4063	INTEGER iflag1(len)
	4064	REAL precip1(len)
	4065	REAL sig1(len, nd), w01(len, nd)
	4066	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
	4067	REAL ftra1(len, nd, ntra)
	4068	REAL ma1(len, nd)
	4069	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
	4070	REAL qcondc1(nloc, nd)
	4071	REAL wd1(nloc), cape1(nloc)
	4072
[2007]	4073	!local variables:
[1992]	4074	INTEGER i, k, j
	4075
	4076	DO i = 1, ncum
	4077	precip1(idcum(i)) = precip(i)
	4078	iflag1(idcum(i)) = iflag(i)
	4079	wd1(idcum(i)) = wd(i)
	4080	cape1(idcum(i)) = cape(i)
	4081	END DO
	4082
	4083	DO k = 1, nl
	4084	DO i = 1, ncum
	4085	sig1(idcum(i), k) = sig(i, k)
	4086	w01(idcum(i), k) = w0(i, k)
	4087	ft1(idcum(i), k) = ft(i, k)
	4088	fq1(idcum(i), k) = fq(i, k)
	4089	fu1(idcum(i), k) = fu(i, k)
	4090	fv1(idcum(i), k) = fv(i, k)
	4091	ma1(idcum(i), k) = ma(i, k)
	4092	upwd1(idcum(i), k) = upwd(i, k)
	4093	dnwd1(idcum(i), k) = dnwd(i, k)
	4094	dnwd01(idcum(i), k) = dnwd0(i, k)
	4095	qcondc1(idcum(i), k) = qcondc(i, k)
	4096	END DO
	4097	END DO
	4098
	4099	DO i = 1, ncum
	4100	sig1(idcum(i), nd) = sig(i, nd)
	4101	END DO
	4102
	4103
[2007]	4104	!AC! do 2100 j=1,ntra
	4105	!AC!c oct3 do 2110 k=1,nl
	4106	!AC! do 2110 k=1,nd ! oct3
	4107	!AC! do 2120 i=1,ncum
	4108	!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
	4109	!AC! 2120 continue
	4110	!AC! 2110 continue
	4111	!AC! 2100 continue
	4112	!
[1992]	4113	RETURN
	4114	END SUBROUTINE cv3_uncompress

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

Download in other formats:

Original Format