Context Navigation

lmdz_cv.F90 @ 5467

Last change on this file since 5467 was 5158, checked in by abarral, 5 months ago
Add missing klon on strataer_emiss_mod.F90 Correct various missing explicit declarations Replace tabs by spaces (tabs are not part of the fortran charset) Continue cleaning modules Removed unused arguments and variables
Property svn:keywords set to `Id`
File size: 53.1 KB

Rev	Line
[4004]	1	! $Id: lmdz_cv.F90 5158 2024-08-02 12:12:03Z fhourdin $
[3927]	2
[5142]	3	MODULE lmdz_cv
	4	!------------------------------------------------------------
	5	! Parameters for convectL:
	6	! (includes - microphysical parameters,
[5158]	7	! - parameters that control the rate of approach
[5142]	8	! to quasi-equilibrium)
[5158]	9	! - noff & minorig (previously in input of convect1)
[5142]	10	!------------------------------------------------------------
	11
	12	IMPLICIT NONE; PRIVATE
	13	PUBLIC elcrit, tlcrit, entp, sigs, sigd, omtrain, omtsnow, coeffr, coeffs &
	14	, dtmax, cu, betad, alpha, damp, delta, noff, minorig, nl, nlp, nlm, &
	15	cv_param, cv_prelim, cv_feed, cv_undilute1, cv_trigger, cv_compress, &
	16	cv_undilute2, cv_closure, cv_mixing, cv_unsat, cv_yield, cv_uncompress
	17
	18	INTEGER noff, minorig, nl, nlp, nlm
	19	REAL elcrit, tlcrit
	20	REAL entp
	21	REAL sigs, sigd
	22	REAL omtrain, omtsnow, coeffr, coeffs
	23	REAL dtmax
	24	REAL cu
	25	REAL betad
	26	REAL alpha, damp
	27	REAL delta
	28
	29	!$OMP THREADPRIVATE(elcrit, tlcrit, entp, sigs, sigd, omtrain, omtsnow, coeffr, coeffs &
	30	!$OMP , dtmax, cu, betad, alpha, damp, delta, noff, minorig, nl, nlp, nlm)
	31
	32	CONTAINS
	33
[3927]	34	SUBROUTINE cv_param(nd)
	35	IMPLICIT NONE
	36
	37	! ------------------------------------------------------------
	38	! Set parameters for convectL
	39	! (includes microphysical parameters and parameters that
	40	! control the rate of approach to quasi-equilibrium)
	41	! ------------------------------------------------------------
	42
	43	! * ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *
	44	! * TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *
	45	! * CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *
	46	! * (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *
	47	! * BETWEEN 0 C AND TLCRIT) *
	48	! * ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *
	49	! * FORMULATION *
	50	! * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
	51	! * SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
	52	! * OF CLOUD *
	53	! * OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *
	54	! * OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *
	55	! * COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *
	56	! * OF RAIN *
	57	! * COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *
	58	! * OF SNOW *
	59	! * CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *
	60	! * TRANSPORT *
	61	! * DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *
	62	! * A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *
	63	! * ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *
	64	! * APPROACH TO QUASI-EQUILIBRIUM *
	65	! * (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *
	66	! * (DAMP MUST BE LESS THAN 1) *
	67
[5142]	68	INTEGER nd
[5132]	69	CHARACTER (LEN = 20) :: modname = 'cv_routines'
	70	CHARACTER (LEN = 80) :: abort_message
[3927]	71
	72	! noff: integer limit for convection (nd-noff)
	73	! minorig: First level of convection
	74
	75	noff = 2
	76	minorig = 2
	77
	78	nl = nd - noff
	79	nlp = nl + 1
	80	nlm = nl - 1
	81
	82	elcrit = 0.0011
	83	tlcrit = -55.0
	84	entp = 1.5
	85	sigs = 0.12
	86	sigd = 0.05
	87	omtrain = 50.0
	88	omtsnow = 5.5
	89	coeffr = 1.0
	90	coeffs = 0.8
	91	dtmax = 0.9
	92
	93	cu = 0.70
	94
	95	betad = 10.0
	96
	97	damp = 0.1
	98	alpha = 0.2
	99
	100	delta = 0.01 ! cld
	101
	102	END SUBROUTINE cv_param
	103
	104	SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm)
[5141]	105	USE lmdz_cvthermo
	106
[3927]	107	IMPLICIT NONE
	108
	109	! =====================================================================
	110	! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
	111	! =====================================================================
	112
	113	! inputs:
	114	INTEGER len, nd, ndp1
	115	REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
	116
	117	! outputs:
	118	REAL lv(len, nd), cpn(len, nd), tv(len, nd)
	119	REAL gz(len, nd), h(len, nd), hm(len, nd)
	120
	121	! local variables:
	122	INTEGER k, i
	123	REAL cpx(len, nd)
	124
[5142]	125
[3927]	126	DO k = 1, nlp
	127	DO i = 1, len
[5132]	128	lv(i, k) = lv0 - clmcpv * (t(i, k) - t0)
	129	cpn(i, k) = cpd * (1.0 - q(i, k)) + cpv * q(i, k)
	130	cpx(i, k) = cpd * (1.0 - q(i, k)) + cl * q(i, k)
	131	tv(i, k) = t(i, k) * (1.0 + q(i, k) * epsim1)
[3927]	132	END DO
	133	END DO
	134
	135	! gz = phi at the full levels (same as p).
	136
	137	DO i = 1, len
	138	gz(i, 1) = 0.0
	139	END DO
	140	DO k = 2, nlp
	141	DO i = 1, len
[5132]	142	gz(i, k) = gz(i, k - 1) + hrd * (tv(i, k - 1) + tv(i, k)) * (p(i, k - 1) - p(i, k)) / ph(i, &
	143	k)
[3927]	144	END DO
	145	END DO
	146
	147	! h = phi + cpT (dry static energy).
	148	! hm = phi + cp(T-Tbase)+Lq
	149
	150	DO k = 1, nlp
	151	DO i = 1, len
[5132]	152	h(i, k) = gz(i, k) + cpn(i, k) * t(i, k)
	153	hm(i, k) = gz(i, k) + cpx(i, k) * (t(i, k) - t(i, 1)) + lv(i, k) * q(i, k)
[3927]	154	END DO
	155	END DO
	156
	157	END SUBROUTINE cv_prelim
	158
	159	SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, &
[5132]	160	qnk, gznk, plcl)
[3927]	161	IMPLICIT NONE
	162
	163	! ================================================================
	164	! Purpose: CONVECTIVE FEED
	165	! ================================================================
	166
[5142]	167
[3927]	168	! inputs:
	169	INTEGER len, nd
	170	REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
	171	REAL hm(len, nd), gz(len, nd)
	172
	173	! outputs:
	174	INTEGER iflag(len), nk(len), icb(len), icbmax
	175	REAL tnk(len), qnk(len), gznk(len), plcl(len)
	176
	177	! local variables:
	178	INTEGER i, k
	179	INTEGER ihmin(len)
	180	REAL work(len)
	181	REAL pnk(len), qsnk(len), rh(len), chi(len)
	182
	183	! -------------------------------------------------------------------
	184	! --- Find level of minimum moist static energy
	185	! --- If level of minimum moist static energy coincides with
	186	! --- or is lower than minimum allowable parcel origin level,
	187	! --- set iflag to 6.
	188	! -------------------------------------------------------------------
	189
	190	DO i = 1, len
	191	work(i) = 1.0E12
	192	ihmin(i) = nl
	193	END DO
	194	DO k = 2, nlp
	195	DO i = 1, len
[5132]	196	IF ((hm(i, k)<work(i)) .AND. (hm(i, k)<hm(i, k - 1))) THEN
[3927]	197	work(i) = hm(i, k)
	198	ihmin(i) = k
	199	END IF
	200	END DO
	201	END DO
	202	DO i = 1, len
	203	ihmin(i) = min(ihmin(i), nlm)
	204	IF (ihmin(i)<=minorig) THEN
	205	iflag(i) = 6
	206	END IF
	207	END DO
	208
	209	! -------------------------------------------------------------------
	210	! --- Find that model level below the level of minimum moist static
	211	! --- energy that has the maximum value of moist static energy
	212	! -------------------------------------------------------------------
	213
	214	DO i = 1, len
	215	work(i) = hm(i, minorig)
	216	nk(i) = minorig
	217	END DO
	218	DO k = minorig + 1, nl
	219	DO i = 1, len
[5132]	220	IF ((hm(i, k)>work(i)) .AND. (k<=ihmin(i))) THEN
[3927]	221	work(i) = hm(i, k)
	222	nk(i) = k
	223	END IF
	224	END DO
	225	END DO
	226	! -------------------------------------------------------------------
	227	! --- Check whether parcel level temperature and specific humidity
	228	! --- are reasonable
	229	! -------------------------------------------------------------------
	230	DO i = 1, len
