Context Navigation

lmdz_cv.f90 @ 5411

Last change on this file since 5411 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 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: 55.9 KB

Rev	Line
[1403]	1	! $Id: lmdz_cv.f90 5158 2024-08-02 12:12:03Z evignon $
[524]	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	!------------------------------------------------------------
[524]	11
[5142]	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
[524]	17
[5142]	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
[524]	28
[5142]	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)
[524]	31
[5142]	32	CONTAINS
[524]	33
[5142]	34	SUBROUTINE cv_param(nd)
	35	IMPLICIT NONE
[524]	36
[5142]	37	! ------------------------------------------------------------
	38	! Set parameters for convectL
	39	! (includes microphysical parameters and parameters that
	40	! control the rate of approach to quasi-equilibrium)
	41	! ------------------------------------------------------------
[524]	42
[5142]	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) *
[524]	67
[5142]	68	INTEGER nd
	69	CHARACTER (LEN = 20) :: modname = 'cv_routines'
	70	CHARACTER (LEN = 80) :: abort_message
[524]	71
[5142]	72	! noff: integer limit for convection (nd-noff)
	73	! minorig: First level of convection
[524]	74
[5142]	75	noff = 2
	76	minorig = 2
[524]	77
[5142]	78	nl = nd - noff
	79	nlp = nl + 1
	80	nlm = nl - 1
[524]	81
[5142]	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
[524]	92
[5142]	93	cu = 0.70
[5141]	94
[5142]	95	betad = 10.0
[524]	96
[5142]	97	damp = 0.1
	98	alpha = 0.2
[524]	99
[5142]	100	delta = 0.01 ! cld
[524]	101
[5142]	102	END SUBROUTINE cv_param
[524]	103
[5142]	104	SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm)
	105	USE lmdz_cvthermo
[524]	106
[5142]	107	IMPLICIT NONE
[524]	108
[5142]	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
	125
	126	DO k = 1, nlp
	127	DO i = 1, len
	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)
	132	END DO
[1992]	133	END DO
[524]	134
[5142]	135	! gz = phi at the full levels (same as p).
[524]	136
[1992]	137	DO i = 1, len
[5142]	138	gz(i, 1) = 0.0
[1992]	139	END DO
[5142]	140	DO k = 2, nlp
	141	DO i = 1, len
	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)
	144	END DO
	145	END DO
[524]	146
[5142]	147	! h = phi + cpT (dry static energy).
	148	! hm = phi + cp(T-Tbase)+Lq
[524]	149
[5142]	150	DO k = 1, nlp
	151	DO i = 1, len
	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)
	154	END DO
[1992]	155	END DO
[524]	156
[5142]	157	END SUBROUTINE cv_prelim
[524]	158
[5142]	159	SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, &
	160	qnk, gznk, plcl)
	161	IMPLICIT NONE
[524]	162
[5142]	163	! ================================================================
	164	! Purpose: CONVECTIVE FEED
	165	! ================================================================
[524]	166
	167
[5142]	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)
[524]	172
[5142]	173	! outputs:
	174	INTEGER iflag(len), nk(len), icb(len), icbmax
	175	REAL tnk(len), qnk(len), gznk(len), plcl(len)
[524]	176
[5142]	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)
[524]	182
[5142]	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	! -------------------------------------------------------------------
[524]	189
[1992]	190	DO i = 1, len
[5142]	191	work(i) = 1.0E12
	192	ihmin(i) = nl
	193	END DO
	194	DO k = 2, nlp
	195	DO i = 1, len
	196	IF ((hm(i, k)<work(i)) .AND. (hm(i, k)<hm(i, k - 1))) THEN
	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
[1992]	206	END IF
	207	END DO
[524]	208
[5142]	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	! -------------------------------------------------------------------
[524]	213
[1992]	214	DO i = 1, len
[5142]	215	work(i) = hm(i, minorig)
	216	nk(i) = minorig
[1992]	217	END DO
[5142]	218	DO k = minorig + 1, nl
	219	DO i = 1, len
	220	IF ((hm(i, k)>work(i)) .AND. (k<=ihmin(i))) THEN
	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
	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
	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))
[524]	244
[5142]	245	rh(i) = qnk(i) / qsnk(i)
	246	rh(i) = min(1.0, rh(i))
	247	chi(i) = tnk(i) / (1669.0 - 122.0 * rh(i) - tnk(i))
	248	plcl(i) = pnk(i) * (rh(i)**chi(i))
	249	IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
	250	) = 8
	251	END DO
	252	! -------------------------------------------------------------------
	253	! --- Calculate first level above lcl (=icb)
	254	! -------------------------------------------------------------------
[1992]	255	DO i = 1, len
[5142]	256	icb(i) = nlm
[1992]	257	END DO
[524]	258
[5142]	259	DO k = minorig, nl
	260	DO i = 1, len
	261	IF ((k>=(nk(i) + 1)) .AND. (p(i, k)<plcl(i))) icb(i) = min(icb(i), k)
	262	END DO
	263	END DO
