Context Navigation

lmdz_cv.f90 @ 5157

Last change on this file since 5157 was 5142, checked in by abarral, 3 months ago
Put cvparam.h, fcg_gcssold.h, planete.h, tsoilnudge.h, YOECUMF.h into modules
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

Line
1	! $Id: lmdz_cv.f90 5142 2024-07-29 13:07:34Z abarral $
2
3	MODULE lmdz_cv
4	!------------------------------------------------------------
5	! Parameters for convectL:
6	! (includes - microphysical parameters,
7	! - parameters that control the rate of approach
8	! to quasi-equilibrium)
9	! - noff & minorig (previously in input of convect1)
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
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
68	INTEGER nd
69	CHARACTER (LEN = 20) :: modname = 'cv_routines'
70	CHARACTER (LEN = 80) :: abort_message
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)
105	USE lmdz_cvthermo
106
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
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
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
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
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
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
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, &
160	qnk, gznk, plcl)
161	IMPLICIT NONE
162
163	! ================================================================
164	! Purpose: CONVECTIVE FEED
165	! ================================================================
166
167
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
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
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
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))
244
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	! -------------------------------------------------------------------
255	DO i = 1, len
256	icb(i) = nlm
257	END DO
258
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
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, &
279	clw)
280	USE lmdz_cvthermo
281
282	IMPLICIT NONE
283
284
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
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
321	END DO
322
323	! * Calculate lifted parcel quantities below cloud base *
324
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
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))
337	alv = lv0 - clmcpv * (ticb(i) - t0)
338
339	! First iteration.
340
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))
354
355	! Second iteration.
356
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)
377	END DO
378
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
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
396
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, &
409	icb(i))<=(tv(i, icb(i)) - dtmax))) iflag(i) = 4
410	END DO
411
412	END SUBROUTINE cv_trigger
413
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
421
422
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)
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)
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)
440
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
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
493	dph(i, k) = ph(i, k) - ph(i, k + 1)
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, &
500	gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
501	USE lmdz_cvthermo
502
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
515
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
558	ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
559	t0)) + gznk(i)
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
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)
571
572	! First iteration.
573
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))
587
588	! Second iteration.
589
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))
603
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)
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
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
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
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)
656	END IF
657	END DO
658	END DO
659	DO i = 1, ncum
660	tvp(i, nlp) = tvp(i, nl) - (gz(i, nlp) - gz(i, nl)) / cpd
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
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
701
702	! K Emanuel fix
703
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
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.
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
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)
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
758
759	! =====================================================================
760	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
761	! =====================================================================
762
763	! initialization:
764	DO i = 1, ncum * nlp
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
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
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, &
780	cpn, iflag, cbmf)
781	USE lmdz_cvthermo
782
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)
789	REAL ph(nloc, nd + 1) ! caution nd instead ndp1 to be consistent...
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
801
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
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
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
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
841	DO i = 1, ncum
842	dtpbl(i) = dtpbl(i) / (ph(i, nk(i)) - ph(i, icb(i)))
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)
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
857
858	END SUBROUTINE cv_closure
859
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
864
865	IMPLICIT NONE
866
867
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)
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
896	DO i = 1, ncum * nlp
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
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
930	END DO
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
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
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
993	DO ij = 1, ncum
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)
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
1014	CALL zilch(bsum, ncum * nlp)
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
1022	IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij))) num1 = num1 + 1
1023	END DO
1024	IF (num1<=0) GO TO 789
1025
1026	DO ij = 1, ncum
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
1039	END DO
1040	DO j = minorig, nl
1041
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))
1069	END IF
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))
1074	END IF
1075	END IF
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
1083	END IF
1084	END DO
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
1094	DO ij = 1, ncum
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
1104	END IF
1105	END DO
1106	789 END DO
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, &
1111	ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
1112	USE lmdz_cvthermo
1113
1114	IMPLICIT NONE
1115
1116
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)
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
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
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
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
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
1202	wdtrain(ij) = g * ep(ij, i) * m(ij, i) * clw(ij, i)
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
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
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
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
1248
1249	! * Calculate precipitating downdraft mass flux under *
1250	! * hydrostatic approximation *
1251
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)
1257
1258	! * Add small amount of inertia to downdraft *
1259
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)
1262
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
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
1281	qstm = qs(ij, i - 1)
1282	END IF
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)
1300	END IF
1301	END IF
1302	END DO
1303	899 END DO
1304
1305	END SUBROUTINE cv_unsat
1306
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
1312
1313	IMPLICIT NONE
1314
1315
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)
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
1353	delti = 1.0 / delt
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
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
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.
1380	precip(i) = wt(i, 1) * sigd * water(i, 1) * 86400 / g
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
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
1393
1394	! * Calculate tendencies of lowest level potential temperature *
1395	! * and mixing ratio *
1396
1397	DO i = 1, ncum
1398	work(i) = 0.01 / (ph(i, 1) - ph(i, 2))
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
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)))
1422	END DO
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
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
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
1477	DO ij = 1, ncum
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)))
1495	END IF
1496	END DO
1497	DO k = 1, i - 1
1498	DO ij = 1, ncum
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
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
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	))
1523	END IF
1524	END DO
1525	END DO
1526	DO ij = 1, ncum
1527	IF (i<=inb(ij)) THEN
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
1547	END IF
1548	END DO
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))))
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)
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))
1587	END DO
1588	END DO
1589	DO ij = 1, ncum
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))
1593	END DO
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
1600	END DO
1601
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
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
1633	ma(i, k) = ma(i, k + 1) + m(i, k)
1634	END DO
1635	END DO
1636
1637
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
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
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
1670	END DO ! cld
1671
1672	END SUBROUTINE cv_yield
1673
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
1678
1679
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
1718
1719	END MODULE lmdz_cv

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

Download in other formats:

Original Format