[5132]	231	IF (((t(i, nk(i))<250.0) .OR. (q(i, nk(i))<=0.0) .OR. (p(i, ihmin(i))< &
	232	400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
[3927]	233	END DO
	234	! -------------------------------------------------------------------
	235	! --- Calculate lifted condensation level of air at parcel origin level
	236	! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
	237	! -------------------------------------------------------------------
	238	DO i = 1, len
	239	tnk(i) = t(i, nk(i))
	240	qnk(i) = q(i, nk(i))
	241	gznk(i) = gz(i, nk(i))
	242	pnk(i) = p(i, nk(i))
	243	qsnk(i) = qs(i, nk(i))
	244
[5132]	245	rh(i) = qnk(i) / qsnk(i)
[3927]	246	rh(i) = min(1.0, rh(i))
[5132]	247	chi(i) = tnk(i) / (1669.0 - 122.0 * rh(i) - tnk(i))
	248	plcl(i) = pnk(i) * (rh(i)**chi(i))
[3927]	249	IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
[5132]	250	) = 8
[3927]	251	END DO
	252	! -------------------------------------------------------------------
	253	! --- Calculate first level above lcl (=icb)
	254	! -------------------------------------------------------------------
	255	DO i = 1, len
	256	icb(i) = nlm
	257	END DO
	258
	259	DO k = minorig, nl
	260	DO i = 1, len
[5132]	261	IF ((k>=(nk(i) + 1)) .AND. (p(i, k)<plcl(i))) icb(i) = min(icb(i), k)
[3927]	262	END DO
	263	END DO
	264
	265	DO i = 1, len
	266	IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
	267	END DO
	268
	269	! Compute icbmax.
	270
	271	icbmax = 2
	272	DO i = 1, len
	273	icbmax = max(icbmax, icb(i))
	274	END DO
	275
	276	END SUBROUTINE cv_feed
	277
	278	SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, &
[5132]	279	clw)
[5141]	280	USE lmdz_cvthermo
	281
[3927]	282	IMPLICIT NONE
	283
[5142]	284
[3927]	285	! inputs:
	286	INTEGER len, nd
	287	INTEGER nk(len), icb(len), icbmax
	288	REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
	289	REAL p(len, nd)
	290
	291	! outputs:
	292	REAL tp(len, nd), tvp(len, nd), clw(len, nd)
	293
	294	! local variables:
	295	INTEGER i, k
	296	REAL tg, qg, alv, s, ahg, tc, denom, es, rg
	297	REAL ah0(len), cpp(len)
	298	REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
	299
	300	! -------------------------------------------------------------------
	301	! --- Calculates the lifted parcel virtual temperature at nk,
	302	! --- the actual temperature, and the adiabatic
	303	! --- liquid water content. The procedure is to solve the equation.
	304	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
	305	! -------------------------------------------------------------------
	306
	307	DO i = 1, len
	308	tnk(i) = t(i, nk(i))
	309	qnk(i) = q(i, nk(i))
	310	gznk(i) = gz(i, nk(i))
	311	ticb(i) = t(i, icb(i))
	312	gzicb(i) = gz(i, icb(i))
	313	END DO
	314
	315	! * Calculate certain parcel quantities, including static energy *
	316
	317	DO i = 1, len
[5132]	318	ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
	319	273.15)) + gznk(i)
	320	cpp(i) = cpd * (1. - qnk(i)) + qnk(i) * cpv
[3927]	321	END DO
	322
	323	! * Calculate lifted parcel quantities below cloud base *
	324
	325	DO k = minorig, icbmax - 1
	326	DO i = 1, len
[5132]	327	tp(i, k) = tnk(i) - (gz(i, k) - gznk(i)) / cpp(i)
	328	tvp(i, k) = tp(i, k) * (1. + qnk(i) * epsi)
[3927]	329	END DO
	330	END DO
	331
	332	! * Find lifted parcel quantities above cloud base *
	333
	334	DO i = 1, len
	335	tg = ticb(i)
	336	qg = qs(i, icb(i))
[5132]	337	alv = lv0 - clmcpv * (ticb(i) - t0)
[3927]	338
	339	! First iteration.
	340
[5132]	341	s = cpd + alv * alv * qg / (rrv * ticb(i) * ticb(i))
	342	s = 1. / s
	343	ahg = cpd * tg + (cl - cpd) * qnk(i) * ticb(i) + alv * qg + gzicb(i)
	344	tg = tg + s * (ah0(i) - ahg)
[3927]	345	tg = max(tg, 35.0)
	346	tc = tg - t0
	347	denom = 243.5 + tc
	348	IF (tc>=0.0) THEN
[5132]	349	es = 6.112 * exp(17.67 * tc / denom)
[3927]	350	ELSE
[5132]	351	es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
[3927]	352	END IF
[5132]	353	qg = eps * es / (p(i, icb(i)) - es * (1. - eps))
[3927]	354
	355	! Second iteration.
	356
[5132]	357	s = cpd + alv * alv * qg / (rrv * ticb(i) * ticb(i))
	358	s = 1. / s
	359	ahg = cpd * tg + (cl - cpd) * qnk(i) * ticb(i) + alv * qg + gzicb(i)
	360	tg = tg + s * (ah0(i) - ahg)
[3927]	361	tg = max(tg, 35.0)
	362	tc = tg - t0
	363	denom = 243.5 + tc
	364	IF (tc>=0.0) THEN
[5132]	365	es = 6.112 * exp(17.67 * tc / denom)
[3927]	366	ELSE
[5132]	367	es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
[3927]	368	END IF
[5132]	369	qg = eps * es / (p(i, icb(i)) - es * (1. - eps))
[3927]	370
[5132]	371	alv = lv0 - clmcpv * (ticb(i) - 273.15)
	372	tp(i, icb(i)) = (ah0(i) - (cl - cpd) * qnk(i) * ticb(i) - gz(i, icb(i)) - alv * qg) / cpd
[3927]	373	clw(i, icb(i)) = qnk(i) - qg
[5132]	374	clw(i, icb(i)) = max(0.0, clw(i, icb(i)))
	375	rg = qg / (1. - qnk(i))
	376	tvp(i, icb(i)) = tp(i, icb(i)) * (1. + rg * epsi)
[3927]	377	END DO
	378
	379	DO k = minorig, icbmax
	380	DO i = 1, len
[5132]	381	tvp(i, k) = tvp(i, k) - tp(i, k) * qnk(i)
[3927]	382	END DO
	383	END DO
	384
	385	END SUBROUTINE cv_undilute1
	386
	387	SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
	388	IMPLICIT NONE
	389
	390	! -------------------------------------------------------------------
	391	! --- Test for instability.
	392	! --- If there was no convection at last time step and parcel
	393	! --- is stable at icb, then set iflag to 4.
	394	! -------------------------------------------------------------------
	395
[5142]	396
[3927]	397	! inputs:
	398	INTEGER len, nd, icb(len)
	399	REAL cbmf(len), tv(len, nd), tvp(len, nd)
	400
	401	! outputs:
	402	INTEGER iflag(len) ! also an input
	403
	404	! local variables:
	405	INTEGER i
	406
	407	DO i = 1, len
	408	IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
[5132]	409	icb(i))<=(tv(i, icb(i)) - dtmax))) iflag(i) = 4
[3927]	410	END DO
	411
	412	END SUBROUTINE cv_trigger
	413
	414	SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
[5132]	415	tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
	416	tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
	417	v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
[5112]	418	USE lmdz_print_control, ONLY: lunout
[5132]	419	USE lmdz_abort_physic, ONLY: abort_physic
[3927]	420	IMPLICIT NONE
	421
[5142]	422
[3927]	423	! inputs:
	424	INTEGER len, ncum, nd, nloc
	425	INTEGER iflag1(len), nk1(len), icb1(len)
	426	REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len)
	427	REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
	428	REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
[5132]	429	REAL p1(len, nd), ph1(len, nd + 1), tv1(len, nd), tp1(len, nd)
[3927]	430	REAL tvp1(len, nd), clw1(len, nd)
	431
	432	! outputs:
	433	INTEGER iflag(nloc), nk(nloc), icb(nloc)
	434	REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
	435	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
	436	REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
[5132]	437	REAL p(nloc, nd), ph(nloc, nd + 1), tv(nloc, nd), tp(nloc, nd)
[3927]	438	REAL tvp(nloc, nd), clw(nloc, nd)
	439	REAL dph(nloc, nd)
	440
	441	! local variables:
	442	INTEGER i, k, nn
[5132]	443	CHARACTER (LEN = 20) :: modname = 'cv_compress'
	444	CHARACTER (LEN = 80) :: abort_message
[3927]	445
	446	DO k = 1, nl + 1
	447	nn = 0
	448	DO i = 1, len
	449	IF (iflag1(i)==0) THEN
	450	nn = nn + 1
	451	t(nn, k) = t1(i, k)
	452	q(nn, k) = q1(i, k)
	453	qs(nn, k) = qs1(i, k)
	454	u(nn, k) = u1(i, k)
	455	v(nn, k) = v1(i, k)
	456	gz(nn, k) = gz1(i, k)
	457	h(nn, k) = h1(i, k)
	458	lv(nn, k) = lv1(i, k)
	459	cpn(nn, k) = cpn1(i, k)
	460	p(nn, k) = p1(i, k)
	461	ph(nn, k) = ph1(i, k)
	462	tv(nn, k) = tv1(i, k)
	463	tp(nn, k) = tp1(i, k)
	464	tvp(nn, k) = tvp1(i, k)
	465	clw(nn, k) = clw1(i, k)
	466	END IF
	467	END DO
	468	END DO
	469
	470	IF (nn/=ncum) THEN
	471	WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
	472	abort_message = ''
	473	CALL abort_physic(modname, abort_message, 1)
	474	END IF
	475
	476	nn = 0
	477	DO i = 1, len
	478	IF (iflag1(i)==0) THEN
	479	nn = nn + 1
	480	cbmf(nn) = cbmf1(i)
	481	plcl(nn) = plcl1(i)
	482	tnk(nn) = tnk1(i)
	483	qnk(nn) = qnk1(i)
	484	gznk(nn) = gznk1(i)
	485	nk(nn) = nk1(i)
	486	icb(nn) = icb1(i)
	487	iflag(nn) = iflag1(i)
	488	END IF
	489	END DO
	490
	491	DO k = 1, nl
	492	DO i = 1, ncum