[524]	264
[5142]	265	DO i = 1, len
	266	IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
	267	END DO
[524]	268
[5142]	269	! Compute icbmax.
[524]	270
[5142]	271	icbmax = 2
	272	DO i = 1, len
	273	icbmax = max(icbmax, icb(i))
	274	END DO
[524]	275
[5142]	276	END SUBROUTINE cv_feed
[5141]	277
[5142]	278	SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, &
	279	clw)
	280	USE lmdz_cvthermo
[524]	281
[5142]	282	IMPLICIT NONE
[524]	283
	284
[5142]	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)
[1403]	290
[5142]	291	! outputs:
	292	REAL tp(len, nd), tvp(len, nd), clw(len, nd)
[524]	293
[5142]	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)
[524]	299
[5142]	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	! -------------------------------------------------------------------
[524]	306
[5142]	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
[524]	314
[5142]	315	! * Calculate certain parcel quantities, including static energy *
[524]	316
[1992]	317	DO i = 1, len
[5142]	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
[1992]	321	END DO
[524]	322
[5142]	323	! * Calculate lifted parcel quantities below cloud base *
[524]	324
[5142]	325	DO k = minorig, icbmax - 1
	326	DO i = 1, len
	327	tp(i, k) = tnk(i) - (gz(i, k) - gznk(i)) / cpp(i)
	328	tvp(i, k) = tp(i, k) * (1. + qnk(i) * epsi)
	329	END DO
	330	END DO
[524]	331
[5142]	332	! * Find lifted parcel quantities above cloud base *
[524]	333
[5142]	334	DO i = 1, len
	335	tg = ticb(i)
	336	qg = qs(i, icb(i))
	337	alv = lv0 - clmcpv * (ticb(i) - t0)
[524]	338
[5142]	339	! First iteration.
[524]	340
[5142]	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)
	345	tg = max(tg, 35.0)
	346	tc = tg - t0
	347	denom = 243.5 + tc
	348	IF (tc>=0.0) THEN
	349	es = 6.112 * exp(17.67 * tc / denom)
	350	ELSE
	351	es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
	352	END IF
	353	qg = eps * es / (p(i, icb(i)) - es * (1. - eps))
[524]	354
[5142]	355	! Second iteration.
[524]	356
[5142]	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)
	361	tg = max(tg, 35.0)
	362	tc = tg - t0
	363	denom = 243.5 + tc
	364	IF (tc>=0.0) THEN
	365	es = 6.112 * exp(17.67 * tc / denom)
	366	ELSE
	367	es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
	368	END IF
	369	qg = eps * es / (p(i, icb(i)) - es * (1. - eps))
	370
	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
	373	clw(i, icb(i)) = qnk(i) - qg
	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)
[1992]	377	END DO
[524]	378
[5142]	379	DO k = minorig, icbmax
	380	DO i = 1, len
	381	tvp(i, k) = tvp(i, k) - tp(i, k) * qnk(i)
	382	END DO
	383	END DO
[524]	384
[5142]	385	END SUBROUTINE cv_undilute1
[524]	386
[5142]	387	SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
	388	IMPLICIT NONE
[524]	389
[5142]	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	! -------------------------------------------------------------------
[524]	395
	396
[5142]	397	! inputs:
	398	INTEGER len, nd, icb(len)
	399	REAL cbmf(len), tv(len, nd), tvp(len, nd)
[524]	400
[5142]	401	! outputs:
	402	INTEGER iflag(len) ! also an input
[524]	403
[5142]	404	! local variables:
	405	INTEGER i
[524]	406
[5142]	407	DO i = 1, len
	408	IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
	409	icb(i))<=(tv(i, icb(i)) - dtmax))) iflag(i) = 4
	410	END DO
[524]	411
[5142]	412	END SUBROUTINE cv_trigger
[524]	413
[5142]	414	SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
	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)
	418	USE lmdz_print_control, ONLY: lunout
	419	USE lmdz_abort_physic, ONLY: abort_physic
	420	IMPLICIT NONE
[524]	421
	422
[5142]	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)
	429	REAL p1(len, nd), ph1(len, nd + 1), tv1(len, nd), tp1(len, nd)
	430	REAL tvp1(len, nd), clw1(len, nd)
[524]	431
[5142]	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)
	437	REAL p(nloc, nd), ph(nloc, nd + 1), tv(nloc, nd), tp(nloc, nd)
	438	REAL tvp(nloc, nd), clw(nloc, nd)
	439	REAL dph(nloc, nd)
[524]	440
[5142]	441	! local variables:
	442	INTEGER i, k, nn
	443	CHARACTER (LEN = 20) :: modname = 'cv_compress'
	444	CHARACTER (LEN = 80) :: abort_message
	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
[1992]	476	nn = 0
	477	DO i = 1, len
	478	IF (iflag1(i)==0) THEN
	479	nn = nn + 1