[5132]	493	dph(i, k) = ph(i, k) - ph(i, k + 1)
[3927]	494	END DO
	495	END DO
	496
	497	END SUBROUTINE cv_compress
	498
	499	SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
[5132]	500	gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
[5141]	501	USE lmdz_cvthermo
	502
[3927]	503	IMPLICIT NONE
	504
	505	! ---------------------------------------------------------------------
	506	! Purpose:
	507	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	508	! &
	509	! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
	510	! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
	511	! &
	512	! FIND THE LEVEL OF NEUTRAL BUOYANCY
	513	! ---------------------------------------------------------------------
	514
[5142]	515
[3927]	516	! inputs:
	517	INTEGER ncum, nd, nloc
	518	INTEGER icb(nloc), nk(nloc)
	519	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
	520	REAL p(nloc, nd), dph(nloc, nd)
	521	REAL tnk(nloc), qnk(nloc), gznk(nloc)
	522	REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
	523
	524	! outputs:
	525	INTEGER inb(nloc), inb1(nloc)
	526	REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
	527	REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
	528	REAL frac(nloc)
	529
	530	! local variables:
	531	INTEGER i, k
	532	REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
	533	REAL by, defrac
	534	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
	535	LOGICAL lcape(nloc)
	536
	537	! =====================================================================
	538	! --- SOME INITIALIZATIONS
	539	! =====================================================================
	540
	541	DO k = 1, nl
	542	DO i = 1, ncum
	543	ep(i, k) = 0.0
	544	sigp(i, k) = sigs
	545	END DO
	546	END DO
	547
	548	! =====================================================================
	549	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	550	! =====================================================================
	551
	552	! --- The procedure is to solve the equation.
	553	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
	554
	555	! * Calculate certain parcel quantities, including static energy *
	556
	557	DO i = 1, ncum
[5132]	558	ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
	559	t0)) + gznk(i)
[3927]	560	END DO
	561
	562
	563	! * Find lifted parcel quantities above cloud base *
	564
	565	DO k = minorig + 1, nl
	566	DO i = 1, ncum
[5132]	567	IF (k>=(icb(i) + 1)) THEN
[3927]	568	tg = t(i, k)
	569	qg = qs(i, k)
[5132]	570	alv = lv0 - clmcpv * (t(i, k) - t0)
[3927]	571
	572	! First iteration.
	573
[5132]	574	s = cpd + alv * alv * qg / (rrv * t(i, k) * t(i, k))
	575	s = 1. / s
	576	ahg = cpd * tg + (cl - cpd) * qnk(i) * t(i, k) + alv * qg + gz(i, k)
	577	tg = tg + s * (ah0(i) - ahg)
[3927]	578	tg = max(tg, 35.0)
	579	tc = tg - t0
	580	denom = 243.5 + tc
	581	IF (tc>=0.0) THEN
[5132]	582	es = 6.112 * exp(17.67 * tc / denom)
[3927]	583	ELSE
[5132]	584	es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
[3927]	585	END IF
[5132]	586	qg = eps * es / (p(i, k) - es * (1. - eps))
[3927]	587
	588	! Second iteration.
	589
[5132]	590	s = cpd + alv * alv * qg / (rrv * t(i, k) * t(i, k))
	591	s = 1. / s
	592	ahg = cpd * tg + (cl - cpd) * qnk(i) * t(i, k) + alv * qg + gz(i, k)
	593	tg = tg + s * (ah0(i) - ahg)
[3927]	594	tg = max(tg, 35.0)
	595	tc = tg - t0
	596	denom = 243.5 + tc
	597	IF (tc>=0.0) THEN
[5132]	598	es = 6.112 * exp(17.67 * tc / denom)
[3927]	599	ELSE
[5132]	600	es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
[3927]	601	END IF
[5132]	602	qg = eps * es / (p(i, k) - es * (1. - eps))
[3927]	603
[5132]	604	alv = lv0 - clmcpv * (t(i, k) - t0)
[5103]	605	! PRINT*,'cpd dans convect2 ',cpd
	606	! PRINT*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
	607	! PRINT*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
[5132]	608	tp(i, k) = (ah0(i) - (cl - cpd) * qnk(i) * t(i, k) - gz(i, k) - alv * qg) / cpd
[5117]	609	! if (.NOT.cpd.gt.1000.) THEN
[5103]	610	! PRINT*,'CPD=',cpd
[3927]	611	! stop
[5117]	612	! END IF
[3927]	613	clw(i, k) = qnk(i) - qg
[5132]	614	clw(i, k) = max(0.0, clw(i, k))
	615	rg = qg / (1. - qnk(i))
	616	tvp(i, k) = tp(i, k) * (1. + rg * epsi)
[3927]	617	END IF
	618	END DO
	619	END DO
	620
	621	! =====================================================================
	622	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
	623	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
	624	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
	625	! =====================================================================
	626
	627	DO k = minorig + 1, nl
	628	DO i = 1, ncum
[5132]	629	IF (k>=(nk(i) + 1)) THEN
[3927]	630	tca = tp(i, k) - t0
	631	IF (tca>=0.0) THEN
	632	elacrit = elcrit
	633	ELSE
[5132]	634	elacrit = elcrit * (1.0 - tca / tlcrit)
[3927]	635	END IF
	636	elacrit = max(elacrit, 0.0)
[5132]	637	ep(i, k) = 1.0 - elacrit / max(clw(i, k), 1.0E-8)
	638	ep(i, k) = max(ep(i, k), 0.0)
	639	ep(i, k) = min(ep(i, k), 1.0)
[3927]	640	sigp(i, k) = sigs
	641	END IF
	642	END DO
	643	END DO
	644
	645	! =====================================================================
	646	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
	647	! --- VIRTUAL TEMPERATURE
	648	! =====================================================================
	649
	650	DO k = minorig + 1, nl
	651	DO i = 1, ncum
[5132]	652	IF (k>=(icb(i) + 1)) THEN
	653	tvp(i, k) = tvp(i, k) * (1.0 - qnk(i) + ep(i, k) * clw(i, k))
[5103]	654	! PRINT*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
	655	! PRINT*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
[3927]	656	END IF
	657	END DO
	658	END DO
	659	DO i = 1, ncum
[5132]	660	tvp(i, nlp) = tvp(i, nl) - (gz(i, nlp) - gz(i, nl)) / cpd
[3927]	661	END DO
	662
	663	! =====================================================================
	664	! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
	665	! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
	666	! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
	667	! =====================================================================
	668
	669	DO i = 1, ncum
	670	cape(i) = 0.0
	671	capem(i) = 0.0
	672	inb(i) = icb(i) + 1
	673	inb1(i) = inb(i)
	674	END DO
	675
	676	! Originial Code
	677
	678	! do 530 k=minorig+1,nl-1
	679	! do 520 i=1,ncum
[5116]	680	! IF(k.ge.(icb(i)+1))THEN
[3927]	681	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	682	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	683	! cape(i)=cape(i)+by
[5116]	684	! IF(by.ge.0.0)inb1(i)=k+1
	685	! IF(cape(i).gt.0.0)THEN
[3927]	686	! inb(i)=k+1
	687	! capem(i)=cape(i)
[5117]	688	! END IF
	689	! END IF
[3927]	690	! 520 continue
	691	! 530 continue
	692	! do 540 i=1,ncum
	693	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
	694	! cape(i)=capem(i)+byp
	695	! defrac=capem(i)-cape(i)
	696	! defrac=max(defrac,0.001)
	697	! frac(i)=-cape(i)/defrac
	698	! frac(i)=min(frac(i),1.0)
	699	! frac(i)=max(frac(i),0.0)
	700	! 540 continue
	701
	702	! K Emanuel fix
	703
[5103]	704	! CALL zilch(byp,ncum)
[3927]	705	! do 530 k=minorig+1,nl-1
	706	! do 520 i=1,ncum
[5116]	707	! IF(k.ge.(icb(i)+1))THEN
[3927]	708	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	709	! cape(i)=cape(i)+by
[5116]	710	! IF(by.ge.0.0)inb1(i)=k+1
	711	! IF(cape(i).gt.0.0)THEN
[3927]	712	! inb(i)=k+1
	713	! capem(i)=cape(i)
	714	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
[5117]	715	! END IF
	716	! END IF
[3927]	717	! 520 continue
	718	! 530 continue
	719	! do 540 i=1,ncum
	720	! inb(i)=max(inb(i),inb1(i))
	721	! cape(i)=capem(i)+byp(i)
	722	! defrac=capem(i)-cape(i)
	723	! defrac=max(defrac,0.001)
	724	! frac(i)=-cape(i)/defrac
	725	! frac(i)=min(frac(i),1.0)
	726	! frac(i)=max(frac(i),0.0)
	727	! 540 continue
	728
	729	! J Teixeira fix
	730
	731	CALL zilch(byp, ncum)
	732	DO i = 1, ncum
	733	lcape(i) = .TRUE.
	734	END DO
	735	DO k = minorig + 1, nl - 1
	736	DO i = 1, ncum
	737	IF (cape(i)<0.0) lcape(i) = .FALSE.