[5142]	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)
[1992]	488	END IF
	489	END DO
[524]	490
[5142]	491	DO k = 1, nl
	492	DO i = 1, ncum
	493	dph(i, k) = ph(i, k) - ph(i, k + 1)
	494	END DO
	495	END DO
[524]	496
[5142]	497	END SUBROUTINE cv_compress
[524]	498
[5142]	499	SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
	500	gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
	501	USE lmdz_cvthermo
[524]	502
[5142]	503	IMPLICIT NONE
[524]	504
[5142]	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	! ---------------------------------------------------------------------
[5141]	514
[524]	515
[5142]	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)
[524]	523
[5142]	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)
[524]	529
[5142]	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)
[524]	536
[5142]	537	! =====================================================================
	538	! --- SOME INITIALIZATIONS
	539	! =====================================================================
[524]	540
[5142]	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
[524]	547
[5142]	548	! =====================================================================
	549	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	550	! =====================================================================
[524]	551
[5142]	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
[1992]	557	DO i = 1, ncum
[5142]	558	ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
	559	t0)) + gznk(i)
[1992]	560	END DO
[524]	561
	562
[5142]	563	! * Find lifted parcel quantities above cloud base *
[524]	564
[5142]	565	DO k = minorig + 1, nl
	566	DO i = 1, ncum
	567	IF (k>=(icb(i) + 1)) THEN
	568	tg = t(i, k)
	569	qg = qs(i, k)
	570	alv = lv0 - clmcpv * (t(i, k) - t0)
[524]	571
[5142]	572	! First iteration.
[524]	573
[5142]	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)
	578	tg = max(tg, 35.0)
	579	tc = tg - t0
	580	denom = 243.5 + tc
	581	IF (tc>=0.0) THEN
	582	es = 6.112 * exp(17.67 * tc / denom)
	583	ELSE
	584	es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
	585	END IF
	586	qg = eps * es / (p(i, k) - es * (1. - eps))
[524]	587
[5142]	588	! Second iteration.
[524]	589
[5142]	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)
	594	tg = max(tg, 35.0)
	595	tc = tg - t0
	596	denom = 243.5 + tc
	597	IF (tc>=0.0) THEN
	598	es = 6.112 * exp(17.67 * tc / denom)
	599	ELSE
	600	es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
	601	END IF
	602	qg = eps * es / (p(i, k) - es * (1. - eps))
[524]	603
[5142]	604	alv = lv0 - clmcpv * (t(i, k) - t0)
	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
	608	tp(i, k) = (ah0(i) - (cl - cpd) * qnk(i) * t(i, k) - gz(i, k) - alv * qg) / cpd
	609	! if (.NOT.cpd.gt.1000.) THEN
	610	! PRINT*,'CPD=',cpd
	611	! stop
	612	! END IF
	613	clw(i, k) = qnk(i) - qg
	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)
[1992]	617	END IF
[5142]	618	END DO
	619	END DO
[524]	620
[5142]	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	! =====================================================================
[524]	626
[5142]	627	DO k = minorig + 1, nl
	628	DO i = 1, ncum
	629	IF (k>=(nk(i) + 1)) THEN
	630	tca = tp(i, k) - t0
	631	IF (tca>=0.0) THEN
	632	elacrit = elcrit
	633	ELSE
	634	elacrit = elcrit * (1.0 - tca / tlcrit)
	635	END IF
	636	elacrit = max(elacrit, 0.0)
	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)
	640	sigp(i, k) = sigs
[1992]	641	END IF
[5142]	642	END DO
[1992]	643	END DO
[524]	644
[5142]	645	! =====================================================================
	646	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
	647	! --- VIRTUAL TEMPERATURE
	648	! =====================================================================
[524]	649
[5142]	650	DO k = minorig + 1, nl
	651	DO i = 1, ncum
	652	IF (k>=(icb(i) + 1)) THEN
	653	tvp(i, k) = tvp(i, k) * (1.0 - qnk(i) + ep(i, k) * clw(i, k))
	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)
[1992]	656	END IF
[5142]	657	END DO
[1992]	658	END DO
[5142]	659	DO i = 1, ncum
	660	tvp(i, nlp) = tvp(i, nl) - (gz(i, nlp) - gz(i, nl)) / cpd
	661	END DO
[524]	662
[5142]	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	! =====================================================================
[1992]	668
	669	DO i = 1, ncum
[5142]	670	cape(i) = 0.0
	671	capem(i) = 0.0
	672	inb(i) = icb(i) + 1
	673	inb1(i) = inb(i)
[1992]	674	END DO
	675
[5142]	676	! Originial Code
[1992]	677
[5142]	678	! do 530 k=minorig+1,nl-1
	679	! do 520 i=1,ncum
	680	! IF(k.ge.(icb(i)+1))THEN
	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
	684	! IF(by.ge.0.0)inb1(i)=k+1
	685	! IF(cape(i).gt.0.0)THEN
	686	! inb(i)=k+1
	687	! capem(i)=cape(i)
	688	! END IF
	689	! END IF
	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
[1992]	701
[5142]	702	! K Emanuel fix
[1992]	703
[5142]	704	! CALL zilch(byp,ncum)
	705	! do 530 k=minorig+1,nl-1
	706	! do 520 i=1,ncum
	707	! IF(k.ge.(icb(i)+1))THEN
	708	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	709	! cape(i)=cape(i)+by
	710	! IF(by.ge.0.0)inb1(i)=k+1
	711	! IF(cape(i).gt.0.0)THEN
	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)
	715	! END IF
	716	! END IF
	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
[1992]	728
[5142]	729	! J Teixeira fix
[1992]	730
[5142]	731	CALL zilch(byp, ncum)
[1992]	732	DO i = 1, ncum
[5142]	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.
	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)
	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
[1992]	747	END IF
[5142]	748	END DO
[1992]	749	END DO
[5142]	750	DO i = 1, ncum
	751	cape(i) = capem(i) + byp(i)
	752	defrac = capem(i) - cape(i)
	753	defrac = max(defrac, 0.001)
	754	frac(i) = -cape(i) / defrac
	755	frac(i) = min(frac(i), 1.0)
	756	frac(i) = max(frac(i), 0.0)
	757	END DO
[1992]	758
[5142]	759	! =====================================================================
	760	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
	761	! =====================================================================
[1992]	762
[5142]	763	! initialization:
	764	DO i = 1, ncum * nlp
	765	hp(i, 1) = h(i, 1)
	766	END DO
[1992]	767
[5142]	768	DO k = minorig + 1, nl
	769	DO i = 1, ncum
	770	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
	771	hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k &
	772	)
	773	END IF
	774	END DO
[1992]	775	END DO
	776
[5142]	777	END SUBROUTINE cv_undilute2
[1992]	778
[5142]	779	SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
	780	cpn, iflag, cbmf)
	781	USE lmdz_cvthermo
[5141]	782
[5142]	783	IMPLICIT NONE
[1992]	784
[5142]	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)
	789	REAL ph(nloc, nd + 1) ! caution nd instead ndp1 to be consistent...
	790	REAL plcl(nloc), cpn(nloc, nd)
[1992]	791
[5142]	792	! outputs:
	793	INTEGER iflag(nloc)
	794	REAL cbmf(nloc) ! also an input
[1992]	795
[5142]	796	! local variables:
	797	INTEGER i, k, icbmax
	798	REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
	799	REAL work(nloc)
[1992]	800
	801
[5142]	802	! -------------------------------------------------------------------
	803	! Compute icbmax.
	804	! -------------------------------------------------------------------
[1992]	805
[5142]	806	icbmax = 2
	807	DO i = 1, ncum
	808	icbmax = max(icbmax, icb(i))
	809	END DO
[1992]	810
[5142]	811	! =====================================================================
	812	! --- CALCULATE CLOUD BASE MASS FLUX
	813	! =====================================================================
[1992]	814
[5142]	815	! tvpplcl = parcel temperature lifted adiabatically from level
	816	! icb-1 to the LCL.
	817	! tvaplcl = virtual temperature at the LCL.
[1992]	818
[5142]	819	DO i = 1, ncum
	820	dtpbl(i) = 0.0
	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))
	825	END DO
[1992]	826
[5142]	827	! -------------------------------------------------------------------
	828	! --- Interpolate difference between lifted parcel and
	829	! --- environmental temperatures to lifted condensation level
	830	! -------------------------------------------------------------------
[1992]	831
[5142]	832	! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
[1992]	833
[5142]	834	DO k = minorig, icbmax
	835	DO i = 1, ncum
	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)
	838	END IF
	839	END DO
	840	END DO
[1992]	841	DO i = 1, ncum
[5142]	842	dtpbl(i) = dtpbl(i) / (ph(i, nk(i)) - ph(i, icb(i)))
	843	dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
[1992]	844	END DO
	845
[5142]	846	! -------------------------------------------------------------------
	847	! --- Adjust cloud base mass flux
	848	! -------------------------------------------------------------------
[1992]	849
[5142]	850	DO i = 1, ncum
	851	work(i) = cbmf(i)
	852	cbmf(i) = max(0.0, (1.0 - damp) * cbmf(i) + 0.1 * alpha * dtmin(i))
	853	IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
	854	iflag(i) = 3
	855	END IF
	856	END DO
[1992]	857
[5142]	858	END SUBROUTINE cv_closure
[1992]	859
[5142]	860	SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
	861	h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
	862	sij, elij)
	863	USE lmdz_cvthermo
[5141]	864
[5142]	865	IMPLICIT NONE
[1992]	866
	867
[5142]	868	! inputs:
	869	INTEGER ncum, nd, nloc
	870	INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
	871	REAL cbmf(nloc), qnk(nloc)
	872	REAL ph(nloc, nd + 1)
	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)