[5132]	738	IF ((k>=(icb(i) + 1)) .AND. lcape(i)) THEN
	739	by = (tvp(i, k) - tv(i, k)) * dph(i, k) / p(i, k)
	740	byp(i) = (tvp(i, k + 1) - tv(i, k + 1)) * dph(i, k + 1) / p(i, k + 1)
[3927]	741	cape(i) = cape(i) + by
	742	IF (by>=0.0) inb1(i) = k + 1
	743	IF (cape(i)>0.0) THEN
	744	inb(i) = k + 1
	745	capem(i) = cape(i)
	746	END IF
	747	END IF
	748	END DO
	749	END DO
	750	DO i = 1, ncum
	751	cape(i) = capem(i) + byp(i)
	752	defrac = capem(i) - cape(i)
	753	defrac = max(defrac, 0.001)
[5132]	754	frac(i) = -cape(i) / defrac
[3927]	755	frac(i) = min(frac(i), 1.0)
	756	frac(i) = max(frac(i), 0.0)
	757	END DO
	758
	759	! =====================================================================
	760	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
	761	! =====================================================================
	762
	763	! initialization:
[5132]	764	DO i = 1, ncum * nlp
[3927]	765	hp(i, 1) = h(i, 1)
	766	END DO
	767
	768	DO k = minorig + 1, nl
	769	DO i = 1, ncum
	770	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
[5132]	771	hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k &
	772	)
[3927]	773	END IF
	774	END DO
	775	END DO
	776
	777	END SUBROUTINE cv_undilute2
	778
	779	SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
[5132]	780	cpn, iflag, cbmf)
[5141]	781	USE lmdz_cvthermo
	782
[3927]	783	IMPLICIT NONE
	784
	785	! inputs:
	786	INTEGER ncum, nd, nloc
	787	INTEGER nk(nloc), icb(nloc)
	788	REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
[5132]	789	REAL ph(nloc, nd + 1) ! caution nd instead ndp1 to be consistent...
[3927]	790	REAL plcl(nloc), cpn(nloc, nd)
	791
	792	! outputs:
	793	INTEGER iflag(nloc)
	794	REAL cbmf(nloc) ! also an input
	795
	796	! local variables:
	797	INTEGER i, k, icbmax
	798	REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
	799	REAL work(nloc)
	800
[5142]	801
[3927]	802	! -------------------------------------------------------------------
	803	! Compute icbmax.
	804	! -------------------------------------------------------------------
	805
	806	icbmax = 2
	807	DO i = 1, ncum
	808	icbmax = max(icbmax, icb(i))
	809	END DO
	810
	811	! =====================================================================
	812	! --- CALCULATE CLOUD BASE MASS FLUX
	813	! =====================================================================
	814
	815	! tvpplcl = parcel temperature lifted adiabatically from level
	816	! icb-1 to the LCL.
	817	! tvaplcl = virtual temperature at the LCL.
	818
	819	DO i = 1, ncum
	820	dtpbl(i) = 0.0
[5132]	821	tvpplcl(i) = tvp(i, icb(i) - 1) - rrd * tvp(i, icb(i) - 1) * (p(i, icb(i) - 1) - plcl(&
	822	i)) / (cpn(i, icb(i) - 1) * p(i, icb(i) - 1))
	823	tvaplcl(i) = tv(i, icb(i)) + (tvp(i, icb(i)) - tvp(i, icb(i) + 1)) * (plcl(i) - p(i &
	824	, icb(i))) / (p(i, icb(i)) - p(i, icb(i) + 1))
[3927]	825	END DO
	826
	827	! -------------------------------------------------------------------
	828	! --- Interpolate difference between lifted parcel and
	829	! --- environmental temperatures to lifted condensation level
	830	! -------------------------------------------------------------------
	831
	832	! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
	833
	834	DO k = minorig, icbmax
	835	DO i = 1, ncum
[5132]	836	IF ((k>=nk(i)) .AND. (k<=(icb(i) - 1))) THEN
	837	dtpbl(i) = dtpbl(i) + (tvp(i, k) - tv(i, k)) * dph(i, k)
[3927]	838	END IF
	839	END DO
	840	END DO
	841	DO i = 1, ncum
[5132]	842	dtpbl(i) = dtpbl(i) / (ph(i, nk(i)) - ph(i, icb(i)))
[3927]	843	dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
	844	END DO
	845
	846	! -------------------------------------------------------------------
	847	! --- Adjust cloud base mass flux
	848	! -------------------------------------------------------------------
	849
	850	DO i = 1, ncum
	851	work(i) = cbmf(i)
[5132]	852	cbmf(i) = max(0.0, (1.0 - damp) * cbmf(i) + 0.1 * alpha * dtmin(i))
[3927]	853	IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
	854	iflag(i) = 3
	855	END IF
	856	END DO
	857
	858	END SUBROUTINE cv_closure
	859
	860	SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
[5132]	861	h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
	862	sij, elij)
[5141]	863	USE lmdz_cvthermo
	864
[3927]	865	IMPLICIT NONE
	866
[5142]	867
[3927]	868	! inputs:
	869	INTEGER ncum, nd, nloc
	870	INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
	871	REAL cbmf(nloc), qnk(nloc)
[5132]	872	REAL ph(nloc, nd + 1)
[3927]	873	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
	874	REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
	875	REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd)
	876
	877	! outputs:
	878	INTEGER nent(nloc, nd)
	879	REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd)
	880	REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
	881	REAL sij(nloc, nd, nd), elij(nloc, nd, nd)
	882
	883	! local variables:
	884	INTEGER i, j, k, ij
	885	INTEGER num1, num2
	886	REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
	887	REAL alt, qp1, smid, sjmin, sjmax, delp, delm
	888	REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc)
	889	REAL bsum(nloc, nd)
	890	LOGICAL lwork(nloc)
	891
	892	! =====================================================================
	893	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	894	! =====================================================================
	895
[5132]	896	DO i = 1, ncum * nlp
[3927]	897	nent(i, 1) = 0
	898	m(i, 1) = 0.0
	899	END DO
	900
	901	DO k = 1, nlp
	902	DO j = 1, nlp
	903	DO i = 1, ncum
	904	qent(i, k, j) = q(i, j)
	905	uent(i, k, j) = u(i, j)
	906	vent(i, k, j) = v(i, j)
	907	elij(i, k, j) = 0.0
	908	ment(i, k, j) = 0.0
	909	sij(i, k, j) = 0.0
	910	END DO
	911	END DO
	912	END DO
	913
	914	! -------------------------------------------------------------------
	915	! --- Calculate rates of mixing, m(i)
	916	! -------------------------------------------------------------------
	917
	918	CALL zilch(work, ncum)
	919
	920	DO j = minorig + 1, nl
	921	DO i = 1, ncum
[5132]	922	IF ((j>=(icb(i) + 1)) .AND. (j<=inb(i))) THEN
[3927]	923	k = min(j, inb1(i))
[5132]	924	dbo = abs(tv(i, k + 1) - tvp(i, k + 1) - tv(i, k - 1) + tvp(i, k - 1)) + &
	925	entp * 0.04 * (ph(i, k) - ph(i, k + 1))
[3927]	926	work(i) = work(i) + dbo
[5132]	927	m(i, j) = cbmf(i) * dbo
[3927]	928	END IF
	929	END DO
	930	END DO
	931	DO k = minorig + 1, nl
	932	DO i = 1, ncum
[5132]	933	IF ((k>=(icb(i) + 1)) .AND. (k<=inb(i))) THEN
	934	m(i, k) = m(i, k) / work(i)
[3927]	935	END IF
	936	END DO
	937	END DO
	938
	939
	940	! =====================================================================
	941	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	942	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	943	! --- FRACTION (sij)
	944	! =====================================================================
	945
	946	DO i = minorig + 1, nl
	947	DO j = minorig + 1, nl
	948	DO ij = 1, ncum