[1992]	876
[5142]	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)
[1992]	882
[5142]	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)
[1992]	891
[5142]	892	! =====================================================================
	893	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	894	! =====================================================================
[1992]	895
[5142]	896	DO i = 1, ncum * nlp
	897	nent(i, 1) = 0
	898	m(i, 1) = 0.0
	899	END DO
[1992]	900
[5142]	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
[1992]	911	END DO
	912	END DO
	913
[5142]	914	! -------------------------------------------------------------------
	915	! --- Calculate rates of mixing, m(i)
	916	! -------------------------------------------------------------------
[1992]	917
[5142]	918	CALL zilch(work, ncum)
[1992]	919
[5142]	920	DO j = minorig + 1, nl
	921	DO i = 1, ncum
	922	IF ((j>=(icb(i) + 1)) .AND. (j<=inb(i))) THEN
	923	k = min(j, inb1(i))
	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))
	926	work(i) = work(i) + dbo
	927	m(i, j) = cbmf(i) * dbo
	928	END IF
	929	END DO
[1992]	930	END DO
[5142]	931	DO k = minorig + 1, nl
	932	DO i = 1, ncum
	933	IF ((k>=(icb(i) + 1)) .AND. (k<=inb(i))) THEN
	934	m(i, k) = m(i, k) / work(i)
	935	END IF
	936	END DO
[1992]	937	END DO
	938
	939
[5142]	940	! =====================================================================
	941	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	942	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	943	! --- FRACTION (sij)
	944	! =====================================================================
[1992]	945
[5142]	946	DO i = minorig + 1, nl
	947	DO j = minorig + 1, nl
	948	DO ij = 1, ncum
	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)
	955	dei = denom
	956	IF (abs(dei)<0.01) dei = 0.01
	957	sij(ij, i, j) = anum / dei
	958	sij(ij, i, i) = 1.0
	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))
	962	stemp = sij(ij, i, j)
	963	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
	964	anum = anum - lv(ij, j) * (qti - qs(ij, j) - cwat * bf2)
	965	denom = denom + lv(ij, j) * (q(ij, i) - qti)
	966	IF (abs(denom)<0.01) denom = 0.01
	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
	970	END IF
	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))
	977	elij(ij, i, j) = altem
	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))
	980	nent(ij, i) = nent(ij, i) + 1
	981	END IF
	982	sij(ij, i, j) = max(0.0, sij(ij, i, j))
	983	sij(ij, i, j) = min(1.0, sij(ij, i, j))
	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
[1992]	993	DO ij = 1, ncum
[5142]	994	IF ((i>=(icb(ij) + 1)) .AND. (i<=inb(ij)) .AND. (nent(ij, i)==0)) THEN
	995	ment(ij, i, i) = m(ij, i)
	996	qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
	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)
[1992]	1000	sij(ij, i, i) = 1.0
	1001	END IF
	1002	END DO
	1003	END DO
	1004
[5142]	1005	DO i = 1, ncum
	1006	sij(i, inb(i), inb(i)) = 1.0
[1992]	1007	END DO
	1008
[5142]	1009	! =====================================================================
	1010	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	1011	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	1012	! =====================================================================
[1992]	1013
[5142]	1014	CALL zilch(bsum, ncum * nlp)
[1992]	1015	DO ij = 1, ncum
[5142]	1016	lwork(ij) = .FALSE.
[1992]	1017	END DO
[5142]	1018	DO i = minorig + 1, nl
[1992]	1019
[5142]	1020	num1 = 0
	1021	DO ij = 1, ncum
	1022	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij))) num1 = num1 + 1
	1023	END DO
	1024	IF (num1<=0) GO TO 789
[1992]	1025
	1026	DO ij = 1, ncum
[5142]	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)
	1032	IF (abs(denom)<0.01) denom = 0.01
	1033	scrit(ij) = anum / denom
	1034	alt = qp1 - qs(ij, i) + scrit(ij) * (q(ij, i) - qp1)
	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
[1992]	1039	END DO
[5142]	1040	DO j = minorig, nl
[1992]	1041
[5142]	1042	num2 = 0
	1043	DO ij = 1, ncum
	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
	1046	END DO
	1047	IF (num2<=0) GO TO 783
	1048
	1049	DO ij = 1, ncum
	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
	1053	IF (j>i) THEN
	1054	smid = min(sij(ij, i, j), scrit(ij))
	1055	sjmax = smid
	1056	sjmin = smid
	1057	IF (smid<smin(ij) .AND. sij(ij, i, j + 1)<smid) THEN
	1058	smin(ij) = smid
	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))
	1061	sjmin = min(sjmin, scrit(ij))
	1062	END IF
	1063	ELSE
	1064	sjmax = max(sij(ij, i, j + 1), scrit(ij))
	1065	smid = max(sij(ij, i, j), scrit(ij))
	1066	sjmin = 0.0
	1067	IF (j>1) sjmin = sij(ij, i, j - 1)
	1068	sjmin = max(sjmin, scrit(ij))
[1992]	1069	END IF
[5142]	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))
[1992]	1074	END IF
	1075	END IF
[5142]	1076	END DO
	1077	783 END DO
	1078	DO ij = 1, ncum
	1079	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
	1080	asij(ij) = max(1.0E-21, asij(ij))
	1081	asij(ij) = 1.0 / asij(ij)
	1082	bsum(ij, i) = 0.0
[1992]	1083	END IF
	1084	END DO
[5142]	1085	DO j = minorig, nl + 1
	1086	DO ij = 1, ncum
	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)
	1090	bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
	1091	END IF
	1092	END DO
	1093	END DO
[1992]	1094	DO ij = 1, ncum
[5142]	1095	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
	1096	i)<1.0E-18) .AND. lwork(ij)) THEN
	1097	nent(ij, i) = 0
	1098	ment(ij, i, i) = m(ij, i)
	1099	qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
	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
[1992]	1104	END IF
	1105	END DO
[5142]	1106	789 END DO
[1992]	1107
[5142]	1108	END SUBROUTINE cv_mixing
[1992]	1109
[5142]	1110	SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
	1111	ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
	1112	USE lmdz_cvthermo
[5141]	1113
[5142]	1114	IMPLICIT NONE
[1992]	1115
	1116
[5142]	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)
	1122	REAL p(nloc, nd), ph(nloc, nd + 1), h(nloc, nd)
	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)
[1992]	1125
[5142]	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)
[1992]	1130
[5142]	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)
[1992]	1138
[5142]	1139	! =====================================================================
	1140	! --- PRECIPITATING DOWNDRAFT CALCULATION
	1141	! =====================================================================
[1992]	1142
[5142]	1143	! Initializations:
[1992]	1144
[5142]	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
[1992]	1152	END DO
	1153
	1154	DO i = 1, ncum
[5142]	1155	qp(i, 1) = q(i, 1)
	1156	up(i, 1) = u(i, 1)
	1157	vp(i, 1) = v(i, 1)
[1992]	1158	END DO
	1159
[5142]	1160	DO k = 2, nl + 1
	1161	DO i = 1, ncum
	1162	qp(i, k) = q(i, k - 1)
	1163	up(i, k) = u(i, k - 1)
	1164	vp(i, k) = v(i, k - 1)
	1165	END DO
	1166	END DO
[1992]	1167
	1168
[5142]	1169	! * Check whether ep(inb)=0, if so, skip precipitating *
	1170	! * downdraft calculation *
[1992]	1171
	1172
[5142]	1173	! * Integrate liquid water equation to find condensed water *
	1174	! * and condensed water flux *
[1992]	1175
[5142]	1176	DO i = 1, ncum
	1177	jtt(i) = 2
	1178	IF (ep(i, inb(i))<=0.0001) iflag(i) = 2
	1179	IF (iflag(i)==0) THEN
	1180	lwork(i) = .TRUE.
	1181	ELSE
	1182	lwork(i) = .FALSE.
	1183	END IF
	1184	END DO
[1992]	1185
[5142]	1186	! * Begin downdraft loop *
[1992]	1187
[5142]	1188	CALL zilch(wdtrain, ncum)
	1189	DO i = nl + 1, 1, -1
[1992]	1190
[5142]	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
[1992]	1196
	1197
[5142]	1198	! * Calculate detrained precipitation *
[1992]	1199
[5142]	1200	DO ij = 1, ncum
	1201	IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
	1202	wdtrain(ij) = g * ep(ij, i) * m(ij, i) * clw(ij, i)
	1203	END IF
[1992]	1204	END DO
	1205
[5142]	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
	1210	awat = elij(ij, j, i) - (1. - ep(ij, i)) * clw(ij, i)
	1211	awat = max(0.0, awat)
	1212	wdtrain(ij) = wdtrain(ij) + g * awat * ment(ij, j, i)
	1213	END IF
	1214	END DO
	1215	END DO
	1216	END IF
[1992]	1217
[5142]	1218	! * Find rain water and evaporation using provisional *
	1219	! * estimates of qp(i)and qp(i-1) *
[1992]	1220
	1221
[5142]	1222	! *** Value of terminal velocity and coeffecient of evaporation for snow
	1223	! ***
[1992]	1224
[5142]	1225	DO ij = 1, ncum
	1226	IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
	1227	coeff = coeffs
	1228	wt(ij, i) = omtsnow
[1992]	1229
[5142]	1230	! *** Value of terminal velocity and coeffecient of evaporation for
	1231	! rain ***
[1992]	1232
[5142]	1233	IF (t(ij, i)>273.0) THEN
	1234	coeff = coeffr
	1235	wt(ij, i) = omtrain
	1236	END IF
	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))
	1239	afac = max(afac, 0.0)
	1240	sigt = sigp(ij, i)
	1241	sigt = max(0.0, sigt)
	1242	sigt = min(1.0, sigt)
	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