[5132]	949	IF ((i>=(icb(ij) + 1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
	950	inb(ij))) THEN
	951	qti = qnk(ij) - ep(ij, i) * clw(ij, i)
	952	bf2 = 1. + lv(ij, j) * lv(ij, j) * qs(ij, j) / (rrv * t(ij, j) * t(ij, j) * cpd)
	953	anum = h(ij, j) - hp(ij, i) + (cpv - cpd) * t(ij, j) * (qti - q(ij, j))
	954	denom = h(ij, i) - hp(ij, i) + (cpd - cpv) * (q(ij, i) - qti) * t(ij, j)
[3927]	955	dei = denom
	956	IF (abs(dei)<0.01) dei = 0.01
[5132]	957	sij(ij, i, j) = anum / dei
[3927]	958	sij(ij, i, i) = 1.0
[5132]	959	altem = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti - qs(ij, j)
	960	altem = altem / bf2
	961	cwat = clw(ij, j) * (1. - ep(ij, j))
[3927]	962	stemp = sij(ij, i, j)
	963	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
[5132]	964	anum = anum - lv(ij, j) * (qti - qs(ij, j) - cwat * bf2)
	965	denom = denom + lv(ij, j) * (q(ij, i) - qti)
[3927]	966	IF (abs(denom)<0.01) denom = 0.01
[5132]	967	sij(ij, i, j) = anum / denom
	968	altem = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti - qs(ij, j)
	969	altem = altem - (bf2 - 1.) * cwat
[3927]	970	END IF
[5132]	971	IF (sij(ij, i, j)>0.0 .AND. sij(ij, i, j)<0.9) THEN
	972	qent(ij, i, j) = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti
	973	uent(ij, i, j) = sij(ij, i, j) * u(ij, i) + &
	974	(1. - sij(ij, i, j)) * u(ij, nk(ij))
	975	vent(ij, i, j) = sij(ij, i, j) * v(ij, i) + &
	976	(1. - sij(ij, i, j)) * v(ij, nk(ij))
[3927]	977	elij(ij, i, j) = altem
[5132]	978	elij(ij, i, j) = max(0.0, elij(ij, i, j))
	979	ment(ij, i, j) = m(ij, i) / (1. - sij(ij, i, j))
[3927]	980	nent(ij, i) = nent(ij, i) + 1
	981	END IF
[5132]	982	sij(ij, i, j) = max(0.0, sij(ij, i, j))
	983	sij(ij, i, j) = min(1.0, sij(ij, i, j))
[3927]	984	END IF
	985	END DO
	986	END DO
	987
	988	! *** If no air can entrain at level i assume that updraft detrains
	989	! ***
	990	! *** at that level and calculate detrained air flux and properties
	991	! ***
	992
	993	DO ij = 1, ncum
[5132]	994	IF ((i>=(icb(ij) + 1)) .AND. (i<=inb(ij)) .AND. (nent(ij, i)==0)) THEN
[3927]	995	ment(ij, i, i) = m(ij, i)
[5132]	996	qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
[3927]	997	uent(ij, i, i) = u(ij, nk(ij))
	998	vent(ij, i, i) = v(ij, nk(ij))
	999	elij(ij, i, i) = clw(ij, i)
	1000	sij(ij, i, i) = 1.0
	1001	END IF
	1002	END DO
	1003	END DO
	1004
	1005	DO i = 1, ncum
	1006	sij(i, inb(i), inb(i)) = 1.0
	1007	END DO
	1008
	1009	! =====================================================================
	1010	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	1011	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	1012	! =====================================================================
	1013
[5132]	1014	CALL zilch(bsum, ncum * nlp)
[3927]	1015	DO ij = 1, ncum
	1016	lwork(ij) = .FALSE.
	1017	END DO
	1018	DO i = minorig + 1, nl
	1019
	1020	num1 = 0
	1021	DO ij = 1, ncum
[5132]	1022	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij))) num1 = num1 + 1
[3927]	1023	END DO
	1024	IF (num1<=0) GO TO 789
	1025
	1026	DO ij = 1, ncum
[5132]	1027	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij))) THEN
	1028	lwork(ij) = (nent(ij, i)/=0)
	1029	qp1 = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
	1030	anum = h(ij, i) - hp(ij, i) - lv(ij, i) * (qp1 - qs(ij, i))
	1031	denom = h(ij, i) - hp(ij, i) + lv(ij, i) * (q(ij, i) - qp1)
[3927]	1032	IF (abs(denom)<0.01) denom = 0.01
[5132]	1033	scrit(ij) = anum / denom
	1034	alt = qp1 - qs(ij, i) + scrit(ij) * (q(ij, i) - qp1)
[3927]	1035	IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
	1036	asij(ij) = 0.0
	1037	smin(ij) = 1.0
	1038	END IF
	1039	END DO
	1040	DO j = minorig, nl
	1041
	1042	num2 = 0
	1043	DO ij = 1, ncum
[5132]	1044	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
	1045	ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
[3927]	1046	END DO
	1047	IF (num2<=0) GO TO 783
	1048
	1049	DO ij = 1, ncum
[5132]	1050	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
	1051	ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
	1052	IF (sij(ij, i, j)>0.0 .AND. sij(ij, i, j)<0.9) THEN
[3927]	1053	IF (j>i) THEN
[5132]	1054	smid = min(sij(ij, i, j), scrit(ij))
[3927]	1055	sjmax = smid
	1056	sjmin = smid
[5132]	1057	IF (smid<smin(ij) .AND. sij(ij, i, j + 1)<smid) THEN
[3927]	1058	smin(ij) = smid
[5132]	1059	sjmax = min(sij(ij, i, j + 1), sij(ij, i, j), scrit(ij))
	1060	sjmin = max(sij(ij, i, j - 1), sij(ij, i, j))
[3927]	1061	sjmin = min(sjmin, scrit(ij))
	1062	END IF
	1063	ELSE
[5132]	1064	sjmax = max(sij(ij, i, j + 1), scrit(ij))
	1065	smid = max(sij(ij, i, j), scrit(ij))
[3927]	1066	sjmin = 0.0
[5132]	1067	IF (j>1) sjmin = sij(ij, i, j - 1)
[3927]	1068	sjmin = max(sjmin, scrit(ij))
	1069	END IF
[5132]	1070	delp = abs(sjmax - smid)
	1071	delm = abs(sjmin - smid)
	1072	asij(ij) = asij(ij) + (delp + delm) * (ph(ij, j) - ph(ij, j + 1))
	1073	ment(ij, i, j) = ment(ij, i, j) * (delp + delm) * (ph(ij, j) - ph(ij, j + 1))
[3927]	1074	END IF
	1075	END IF
	1076	END DO
[5132]	1077	783 END DO
[3927]	1078	DO ij = 1, ncum
[5132]	1079	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
[3927]	1080	asij(ij) = max(1.0E-21, asij(ij))
[5132]	1081	asij(ij) = 1.0 / asij(ij)
[3927]	1082	bsum(ij, i) = 0.0
	1083	END IF
	1084	END DO
	1085	DO j = minorig, nl + 1
	1086	DO ij = 1, ncum
[5132]	1087	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
	1088	ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
	1089	ment(ij, i, j) = ment(ij, i, j) * asij(ij)
[3927]	1090	bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
	1091	END IF
	1092	END DO
	1093	END DO
	1094	DO ij = 1, ncum
[5132]	1095	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
	1096	i)<1.0E-18) .AND. lwork(ij)) THEN
[3927]	1097	nent(ij, i) = 0
	1098	ment(ij, i, i) = m(ij, i)
[5132]	1099	qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
[3927]	1100	uent(ij, i, i) = u(ij, nk(ij))
	1101	vent(ij, i, i) = v(ij, nk(ij))
	1102	elij(ij, i, i) = clw(ij, i)
	1103	sij(ij, i, i) = 1.0
	1104	END IF
	1105	END DO
[5132]	1106	789 END DO
[3927]	1107
	1108	END SUBROUTINE cv_mixing
	1109
	1110	SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
[5132]	1111	ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
[5141]	1112	USE lmdz_cvthermo
	1113
[3927]	1114	IMPLICIT NONE
	1115
[5142]	1116
[3927]	1117	! inputs:
	1118	INTEGER ncum, nd, nloc
	1119	INTEGER inb(nloc)
	1120	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd)
	1121	REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd)
[5132]	1122	REAL p(nloc, nd), ph(nloc, nd + 1), h(nloc, nd)
[3927]	1123	REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd)
	1124	REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd)
	1125
	1126	! outputs:
	1127	INTEGER iflag(nloc) ! also an input
	1128	REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd)
	1129	REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd)
	1130
	1131	! local variables:
	1132	INTEGER i, j, k, ij, num1
	1133	INTEGER jtt(nloc)
	1134	REAL awat, coeff, qsm, afac, sigt, b6, c6, revap
	1135	REAL dhdp, fac, qstm, rat
	1136	REAL wdtrain(nloc)
	1137	LOGICAL lwork(nloc)
	1138
	1139	! =====================================================================
	1140	! --- PRECIPITATING DOWNDRAFT CALCULATION
	1141	! =====================================================================
	1142
	1143	! Initializations:
	1144
	1145	DO i = 1, ncum
	1146	DO k = 1, nl + 1
	1147	wt(i, k) = omtsnow
	1148	mp(i, k) = 0.0
	1149	evap(i, k) = 0.0
	1150	water(i, k) = 0.0
	1151	END DO
	1152	END DO
	1153
	1154	DO i = 1, ncum
	1155	qp(i, 1) = q(i, 1)
	1156	up(i, 1) = u(i, 1)
	1157	vp(i, 1) = v(i, 1)
	1158	END DO
	1159
	1160	DO k = 2, nl + 1
	1161	DO i = 1, ncum
[5132]	1162	qp(i, k) = q(i, k - 1)
	1163	up(i, k) = u(i, k - 1)
	1164	vp(i, k) = v(i, k - 1)
[3927]	1165	END DO
	1166	END DO
	1167
	1168
	1169	! * Check whether ep(inb)=0, if so, skip precipitating *
	1170	! * downdraft calculation *
	1171
	1172
	1173	! * Integrate liquid water equation to find condensed water *
	1174	! * and condensed water flux *
	1175
	1176	DO i = 1, ncum
	1177	jtt(i) = 2