[1992]	1248
[5142]	1249	! * Calculate precipitating downdraft mass flux under *
	1250	! * hydrostatic approximation *
[1992]	1251
[5142]	1252	IF (i>1) THEN
	1253	dhdp = (h(ij, i) - h(ij, i - 1)) / (p(ij, i - 1) - p(ij, i))
	1254	dhdp = max(dhdp, 10.0)
	1255	mp(ij, i) = 100. * ginv * lv(ij, i) * sigd * evap(ij, i) / dhdp
	1256	mp(ij, i) = max(mp(ij, i), 0.0)
[1992]	1257
[5142]	1258	! * Add small amount of inertia to downdraft *
[1992]	1259
[5142]	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)
[1992]	1262
[5142]	1263	! *** Force mp to decrease linearly to zero
	1264	! ***
	1265	! *** between about 950 mb and the surface
	1266	! ***
	1267
	1268	IF (p(ij, i)>(0.949 * p(ij, 1))) THEN
	1269	jtt(ij) = max(jtt(ij), i)
	1270	mp(ij, i) = mp(ij, jtt(ij)) * (p(ij, 1) - p(ij, i)) / &
	1271	(p(ij, 1) - p(ij, jtt(ij)))
	1272	END IF
[1992]	1273	END IF
	1274
[5142]	1275	! * Find mixing ratio of precipitating downdraft *
[1992]	1276
[5142]	1277	IF (i/=inb(ij)) THEN
	1278	IF (i==1) THEN
	1279	qstm = qs(ij, 1)
	1280	ELSE
	1281	qstm = qs(ij, i - 1)
[1992]	1282	END IF
[5142]	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)
	1289	ELSE
	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)
	1296	END IF
	1297	END IF
	1298	qp(ij, i) = min(qp(ij, i), qstm)
	1299	qp(ij, i) = max(qp(ij, i), 0.0)
[1992]	1300	END IF
	1301	END IF
[5142]	1302	END DO
	1303	899 END DO
[1992]	1304
[5142]	1305	END SUBROUTINE cv_unsat
[1992]	1306
[5142]	1307	SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
	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)
	1311	USE lmdz_cvthermo
[5141]	1312
[5142]	1313	IMPLICIT NONE
[1992]	1314
	1315
[5142]	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)
	1323	REAL p(nloc, nd), ph(nloc, nd + 1), h(nloc, nd)
	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)
[1992]	1332
[5142]	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)
[1992]	1341
[5142]	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
[1992]	1349
	1350
[5142]	1351	! -- initializations:
[1992]	1352
[5142]	1353	delti = 1.0 / delt
[1992]	1354
[5142]	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
	1365	lvcp(i, k) = lv(i, k) / cpn(i, k)
	1366	qcondc(i, k) = 0.0 ! cld
	1367	qcond(i, k) = 0.0 ! cld
	1368	nqcond(i, k) = 0.0 ! cld
	1369	END DO
[1992]	1370	END DO
	1371
	1372
[5142]	1373	! * Calculate surface precipitation in mm/day *
[1992]	1374
[5142]	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.
	1380	precip(i) = wt(i, 1) * sigd * water(i, 1) * 86400 / g
	1381	END IF
	1382	END DO
[1992]	1383
	1384
[5142]	1385	! * Calculate downdraft velocity scale and surface temperature and *
	1386	! * water vapor fluctuations *
[1992]	1387
[5142]	1388	DO i = 1, ncum
	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
	1392	END DO
[1992]	1393
[5142]	1394	! * Calculate tendencies of lowest level potential temperature *
	1395	! * and mixing ratio *
[1992]	1396
	1397	DO i = 1, ncum
[5142]	1398	work(i) = 0.01 / (ph(i, 1) - ph(i, 2))
	1399	am(i) = 0.0
[1992]	1400	END DO
[5142]	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
[1992]	1408	DO i = 1, ncum
[5142]	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)))
[1992]	1422	END DO
[5142]	1423	DO j = 2, nl
	1424	DO i = 1, ncum
	1425	IF (j<=inb(i)) THEN
	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))
	1429	END IF
	1430	END DO
	1431	END DO
[1992]	1432
[5142]	1433	! * Calculate tendencies of potential temperature and mixing ratio *
	1434	! * at levels above the lowest level *
[1992]	1435
[5142]	1436	! * First find the net saturated updraft and downdraft mass fluxes *
	1437	! * through each level *
[1992]	1438
[5142]	1439	DO i = 2, nl + 1
[1992]	1440
[5142]	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
[1992]	1446
[5142]	1447	CALL zilch(amp1, ncum)
	1448	CALL zilch(ad, ncum)
[1992]	1449
[5142]	1450	DO k = i + 1, nl + 1
	1451	DO ij = 1, ncum
	1452	IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij) + 1))) THEN
	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
	1461	IF ((j<=(inb(ij) + 1)) .AND. (i<=inb(ij))) THEN
	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