[5132]	1178	IF (ep(i, inb(i))<=0.0001) iflag(i) = 2
[3927]	1179	IF (iflag(i)==0) THEN
	1180	lwork(i) = .TRUE.
	1181	ELSE
	1182	lwork(i) = .FALSE.
	1183	END IF
	1184	END DO
	1185
	1186	! * Begin downdraft loop *
	1187
	1188	CALL zilch(wdtrain, ncum)
	1189	DO i = nl + 1, 1, -1
	1190
	1191	num1 = 0
	1192	DO ij = 1, ncum
	1193	IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1
	1194	END DO
	1195	IF (num1<=0) GO TO 899
	1196
	1197
	1198	! * Calculate detrained precipitation *
	1199
	1200	DO ij = 1, ncum
	1201	IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
[5132]	1202	wdtrain(ij) = g * ep(ij, i) * m(ij, i) * clw(ij, i)
[3927]	1203	END IF
	1204	END DO
	1205
	1206	IF (i>1) THEN
	1207	DO j = 1, i - 1
	1208	DO ij = 1, ncum
	1209	IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
[5132]	1210	awat = elij(ij, j, i) - (1. - ep(ij, i)) * clw(ij, i)
[3927]	1211	awat = max(0.0, awat)
[5132]	1212	wdtrain(ij) = wdtrain(ij) + g * awat * ment(ij, j, i)
[3927]	1213	END IF
	1214	END DO
	1215	END DO
	1216	END IF
	1217
	1218	! * Find rain water and evaporation using provisional *
	1219	! * estimates of qp(i)and qp(i-1) *
	1220
	1221
	1222	! *** Value of terminal velocity and coeffecient of evaporation for snow
	1223	! ***
	1224
	1225	DO ij = 1, ncum
	1226	IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
	1227	coeff = coeffs
	1228	wt(ij, i) = omtsnow
	1229
	1230	! *** Value of terminal velocity and coeffecient of evaporation for
	1231	! rain ***
	1232
[5132]	1233	IF (t(ij, i)>273.0) THEN
[3927]	1234	coeff = coeffr
	1235	wt(ij, i) = omtrain
	1236	END IF
[5132]	1237	qsm = 0.5 * (q(ij, i) + qp(ij, i + 1))
	1238	afac = coeff * ph(ij, i) * (qs(ij, i) - qsm) / (1.0E4 + 2.0E3 * ph(ij, i) * qs(ij, i))
[3927]	1239	afac = max(afac, 0.0)
	1240	sigt = sigp(ij, i)
	1241	sigt = max(0.0, sigt)
	1242	sigt = min(1.0, sigt)
[5132]	1243	b6 = 100. * (ph(ij, i) - ph(ij, i + 1)) * sigt * afac / wt(ij, i)
	1244	c6 = (water(ij, i + 1) * wt(ij, i + 1) + wdtrain(ij) / sigd) / wt(ij, i)
	1245	revap = 0.5 * (-b6 + sqrt(b6 * b6 + 4. * c6))
	1246	evap(ij, i) = sigt * afac * revap
	1247	water(ij, i) = revap * revap
[3927]	1248
	1249	! * Calculate precipitating downdraft mass flux under *
	1250	! * hydrostatic approximation *
	1251
	1252	IF (i>1) THEN
[5132]	1253	dhdp = (h(ij, i) - h(ij, i - 1)) / (p(ij, i - 1) - p(ij, i))
[3927]	1254	dhdp = max(dhdp, 10.0)
[5132]	1255	mp(ij, i) = 100. * ginv * lv(ij, i) * sigd * evap(ij, i) / dhdp
	1256	mp(ij, i) = max(mp(ij, i), 0.0)
[3927]	1257
	1258	! * Add small amount of inertia to downdraft *
	1259
[5132]	1260	fac = 20.0 / (ph(ij, i - 1) - ph(ij, i))
	1261	mp(ij, i) = (fac * mp(ij, i + 1) + mp(ij, i)) / (1. + fac)
[3927]	1262
	1263	! *** Force mp to decrease linearly to zero
	1264	! ***
	1265	! *** between about 950 mb and the surface
	1266	! ***
	1267
[5132]	1268	IF (p(ij, i)>(0.949 * p(ij, 1))) THEN
[3927]	1269	jtt(ij) = max(jtt(ij), i)
[5132]	1270	mp(ij, i) = mp(ij, jtt(ij)) * (p(ij, 1) - p(ij, i)) / &
	1271	(p(ij, 1) - p(ij, jtt(ij)))
[3927]	1272	END IF
	1273	END IF
	1274
	1275	! * Find mixing ratio of precipitating downdraft *
	1276
	1277	IF (i/=inb(ij)) THEN
	1278	IF (i==1) THEN
	1279	qstm = qs(ij, 1)
	1280	ELSE
[5132]	1281	qstm = qs(ij, i - 1)
[3927]	1282	END IF
[5132]	1283	IF (mp(ij, i)>mp(ij, i + 1)) THEN
	1284	rat = mp(ij, i + 1) / mp(ij, i)
	1285	qp(ij, i) = qp(ij, i + 1) * rat + q(ij, i) * (1.0 - rat) + &
	1286	100. * ginv * sigd * (ph(ij, i) - ph(ij, i + 1)) * (evap(ij, i) / mp(ij, i))
	1287	up(ij, i) = up(ij, i + 1) * rat + u(ij, i) * (1. - rat)
	1288	vp(ij, i) = vp(ij, i + 1) * rat + v(ij, i) * (1. - rat)
[3927]	1289	ELSE
[5132]	1290	IF (mp(ij, i + 1)>0.0) THEN
	1291	qp(ij, i) = (gz(ij, i + 1) - gz(ij, i) + qp(ij, i + 1) * (lv(ij, i + 1) + t(ij, &
	1292	i + 1) * (cl - cpd)) + cpd * (t(ij, i + 1) - t(ij, &
	1293	i))) / (lv(ij, i) + t(ij, i) * (cl - cpd))
	1294	up(ij, i) = up(ij, i + 1)
	1295	vp(ij, i) = vp(ij, i + 1)
[3927]	1296	END IF
	1297	END IF
[5132]	1298	qp(ij, i) = min(qp(ij, i), qstm)
	1299	qp(ij, i) = max(qp(ij, i), 0.0)
[3927]	1300	END IF
	1301	END IF
	1302	END DO
[5132]	1303	899 END DO
[3927]	1304
	1305	END SUBROUTINE cv_unsat
	1306
	1307	SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
[5132]	1308	ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
	1309	ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
	1310	precip, cbmf, ft, fq, fu, fv, ma, qcondc)
[5141]	1311	USE lmdz_cvthermo
	1312
[3927]	1313	IMPLICIT NONE
	1314
[5142]	1315
[3927]	1316	! inputs
	1317	INTEGER ncum, nd, nloc
	1318	INTEGER nk(nloc), icb(nloc), inb(nloc)
	1319	INTEGER nent(nloc, nd)
	1320	REAL delt
	1321	REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd)
	1322	REAL gz(nloc, nd)
[5132]	1323	REAL p(nloc, nd), ph(nloc, nd + 1), h(nloc, nd)
[3927]	1324	REAL hp(nloc, nd), lv(nloc, nd)
	1325	REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc)
	1326	REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd)
	1327	REAL up(nloc, nd), vp(nloc, nd)
	1328	REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd)
	1329	REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd)
	1330	REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
	1331	REAL tv(nloc, nd), tvp(nloc, nd)
	1332
	1333	! outputs
	1334	INTEGER iflag(nloc) ! also an input
	1335	REAL cbmf(nloc) ! also an input
	1336	REAL wd(nloc), tprime(nloc), qprime(nloc)
	1337	REAL precip(nloc)
	1338	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	1339	REAL ma(nloc, nd)
	1340	REAL qcondc(nloc, nd)
	1341
	1342	! local variables
	1343	INTEGER i, j, ij, k, num1
	1344	REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti
	1345	REAL work(nloc), am(nloc), amp1(nloc), ad(nloc)
	1346	REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd)
	1347	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
	1348	REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld
	1349
	1350
	1351	! -- initializations:
	1352
[5132]	1353	delti = 1.0 / delt
[3927]	1354
	1355	DO i = 1, ncum
	1356	precip(i) = 0.0
	1357	wd(i) = 0.0
	1358	tprime(i) = 0.0
	1359	qprime(i) = 0.0
	1360	DO k = 1, nl + 1
	1361	ft(i, k) = 0.0
	1362	fu(i, k) = 0.0
	1363	fv(i, k) = 0.0
	1364	fq(i, k) = 0.0
[5132]	1365	lvcp(i, k) = lv(i, k) / cpn(i, k)
[3927]	1366	qcondc(i, k) = 0.0 ! cld
	1367	qcond(i, k) = 0.0 ! cld
	1368	nqcond(i, k) = 0.0 ! cld
	1369	END DO
	1370	END DO
	1371
	1372
	1373	! * Calculate surface precipitation in mm/day *
	1374
	1375	DO i = 1, ncum
	1376	IF (iflag(i)<=1) THEN
	1377	! c precip(i)=precip(i)+wt(i,1)sigdwater(i,1)3600.24000.
	1378	! c & /(rowl*g)
	1379	! c precip(i)=precip(i)*delt/86400.