[1992]	1477	DO ij = 1, ncum
[5142]	1478	IF (i<=inb(ij)) THEN
	1479	dpinv = 0.01 / (ph(ij, i) - ph(ij, i + 1))
	1480	cpinv = 1.0 / cpn(ij, i)
	1481
	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)))
[1992]	1495	END IF
	1496	END DO
[5142]	1497	DO k = 1, i - 1
[1992]	1498	DO ij = 1, ncum
[5142]	1499	IF (i<=inb(ij)) THEN
	1500	awat = elij(ij, k, i) - (1. - ep(ij, i)) * clw(ij, i)
	1501	awat = max(awat, 0.0)
	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	))
	1508	! (saturated updrafts resulting from mixing) ! cld
	1509	qcond(ij, i) = qcond(ij, i) + (elij(ij, k, i) - awat) ! cld
	1510	nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
[1992]	1511	END IF
	1512	END DO
	1513	END DO
[5142]	1514	DO k = i, nl + 1
[1992]	1515	DO ij = 1, ncum
[5142]	1516	IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
	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	))
[1992]	1523	END IF
	1524	END DO
	1525	END DO
	1526	DO ij = 1, ncum
	1527	IF (i<=inb(ij)) THEN
[5142]	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
	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
	1540	IF (nent(ij, i)==0) THEN ! cld
	1541	qcond(ij, i) = qcond(ij, i) + (1. - ep(ij, i)) * clw(ij, i) ! cld
	1542	nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
	1543	END IF ! cld
	1544	IF (nqcond(ij, i)/=0.) THEN ! cld
	1545	qcond(ij, i) = qcond(ij, i) / nqcond(ij, i) ! cld
	1546	END IF ! cld
[1992]	1547	END IF
	1548	END DO
[5142]	1549	1500 END DO
	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))
	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)
	1560	ftold = ft(ij, inb(ij))
	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)
	1565	fuold = fu(ij, inb(ij))
	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))))
	1569	fvold = fv(ij, inb(ij))
	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))))
[1992]	1573	END DO
[5142]	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)
	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))
[1992]	1587	END DO
	1588	END DO
	1589	DO ij = 1, ncum
[5142]	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))
[1992]	1593	END DO
[5142]	1594	DO ij = 1, ncum
	1595	DO i = 1, inb(ij)
	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))
	1599	END DO
[1992]	1600	END DO
	1601
[5142]	1602	DO k = 1, nl + 1
	1603	DO i = 1, ncum
	1604	IF ((q(i, k) + delt * fq(i, k))<0.0) iflag(i) = 10
	1605	END DO
[1992]	1606	END DO
	1607
	1608	DO i = 1, ncum
	1609	IF (iflag(i)>2) THEN
[5142]	1610	precip(i) = 0.0
	1611	cbmf(i) = 0.0
[1992]	1612	END IF
	1613	END DO
[5142]	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
[1992]	1625
[5142]	1626	DO k = 1, nl + 1
	1627	DO i = 1, ncum
	1628	ma(i, k) = 0.
	1629	END DO
[1992]	1630	END DO
[5142]	1631	DO k = nl, 1, -1
	1632	DO i = 1, ncum
	1633	ma(i, k) = ma(i, k + 1) + m(i, k)
	1634	END DO
[1992]	1635	END DO
	1636
	1637
[5142]	1638	! * diagnose the in-cloud mixing ratio * ! cld
	1639	! * of condensed water * ! cld
	1640	! cld
	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
[1992]	1646	END DO ! cld
[5142]	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
[1992]	1651	END DO ! cld
[5142]	1652	DO i = icb(ij), inb(ij) - 1 ! cld
	1653	ax(ij, i) = 0. ! cld
	1654	DO j = icb(ij), i ! cld
	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
	1657	END DO ! cld
	1658	IF (ax(ij, i)>0.0) THEN ! cld
	1659	wa(ij, i) = sqrt(2. * ax(ij, i)) ! cld
	1660	END IF ! cld
	1661	END DO ! cld
	1662	DO i = 1, nl ! cld
	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
	1669	END DO ! cld
[1992]	1670	END DO ! cld
	1671
[5142]	1672	END SUBROUTINE cv_yield
[1992]	1673
[5142]	1674	SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
	1675	fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
	1676	qcondc1)
	1677	IMPLICIT NONE
[1992]	1678
	1679
[5142]	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
[1992]	1688
[5142]	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
[1992]	1695
[5142]	1696	! local variables:
	1697	INTEGER i, k
[1992]	1698
	1699	DO i = 1, ncum
[5142]	1700	precip1(idcum(i)) = precip(i)
	1701	cbmf1(idcum(i)) = cbmf(i)
	1702	iflag1(idcum(i)) = iflag(i)
[1992]	1703	END DO
	1704
[5142]	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
[1992]	1715
[5142]	1716	END SUBROUTINE cv_uncompress
	1717
	1718
	1719	END MODULE lmdz_cv

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

Download in other formats:

Original Format