[5132]	1380	precip(i) = wt(i, 1) * sigd * water(i, 1) * 86400 / g
[3927]	1381	END IF
	1382	END DO
	1383
	1384
	1385	! * Calculate downdraft velocity scale and surface temperature and *
	1386	! * water vapor fluctuations *
	1387
	1388	DO i = 1, ncum
[5132]	1389	wd(i) = betad * abs(mp(i, icb(i))) * 0.01 * rrd * t(i, icb(i)) / (sigd * p(i, icb(i)))
	1390	qprime(i) = 0.5 * (qp(i, 1) - q(i, 1))
	1391	tprime(i) = lv0 * qprime(i) / cpd
[3927]	1392	END DO
	1393
	1394	! * Calculate tendencies of lowest level potential temperature *
	1395	! * and mixing ratio *
	1396
	1397	DO i = 1, ncum
[5132]	1398	work(i) = 0.01 / (ph(i, 1) - ph(i, 2))
[3927]	1399	am(i) = 0.0
	1400	END DO
	1401	DO k = 2, nl
	1402	DO i = 1, ncum
	1403	IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN
	1404	am(i) = am(i) + m(i, k)
	1405	END IF
	1406	END DO
	1407	END DO
	1408	DO i = 1, ncum
[5132]	1409	IF ((g * work(i) * am(i))>=delti) iflag(i) = 1
	1410	ft(i, 1) = ft(i, 1) + g * work(i) * am(i) * (t(i, 2) - t(i, 1) + (gz(i, 2) - gz(i, &
	1411	1)) / cpn(i, 1))
	1412	ft(i, 1) = ft(i, 1) - lvcp(i, 1) * sigd * evap(i, 1)
	1413	ft(i, 1) = ft(i, 1) + sigd * wt(i, 2) * (cl - cpd) * water(i, 2) * (t(i, 2) - t(i, 1)) * &
	1414	work(i) / cpn(i, 1)
	1415	fq(i, 1) = fq(i, 1) + g * mp(i, 2) * (qp(i, 2) - q(i, 1)) * work(i) + &
	1416	sigd * evap(i, 1)
	1417	fq(i, 1) = fq(i, 1) + g * am(i) * (q(i, 2) - q(i, 1)) * work(i)
	1418	fu(i, 1) = fu(i, 1) + g * work(i) * (mp(i, 2) * (up(i, 2) - u(i, 1)) + am(i) * (u(i, &
	1419	2) - u(i, 1)))
	1420	fv(i, 1) = fv(i, 1) + g * work(i) * (mp(i, 2) * (vp(i, 2) - v(i, 1)) + am(i) * (v(i, &
	1421	2) - v(i, 1)))
[3927]	1422	END DO
	1423	DO j = 2, nl
	1424	DO i = 1, ncum
	1425	IF (j<=inb(i)) THEN
[5132]	1426	fq(i, 1) = fq(i, 1) + g * work(i) * ment(i, j, 1) * (qent(i, j, 1) - q(i, 1))
	1427	fu(i, 1) = fu(i, 1) + g * work(i) * ment(i, j, 1) * (uent(i, j, 1) - u(i, 1))
	1428	fv(i, 1) = fv(i, 1) + g * work(i) * ment(i, j, 1) * (vent(i, j, 1) - v(i, 1))
[3927]	1429	END IF
	1430	END DO
	1431	END DO
	1432
	1433	! * Calculate tendencies of potential temperature and mixing ratio *
	1434	! * at levels above the lowest level *
	1435
	1436	! * First find the net saturated updraft and downdraft mass fluxes *
	1437	! * through each level *
	1438
	1439	DO i = 2, nl + 1
	1440
	1441	num1 = 0
	1442	DO ij = 1, ncum
	1443	IF (i<=inb(ij)) num1 = num1 + 1
	1444	END DO
	1445	IF (num1<=0) GO TO 1500
	1446
	1447	CALL zilch(amp1, ncum)
	1448	CALL zilch(ad, ncum)
	1449
	1450	DO k = i + 1, nl + 1
	1451	DO ij = 1, ncum
[5132]	1452	IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij) + 1))) THEN
[3927]	1453	amp1(ij) = amp1(ij) + m(ij, k)
	1454	END IF
	1455	END DO
	1456	END DO
	1457
	1458	DO k = 1, i
	1459	DO j = i + 1, nl + 1
	1460	DO ij = 1, ncum
[5132]	1461	IF ((j<=(inb(ij) + 1)) .AND. (i<=inb(ij))) THEN
[3927]	1462	amp1(ij) = amp1(ij) + ment(ij, k, j)
	1463	END IF
	1464	END DO
	1465	END DO
	1466	END DO
	1467	DO k = 1, i - 1
	1468	DO j = i, nl + 1
	1469	DO ij = 1, ncum
	1470	IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN
	1471	ad(ij) = ad(ij) + ment(ij, j, k)
	1472	END IF
	1473	END DO
	1474	END DO
	1475	END DO
	1476
	1477	DO ij = 1, ncum
	1478	IF (i<=inb(ij)) THEN
[5132]	1479	dpinv = 0.01 / (ph(ij, i) - ph(ij, i + 1))
	1480	cpinv = 1.0 / cpn(ij, i)
[3927]	1481
[5132]	1482	ft(ij, i) = ft(ij, i) + g * dpinv * (amp1(ij) * (t(ij, i + 1) - t(ij, &
	1483	i) + (gz(ij, i + 1) - gz(ij, i)) * cpinv) - ad(ij) * (t(ij, i) - t(ij, &
	1484	i - 1) + (gz(ij, i) - gz(ij, i - 1)) * cpinv)) - sigd * lvcp(ij, i) * evap(ij, i)
	1485	ft(ij, i) = ft(ij, i) + g * dpinv * ment(ij, i, i) * (hp(ij, i) - h(ij, i) + t(ij &
	1486	, i) * (cpv - cpd) * (q(ij, i) - qent(ij, i, i))) * cpinv
	1487	ft(ij, i) = ft(ij, i) + sigd * wt(ij, i + 1) * (cl - cpd) * water(ij, i + 1) * (t(&
	1488	ij, i + 1) - t(ij, i)) * dpinv * cpinv
	1489	fq(ij, i) = fq(ij, i) + g * dpinv * (amp1(ij) * (q(ij, i + 1) - q(ij, &
	1490	i)) - ad(ij) * (q(ij, i) - q(ij, i - 1)))
	1491	fu(ij, i) = fu(ij, i) + g * dpinv * (amp1(ij) * (u(ij, i + 1) - u(ij, &
	1492	i)) - ad(ij) * (u(ij, i) - u(ij, i - 1)))
	1493	fv(ij, i) = fv(ij, i) + g * dpinv * (amp1(ij) * (v(ij, i + 1) - v(ij, &
	1494	i)) - ad(ij) * (v(ij, i) - v(ij, i - 1)))
[3927]	1495	END IF
	1496	END DO
	1497	DO k = 1, i - 1
	1498	DO ij = 1, ncum
	1499	IF (i<=inb(ij)) THEN
[5132]	1500	awat = elij(ij, k, i) - (1. - ep(ij, i)) * clw(ij, i)
[3927]	1501	awat = max(awat, 0.0)
[5132]	1502	fq(ij, i) = fq(ij, i) + g * dpinv * ment(ij, k, i) * (qent(ij, k, i) - awat - q &
	1503	(ij, i))
	1504	fu(ij, i) = fu(ij, i) + g * dpinv * ment(ij, k, i) * (uent(ij, k, i) - u(ij, i &
	1505	))
	1506	fv(ij, i) = fv(ij, i) + g * dpinv * ment(ij, k, i) * (vent(ij, k, i) - v(ij, i &
	1507	))
[3927]	1508	! (saturated updrafts resulting from mixing) ! cld
[5132]	1509	qcond(ij, i) = qcond(ij, i) + (elij(ij, k, i) - awat) ! cld
[3927]	1510	nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
	1511	END IF
	1512	END DO
	1513	END DO
	1514	DO k = i, nl + 1
	1515	DO ij = 1, ncum
	1516	IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
[5132]	1517	fq(ij, i) = fq(ij, i) + g * dpinv * ment(ij, k, i) * (qent(ij, k, i) - q(ij, i &
	1518	))
	1519	fu(ij, i) = fu(ij, i) + g * dpinv * ment(ij, k, i) * (uent(ij, k, i) - u(ij, i &
	1520	))
	1521	fv(ij, i) = fv(ij, i) + g * dpinv * ment(ij, k, i) * (vent(ij, k, i) - v(ij, i &
	1522	))
[3927]	1523	END IF
	1524	END DO
	1525	END DO
	1526	DO ij = 1, ncum
	1527	IF (i<=inb(ij)) THEN
[5132]	1528	fq(ij, i) = fq(ij, i) + sigd * evap(ij, i) + g * (mp(ij, i + 1) * (qp(ij, &
	1529	i + 1) - q(ij, i)) - mp(ij, i) * (qp(ij, i) - q(ij, i - 1))) * dpinv
	1530	fu(ij, i) = fu(ij, i) + g * (mp(ij, i + 1) * (up(ij, i + 1) - u(ij, &
	1531	i)) - mp(ij, i) * (up(ij, i) - u(ij, i - 1))) * dpinv
	1532	fv(ij, i) = fv(ij, i) + g * (mp(ij, i + 1) * (vp(ij, i + 1) - v(ij, &
	1533	i)) - mp(ij, i) * (vp(ij, i) - v(ij, i - 1))) * dpinv
[3927]	1534	! (saturated downdrafts resulting from mixing) ! cld
	1535	DO k = i + 1, inb(ij) ! cld
	1536	qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld
	1537	nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
	1538	END DO ! cld
	1539	! (particular case: no detraining level is found) ! cld
[5132]	1540	IF (nent(ij, i)==0) THEN ! cld
	1541	qcond(ij, i) = qcond(ij, i) + (1. - ep(ij, i)) * clw(ij, i) ! cld
[3927]	1542	nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
	1543	END IF ! cld
[5132]	1544	IF (nqcond(ij, i)/=0.) THEN ! cld
	1545	qcond(ij, i) = qcond(ij, i) / nqcond(ij, i) ! cld
[3927]	1546	END IF ! cld
	1547	END IF
	1548	END DO
[5132]	1549	1500 END DO
[3927]	1550
	1551	! * Adjust tendencies at top of convection layer to reflect *
	1552	! * actual position of the level zero cape *
	1553
	1554	DO ij = 1, ncum
	1555	fqold = fq(ij, inb(ij))
[5132]	1556	fq(ij, inb(ij)) = fq(ij, inb(ij)) * (1. - frac(ij))
	1557	fq(ij, inb(ij) - 1) = fq(ij, inb(ij) - 1) + frac(ij) * fqold * ((ph(ij, &
	1558	inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, &
	1559	inb(ij)))) * lv(ij, inb(ij)) / lv(ij, inb(ij) - 1)
[3927]	1560	ftold = ft(ij, inb(ij))
[5132]	1561	ft(ij, inb(ij)) = ft(ij, inb(ij)) * (1. - frac(ij))
	1562	ft(ij, inb(ij) - 1) = ft(ij, inb(ij) - 1) + frac(ij) * ftold * ((ph(ij, &
	1563	inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, &
	1564	inb(ij)))) * cpn(ij, inb(ij)) / cpn(ij, inb(ij) - 1)
[3927]	1565	fuold = fu(ij, inb(ij))
[5132]	1566	fu(ij, inb(ij)) = fu(ij, inb(ij)) * (1. - frac(ij))
	1567	fu(ij, inb(ij) - 1) = fu(ij, inb(ij) - 1) + frac(ij) * fuold * ((ph(ij, &
	1568	inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, inb(ij))))
[3927]	1569	fvold = fv(ij, inb(ij))
[5132]	1570	fv(ij, inb(ij)) = fv(ij, inb(ij)) * (1. - frac(ij))
	1571	fv(ij, inb(ij) - 1) = fv(ij, inb(ij) - 1) + frac(ij) * fvold * ((ph(ij, &
	1572	inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, inb(ij))))
[3927]	1573	END DO
	1574
	1575	! * Very slightly adjust tendencies to force exact *
	1576	! * enthalpy, momentum and tracer conservation *
	1577
	1578	DO ij = 1, ncum
	1579	ents(ij) = 0.0
	1580	uav(ij) = 0.0
	1581	vav(ij) = 0.0
	1582	DO i = 1, inb(ij)
[5132]	1583	ents(ij) = ents(ij) + (cpn(ij, i) * ft(ij, i) + lv(ij, i) * fq(ij, i)) * (ph(ij, i) - &
	1584	ph(ij, i + 1))
	1585	uav(ij) = uav(ij) + fu(ij, i) * (ph(ij, i) - ph(ij, i + 1))
	1586	vav(ij) = vav(ij) + fv(ij, i) * (ph(ij, i) - ph(ij, i + 1))
[3927]	1587	END DO
	1588	END DO
	1589	DO ij = 1, ncum
[5132]	1590	ents(ij) = ents(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
	1591	uav(ij) = uav(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
	1592	vav(ij) = vav(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
[3927]	1593	END DO
	1594	DO ij = 1, ncum
	1595	DO i = 1, inb(ij)
[5132]	1596	ft(ij, i) = ft(ij, i) - ents(ij) / cpn(ij, i)
	1597	fu(ij, i) = (1. - cu) * (fu(ij, i) - uav(ij))
	1598	fv(ij, i) = (1. - cu) * (fv(ij, i) - vav(ij))
[3927]	1599	END DO
	1600	END DO
	1601
	1602	DO k = 1, nl + 1
	1603	DO i = 1, ncum
[5132]	1604	IF ((q(i, k) + delt * fq(i, k))<0.0) iflag(i) = 10
[3927]	1605	END DO
	1606	END DO
	1607
	1608	DO i = 1, ncum
	1609	IF (iflag(i)>2) THEN
	1610	precip(i) = 0.0
	1611	cbmf(i) = 0.0
	1612	END IF
	1613	END DO
	1614	DO k = 1, nl
	1615	DO i = 1, ncum
	1616	IF (iflag(i)>2) THEN
	1617	ft(i, k) = 0.0
	1618	fq(i, k) = 0.0
	1619	fu(i, k) = 0.0
	1620	fv(i, k) = 0.0
	1621	qcondc(i, k) = 0.0 ! cld
	1622	END IF
	1623	END DO
	1624	END DO
	1625
	1626	DO k = 1, nl + 1
	1627	DO i = 1, ncum
	1628	ma(i, k) = 0.
	1629	END DO
	1630	END DO
	1631	DO k = nl, 1, -1
	1632	DO i = 1, ncum
[5132]	1633	ma(i, k) = ma(i, k + 1) + m(i, k)
[3927]	1634	END DO
	1635	END DO
	1636
	1637
	1638	! * diagnose the in-cloud mixing ratio * ! cld
	1639	! * of condensed water * ! cld
[5113]	1640	! cld
[3927]	1641	DO ij = 1, ncum ! cld
	1642	DO i = 1, nd ! cld
	1643	mac(ij, i) = 0.0 ! cld
	1644	wa(ij, i) = 0.0 ! cld
	1645	siga(ij, i) = 0.0 ! cld
	1646	END DO ! cld
	1647	DO i = nk(ij), inb(ij) ! cld
	1648	DO k = i + 1, inb(ij) + 1 ! cld
	1649	mac(ij, i) = mac(ij, i) + m(ij, k) ! cld
	1650	END DO ! cld
	1651	END DO ! cld
	1652	DO i = icb(ij), inb(ij) - 1 ! cld
	1653	ax(ij, i) = 0. ! cld
	1654	DO j = icb(ij), i ! cld
[5132]	1655	ax(ij, i) = ax(ij, i) + rrd * (tvp(ij, j) - tv(ij, j)) & ! cld
	1656	* (ph(ij, j) - ph(ij, j + 1)) / p(ij, j) ! cld
[3927]	1657	END DO ! cld
[5132]	1658	IF (ax(ij, i)>0.0) THEN ! cld
	1659	wa(ij, i) = sqrt(2. * ax(ij, i)) ! cld
[3927]	1660	END IF ! cld
	1661	END DO ! cld
	1662	DO i = 1, nl ! cld
[5132]	1663	IF (wa(ij, i)>0.0) & ! cld
	1664	siga(ij, i) = mac(ij, i) / wa(ij, i) & ! cld
	1665	* rrd * tvp(ij, i) / p(ij, i) / 100. / delta ! cld
	1666	siga(ij, i) = min(siga(ij, i), 1.0) ! cld
	1667	qcondc(ij, i) = siga(ij, i) * clw(ij, i) * (1. - ep(ij, i)) & ! cld
	1668	+ (1. - siga(ij, i)) * qcond(ij, i) ! cld
[3927]	1669	END DO ! cld
	1670	END DO ! cld
	1671
	1672	END SUBROUTINE cv_yield
	1673
	1674	SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
[5132]	1675	fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
	1676	qcondc1)
[3927]	1677	IMPLICIT NONE
	1678
[5142]	1679
[3927]	1680	! inputs:
	1681	INTEGER len, ncum, nd, nloc
	1682	INTEGER idcum(nloc)
	1683	INTEGER iflag(nloc)
	1684	REAL precip(nloc), cbmf(nloc)
	1685	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
	1686	REAL ma(nloc, nd)
	1687	REAL qcondc(nloc, nd) !cld
	1688
	1689	! outputs:
	1690	INTEGER iflag1(len)
	1691	REAL precip1(len), cbmf1(len)
	1692	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
	1693	REAL ma1(len, nd)
	1694	REAL qcondc1(len, nd) !cld
	1695
	1696	! local variables:
	1697	INTEGER i, k
	1698
	1699	DO i = 1, ncum
	1700	precip1(idcum(i)) = precip(i)
	1701	cbmf1(idcum(i)) = cbmf(i)
	1702	iflag1(idcum(i)) = iflag(i)
	1703	END DO
	1704
	1705	DO k = 1, nl
	1706	DO i = 1, ncum
	1707	ft1(idcum(i), k) = ft(i, k)
	1708	fq1(idcum(i), k) = fq(i, k)
	1709	fu1(idcum(i), k) = fu(i, k)
	1710	fv1(idcum(i), k) = fv(i, k)
	1711	ma1(idcum(i), k) = ma(i, k)
	1712	qcondc1(idcum(i), k) = qcondc(i, k)
	1713	END DO
	1714	END DO
	1715
	1716	END SUBROUTINE cv_uncompress
	1717
[5142]	1718
	1719	END MODULE lmdz_cv

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

Download in other formats:

Original Format