Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

cv30_routines.F90

Last change on this file was 2520, checked in by lguez, 9 years ago
Bug fix: icb1 was not made >= 2.
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: 98.3 KB

Line
1
2	! $Id: cv30_routines.F90 2520 2016-05-23 08:32:55Z aclsce $
3
4
5
6	SUBROUTINE cv30_param(nd, delt)
7	IMPLICIT NONE
8
9	! ------------------------------------------------------------
10	! Set parameters for convectL for iflag_con = 3
11	! ------------------------------------------------------------
12
13
14	! * PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *
15	! * PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *
16	! * PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *
17	! * EFFICIENCY IS ASSUMED TO BE UNITY *
18	! * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
19	! * SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
20	! * OF CLOUD *
21
22	! [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
23	! * ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *
24	! * APPROACH TO QUASI-EQUILIBRIUM *
25	! * (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *
26	! * (BETA MUST BE LESS THAN OR EQUAL TO 1) *
27
28	! * DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *
29	! * APPROACH TO QUASI-EQUILIBRIUM *
30	! * IT MUST BE LESS THAN 0 *
31
32	include "cv30param.h"
33	include "conema3.h"
34
35	INTEGER nd
36	REAL delt ! timestep (seconds)
37
38	! noff: integer limit for convection (nd-noff)
39	! minorig: First level of convection
40
41	! -- limit levels for convection:
42
43	noff = 1
44	minorig = 1
45	nl = nd - noff
46	nlp = nl + 1
47	nlm = nl - 1
48
49	! -- "microphysical" parameters:
50
51	sigd = 0.01
52	spfac = 0.15
53	pbcrit = 150.0
54	ptcrit = 500.0
55	! IM cf. FH epmax = 0.993
56
57	omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
58
59	! -- misc:
60
61	dtovsh = -0.2 ! dT for overshoot
62	dpbase = -40. ! definition cloud base (400m above LCL)
63	dttrig = 5. ! (loose) condition for triggering
64
65	! -- rate of approach to quasi-equilibrium:
66
67	dtcrit = -2.0
68	tau = 8000.
69	beta = 1.0 - delt/tau
70	alpha = 1.5E-3*delt/tau
71	! increase alpha to compensate W decrease:
72	alpha = alpha*1.5
73
74	! -- interface cloud parameterization:
75
76	delta = 0.01 ! cld
77
78	! -- interface with boundary-layer (gust factor): (sb)
79
80	betad = 10.0 ! original value (from convect 4.3)
81
82	RETURN
83	END SUBROUTINE cv30_param
84
85	SUBROUTINE cv30_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm, &
86	th)
87	IMPLICIT NONE
88
89	! =====================================================================
90	! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
91	! "ori": from convect4.3 (vectorized)
92	! "convect3": to be exactly consistent with convect3
93	! =====================================================================
94
95	! inputs:
96	INTEGER len, nd, ndp1
97	REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
98
99	! outputs:
100	REAL lv(len, nd), cpn(len, nd), tv(len, nd)
101	REAL gz(len, nd), h(len, nd), hm(len, nd)
102	REAL th(len, nd)
103
104	! local variables:
105	INTEGER k, i
106	REAL rdcp
107	REAL tvx, tvy ! convect3
108	REAL cpx(len, nd)
109
110	include "cvthermo.h"
111	include "cv30param.h"
112
113
114	! ori do 110 k=1,nlp
115	DO k = 1, nl ! convect3
116	DO i = 1, len
117	! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
118	lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
119	cpn(i, k) = cpd(1.0-q(i,k)) + cpvq(i, k)
120	cpx(i, k) = cpd(1.0-q(i,k)) + clq(i, k)
121	! ori tv(i,k)=t(i,k)(1.0+q(i,k)epsim1)
122	tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
123	rdcp = (rrd(1.-q(i,k))+q(i,k)rrv)/cpn(i, k)
124	th(i, k) = t(i, k)(1000.0/p(i,k))*rdcp
125	END DO
126	END DO
127
128	! gz = phi at the full levels (same as p).
129
130	DO i = 1, len
131	gz(i, 1) = 0.0
132	END DO
133	! ori do 140 k=2,nlp
134	DO k = 2, nl ! convect3
135	DO i = 1, len
136	tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3
137	tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
138	gz(i, k) = gz(i, k-1) + 0.5rrd(tvx+tvy) & !convect3
139	*(p(i,k-1)-p(i,k))/ph(i, k) !convect3
140
141	! ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
142	! ori & *(p(i,k-1)-p(i,k))/ph(i,k)
143	END DO
144	END DO
145
146	! h = phi + cpT (dry static energy).
147	! hm = phi + cp(T-Tbase)+Lq
148
149	! ori do 170 k=1,nlp
150	DO k = 1, nl ! convect3
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	RETURN
158	END SUBROUTINE cv30_prelim
159
160	SUBROUTINE cv30_feed(len, nd, t, q, qs, p, ph, hm, gz, nk, icb, icbmax, &
161	iflag, tnk, qnk, gznk, plcl)
162	IMPLICIT NONE
163
164	! ================================================================
165	! Purpose: CONVECTIVE FEED
166
167	! Main differences with cv_feed:
168	! - ph added in input
169	! - here, nk(i)=minorig
170	! - icb defined differently (plcl compared with ph instead of p)
171
172	! Main differences with convect3:
173	! - we do not compute dplcldt and dplcldr of CLIFT anymore
174	! - values iflag different (but tests identical)
175	! - A,B explicitely defined (!...)
176	! ================================================================
177
178	include "cv30param.h"
179
180	! inputs:
181	INTEGER len, nd
182	REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
183	REAL hm(len, nd), gz(len, nd)
184	REAL ph(len, nd+1)
185
186	! outputs:
187	INTEGER iflag(len), nk(len), icb(len), icbmax
188	REAL tnk(len), qnk(len), gznk(len), plcl(len)
189
190	! local variables:
191	INTEGER i, k
192	INTEGER ihmin(len)
193	REAL work(len)
194	REAL pnk(len), qsnk(len), rh(len), chi(len)
195	REAL a, b ! convect3
196	! ym
197	plcl = 0.0
198	! @ !-------------------------------------------------------------------
199	! @ ! --- Find level of minimum moist static energy
200	! @ ! --- If level of minimum moist static energy coincides with
201	! @ ! --- or is lower than minimum allowable parcel origin level,
202	! @ ! --- set iflag to 6.
203	! @ !-------------------------------------------------------------------
204	! @
205	! @ do 180 i=1,len
206	! @ work(i)=1.0e12
207	! @ ihmin(i)=nl
208	! @ 180 continue
209	! @ do 200 k=2,nlp
210	! @ do 190 i=1,len
211	! @ if((hm(i,k).lt.work(i)).and.
212	! @ & (hm(i,k).lt.hm(i,k-1)))then
213	! @ work(i)=hm(i,k)
214	! @ ihmin(i)=k
215	! @ endif
216	! @ 190 continue
217	! @ 200 continue
218	! @ do 210 i=1,len
219	! @ ihmin(i)=min(ihmin(i),nlm)
220	! @ if(ihmin(i).le.minorig)then
221	! @ iflag(i)=6
222	! @ endif
223	! @ 210 continue
224	! @ c
225	! @ !-------------------------------------------------------------------
226	! @ ! --- Find that model level below the level of minimum moist static
227	! @ ! --- energy that has the maximum value of moist static energy
228	! @ !-------------------------------------------------------------------
229	! @
230	! @ do 220 i=1,len
231	! @ work(i)=hm(i,minorig)
232	! @ nk(i)=minorig
233	! @ 220 continue
234	! @ do 240 k=minorig+1,nl
235	! @ do 230 i=1,len
236	! @ if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
237	! @ work(i)=hm(i,k)
238	! @ nk(i)=k
239	! @ endif
240	! @ 230 continue
241	! @ 240 continue
242
243	! -------------------------------------------------------------------
244	! --- Origin level of ascending parcels for convect3:
245	! -------------------------------------------------------------------
246
247	DO i = 1, len
248	nk(i) = minorig
249	END DO
250
251	! -------------------------------------------------------------------
252	! --- Check whether parcel level temperature and specific humidity
253	! --- are reasonable
254	! -------------------------------------------------------------------
255	DO i = 1, len
256	IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0)) & ! @ & .or.(
257	! p(i,ihmin(i)).lt.400.0
258	! ) )
259	.AND. (iflag(i)==0)) iflag(i) = 7
260	END DO
261	! -------------------------------------------------------------------
262	! --- Calculate lifted condensation level of air at parcel origin level
263	! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
264	! -------------------------------------------------------------------
265
266	a = 1669.0 ! convect3
267	b = 122.0 ! convect3
268
269	DO i = 1, len
270
271	IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002
272
273	tnk(i) = t(i, nk(i))
274	qnk(i) = q(i, nk(i))
275	gznk(i) = gz(i, nk(i))
276	pnk(i) = p(i, nk(i))
277	qsnk(i) = qs(i, nk(i))
278
279	rh(i) = qnk(i)/qsnk(i)
280	! ori rh(i)=min(1.0,rh(i)) ! removed for convect3
281	! ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
282	chi(i) = tnk(i)/(a-b*rh(i)-tnk(i)) ! convect3
283	plcl(i) = pnk(i)(rh(i)*chi(i))
284	IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag &
285	(i) = 8
286
287	END IF ! iflag=7
288
289	END DO
290
291	! -------------------------------------------------------------------
292	! --- Calculate first level above lcl (=icb)
293	! -------------------------------------------------------------------
294
295	! @ do 270 i=1,len
296	! @ icb(i)=nlm
297	! @ 270 continue
298	! @c
299	! @ do 290 k=minorig,nl
300	! @ do 280 i=1,len
301	! @ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
302	! @ & icb(i)=min(icb(i),k)
303	! @ 280 continue
304	! @ 290 continue
305	! @c
306	! @ do 300 i=1,len
307	! @ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
308	! @ 300 continue
309
310	DO i = 1, len
311	icb(i) = nlm
312	END DO
313
314	! la modification consiste a comparer plcl a ph et non a p:
315	! icb est defini par : ph(icb)<plcl<ph(icb-1)
316	! @ do 290 k=minorig,nl
317	DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
318	DO i = 1, len
319	IF (ph(i,k)<plcl(i)) icb(i) = min(icb(i), k)
320	END DO
321	END DO
322
323	DO i = 1, len
324	! @ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
325	IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
326	END DO
327
328	DO i = 1, len
329	icb(i) = icb(i) - 1 ! icb sup ou egal a 2
330	END DO
331
332	! Compute icbmax.
333
334	icbmax = 2
335	DO i = 1, len
336	! ! icbmax=max(icbmax,icb(i))
337	IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
338	END DO
339
340	RETURN
341	END SUBROUTINE cv30_feed
342
343	SUBROUTINE cv30_undilute1(len, nd, t, q, qs, gz, plcl, p, nk, icb, tp, tvp, &
344	clw, icbs)
345	IMPLICIT NONE
346
347	! ----------------------------------------------------------------
348	! Equivalent de TLIFT entre NK et ICB+1 inclus
349
350	! Differences with convect4:
351	! - specify plcl in input
352	! - icbs is the first level above LCL (may differ from icb)
353	! - in the iterations, used x(icbs) instead x(icb)
354	! - many minor differences in the iterations
355	! - tvp is computed in only one time
356	! - icbs: first level above Plcl (IMIN de TLIFT) in output
357	! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
358	! ----------------------------------------------------------------
359
360	include "cvthermo.h"
361	include "cv30param.h"
362
363	! inputs:
364	INTEGER len, nd
365	INTEGER nk(len), icb(len)
366	REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
367	REAL p(len, nd)
368	REAL plcl(len) ! convect3
369
370	! outputs:
371	REAL tp(len, nd), tvp(len, nd), clw(len, nd)
372
373	! local variables:
374	INTEGER i, k
375	INTEGER icb1(len), icbs(len), icbsmax2 ! convect3
376	REAL tg, qg, alv, s, ahg, tc, denom, es, rg
377	REAL ah0(len), cpp(len)
378	REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
379	REAL qsicb(len) ! convect3
380	REAL cpinv(len) ! convect3
381
382	! -------------------------------------------------------------------
383	! --- Calculates the lifted parcel virtual temperature at nk,
384	! --- the actual temperature, and the adiabatic
385	! --- liquid water content. The procedure is to solve the equation.
386	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
387	! -------------------------------------------------------------------
388
389	DO i = 1, len
390	tnk(i) = t(i, nk(i))
391	qnk(i) = q(i, nk(i))
392	gznk(i) = gz(i, nk(i))
393	! ori ticb(i)=t(i,icb(i))
394	! ori gzicb(i)=gz(i,icb(i))
395	END DO
396
397	! * Calculate certain parcel quantities, including static energy *
398
399	DO i = 1, len
400	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))tnk(i) + qnk(i)(lv0-clmcpv*(tnk(i)- &
401	273.15)) + gznk(i)
402	cpp(i) = cpd(1.-qnk(i)) + qnk(i)cpv
403	cpinv(i) = 1./cpp(i)
404	END DO
405
406	! * Calculate lifted parcel quantities below cloud base *
407
408	DO i = 1, len !convect3
409	icb1(i) = min(max(icb(i), 2), nl)
410	! if icb is below LCL, start loop at ICB+1:
411	! (icbs est le premier niveau au-dessus du LCL)
412	icbs(i) = icb1(i) !convect3
413	IF (plcl(i)<p(i,icb1(i))) THEN
414	icbs(i) = min(icbs(i)+1, nl) !convect3
415	END IF
416	END DO !convect3
417
418	DO i = 1, len !convect3
419	ticb(i) = t(i, icbs(i)) !convect3
420	gzicb(i) = gz(i, icbs(i)) !convect3
421	qsicb(i) = qs(i, icbs(i)) !convect3
422	END DO !convect3
423
424
425	! Re-compute icbsmax (icbsmax2): !convect3
426	! !convect3
427	icbsmax2 = 2 !convect3
428	DO i = 1, len !convect3
429	icbsmax2 = max(icbsmax2, icbs(i)) !convect3
430	END DO !convect3
431
432	! initialization outputs:
433
434	DO k = 1, icbsmax2 ! convect3
435	DO i = 1, len ! convect3
436	tp(i, k) = 0.0 ! convect3
437	tvp(i, k) = 0.0 ! convect3
438	clw(i, k) = 0.0 ! convect3
439	END DO ! convect3
440	END DO ! convect3
441
442	! tp and tvp below cloud base:
443
444	DO k = minorig, icbsmax2 - 1
445	DO i = 1, len
446	tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
447	tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
448	END DO
449	END DO
450
451	! * Find lifted parcel quantities above cloud base *
452
453	DO i = 1, len
454	tg = ticb(i)
455	! ori qg=qs(i,icb(i))
456	qg = qsicb(i) ! convect3
457	! debug alv=lv0-clmcpv*(ticb(i)-t0)
458	alv = lv0 - clmcpv*(ticb(i)-273.15)
459
460	! First iteration.
461
462	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
463	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
464	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
465	s = 1./s
466	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
467	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
468	tg = tg + s*(ah0(i)-ahg)
469	! ori tg=max(tg,35.0)
470	! debug tc=tg-t0
471	tc = tg - 273.15
472	denom = 243.5 + tc
473	denom = max(denom, 1.0) ! convect3
474	! ori if(tc.ge.0.0)then
475	es = 6.112exp(17.67tc/denom)
476	! ori else
477	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
478	! ori endif
479	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
480	qg = epses/(p(i,icbs(i))-es(1.-eps))
481
482	! Second iteration.
483
484
485	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
486	! ori s=1./s
487	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
488	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
489	tg = tg + s*(ah0(i)-ahg)
490	! ori tg=max(tg,35.0)
491	! debug tc=tg-t0
492	tc = tg - 273.15
493	denom = 243.5 + tc
494	denom = max(denom, 1.0) ! convect3
495	! ori if(tc.ge.0.0)then
496	es = 6.112exp(17.67tc/denom)
497	! ori else
498	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
499	! ori end if
500	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
501	qg = epses/(p(i,icbs(i))-es(1.-eps))
502
503	alv = lv0 - clmcpv*(ticb(i)-273.15)
504
505	! ori c approximation here:
506	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
507	! ori & -gz(i,icb(i))-alv*qg)/cpd
508
509	! convect3: no approximation:
510	tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alvqg)/(cpd+(cl-cpd)qnk(i))
511
512	! ori clw(i,icb(i))=qnk(i)-qg
513	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
514	clw(i, icbs(i)) = qnk(i) - qg
515	clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
516
517	rg = qg/(1.-qnk(i))
518	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
519	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
520	tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
521
522	END DO
523
524	! ori do 380 k=minorig,icbsmax2
525	! ori do 370 i=1,len
526	! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
527	! ori 370 continue
528	! ori 380 continue
529
530
531	! -- The following is only for convect3:
532
533	! * icbs is the first level above the LCL:
534	! if plcl<p(icb), then icbs=icb+1
535	! if plcl>p(icb), then icbs=icb
536
537	! * the routine above computes tvp from minorig to icbs (included).
538
539	! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
540	! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
541
542	! * therefore, in the case icbs=icb, we compute tvp at level icb+1
543	! (tvp at other levels will be computed in cv3_undilute2.F)
544
545
546	DO i = 1, len
547	ticb(i) = t(i, icb(i)+1)
548	gzicb(i) = gz(i, icb(i)+1)
549	qsicb(i) = qs(i, icb(i)+1)
550	END DO
551
552	DO i = 1, len
553	tg = ticb(i)
554	qg = qsicb(i) ! convect3
555	! debug alv=lv0-clmcpv*(ticb(i)-t0)
556	alv = lv0 - clmcpv*(ticb(i)-273.15)
557
558	! First iteration.
559
560	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
561	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
562	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
563	s = 1./s
564	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
565	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
566	tg = tg + s*(ah0(i)-ahg)
567	! ori tg=max(tg,35.0)
568	! debug tc=tg-t0
569	tc = tg - 273.15
570	denom = 243.5 + tc
571	denom = max(denom, 1.0) ! convect3
572	! ori if(tc.ge.0.0)then
573	es = 6.112exp(17.67tc/denom)
574	! ori else
575	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
576	! ori endif
577	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
578	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
579
580	! Second iteration.
581
582
583	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
584	! ori s=1./s
585	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
586	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
587	tg = tg + s*(ah0(i)-ahg)
588	! ori tg=max(tg,35.0)
589	! debug tc=tg-t0
590	tc = tg - 273.15
591	denom = 243.5 + tc
592	denom = max(denom, 1.0) ! convect3
593	! ori if(tc.ge.0.0)then
594	es = 6.112exp(17.67tc/denom)
595	! ori else
596	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
597	! ori end if
598	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
599	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
600
601	alv = lv0 - clmcpv*(ticb(i)-273.15)
602
603	! ori c approximation here:
604	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
605	! ori & -gz(i,icb(i))-alv*qg)/cpd
606
607	! convect3: no approximation:
608	tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alvqg)/(cpd+(cl-cpd)qnk(i))
609
610	! ori clw(i,icb(i))=qnk(i)-qg
611	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
612	clw(i, icb(i)+1) = qnk(i) - qg
613	clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
614
615	rg = qg/(1.-qnk(i))
616	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
617	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
618	tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
619
620	END DO
621
622	RETURN
623	END SUBROUTINE cv30_undilute1
624
625	SUBROUTINE cv30_trigger(len, nd, icb, plcl, p, th, tv, tvp, pbase, buoybase, &
626	iflag, sig, w0)
627	IMPLICIT NONE
628
629	! -------------------------------------------------------------------
630	! --- TRIGGERING
631
632	! - computes the cloud base
633	! - triggering (crude in this version)
634	! - relaxation of sig and w0 when no convection
635
636	! Caution1: if no convection, we set iflag=4
637	! (it used to be 0 in convect3)
638
639	! Caution2: at this stage, tvp (and thus buoy) are know up
640	! through icb only!
641	! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
642	! -------------------------------------------------------------------
643
644	include "cv30param.h"
645
646	! input:
647	INTEGER len, nd
648	INTEGER icb(len)
649	REAL plcl(len), p(len, nd)
650	REAL th(len, nd), tv(len, nd), tvp(len, nd)
651
652	! output:
653	REAL pbase(len), buoybase(len)
654
655	! input AND output:
656	INTEGER iflag(len)
657	REAL sig(len, nd), w0(len, nd)
658
659	! local variables:
660	INTEGER i, k
661	REAL tvpbase, tvbase, tdif, ath, ath1
662
663
664	! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
665
666	DO i = 1, len
667	pbase(i) = plcl(i) + dpbase
668	tvpbase = tvp(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
669	(p(i,icb(i))-p(i,icb(i)+1)) + tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/( &
670	p(i,icb(i))-p(i,icb(i)+1))
671	tvbase = tv(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
672	(p(i,icb(i))-p(i,icb(i)+1)) + tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/(p &
673	(i,icb(i))-p(i,icb(i)+1))
674	buoybase(i) = tvpbase - tvbase
675	END DO
676
677
678	! * make sure that column is dry adiabatic between the surface *
679	! * and cloud base, and that lifted air is positively buoyant *
680	! * at cloud base *
681	! * if not, return to calling program after resetting *
682	! * sig(i) and w0(i) *
683
684
685	! oct3 do 200 i=1,len
686	! oct3
687	! oct3 tdif = buoybase(i)
688	! oct3 ath1 = th(i,1)
689	! oct3 ath = th(i,icb(i)-1) - dttrig
690	! oct3
691	! oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
692	! oct3 do 60 k=1,nl
693	! oct3 sig(i,k) = betasig(i,k) - 2.alphatdiftdif
694	! oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
695	! oct3 w0(i,k) = beta*w0(i,k)
696	! oct3 60 continue
697	! oct3 iflag(i)=4 ! pour version vectorisee
698	! oct3c convect3 iflag(i)=0
699	! oct3cccc return
700	! oct3 endif
701	! oct3
702	! oct3200 continue
703
704	! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
705
706	DO k = 1, nl
707	DO i = 1, len
708
709	tdif = buoybase(i)
710	ath1 = th(i, 1)
711	ath = th(i, icb(i)-1) - dttrig
712
713	IF (tdif<dtcrit .OR. ath>ath1) THEN
714	sig(i, k) = betasig(i, k) - 2.alphatdiftdif
715	sig(i, k) = amax1(sig(i,k), 0.0)
716	w0(i, k) = beta*w0(i, k)
717	iflag(i) = 4 ! pour version vectorisee
718	! convect3 iflag(i)=0
719	END IF
720
721	END DO
722	END DO
723
724	! fin oct3 --
725
726	RETURN
727	END SUBROUTINE cv30_trigger
728
729	SUBROUTINE cv30_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, &
730	plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, &
731	th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, &
732	iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, &
733	v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0)
734	USE print_control_mod, ONLY: lunout
735	IMPLICIT NONE
736
737	include "cv30param.h"
738
739	! inputs:
740	INTEGER len, ncum, nd, ntra, nloc
741	INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
742	REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
743	REAL pbase1(len), buoybase1(len)
744	REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
745	REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
746	REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
747	REAL tvp1(len, nd), clw1(len, nd)
748	REAL th1(len, nd)
749	REAL sig1(len, nd), w01(len, nd)
750	REAL tra1(len, nd, ntra)
751
752	! outputs:
753	! en fait, on a nloc=len pour l'instant (cf cv_driver)
754	INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
755	REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
756	REAL pbase(nloc), buoybase(nloc)
757	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
758	REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
759	REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
760	REAL tvp(nloc, nd), clw(nloc, nd)
761	REAL th(nloc, nd)
762	REAL sig(nloc, nd), w0(nloc, nd)
763	REAL tra(nloc, nd, ntra)
764
765	! local variables:
766	INTEGER i, k, nn, j
767
768	CHARACTER (LEN=20) :: modname = 'cv30_compress'
769	CHARACTER (LEN=80) :: abort_message
770
771
772	DO k = 1, nl + 1
773	nn = 0
774	DO i = 1, len
775	IF (iflag1(i)==0) THEN
776	nn = nn + 1
777	sig(nn, k) = sig1(i, k)
778	w0(nn, k) = w01(i, k)
779	t(nn, k) = t1(i, k)
780	q(nn, k) = q1(i, k)
781	qs(nn, k) = qs1(i, k)
782	u(nn, k) = u1(i, k)
783	v(nn, k) = v1(i, k)
784	gz(nn, k) = gz1(i, k)
785	h(nn, k) = h1(i, k)
786	lv(nn, k) = lv1(i, k)
787	cpn(nn, k) = cpn1(i, k)
788	p(nn, k) = p1(i, k)
789	ph(nn, k) = ph1(i, k)
790	tv(nn, k) = tv1(i, k)
791	tp(nn, k) = tp1(i, k)
792	tvp(nn, k) = tvp1(i, k)
793	clw(nn, k) = clw1(i, k)
794	th(nn, k) = th1(i, k)
795	END IF
796	END DO
797	END DO
798
799	! do 121 j=1,ntra
800	! do 111 k=1,nd
801	! nn=0
802	! do 101 i=1,len
803	! if(iflag1(i).eq.0)then
804	! nn=nn+1
805	! tra(nn,k,j)=tra1(i,k,j)
806	! endif
807	! 101 continue
808	! 111 continue
809	! 121 continue
810
811	IF (nn/=ncum) THEN
812	WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
813	abort_message = ''
814	CALL abort_physic(modname, abort_message, 1)
815	END IF
816
817	nn = 0
818	DO i = 1, len
819	IF (iflag1(i)==0) THEN
820	nn = nn + 1
821	pbase(nn) = pbase1(i)
822	buoybase(nn) = buoybase1(i)
823	plcl(nn) = plcl1(i)
824	tnk(nn) = tnk1(i)
825	qnk(nn) = qnk1(i)
826	gznk(nn) = gznk1(i)
827	nk(nn) = nk1(i)
828	icb(nn) = icb1(i)
829	icbs(nn) = icbs1(i)
830	iflag(nn) = iflag1(i)
831	END IF
832	END DO
833
834	RETURN
835	END SUBROUTINE cv30_compress
836
837	SUBROUTINE cv30_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
838	q, qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
839	ep, sigp, buoy)
840	! epmax_cape: ajout arguments
841	IMPLICIT NONE
842
843	! ---------------------------------------------------------------------
844	! Purpose:
845	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
846	! &
847	! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
848	! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
849	! &
850	! FIND THE LEVEL OF NEUTRAL BUOYANCY
851
852	! Main differences convect3/convect4:
853	! - icbs (input) is the first level above LCL (may differ from icb)
854	! - many minor differences in the iterations
855	! - condensed water not removed from tvp in convect3
856	! - vertical profile of buoyancy computed here (use of buoybase)
857	! - the determination of inb is different
858	! - no inb1, only inb in output
859	! ---------------------------------------------------------------------
860
861	include "cvthermo.h"
862	include "cv30param.h"
863	include "conema3.h"
864
865	! inputs:
866	INTEGER ncum, nd, nloc
867	INTEGER icb(nloc), icbs(nloc), nk(nloc)
868	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
869	REAL p(nloc, nd)
870	REAL tnk(nloc), qnk(nloc), gznk(nloc)
871	REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
872	REAL pbase(nloc), buoybase(nloc), plcl(nloc)
873
874	! outputs:
875	INTEGER inb(nloc)
876	REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
877	REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
878	REAL buoy(nloc, nd)
879
880	! local variables:
881	INTEGER i, k
882	REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
883	REAL by, defrac, pden
884	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
885	LOGICAL lcape(nloc)
886
887	! =====================================================================
888	! --- SOME INITIALIZATIONS
889	! =====================================================================
890
891	DO k = 1, nl
892	DO i = 1, ncum
893	ep(i, k) = 0.0
894	sigp(i, k) = spfac
895	END DO
896	END DO
897
898	! =====================================================================
899	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
900	! =====================================================================
901
902	! --- The procedure is to solve the equation.
903	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
904
905	! * Calculate certain parcel quantities, including static energy *
906
907
908	DO i = 1, ncum
909	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))*tnk(i) & ! debug &
910	! +qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
911	+qnk(i)(lv0-clmcpv(tnk(i)-273.15)) + gznk(i)
912	END DO
913
914
915	! * Find lifted parcel quantities above cloud base *
916
917
918	DO k = minorig + 1, nl
919	DO i = 1, ncum
920	! ori if(k.ge.(icb(i)+1))then
921	IF (k>=(icbs(i)+1)) THEN ! convect3
922	tg = t(i, k)
923	qg = qs(i, k)
924	! debug alv=lv0-clmcpv*(t(i,k)-t0)
925	alv = lv0 - clmcpv*(t(i,k)-273.15)
926
927	! First iteration.
928
929	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
930	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
931	+alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
932	s = 1./s
933	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
934	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
935	tg = tg + s*(ah0(i)-ahg)
936	! ori tg=max(tg,35.0)
937	! debug tc=tg-t0
938	tc = tg - 273.15
939	denom = 243.5 + tc
940	denom = max(denom, 1.0) ! convect3
941	! ori if(tc.ge.0.0)then
942	es = 6.112exp(17.67tc/denom)
943	! ori else
944	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
945	! ori endif
946	qg = epses/(p(i,k)-es(1.-eps))
947
948	! Second iteration.
949
950	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
951	! ori s=1./s
952	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
953	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
954	tg = tg + s*(ah0(i)-ahg)
955	! ori tg=max(tg,35.0)
956	! debug tc=tg-t0
957	tc = tg - 273.15
958	denom = 243.5 + tc
959	denom = max(denom, 1.0) ! convect3
960	! ori if(tc.ge.0.0)then
961	es = 6.112exp(17.67tc/denom)
962	! ori else
963	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
964	! ori endif
965	qg = epses/(p(i,k)-es(1.-eps))
966
967	! debug alv=lv0-clmcpv*(t(i,k)-t0)
968	alv = lv0 - clmcpv*(t(i,k)-273.15)
969	! print*,'cpd dans convect2 ',cpd
970	! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
971	! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
972
973	! ori c approximation here:
974	! ori
975	! tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
976
977	! convect3: no approximation:
978	tp(i, k) = (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
979
980	clw(i, k) = qnk(i) - qg
981	clw(i, k) = max(0.0, clw(i,k))
982	rg = qg/(1.-qnk(i))
983	! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
984	! convect3: (qg utilise au lieu du vrai mixing ratio rg):
985	tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
986	END IF
987	END DO
988	END DO
989
990	! =====================================================================
991	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
992	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
993	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
994	! =====================================================================
995
996	! ori do 320 k=minorig+1,nl
997	DO k = 1, nl ! convect3
998	DO i = 1, ncum
999	pden = ptcrit - pbcrit
1000	ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1001	ep(i, k) = amax1(ep(i,k), 0.0)
1002	ep(i, k) = amin1(ep(i,k), epmax)
1003	sigp(i, k) = spfac
1004	! ori if(k.ge.(nk(i)+1))then
1005	! ori tca=tp(i,k)-t0
1006	! ori if(tca.ge.0.0)then
1007	! ori elacrit=elcrit
1008	! ori else
1009	! ori elacrit=elcrit*(1.0-tca/tlcrit)
1010	! ori endif
1011	! ori elacrit=max(elacrit,0.0)
1012	! ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
1013	! ori ep(i,k)=max(ep(i,k),0.0 )
1014	! ori ep(i,k)=min(ep(i,k),1.0 )
1015	! ori sigp(i,k)=sigs
1016	! ori endif
1017	END DO
1018	END DO
1019
1020	! =====================================================================
1021	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1022	! --- VIRTUAL TEMPERATURE
1023	! =====================================================================
1024
1025	! dans convect3, tvp est calcule en une seule fois, et sans retirer
1026	! l'eau condensee (~> reversible CAPE)
1027
1028	! ori do 340 k=minorig+1,nl
1029	! ori do 330 i=1,ncum
1030	! ori if(k.ge.(icb(i)+1))then
1031	! ori tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
1032	! oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
1033	! oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1034	! ori endif
1035	! ori 330 continue
1036	! ori 340 continue
1037
1038	! ori do 350 i=1,ncum
1039	! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1040	! ori 350 continue
1041
1042	DO i = 1, ncum ! convect3
1043	tp(i, nlp) = tp(i, nl) ! convect3
1044	END DO ! convect3
1045
1046	! =====================================================================
1047	! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1048	! =====================================================================
1049
1050	! -- this is for convect3 only:
1051
1052	! first estimate of buoyancy:
1053
1054	DO i = 1, ncum
1055	DO k = 1, nl
1056	buoy(i, k) = tvp(i, k) - tv(i, k)
1057	END DO
1058	END DO
1059
1060	! set buoyancy=buoybase for all levels below base
1061	! for safety, set buoy(icb)=buoybase
1062
1063	DO i = 1, ncum
1064	DO k = 1, nl
1065	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1066	buoy(i, k) = buoybase(i)
1067	END IF
1068	END DO
1069	! IM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i)
1070	buoy(i, icb(i)) = buoybase(i)
1071	END DO
1072
1073	! -- end convect3
1074
1075	! =====================================================================
1076	! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1077	! --- LEVEL OF NEUTRAL BUOYANCY
1078	! =====================================================================
1079
1080	! -- this is for convect3 only:
1081
1082	DO i = 1, ncum
1083	inb(i) = nl - 1
1084	END DO
1085
1086	DO i = 1, ncum
1087	DO k = 1, nl - 1
1088	IF ((k>=icb(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1089	inb(i) = min(inb(i), k)
1090	END IF
1091	END DO
1092	END DO
1093
1094	! -- end convect3
1095
1096	! ori do 510 i=1,ncum
1097	! ori cape(i)=0.0
1098	! ori capem(i)=0.0
1099	! ori inb(i)=icb(i)+1
1100	! ori inb1(i)=inb(i)
1101	! ori 510 continue
1102
1103	! Originial Code
1104
1105	! do 530 k=minorig+1,nl-1
1106	! do 520 i=1,ncum
1107	! if(k.ge.(icb(i)+1))then
1108	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1109	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1110	! cape(i)=cape(i)+by
1111	! if(by.ge.0.0)inb1(i)=k+1
1112	! if(cape(i).gt.0.0)then
1113	! inb(i)=k+1
1114	! capem(i)=cape(i)
1115	! endif
1116	! endif
1117	! 520 continue
1118	! 530 continue
1119	! do 540 i=1,ncum
1120	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1121	! cape(i)=capem(i)+byp
1122	! defrac=capem(i)-cape(i)
1123	! defrac=max(defrac,0.001)
1124	! frac(i)=-cape(i)/defrac
1125	! frac(i)=min(frac(i),1.0)
1126	! frac(i)=max(frac(i),0.0)
1127	! 540 continue
1128
1129	! K Emanuel fix
1130
1131	! call zilch(byp,ncum)
1132	! do 530 k=minorig+1,nl-1
1133	! do 520 i=1,ncum
1134	! if(k.ge.(icb(i)+1))then
1135	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1136	! cape(i)=cape(i)+by
1137	! if(by.ge.0.0)inb1(i)=k+1
1138	! if(cape(i).gt.0.0)then
1139	! inb(i)=k+1
1140	! capem(i)=cape(i)
1141	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1142	! endif
1143	! endif
1144	! 520 continue
1145	! 530 continue
1146	! do 540 i=1,ncum
1147	! inb(i)=max(inb(i),inb1(i))
1148	! cape(i)=capem(i)+byp(i)
1149	! defrac=capem(i)-cape(i)
1150	! defrac=max(defrac,0.001)
1151	! frac(i)=-cape(i)/defrac
1152	! frac(i)=min(frac(i),1.0)
1153	! frac(i)=max(frac(i),0.0)
1154	! 540 continue
1155
1156	! J Teixeira fix
1157
1158	! ori call zilch(byp,ncum)
1159	! ori do 515 i=1,ncum
1160	! ori lcape(i)=.true.
1161	! ori 515 continue
1162	! ori do 530 k=minorig+1,nl-1
1163	! ori do 520 i=1,ncum
1164	! ori if(cape(i).lt.0.0)lcape(i)=.false.
1165	! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1166	! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1167	! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1168	! ori cape(i)=cape(i)+by
1169	! ori if(by.ge.0.0)inb1(i)=k+1
1170	! ori if(cape(i).gt.0.0)then
1171	! ori inb(i)=k+1
1172	! ori capem(i)=cape(i)
1173	! ori endif
1174	! ori endif
1175	! ori 520 continue
1176	! ori 530 continue
1177	! ori do 540 i=1,ncum
1178	! ori cape(i)=capem(i)+byp(i)
1179	! ori defrac=capem(i)-cape(i)
1180	! ori defrac=max(defrac,0.001)
1181	! ori frac(i)=-cape(i)/defrac
1182	! ori frac(i)=min(frac(i),1.0)
1183	! ori frac(i)=max(frac(i),0.0)
1184	! ori 540 continue
1185
1186	! =====================================================================
1187	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1188	! =====================================================================
1189
1190	! ym do i=1,ncum*nlp
1191	! ym hp(i,1)=h(i,1)
1192	! ym enddo
1193
1194	DO k = 1, nlp
1195	DO i = 1, ncum
1196	hp(i, k) = h(i, k)
1197	END DO
1198	END DO
1199
1200	DO k = minorig + 1, nl
1201	DO i = 1, ncum
1202	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1203	hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k &
1204	)
1205	END IF
1206	END DO
1207	END DO
1208
1209	RETURN
1210	END SUBROUTINE cv30_undilute2
1211
1212	SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
1213	sig, w0, cape, m)
1214	IMPLICIT NONE
1215
1216	! ===================================================================
1217	! --- CLOSURE OF CONVECT3
1218
1219	! vectorization: S. Bony
1220	! ===================================================================
1221
1222	include "cvthermo.h"
1223	include "cv30param.h"
1224
1225	! input:
1226	INTEGER ncum, nd, nloc
1227	INTEGER icb(nloc), inb(nloc)
1228	REAL pbase(nloc)
1229	REAL p(nloc, nd), ph(nloc, nd+1)
1230	REAL tv(nloc, nd), buoy(nloc, nd)
1231
1232	! input/output:
1233	REAL sig(nloc, nd), w0(nloc, nd)
1234
1235	! output:
1236	REAL cape(nloc)
1237	REAL m(nloc, nd)
1238
1239	! local variables:
1240	INTEGER i, j, k, icbmax
1241	REAL deltap, fac, w, amu
1242	REAL dtmin(nloc, nd), sigold(nloc, nd)
1243
1244	! -------------------------------------------------------
1245	! -- Initialization
1246	! -------------------------------------------------------
1247
1248	DO k = 1, nl
1249	DO i = 1, ncum
1250	m(i, k) = 0.0
1251	END DO
1252	END DO
1253
1254	! -------------------------------------------------------
1255	! -- Reset sig(i) and w0(i) for i>inb and i<icb
1256	! -------------------------------------------------------
1257
1258	! update sig and w0 above LNB:
1259
1260	DO k = 1, nl - 1
1261	DO i = 1, ncum
1262	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
1263	sig(i, k) = betasig(i, k) + 2.alphabuoy(i, inb(i))abs(buoy(i,inb( &
1264	i)))
1265	sig(i, k) = amax1(sig(i,k), 0.0)
1266	w0(i, k) = beta*w0(i, k)
1267	END IF
1268	END DO
1269	END DO
1270
1271	! compute icbmax:
1272
1273	icbmax = 2
1274	DO i = 1, ncum
1275	icbmax = max(icbmax, icb(i))
1276	END DO
1277
1278	! update sig and w0 below cloud base:
1279
1280	DO k = 1, icbmax
1281	DO i = 1, ncum
1282	IF (k<=icb(i)) THEN
1283	sig(i, k) = betasig(i, k) - 2.alphabuoy(i, icb(i))buoy(i, icb(i))
1284	sig(i, k) = amax1(sig(i,k), 0.0)
1285	w0(i, k) = beta*w0(i, k)
1286	END IF
1287	END DO
1288	END DO
1289
1290	! ! if(inb.lt.(nl-1))then
1291	! ! do 85 i=inb+1,nl-1
1292	! ! sig(i)=betasig(i)+2.alphabuoy(inb)
1293	! ! 1 abs(buoy(inb))
1294	! ! sig(i)=amax1(sig(i),0.0)
1295	! ! w0(i)=beta*w0(i)
1296	! ! 85 continue
1297	! ! end if
1298
1299	! ! do 87 i=1,icb
1300	! ! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
1301	! ! sig(i)=amax1(sig(i),0.0)
1302	! ! w0(i)=beta*w0(i)
1303	! ! 87 continue
1304
1305	! -------------------------------------------------------------
1306	! -- Reset fractional areas of updrafts and w0 at initial time
1307	! -- and after 10 time steps of no convection
1308	! -------------------------------------------------------------
1309
1310	DO k = 1, nl - 1
1311	DO i = 1, ncum
1312	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
1313	sig(i, k) = 0.0
1314	w0(i, k) = 0.0
1315	END IF
1316	END DO
1317	END DO
1318
1319	! -------------------------------------------------------------
1320	! -- Calculate convective available potential energy (cape),
1321	! -- vertical velocity (w), fractional area covered by
1322	! -- undilute updraft (sig), and updraft mass flux (m)
1323	! -------------------------------------------------------------
1324
1325	DO i = 1, ncum
1326	cape(i) = 0.0
1327	END DO
1328
1329	! compute dtmin (minimum buoyancy between ICB and given level k):
1330
1331	DO i = 1, ncum
1332	DO k = 1, nl
1333	dtmin(i, k) = 100.0
1334	END DO
1335	END DO
1336
1337	DO i = 1, ncum
1338	DO k = 1, nl
1339	DO j = minorig, nl
1340	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- &
1341	1))) THEN
1342	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
1343	END IF
1344	END DO
1345	END DO
1346	END DO
1347
1348	! the interval on which cape is computed starts at pbase :
1349	DO k = 1, nl
1350	DO i = 1, ncum
1351
1352	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
1353
1354	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
1355	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
1356	cape(i) = amax1(0.0, cape(i))
1357	sigold(i, k) = sig(i, k)
1358
1359	! dtmin(i,k)=100.0
1360	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
1361	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1362	! 97 continue
1363
1364	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
1365	sig(i, k) = amax1(sig(i,k), 0.0)
1366	sig(i, k) = amin1(sig(i,k), 0.01)
1367	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
1368	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
1369	amu = 0.5(sig(i,k)+sigold(i,k))w
1370	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
1371	w0(i, k) = w
1372	END IF
1373
1374	END DO
1375	END DO
1376
1377	DO i = 1, ncum
1378	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
1379	m(i, icb(i)) = 0.5m(i, icb(i)+1)(ph(i,icb(i))-ph(i,icb(i)+1))/ &
1380	(ph(i,icb(i)+1)-ph(i,icb(i)+2))
1381	sig(i, icb(i)) = sig(i, icb(i)+1)
1382	sig(i, icb(i)-1) = sig(i, icb(i))
1383	END DO
1384
1385
1386	! ! cape=0.0
1387	! ! do 98 i=icb+1,inb
1388	! ! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1389	! ! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
1390	! ! dcape=rrdbuoy(i-1)deltap/p(i-1)
1391	! ! dlnp=deltap/p(i-1)
1392	! ! cape=amax1(0.0,cape)
1393	! ! sigold=sig(i)
1394
1395	! ! dtmin=100.0
1396	! ! do 97 j=icb,i-1
1397	! ! dtmin=amin1(dtmin,buoy(j))
1398	! ! 97 continue
1399
1400	! ! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
1401	! ! sig(i)=amax1(sig(i),0.0)
1402	! ! sig(i)=amin1(sig(i),0.01)
1403	! ! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1404	! ! w=(1.-beta)facsqrt(cape)+beta*w0(i)
1405	! ! amu=0.5(sig(i)+sigold)w
1406	! ! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
1407	! ! w0(i)=w
1408	! ! 98 continue
1409	! ! w0(icb)=0.5*w0(icb+1)
1410	! ! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1411	! ! sig(icb)=sig(icb+1)
1412	! ! sig(icb-1)=sig(icb)
1413
1414	RETURN
1415	END SUBROUTINE cv30_closure
1416
1417	SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
1418	u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, &
1419	vent, sij, elij, ments, qents, traent)
1420	IMPLICIT NONE
1421
1422	! ---------------------------------------------------------------------
1423	! a faire:
1424	! - changer rr(il,1) -> qnk(il)
1425	! - vectorisation de la partie normalisation des flux (do 789...)
1426	! ---------------------------------------------------------------------
1427
1428	include "cvthermo.h"
1429	include "cv30param.h"
1430
1431	! inputs:
1432	INTEGER ncum, nd, na, ntra, nloc
1433	INTEGER icb(nloc), inb(nloc), nk(nloc)
1434	REAL sig(nloc, nd)
1435	REAL qnk(nloc)
1436	REAL ph(nloc, nd+1)
1437	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1438	REAL u(nloc, nd), v(nloc, nd)
1439	REAL tra(nloc, nd, ntra) ! input of convect3
1440	REAL lv(nloc, na), h(nloc, na), hp(nloc, na)
1441	REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
1442	REAL m(nloc, na) ! input of convect3
1443
1444	! outputs:
1445	REAL ment(nloc, na, na), qent(nloc, na, na)
1446	REAL uent(nloc, na, na), vent(nloc, na, na)
1447	REAL sij(nloc, na, na), elij(nloc, na, na)
1448	REAL traent(nloc, nd, nd, ntra)
1449	REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
1450	REAL sigij(nloc, nd, nd)
1451
1452	! local variables:
1453	INTEGER i, j, k, il, im, jm
1454	INTEGER num1, num2
1455	INTEGER nent(nloc, na)
1456	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1457	REAL alt, smid, sjmin, sjmax, delp, delm
1458	REAL asij(nloc), smax(nloc), scrit(nloc)
1459	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
1460	REAL wgh
1461	REAL zm(nloc, na)
1462	LOGICAL lwork(nloc)
1463
1464	! =====================================================================
1465	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1466	! =====================================================================
1467
1468	! ori do 360 i=1,ncum*nlp
1469	DO j = 1, nl
1470	DO i = 1, ncum
1471	nent(i, j) = 0
1472	! in convect3, m is computed in cv3_closure
1473	! ori m(i,1)=0.0
1474	END DO
1475	END DO
1476
1477	! ori do 400 k=1,nlp
1478	! ori do 390 j=1,nlp
1479	DO j = 1, nl
1480	DO k = 1, nl
1481	DO i = 1, ncum
1482	qent(i, k, j) = rr(i, j)
1483	uent(i, k, j) = u(i, j)
1484	vent(i, k, j) = v(i, j)
1485	elij(i, k, j) = 0.0
1486	! ym ment(i,k,j)=0.0
1487	! ym sij(i,k,j)=0.0
1488	END DO
1489	END DO
1490	END DO
1491
1492	! ym
1493	ment(1:ncum, 1:nd, 1:nd) = 0.0
1494	sij(1:ncum, 1:nd, 1:nd) = 0.0
1495
1496	! do k=1,ntra
1497	! do j=1,nd ! instead nlp
1498	! do i=1,nd ! instead nlp
1499	! do il=1,ncum
1500	! traent(il,i,j,k)=tra(il,j,k)
1501	! enddo
1502	! enddo
1503	! enddo
1504	! enddo
1505	zm(:, :) = 0.
1506
1507	! =====================================================================
1508	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1509	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1510	! --- FRACTION (sij)
1511	! =====================================================================
1512
1513	DO i = minorig + 1, nl
1514
1515	DO j = minorig, nl
1516	DO il = 1, ncum
1517	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
1518	1)) .AND. (j<=inb(il))) THEN
1519
1520	rti = rr(il, 1) - ep(il, i)*clw(il, i)
1521	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
1522	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
1523	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
1524	dei = denom
1525	IF (abs(dei)<0.01) dei = 0.01
1526	sij(il, i, j) = anum/dei
1527	sij(il, i, i) = 1.0
1528	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
1529	altem = altem/bf2
1530	cwat = clw(il, j)*(1.-ep(il,j))
1531	stemp = sij(il, i, j)
1532	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
1533	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
1534	denom = denom + lv(il, j)*(rr(il,i)-rti)
1535	IF (abs(denom)<0.01) denom = 0.01
1536	sij(il, i, j) = anum/denom
1537	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - &
1538	rs(il, j)
1539	altem = altem - (bf2-1.)*cwat
1540	END IF
1541	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
1542	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
1543	uent(il, i, j) = sij(il, i, j)*u(il, i) + &
1544	(1.-sij(il,i,j))*u(il, nk(il))
1545	vent(il, i, j) = sij(il, i, j)*v(il, i) + &
1546	(1.-sij(il,i,j))*v(il, nk(il))
1547	! !!! do k=1,ntra
1548	! !!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1549	! !!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1550	! !!! end do
1551	elij(il, i, j) = altem
1552	elij(il, i, j) = amax1(0.0, elij(il,i,j))
1553	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
1554	nent(il, i) = nent(il, i) + 1
1555	END IF
1556	sij(il, i, j) = amax1(0.0, sij(il,i,j))
1557	sij(il, i, j) = amin1(1.0, sij(il,i,j))
1558	END IF ! new
1559	END DO
1560	END DO
1561
1562	! do k=1,ntra
1563	! do j=minorig,nl
1564	! do il=1,ncum
1565	! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1566	! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1567	! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1568	! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1569	! endif
1570	! enddo
1571	! enddo
1572	! enddo
1573
1574
1575	! *** if no air can entrain at level i assume that updraft detrains
1576	! ***
1577	! *** at that level and calculate detrained air flux and properties
1578	! ***
1579
1580
1581	! @ do 170 i=icb(il),inb(il)
1582
1583	DO il = 1, ncum
1584	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
1585	! @ if(nent(il,i).eq.0)then
1586	ment(il, i, i) = m(il, i)
1587	qent(il, i, i) = rr(il, nk(il)) - ep(il, i)*clw(il, i)
1588	uent(il, i, i) = u(il, nk(il))
1589	vent(il, i, i) = v(il, nk(il))
1590	elij(il, i, i) = clw(il, i)
1591	! MAF sij(il,i,i)=1.0
1592	sij(il, i, i) = 0.0
1593	END IF
1594	END DO
1595	END DO
1596
1597	! do j=1,ntra
1598	! do i=minorig+1,nl
1599	! do il=1,ncum
1600	! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1601	! traent(il,i,i,j)=tra(il,nk(il),j)
1602	! endif
1603	! enddo
1604	! enddo
1605	! enddo
1606
1607	DO j = minorig, nl
1608	DO i = minorig, nl
1609	DO il = 1, ncum
1610	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
1611	inb(il))) THEN
1612	sigij(il, i, j) = sij(il, i, j)
1613	END IF
1614	END DO
1615	END DO
1616	END DO
1617	! @ enddo
1618
1619	! @170 continue
1620
1621	! =====================================================================
1622	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
1623	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
1624	! =====================================================================
1625
1626	! ym call zilch(asum,ncum*nd)
1627	! ym call zilch(bsum,ncum*nd)
1628	! ym call zilch(csum,ncum*nd)
1629	CALL zilch(asum, nloc*nd)
1630	CALL zilch(csum, nloc*nd)
1631	CALL zilch(csum, nloc*nd)
1632
1633	DO il = 1, ncum
1634	lwork(il) = .FALSE.
1635	END DO
1636
1637	DO i = minorig + 1, nl
1638
1639	num1 = 0
1640	DO il = 1, ncum
1641	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
1642	END DO
1643	IF (num1<=0) GO TO 789
1644
1645
1646	DO il = 1, ncum
1647	IF (i>=icb(il) .AND. i<=inb(il)) THEN
1648	lwork(il) = (nent(il,i)/=0)
1649	qp = rr(il, 1) - ep(il, i)*clw(il, i)
1650	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
1651	(cpv-cpd)t(il, i)(qp-rr(il,i))
1652	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
1653	(cpd-cpv)t(il, i)(rr(il,i)-qp)
1654	IF (abs(denom)<0.01) denom = 0.01
1655	scrit(il) = anum/denom
1656	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
1657	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
1658	smax(il) = 0.0
1659	asij(il) = 0.0
1660	END IF
1661	END DO
1662
1663	DO j = nl, minorig, -1
1664
1665	num2 = 0
1666	DO il = 1, ncum
1667	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
1668	il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
1669	END DO
1670	IF (num2<=0) GO TO 175
1671
1672	DO il = 1, ncum
1673	IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
1674	il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
1675
1676	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
1677	wgh = 1.0
1678	IF (j>i) THEN
1679	sjmax = amax1(sij(il,i,j+1), smax(il))
1680	sjmax = amin1(sjmax, scrit(il))
1681	smax(il) = amax1(sij(il,i,j), smax(il))
1682	sjmin = amax1(sij(il,i,j-1), smax(il))
1683	sjmin = amin1(sjmin, scrit(il))
1684	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
1685	smid = amin1(sij(il,i,j), scrit(il))
1686	ELSE
1687	sjmax = amax1(sij(il,i,j+1), scrit(il))
1688	smid = amax1(sij(il,i,j), scrit(il))
1689	sjmin = 0.0
1690	IF (j>1) sjmin = sij(il, i, j-1)
1691	sjmin = amax1(sjmin, scrit(il))
1692	END IF
1693	delp = abs(sjmax-smid)
1694	delm = abs(sjmin-smid)
1695	asij(il) = asij(il) + wgh*(delp+delm)
1696	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
1697	END IF
1698	END IF
1699	END DO
1700
1701	175 END DO
1702
1703	DO il = 1, ncum
1704	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1705	asij(il) = amax1(1.0E-16, asij(il))
1706	asij(il) = 1.0/asij(il)
1707	asum(il, i) = 0.0
1708	bsum(il, i) = 0.0
1709	csum(il, i) = 0.0
1710	END IF
1711	END DO
1712
1713	DO j = minorig, nl
1714	DO il = 1, ncum
1715	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1716	il)-1) .AND. j<=inb(il)) THEN
1717	ment(il, i, j) = ment(il, i, j)*asij(il)
1718	END IF
1719	END DO
1720	END DO
1721
1722	DO j = minorig, nl
1723	DO il = 1, ncum
1724	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1725	il)-1) .AND. j<=inb(il)) THEN
1726	asum(il, i) = asum(il, i) + ment(il, i, j)
1727	ment(il, i, j) = ment(il, i, j)*sig(il, j)
1728	bsum(il, i) = bsum(il, i) + ment(il, i, j)
1729	END IF
1730	END DO
1731	END DO
1732
1733	DO il = 1, ncum
1734	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1735	bsum(il, i) = amax1(bsum(il,i), 1.0E-16)
1736	bsum(il, i) = 1.0/bsum(il, i)
1737	END IF
1738	END DO
1739
1740	DO j = minorig, nl
1741	DO il = 1, ncum
1742	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1743	il)-1) .AND. j<=inb(il)) THEN
1744	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
1745	END IF
1746	END DO
1747	END DO
1748
1749	DO j = minorig, nl
1750	DO il = 1, ncum
1751	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1752	il)-1) .AND. j<=inb(il)) THEN
1753	csum(il, i) = csum(il, i) + ment(il, i, j)
1754	END IF
1755	END DO
1756	END DO
1757
1758	DO il = 1, ncum
1759	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1760	csum(il,i)<m(il,i)) THEN
1761	nent(il, i) = 0
1762	ment(il, i, i) = m(il, i)
1763	qent(il, i, i) = rr(il, 1) - ep(il, i)*clw(il, i)
1764	uent(il, i, i) = u(il, nk(il))
1765	vent(il, i, i) = v(il, nk(il))
1766	elij(il, i, i) = clw(il, i)
1767	! MAF sij(il,i,i)=1.0
1768	sij(il, i, i) = 0.0
1769	END IF
1770	END DO ! il
1771
1772	! do j=1,ntra
1773	! do il=1,ncum
1774	! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1775	! : .and. csum(il,i).lt.m(il,i) ) then
1776	! traent(il,i,i,j)=tra(il,nk(il),j)
1777	! endif
1778	! enddo
1779	! enddo
1780	789 END DO
1781
1782	! MAF: renormalisation de MENT
1783	DO jm = 1, nd
1784	DO im = 1, nd
1785	DO il = 1, ncum
1786	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
1787	END DO
1788	END DO
1789	END DO
1790
1791	DO jm = 1, nd
1792	DO im = 1, nd
1793	DO il = 1, ncum
1794	IF (zm(il,im)/=0.) THEN
1795	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
1796	END IF
1797	END DO
1798	END DO
1799	END DO
1800
1801	DO jm = 1, nd
1802	DO im = 1, nd
1803	DO il = 1, ncum
1804	qents(il, im, jm) = qent(il, im, jm)
1805	ments(il, im, jm) = ment(il, im, jm)
1806	END DO
1807	END DO
1808	END DO
1809
1810	RETURN
1811	END SUBROUTINE cv30_mixing
1812
1813
1814	SUBROUTINE cv30_unsat(nloc, ncum, nd, na, ntra, icb, inb, t, rr, rs, gz, u, &
1815	v, tra, p, ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, &
1816	mp, rp, up, vp, trap, wt, water, evap, b & ! RomP-jyg
1817	, wdtraina, wdtrainm) ! 26/08/10 RomP-jyg
1818	IMPLICIT NONE
1819
1820
1821	include "cvthermo.h"
1822	include "cv30param.h"
1823	include "cvflag.h"
1824
1825	! inputs:
1826	INTEGER ncum, nd, na, ntra, nloc
1827	INTEGER icb(nloc), inb(nloc)
1828	REAL delt, plcl(nloc)
1829	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1830	REAL u(nloc, nd), v(nloc, nd)
1831	REAL tra(nloc, nd, ntra)
1832	REAL p(nloc, nd), ph(nloc, nd+1)
1833	REAL th(nloc, na), gz(nloc, na)
1834	REAL lv(nloc, na), ep(nloc, na), sigp(nloc, na), clw(nloc, na)
1835	REAL cpn(nloc, na), tv(nloc, na)
1836	REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
1837
1838	! outputs:
1839	REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
1840	REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
1841	REAL trap(nloc, na, ntra)
1842	REAL b(nloc, na)
1843	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
1844	! lascendance adiabatique et des flux melanges Pa et Pm.
1845	! Distinction des wdtrain
1846	! Pa = wdtrainA Pm = wdtrainM
1847	REAL wdtraina(nloc, na), wdtrainm(nloc, na)
1848
1849	! local variables
1850	INTEGER i, j, k, il, num1
1851	REAL tinv, delti
1852	REAL awat, afac, afac1, afac2, bfac
1853	REAL pr1, pr2, sigt, b6, c6, revap, tevap, delth
1854	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1855	REAL ampmax
1856	REAL lvcp(nloc, na)
1857	REAL wdtrain(nloc)
1858	LOGICAL lwork(nloc)
1859
1860
1861	! ------------------------------------------------------
1862
1863	delti = 1./delt
1864	tinv = 1./3.
1865
1866	mp(:, :) = 0.
1867
1868	DO i = 1, nl
1869	DO il = 1, ncum
1870	mp(il, i) = 0.0
1871	rp(il, i) = rr(il, i)
1872	up(il, i) = u(il, i)
1873	vp(il, i) = v(il, i)
1874	wt(il, i) = 0.001
1875	water(il, i) = 0.0
1876	evap(il, i) = 0.0
1877	b(il, i) = 0.0
1878	lvcp(il, i) = lv(il, i)/cpn(il, i)
1879	END DO
1880	END DO
1881
1882	! do k=1,ntra
1883	! do i=1,nd
1884	! do il=1,ncum
1885	! trap(il,i,k)=tra(il,i,k)
1886	! enddo
1887	! enddo
1888	! enddo
1889	! ! RomP >>>
1890	DO i = 1, nd
1891	DO il = 1, ncum
1892	wdtraina(il, i) = 0.0
1893	wdtrainm(il, i) = 0.0
1894	END DO
1895	END DO
1896	! ! RomP <<<
1897
1898	! * check whether ep(inb)=0, if so, skip precipitating *
1899	! * downdraft calculation *
1900
1901
1902	DO il = 1, ncum
1903	lwork(il) = .TRUE.
1904	IF (ep(il,inb(il))<0.0001) lwork(il) = .FALSE.
1905	END DO
1906
1907	CALL zilch(wdtrain, ncum)
1908
1909	DO i = nl + 1, 1, -1
1910
1911	num1 = 0
1912	DO il = 1, ncum
1913	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
1914	END DO
1915	IF (num1<=0) GO TO 400
1916
1917
1918	! * integrate liquid water equation to find condensed water *
1919	! * and condensed water flux *
1920
1921
1922
1923	! * begin downdraft loop *
1924
1925
1926
1927	! * calculate detrained precipitation *
1928
1929	DO il = 1, ncum
1930	IF (i<=inb(il) .AND. lwork(il)) THEN
1931	IF (cvflag_grav) THEN
1932	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
1933	wdtraina(il, i) = wdtrain(il)/grav ! Pa 26/08/10 RomP
1934	ELSE
1935	wdtrain(il) = 10.0ep(il, i)m(il, i)*clw(il, i)
1936	wdtraina(il, i) = wdtrain(il)/10. ! Pa 26/08/10 RomP
1937	END IF
1938	END IF
1939	END DO
1940
1941	IF (i>1) THEN
1942
1943	DO j = 1, i - 1
1944	DO il = 1, ncum
1945	IF (i<=inb(il) .AND. lwork(il)) THEN
1946	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
1947	awat = amax1(awat, 0.0)
1948	IF (cvflag_grav) THEN
1949	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
1950	ELSE
1951	wdtrain(il) = wdtrain(il) + 10.0awatment(il, j, i)
1952	END IF
1953	END IF
1954	END DO
1955	END DO
1956	DO il = 1, ncum
1957	IF (cvflag_grav) THEN
1958	wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) ! Pm 26/08/10 RomP
1959	ELSE
1960	wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) ! Pm 26/08/10 RomP
1961	END IF
1962	END DO
1963
1964	END IF
1965
1966
1967	! * find rain water and evaporation using provisional *
1968	! * estimates of rp(i)and rp(i-1) *
1969
1970
1971	DO il = 1, ncum
1972
1973	IF (i<=inb(il) .AND. lwork(il)) THEN
1974
1975	wt(il, i) = 45.0
1976
1977	IF (i<inb(il)) THEN
1978	rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, &
1979	i))+gz(il,i+1)-gz(il,i))/lv(il, i)
1980	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
1981	END IF
1982	rp(il, i) = amax1(rp(il,i), 0.0)
1983	rp(il, i) = amin1(rp(il,i), rs(il,i))
1984	rp(il, inb(il)) = rr(il, inb(il))
1985
1986	IF (i==1) THEN
1987	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
1988	ELSE
1989	rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, &
1990	i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
1991	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
1992	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
1993	rp(il, i-1) = amax1(rp(il,i-1), 0.0)
1994	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i) &
1995	)
1996	afac2 = p(il, i-1)*(rs(il,i-1)-rp(il,i-1))/ &
1997	(1.0E4+2000.0p(il,i-1)rs(il,i-1))
1998	afac = 0.5*(afac1+afac2)
1999	END IF
2000	IF (i==inb(il)) afac = 0.0
2001	afac = amax1(afac, 0.0)
2002	bfac = 1./(sigd*wt(il,i))
2003
2004	! jyg1
2005	! cc sigt=1.0
2006	! cc if(i.ge.icb)sigt=sigp(i)
2007	! prise en compte de la variation progressive de sigt dans
2008	! les couches icb et icb-1:
2009	! pour plcl<ph(i+1), pr1=0 & pr2=1
2010	! pour plcl>ph(i), pr1=1 & pr2=0
2011	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2012	! sur le nuage, et pr2 est la proportion sous la base du
2013	! nuage.
2014	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2015	pr1 = max(0., min(1.,pr1))
2016	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2017	pr2 = max(0., min(1.,pr2))
2018	sigt = sigp(il, i)*pr1 + pr2
2019	! jyg2
2020
2021	b6 = bfac50.sigd(ph(il,i)-ph(il,i+1))sigt*afac
2022	c6 = water(il, i+1) + bfacwdtrain(il) - 50.sigdbfac(ph(il,i)-ph( &
2023	il,i+1))*evap(il, i+1)
2024	IF (c6>0.0) THEN
2025	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
2026	evap(il, i) = sigtafacrevap
2027	water(il, i) = revap*revap
2028	ELSE
2029	evap(il, i) = -evap(il, i+1) + 0.02(wdtrain(il)+sigdwt(il,i)* &
2030	water(il,i+1))/(sigd*(ph(il,i)-ph(il,i+1)))
2031	END IF
2032
2033	! * calculate precipitating downdraft mass flux under *
2034	! * hydrostatic approximation *
2035
2036	IF (i/=1) THEN
2037
2038	tevap = amax1(0.0, evap(il,i))
2039	delth = amax1(0.001, (th(il,i)-th(il,i-1)))
2040	IF (cvflag_grav) THEN
2041	mp(il, i) = 100.ginvlvcp(il, i)sigdtevap*(p(il,i-1)-p(il,i))/ &
2042	delth
2043	ELSE
2044	mp(il, i) = 10.lvcp(il, i)sigdtevap(p(il,i-1)-p(il,i))/delth
2045	END IF
2046
2047	! * if hydrostatic assumption fails, *
2048	! * solve cubic difference equation for downdraft theta *
2049	! * and mass flux from two simultaneous differential eqns *
2050
2051	amfac = sigdsigd70.0ph(il, i)(p(il,i-1)-p(il,i))* &
2052	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2053	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2054	IF (amp2>(0.1*amfac)) THEN
2055	xf = 100.0sigdsigdsigd(ph(il,i)-ph(il,i+1))
2056	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i)/(lvcp(il,i)* &
2057	sigd*th(il,i))
2058	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
2059	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2060	50.(p(il,i-1)-p(il,i))xf*tevap
2061	fac2 = 1.0
2062	IF (bf<0.0) fac2 = -1.0
2063	bf = abs(bf)
2064	ur = 0.25bfbf - afafaftinvtinv*tinv
2065	IF (ur>=0.0) THEN
2066	sru = sqrt(ur)
2067	fac = 1.0
2068	IF ((0.5*bf-sru)<0.0) fac = -1.0
2069	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
2070	fac(abs(0.5bf-sru))**tinv
2071	ELSE
2072	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
2073	IF (fac2<0.0) d = 3.14159 - d
2074	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
2075	END IF
2076	mp(il, i) = amax1(0.0, mp(il,i))
2077
2078	IF (cvflag_grav) THEN
2079	! jyg : il y a vraisemblablement une erreur dans la ligne 2
2080	! suivante:
2081	! il faut diviser par (mp(il,i)sigdgrav) et non par
2082	! (mp(il,i)+sigd*0.1).
2083	! Et il faut bien revoir les facteurs 100.
2084	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap/(mp(il, &
2085	i)+sigd0.1) - 10.0(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
2086	)sigdth(il,i))
2087	ELSE
2088	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap/(mp(il, &
2089	i)+sigd0.1) - 10.0(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
2090	)sigdth(il,i))
2091	END IF
2092	b(il, i-1) = amax1(b(il,i-1), 0.0)
2093	END IF
2094
2095	! *** limit magnitude of mp(i) to meet cfl condition
2096	! ***
2097
2098	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
2099	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
2100	ampmax = amin1(ampmax, amp2)
2101	mp(il, i) = amin1(mp(il,i), ampmax)
2102
2103	! *** force mp to decrease linearly to zero
2104	! ***
2105	! *** between cloud base and the surface
2106	! ***
2107
2108	IF (p(il,i)>p(il,icb(il))) THEN
2109	mp(il, i) = mp(il, icb(il))*(p(il,1)-p(il,i))/ &
2110	(p(il,1)-p(il,icb(il)))
2111	END IF
2112
2113	END IF ! i.eq.1
2114
2115	! * find mixing ratio of precipitating downdraft *
2116
2117
2118	IF (i/=inb(il)) THEN
2119
2120	rp(il, i) = rr(il, i)
2121
2122	IF (mp(il,i)>mp(il,i+1)) THEN
2123
2124	IF (cvflag_grav) THEN
2125	rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
2126	rr(il, i)(mp(il,i)-mp(il,i+1)) + 100.ginv0.5sigd*(ph(il,i &
2127	)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2128	ELSE
2129	rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
2130	rr(il, i)(mp(il,i)-mp(il,i+1)) + 5.sigd*(ph(il,i)-ph(il,i+1 &
2131	))*(evap(il,i+1)+evap(il,i))
2132	END IF
2133	rp(il, i) = rp(il, i)/mp(il, i)
2134	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+ &
2135	1))
2136	up(il, i) = up(il, i)/mp(il, i)
2137	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+ &
2138	1))
2139	vp(il, i) = vp(il, i)/mp(il, i)
2140
2141	! do j=1,ntra
2142	! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2143	! testmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2144	! : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2145	! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2146	! end do
2147
2148	ELSE
2149
2150	IF (mp(il,i+1)>1.0E-16) THEN
2151	IF (cvflag_grav) THEN
2152	rp(il, i) = rp(il, i+1) + 100.ginv0.5sigd(ph(il,i)-ph(il, &
2153	i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1)
2154	ELSE
2155	rp(il, i) = rp(il, i+1) + 5.sigd(ph(il,i)-ph(il,i+1))*(evap &
2156	(il,i+1)+evap(il,i))/mp(il, i+1)
2157	END IF
2158	up(il, i) = up(il, i+1)
2159	vp(il, i) = vp(il, i+1)
2160
2161	! do j=1,ntra
2162	! trap(il,i,j)=trap(il,i+1,j)
2163	! end do
2164
2165	END IF
2166	END IF
2167	rp(il, i) = amin1(rp(il,i), rs(il,i))
2168	rp(il, i) = amax1(rp(il,i), 0.0)
2169
2170	END IF
2171	END IF
2172	END DO
2173
2174	400 END DO
2175
2176	RETURN
2177	END SUBROUTINE cv30_unsat
2178
2179	SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, &
2180	tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, &
2181	wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, &
2182	tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, &
2183	mike, tls, tps, qcondc, wd)
2184	IMPLICIT NONE
2185
2186	include "cvthermo.h"
2187	include "cv30param.h"
2188	include "cvflag.h"
2189	include "conema3.h"
2190
2191	! inputs:
2192	INTEGER ncum, nd, na, ntra, nloc
2193	INTEGER icb(nloc), inb(nloc)
2194	REAL delt
2195	REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
2196	REAL tra(nloc, nd, ntra), sig(nloc, nd)
2197	REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
2198	REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
2199	REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
2200	REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
2201	REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
2202	REAL water(nloc, na), evap(nloc, na), b(nloc, na)
2203	REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
2204	! ym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2205	REAL vent(nloc, na, na), elij(nloc, na, na)
2206	INTEGER nent(nloc, na)
2207	REAL traent(nloc, na, na, ntra)
2208	REAL tv(nloc, nd), tvp(nloc, nd)
2209
2210	! input/output:
2211	INTEGER iflag(nloc)
2212
2213	! outputs:
2214	REAL precip(nloc)
2215	REAL vprecip(nloc, nd+1)
2216	REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
2217	REAL ftra(nloc, nd, ntra)
2218	REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
2219	REAL dnwd0(nloc, nd), mike(nloc, nd)
2220	REAL tls(nloc, nd), tps(nloc, nd)
2221	REAL qcondc(nloc, nd) ! cld
2222	REAL wd(nloc) ! gust
2223
2224	! local variables:
2225	INTEGER i, k, il, n, j, num1
2226	REAL rat, awat, delti
2227	REAL ax, bx, cx, dx, ex
2228	REAL cpinv, rdcp, dpinv
2229	REAL lvcp(nloc, na), mke(nloc, na)
2230	REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
2231	! !! real up1(nloc), dn1(nloc)
2232	REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
2233	REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2234	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
2235	REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
2236
2237
2238	! -------------------------------------------------------------
2239
2240	! initialization:
2241
2242	delti = 1.0/delt
2243
2244	DO il = 1, ncum
2245	precip(il) = 0.0
2246	wd(il) = 0.0 ! gust
2247	vprecip(il, nd+1) = 0.
2248	END DO
2249
2250	DO i = 1, nd
2251	DO il = 1, ncum
2252	vprecip(il, i) = 0.0
2253	ft(il, i) = 0.0
2254	fr(il, i) = 0.0
2255	fu(il, i) = 0.0
2256	fv(il, i) = 0.0
2257	qcondc(il, i) = 0.0 ! cld
2258	qcond(il, i) = 0.0 ! cld
2259	nqcond(il, i) = 0.0 ! cld
2260	END DO
2261	END DO
2262
2263	! do j=1,ntra
2264	! do i=1,nd
2265	! do il=1,ncum
2266	! ftra(il,i,j)=0.0
2267	! enddo
2268	! enddo
2269	! enddo
2270
2271	DO i = 1, nl
2272	DO il = 1, ncum
2273	lvcp(il, i) = lv(il, i)/cpn(il, i)
2274	END DO
2275	END DO
2276
2277
2278
2279	! * calculate surface precipitation in mm/day *
2280
2281	DO il = 1, ncum
2282	IF (ep(il,inb(il))>=0.0001) THEN
2283	IF (cvflag_grav) THEN
2284	precip(il) = wt(il, 1)sigdwater(il, 1)86400.1000./(rowl*grav)
2285	ELSE
2286	precip(il) = wt(il, 1)sigdwater(il, 1)*8640.
2287	END IF
2288	END IF
2289	END DO
2290
2291	! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
2292
2293	! MAF rajout pour lessivage
2294	DO k = 1, nl
2295	DO il = 1, ncum
2296	IF (k<=inb(il)) THEN
2297	IF (cvflag_grav) THEN
2298	vprecip(il, k) = wt(il, k)sigdwater(il, k)/grav
2299	ELSE
2300	vprecip(il, k) = wt(il, k)sigdwater(il, k)/10.
2301	END IF
2302	END IF
2303	END DO
2304	END DO
2305
2306
2307	! * Calculate downdraft velocity scale *
2308	! * NE PAS UTILISER POUR L'INSTANT *
2309
2310	! ! do il=1,ncum
2311	! ! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
2312	! ! : /(sigd*p(il,icb(il)))
2313	! ! enddo
2314
2315
2316	! * calculate tendencies of lowest level potential temperature *
2317	! * and mixing ratio *
2318
2319	DO il = 1, ncum
2320	work(il) = 1.0/(ph(il,1)-ph(il,2))
2321	am(il) = 0.0
2322	END DO
2323
2324	DO k = 2, nl
2325	DO il = 1, ncum
2326	IF (k<=inb(il)) THEN
2327	am(il) = am(il) + m(il, k)
2328	END IF
2329	END DO
2330	END DO
2331
2332	DO il = 1, ncum
2333
2334	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
2335	IF (cvflag_grav) THEN
2336	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
2337	ft(il, 1) = 0.01gravwork(il)am(il)(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
2338	1))/cpn(il,1))
2339	ELSE
2340	IF ((0.1work(il)am(il))>=delti) iflag(il) = 1 !consistency vect
2341	ft(il, 1) = 0.1work(il)am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
2342	1))/cpn(il,1))
2343	END IF
2344
2345	ft(il, 1) = ft(il, 1) - 0.5lvcp(il, 1)sigd*(evap(il,1)+evap(il,2))
2346
2347	IF (cvflag_grav) THEN
2348	ft(il, 1) = ft(il, 1) - 0.009gravsigdmp(il, 2)t(il, 1)b(il, 1) &
2349	work(il)
2350	ELSE
2351	ft(il, 1) = ft(il, 1) - 0.09sigdmp(il, 2)t(il, 1)b(il, 1)*work(il)
2352	END IF
2353
2354	ft(il, 1) = ft(il, 1) + 0.01sigdwt(il, 1)(cl-cpd)water(il, 2)*(t(il,2 &
2355	)-t(il,1))*work(il)/cpn(il, 1)
2356
2357	IF (cvflag_grav) THEN
2358	! jyg1 Correction pour mieux conserver l'eau (conformite avec
2359	! CONVECT4.3)
2360	! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas
2361	! evap)
2362	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr(il,1))work(il) + &
2363	sigd0.5(evap(il,1)+evap(il,2))
2364	! +tard : +sigd*evap(il,1)
2365
2366	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il)
2367
2368	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il, &
2369	1))+am(il)*(u(il,2)-u(il,1)))
2370	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il, &
2371	1))+am(il)*(v(il,2)-v(il,1)))
2372	ELSE ! cvflag_grav
2373	fr(il, 1) = 0.1mp(il, 2)(rp(il,2)-rr(il,1))*work(il) + &
2374	sigd0.5(evap(il,1)+evap(il,2))
2375	fr(il, 1) = fr(il, 1) + 0.1am(il)(rr(il,2)-rr(il,1))*work(il)
2376	fu(il, 1) = fu(il, 1) + 0.1work(il)(mp(il,2)*(up(il,2)-u(il, &
2377	1))+am(il)*(u(il,2)-u(il,1)))
2378	fv(il, 1) = fv(il, 1) + 0.1work(il)(mp(il,2)*(vp(il,2)-v(il, &
2379	1))+am(il)*(v(il,2)-v(il,1)))
2380	END IF ! cvflag_grav
2381
2382	END DO ! il
2383
2384	! do j=1,ntra
2385	! do il=1,ncum
2386	! if (cvflag_grav) then
2387	! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
2388	! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2389	! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2390	! else
2391	! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2392	! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2393	! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2394	! endif
2395	! enddo
2396	! enddo
2397
2398	DO j = 2, nl
2399	DO il = 1, ncum
2400	IF (j<=inb(il)) THEN
2401	IF (cvflag_grav) THEN
2402	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il, &
2403	j,1)-rr(il,1))
2404	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il, &
2405	j,1)-u(il,1))
2406	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il, &
2407	j,1)-v(il,1))
2408	ELSE ! cvflag_grav
2409	fr(il, 1) = fr(il, 1) + 0.1work(il)ment(il, j, 1)*(qent(il,j,1)- &
2410	rr(il,1))
2411	fu(il, 1) = fu(il, 1) + 0.1work(il)ment(il, j, 1)*(uent(il,j,1)-u &
2412	(il,1))
2413	fv(il, 1) = fv(il, 1) + 0.1work(il)ment(il, j, 1)*(vent(il,j,1)-v &
2414	(il,1))
2415	END IF ! cvflag_grav
2416	END IF ! j
2417	END DO
2418	END DO
2419
2420	! do k=1,ntra
2421	! do j=2,nl
2422	! do il=1,ncum
2423	! if (j.le.inb(il)) then
2424
2425	! if (cvflag_grav) then
2426	! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
2427	! : *(traent(il,j,1,k)-tra(il,1,k))
2428	! else
2429	! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
2430	! : *(traent(il,j,1,k)-tra(il,1,k))
2431	! endif
2432
2433	! endif
2434	! enddo
2435	! enddo
2436	! enddo
2437
2438
2439	! * calculate tendencies of potential temperature and mixing ratio *
2440	! * at levels above the lowest level *
2441
2442	! * first find the net saturated updraft and downdraft mass fluxes *
2443	! * through each level *
2444
2445
2446	DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
2447
2448	num1 = 0
2449	DO il = 1, ncum
2450	IF (i<=inb(il)) num1 = num1 + 1
2451	END DO
2452	IF (num1<=0) GO TO 500
2453
2454	CALL zilch(amp1, ncum)
2455	CALL zilch(ad, ncum)
2456
2457	DO k = i + 1, nl + 1
2458	DO il = 1, ncum
2459	IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN
2460	amp1(il) = amp1(il) + m(il, k)
2461	END IF
2462	END DO
2463	END DO
2464
2465	DO k = 1, i
2466	DO j = i + 1, nl + 1
2467	DO il = 1, ncum
2468	IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
2469	amp1(il) = amp1(il) + ment(il, k, j)
2470	END IF
2471	END DO
2472	END DO
2473	END DO
2474
2475	DO k = 1, i - 1
2476	DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
2477	DO il = 1, ncum
2478	IF (i<=inb(il) .AND. j<=inb(il)) THEN
2479	ad(il) = ad(il) + ment(il, j, k)
2480	END IF
2481	END DO
2482	END DO
2483	END DO
2484
2485	DO il = 1, ncum
2486	IF (i<=inb(il)) THEN
2487	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2488	cpinv = 1.0/cpn(il, i)
2489
2490	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
2491	IF (cvflag_grav) THEN
2492	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
2493	ELSE
2494	IF ((0.1dpinvamp1(il))>=delti) iflag(il) = 1 ! vecto
2495	END IF
2496
2497	IF (cvflag_grav) THEN
2498	ft(il, i) = 0.01gravdpinv(amp1(il)(t(il,i+1)-t(il, &
2499	i)+(gz(il,i+1)-gz(il,i))cpinv)-ad(il)(t(il,i)-t(il, &
2500	i-1)+(gz(il,i)-gz(il,i-1))cpinv)) - 0.5sigdlvcp(il, i)(evap( &
2501	il,i)+evap(il,i+1))
2502	rat = cpn(il, i-1)*cpinv
2503	ft(il, i) = ft(il, i) - 0.009gravsigd(mp(il,i+1)t(il,i)*b(il,i) &
2504	-mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2505	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h( &
2506	il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2507	ELSE ! cvflag_grav
2508	ft(il, i) = 0.1dpinv(amp1(il)*(t(il,i+1)-t(il, &
2509	i)+(gz(il,i+1)-gz(il,i))cpinv)-ad(il)(t(il,i)-t(il, &
2510	i-1)+(gz(il,i)-gz(il,i-1))cpinv)) - 0.5sigdlvcp(il, i)(evap( &
2511	il,i)+evap(il,i+1))
2512	rat = cpn(il, i-1)*cpinv
2513	ft(il, i) = ft(il, i) - 0.09sigd(mp(il,i+1)t(il,i)b(il,i)-mp(il &
2514	,i)t(il,i-1)ratb(il,i-1))dpinv
2515	ft(il, i) = ft(il, i) + 0.1dpinvment(il, i, i)*(hp(il,i)-h(il,i)+ &
2516	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2517	END IF ! cvflag_grav
2518
2519
2520	ft(il, i) = ft(il, i) + 0.01sigdwt(il, i)(cl-cpd)water(il, i+1)*( &
2521	t(il,i+1)-t(il,i))dpinvcpinv
2522
2523	IF (cvflag_grav) THEN
2524	fr(il, i) = 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il, &
2525	i))-ad(il)*(rr(il,i)-rr(il,i-1)))
2526	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il, &
2527	i))-ad(il)*(u(il,i)-u(il,i-1)))
2528	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il, &
2529	i))-ad(il)*(v(il,i)-v(il,i-1)))
2530	ELSE ! cvflag_grav
2531	fr(il, i) = 0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il, &
2532	i))-ad(il)*(rr(il,i)-rr(il,i-1)))
2533	fu(il, i) = fu(il, i) + 0.1dpinv(amp1(il)*(u(il,i+1)-u(il, &
2534	i))-ad(il)*(u(il,i)-u(il,i-1)))
2535	fv(il, i) = fv(il, i) + 0.1dpinv(amp1(il)*(v(il,i+1)-v(il, &
2536	i))-ad(il)*(v(il,i)-v(il,i-1)))
2537	END IF ! cvflag_grav
2538
2539	END IF ! i
2540	END DO
2541
2542	! do k=1,ntra
2543	! do il=1,ncum
2544	! if (i.le.inb(il)) then
2545	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2546	! cpinv=1.0/cpn(il,i)
2547	! if (cvflag_grav) then
2548	! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
2549	! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2550	! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2551	! else
2552	! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2553	! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2554	! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2555	! endif
2556	! endif
2557	! enddo
2558	! enddo
2559
2560	DO k = 1, i - 1
2561	DO il = 1, ncum
2562	IF (i<=inb(il)) THEN
2563	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2564	cpinv = 1.0/cpn(il, i)
2565
2566	awat = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
2567	awat = amax1(awat, 0.0)
2568
2569	IF (cvflag_grav) THEN
2570	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k &
2571	,i)-awat-rr(il,i))
2572	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
2573	,i)-u(il,i))
2574	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k &
2575	,i)-v(il,i))
2576	ELSE ! cvflag_grav
2577	fr(il, i) = fr(il, i) + 0.1dpinvment(il, k, i)*(qent(il,k,i)- &
2578	awat-rr(il,i))
2579	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
2580	,i)-u(il,i))
2581	fv(il, i) = fv(il, i) + 0.1dpinvment(il, k, i)*(vent(il,k,i)-v( &
2582	il,i))
2583	END IF ! cvflag_grav
2584
2585	! (saturated updrafts resulting from mixing) ! cld
2586	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat) ! cld
2587	nqcond(il, i) = nqcond(il, i) + 1. ! cld
2588	END IF ! i
2589	END DO
2590	END DO
2591
2592	! do j=1,ntra
2593	! do k=1,i-1
2594	! do il=1,ncum
2595	! if (i.le.inb(il)) then
2596	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2597	! cpinv=1.0/cpn(il,i)
2598	! if (cvflag_grav) then
2599	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2600	! : *(traent(il,k,i,j)-tra(il,i,j))
2601	! else
2602	! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2603	! : *(traent(il,k,i,j)-tra(il,i,j))
2604	! endif
2605	! endif
2606	! enddo
2607	! enddo
2608	! enddo
2609
2610	DO k = i, nl + 1
2611	DO il = 1, ncum
2612	IF (i<=inb(il) .AND. k<=inb(il)) THEN
2613	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2614	cpinv = 1.0/cpn(il, i)
2615
2616	IF (cvflag_grav) THEN
2617	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k &
2618	,i)-rr(il,i))
2619	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k &
2620	,i)-u(il,i))
2621	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k &
2622	,i)-v(il,i))
2623	ELSE ! cvflag_grav
2624	fr(il, i) = fr(il, i) + 0.1dpinvment(il, k, i)*(qent(il,k,i)-rr &
2625	(il,i))
2626	fu(il, i) = fu(il, i) + 0.1dpinvment(il, k, i)*(uent(il,k,i)-u( &
2627	il,i))
2628	fv(il, i) = fv(il, i) + 0.1dpinvment(il, k, i)*(vent(il,k,i)-v( &
2629	il,i))
2630	END IF ! cvflag_grav
2631	END IF ! i and k
2632	END DO
2633	END DO
2634
2635	! do j=1,ntra
2636	! do k=i,nl+1
2637	! do il=1,ncum
2638	! if (i.le.inb(il) .and. k.le.inb(il)) then
2639	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2640	! cpinv=1.0/cpn(il,i)
2641	! if (cvflag_grav) then
2642	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2643	! : *(traent(il,k,i,j)-tra(il,i,j))
2644	! else
2645	! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2646	! : *(traent(il,k,i,j)-tra(il,i,j))
2647	! endif
2648	! endif ! i and k
2649	! enddo
2650	! enddo
2651	! enddo
2652
2653	DO il = 1, ncum
2654	IF (i<=inb(il)) THEN
2655	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2656	cpinv = 1.0/cpn(il, i)
2657
2658	IF (cvflag_grav) THEN
2659	! sb: on ne fait pas encore la correction permettant de mieux
2660	! conserver l'eau:
2661	fr(il, i) = fr(il, i) + 0.5sigd(evap(il,i)+evap(il,i+1)) + &
2662	0.01grav(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)(rp(il, &
2663	i)-rr(il,i-1)))*dpinv
2664
2665	fu(il, i) = fu(il, i) + 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il, &
2666	i))-mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2667	fv(il, i) = fv(il, i) + 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il, &
2668	i))-mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2669	ELSE ! cvflag_grav
2670	fr(il, i) = fr(il, i) + 0.5sigd(evap(il,i)+evap(il,i+1)) + &
2671	0.1(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il,i)-rr(il, &
2672	i-1)))*dpinv
2673	fu(il, i) = fu(il, i) + 0.1(mp(il,i+1)(up(il,i+1)-u(il, &
2674	i))-mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2675	fv(il, i) = fv(il, i) + 0.1(mp(il,i+1)(vp(il,i+1)-v(il, &
2676	i))-mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2677	END IF ! cvflag_grav
2678
2679	END IF ! i
2680	END DO
2681
2682	! sb: interface with the cloud parameterization: ! cld
2683
2684	DO k = i + 1, nl
2685	DO il = 1, ncum
2686	IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld
2687	! (saturated downdrafts resulting from mixing) ! cld
2688	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
2689	nqcond(il, i) = nqcond(il, i) + 1. ! cld
2690	END IF ! cld
2691	END DO ! cld
2692	END DO ! cld
2693
2694	! (particular case: no detraining level is found) ! cld
2695	DO il = 1, ncum ! cld
2696	IF (i<=inb(il) .AND. nent(il,i)==0) THEN ! cld
2697	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
2698	nqcond(il, i) = nqcond(il, i) + 1. ! cld
2699	END IF ! cld
2700	END DO ! cld
2701
2702	DO il = 1, ncum ! cld
2703	IF (i<=inb(il) .AND. nqcond(il,i)/=0.) THEN ! cld
2704	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
2705	END IF ! cld
2706	END DO
2707
2708	! do j=1,ntra
2709	! do il=1,ncum
2710	! if (i.le.inb(il)) then
2711	! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2712	! cpinv=1.0/cpn(il,i)
2713
2714	! if (cvflag_grav) then
2715	! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
2716	! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2717	! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2718	! else
2719	! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2720	! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2721	! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2722	! endif
2723	! endif ! i
2724	! enddo
2725	! enddo
2726
2727	500 END DO
2728
2729
2730	! * move the detrainment at level inb down to level inb-1 *
2731	! * in such a way as to preserve the vertically *
2732	! * integrated enthalpy and water tendencies *
2733
2734	DO il = 1, ncum
2735
2736	ax = 0.1ment(il, inb(il), inb(il))(hp(il,inb(il))-h(il,inb(il))+t(il, &
2737	inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il), &
2738	inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2739	ft(il, inb(il)) = ft(il, inb(il)) - ax
2740	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il &
2741	))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il, &
2742	inb(il))))
2743
2744	bx = 0.1ment(il, inb(il), inb(il))(qent(il,inb(il),inb(il))-rr(il,inb( &
2745	il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2746	fr(il, inb(il)) = fr(il, inb(il)) - bx
2747	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+ &
2748	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2749
2750	cx = 0.1ment(il, inb(il), inb(il))(uent(il,inb(il),inb(il))-u(il,inb(il &
2751	)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2752	fu(il, inb(il)) = fu(il, inb(il)) - cx
2753	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+ &
2754	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2755
2756	dx = 0.1ment(il, inb(il), inb(il))(vent(il,inb(il),inb(il))-v(il,inb(il &
2757	)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2758	fv(il, inb(il)) = fv(il, inb(il)) - dx
2759	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+ &
2760	1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2761
2762	END DO
2763
2764	! do j=1,ntra
2765	! do il=1,ncum
2766	! ex=0.1*ment(il,inb(il),inb(il))
2767	! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2768	! : /(ph(il,inb(il))-ph(il,inb(il)+1))
2769	! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2770	! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2771	! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2772	! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2773	! enddo
2774	! enddo
2775
2776
2777	! * homoginize tendencies below cloud base *
2778
2779
2780	DO il = 1, ncum
2781	asum(il) = 0.0
2782	bsum(il) = 0.0
2783	csum(il) = 0.0
2784	dsum(il) = 0.0
2785	END DO
2786
2787	DO i = 1, nl
2788	DO il = 1, ncum
2789	IF (i<=(icb(il)-1)) THEN
2790	asum(il) = asum(il) + ft(il, i)*(ph(il,i)-ph(il,i+1))
2791	bsum(il) = bsum(il) + fr(il, i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il, &
2792	1)))*(ph(il,i)-ph(il,i+1))
2793	csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
2794	1)))*(ph(il,i)-ph(il,i+1))
2795	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
2796	END IF
2797	END DO
2798	END DO
2799
2800	! !!! do 700 i=1,icb(il)-1
2801	DO i = 1, nl
2802	DO il = 1, ncum
2803	IF (i<=(icb(il)-1)) THEN
2804	ft(il, i) = asum(il)t(il, i)/(th(il,i)dsum(il))
2805	fr(il, i) = bsum(il)/csum(il)
2806	END IF
2807	END DO
2808	END DO
2809
2810
2811	! * reset counter and return *
2812
2813	DO il = 1, ncum
2814	sig(il, nd) = 2.0
2815	END DO
2816
2817
2818	DO i = 1, nd
2819	DO il = 1, ncum
2820	upwd(il, i) = 0.0
2821	dnwd(il, i) = 0.0
2822	END DO
2823	END DO
2824
2825	DO i = 1, nl
2826	DO il = 1, ncum
2827	dnwd0(il, i) = -mp(il, i)
2828	END DO
2829	END DO
2830	DO i = nl + 1, nd
2831	DO il = 1, ncum
2832	dnwd0(il, i) = 0.
2833	END DO
2834	END DO
2835
2836
2837	DO i = 1, nl
2838	DO il = 1, ncum
2839	IF (i>=icb(il) .AND. i<=inb(il)) THEN
2840	upwd(il, i) = 0.0
2841	dnwd(il, i) = 0.0
2842	END IF
2843	END DO
2844	END DO
2845
2846	DO i = 1, nl
2847	DO k = 1, nl
2848	DO il = 1, ncum
2849	up1(il, k, i) = 0.0
2850	dn1(il, k, i) = 0.0
2851	END DO
2852	END DO
2853	END DO
2854
2855	DO i = 1, nl
2856	DO k = i, nl
2857	DO n = 1, i - 1
2858	DO il = 1, ncum
2859	IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
2860	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
2861	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
2862	END IF
2863	END DO
2864	END DO
2865	END DO
2866	END DO
2867
2868	DO i = 2, nl
2869	DO k = i, nl
2870	DO il = 1, ncum
2871	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
2872	! then
2873	IF (i<=inb(il) .AND. k<=inb(il)) THEN
2874	upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
2875	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
2876	END IF
2877	END DO
2878	END DO
2879	END DO
2880
2881
2882	! !!! DO il=1,ncum
2883	! !!! do i=icb(il),inb(il)
2884	! !!!
2885	! !!! upwd(il,i)=0.0
2886	! !!! dnwd(il,i)=0.0
2887	! !!! do k=i,inb(il)
2888	! !!! up1=0.0
2889	! !!! dn1=0.0
2890	! !!! do n=1,i-1
2891	! !!! up1=up1+ment(il,n,k)
2892	! !!! dn1=dn1-ment(il,k,n)
2893	! !!! enddo
2894	! !!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
2895	! !!! dnwd(il,i)=dnwd(il,i)+dn1
2896	! !!! enddo
2897	! !!! enddo
2898	! !!!
2899	! !!! ENDDO
2900
2901	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2902	! determination de la variation de flux ascendant entre
2903	! deux niveau non dilue mike
2904	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2905
2906	DO i = 1, nl
2907	DO il = 1, ncum
2908	mike(il, i) = m(il, i)
2909	END DO
2910	END DO
2911
2912	DO i = nl + 1, nd
2913	DO il = 1, ncum
2914	mike(il, i) = 0.
2915	END DO
2916	END DO
2917
2918	DO i = 1, nd
2919	DO il = 1, ncum
2920	ma(il, i) = 0
2921	END DO
2922	END DO
2923
2924	DO i = 1, nl
2925	DO j = i, nl
2926	DO il = 1, ncum
2927	ma(il, i) = ma(il, i) + m(il, j)
2928	END DO
2929	END DO
2930	END DO
2931
2932	DO i = nl + 1, nd
2933	DO il = 1, ncum
2934	ma(il, i) = 0.
2935	END DO
2936	END DO
2937
2938	DO i = 1, nl
2939	DO il = 1, ncum
2940	IF (i<=(icb(il)-1)) THEN
2941	ma(il, i) = 0
2942	END IF
2943	END DO
2944	END DO
2945
2946	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2947	! icb represente de niveau ou se trouve la
2948	! base du nuage , et inb le top du nuage
2949	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2950
2951	DO i = 1, nd
2952	DO il = 1, ncum
2953	mke(il, i) = upwd(il, i) + dnwd(il, i)
2954	END DO
2955	END DO
2956
2957	DO i = 1, nd
2958	DO il = 1, ncum
2959	rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd*(1.-rr(il, &
2960	i))+rr(il,i)*cpv)
2961	tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp
2962	tps(il, i) = tp(il, i)
2963	END DO
2964	END DO
2965
2966
2967	! * diagnose the in-cloud mixing ratio * ! cld
2968	! * of condensed water * ! cld
2969	! ! cld
2970
2971	DO i = 1, nd ! cld
2972	DO il = 1, ncum ! cld
2973	mac(il, i) = 0.0 ! cld
2974	wa(il, i) = 0.0 ! cld
2975	siga(il, i) = 0.0 ! cld
2976	sax(il, i) = 0.0 ! cld
2977	END DO ! cld
2978	END DO ! cld
2979
2980	DO i = minorig, nl ! cld
2981	DO k = i + 1, nl + 1 ! cld
2982	DO il = 1, ncum ! cld
2983	IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN ! cld
2984	mac(il, i) = mac(il, i) + m(il, k) ! cld
2985	END IF ! cld
2986	END DO ! cld
2987	END DO ! cld
2988	END DO ! cld
2989
2990	DO i = 1, nl ! cld
2991	DO j = 1, i ! cld
2992	DO il = 1, ncum ! cld
2993	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
2994	.AND. j>=icb(il)) THEN ! cld
2995	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
2996	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
2997	END IF ! cld
2998	END DO ! cld
2999	END DO ! cld
3000	END DO ! cld
3001
3002	DO i = 1, nl ! cld
3003	DO il = 1, ncum ! cld
3004	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
3005	.AND. sax(il,i)>0.0) THEN ! cld
3006	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
3007	END IF ! cld
3008	END DO ! cld
3009	END DO ! cld
3010
3011	DO i = 1, nl ! cld
3012	DO il = 1, ncum ! cld
3013	IF (wa(il,i)>0.0) & ! cld
3014	siga(il, i) = mac(il, i)/wa(il, i) & ! cld
3015	rrdtvp(il, i)/p(il, i)/100./delta ! cld
3016	siga(il, i) = min(siga(il,i), 1.0) ! cld
3017	! IM cf. FH
3018	IF (iflag_clw==0) THEN
3019	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
3020	+(1.-siga(il,i))*qcond(il, i) ! cld
3021	ELSE IF (iflag_clw==1) THEN
3022	qcondc(il, i) = qcond(il, i) ! cld
3023	END IF
3024
3025	END DO ! cld
3026	END DO ! cld
3027
3028	RETURN
3029	END SUBROUTINE cv30_yield
3030
3031	! !RomP >>>
3032	SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, &
3033	d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
3034	IMPLICIT NONE
3035
3036	include "cv30param.h"
3037
3038	! inputs:
3039	INTEGER ncum, nd, na, nloc, len
3040	REAL ment(nloc, na, na), sij(nloc, na, na)
3041	REAL clw(nloc, nd), elij(nloc, na, na)
3042	REAL ep(nloc, na)
3043	INTEGER icb(nloc), inb(nloc)
3044	REAL vprecip(nloc, nd+1)
3045	! ouputs:
3046	REAL da(nloc, na), phi(nloc, na, na)
3047	REAL phi2(nloc, na, na)
3048	REAL d1a(nloc, na), dam(nloc, na)
3049	REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
3050	! variables pour tracer dans precip de l'AA et des mel
3051	! local variables:
3052	INTEGER i, j, k
3053	REAL epm(nloc, na, na)
3054
3055	! variables d'Emanuel : du second indice au troisieme
3056	! ---> tab(i,k,j) -> de l origine k a l arrivee j
3057	! ment, sij, elij
3058	! variables personnelles : du troisieme au second indice
3059	! ---> tab(i,j,k) -> de k a j
3060	! phi, phi2
3061
3062	! initialisations
3063	DO j = 1, na
3064	DO i = 1, ncum
3065	da(i, j) = 0.
3066	d1a(i, j) = 0.
3067	dam(i, j) = 0.
3068	eplamm(i, j) = 0.
3069	END DO
3070	END DO
3071	DO k = 1, na
3072	DO j = 1, na
3073	DO i = 1, ncum
3074	epm(i, j, k) = 0.
3075	epmlmmm(i, j, k) = 0.
3076	phi(i, j, k) = 0.
3077	phi2(i, j, k) = 0.
3078	END DO
3079	END DO
3080	END DO
3081
3082	! fraction deau condensee dans les melanges convertie en precip : epm
3083	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3084	DO j = 1, na
3085	DO k = 1, j - 1
3086	DO i = 1, ncum
3087	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3088	! !jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3089	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
3090	! !
3091	epm(i, j, k) = max(epm(i,j,k), 0.0)
3092	END IF
3093	END DO
3094	END DO
3095	END DO
3096
3097	DO j = 1, na
3098	DO k = 1, na
3099	DO i = 1, ncum
3100	IF (k>=icb(i) .AND. k<=inb(i)) THEN
3101	eplamm(i, j) = eplamm(i, j) + ep(i, j)clw(i, j)ment(i, j, k)*(1.- &
3102	sij(i,j,k))
3103	END IF
3104	END DO
3105	END DO
3106	END DO
3107
3108	DO j = 1, na
3109	DO k = 1, j - 1
3110	DO i = 1, ncum
3111	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3112	epmlmmm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
3113	END IF
3114	END DO
3115	END DO
3116	END DO
3117
3118	! matrices pour calculer la tendance des concentrations dans cvltr.F90
3119	DO j = 1, na
3120	DO k = 1, na
3121	DO i = 1, ncum
3122	da(i, j) = da(i, j) + (1.-sij(i,k,j))*ment(i, k, j)
3123	phi(i, j, k) = sij(i, k, j)*ment(i, k, j)
3124	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sij(i,k,j))
3125	END DO
3126	END DO
3127	END DO
3128
3129	DO j = 1, na
3130	DO k = 1, j - 1
3131	DO i = 1, ncum
3132	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, j, k)(1.-ep(i,k))*(1.- &
3133	sij(i,k,j))
3134	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
3135	END DO
3136	END DO
3137	END DO
3138
3139	RETURN
3140	END SUBROUTINE cv30_tracer
3141	! RomP <<<
3142
3143	SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
3144	vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
3145	dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
3146	epmlmmm, eplamm, wdtraina, wdtrainm,epmax_diag, iflag1, precip1, vprecip1, evap1, &
3147	ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
3148	dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
3149	elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1,epmax_diag1) ! epmax_cape
3150	IMPLICIT NONE
3151
3152	include "cv30param.h"
3153
3154	! inputs:
3155	INTEGER len, ncum, nd, ntra, nloc
3156	INTEGER idcum(nloc)
3157	INTEGER iflag(nloc)
3158	INTEGER inb(nloc)
3159	REAL precip(nloc)
3160	REAL vprecip(nloc, nd+1), evap(nloc, nd)
3161	REAL ep(nloc, nd)
3162	REAL sig(nloc, nd), w0(nloc, nd)
3163	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
3164	REAL ftra(nloc, nd, ntra)
3165	REAL ma(nloc, nd)
3166	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
3167	REAL qcondc(nloc, nd)
3168	REAL wd(nloc), cape(nloc)
3169	REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
3170	REAL epmax_diag(nloc) ! epmax_cape
3171	! RomP >>>
3172	REAL phi2(nloc, nd, nd)
3173	REAL d1a(nloc, nd), dam(nloc, nd)
3174	REAL wdtraina(nloc, nd), wdtrainm(nloc, nd)
3175	REAL sij(nloc, nd, nd)
3176	REAL elij(nloc, nd, nd), clw(nloc, nd)
3177	REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd)
3178	! RomP <<<
3179
3180	! outputs:
3181	INTEGER iflag1(len)
3182	INTEGER inb1(len)
3183	REAL precip1(len)
3184	REAL vprecip1(len, nd+1), evap1(len, nd) !<<< RomP
3185	REAL ep1(len, nd) !<<< RomP
3186	REAL sig1(len, nd), w01(len, nd)
3187	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
3188	REAL ftra1(len, nd, ntra)
3189	REAL ma1(len, nd)
3190	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
3191	REAL qcondc1(nloc, nd)
3192	REAL wd1(nloc), cape1(nloc)
3193	REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd)
3194	REAL epmax_diag1(len) ! epmax_cape
3195	! RomP >>>
3196	REAL phi21(len, nd, nd)
3197	REAL d1a1(len, nd), dam1(len, nd)
3198	REAL wdtraina1(len, nd), wdtrainm1(len, nd)
3199	REAL sij1(len, nd, nd)
3200	REAL elij1(len, nd, nd), clw1(len, nd)
3201	REAL epmlmmm1(len, nd, nd), eplamm1(len, nd)
3202	! RomP <<<
3203
3204	! local variables:
3205	INTEGER i, k, j
3206
3207	DO i = 1, ncum
3208	precip1(idcum(i)) = precip(i)
3209	iflag1(idcum(i)) = iflag(i)
3210	wd1(idcum(i)) = wd(i)
3211	inb1(idcum(i)) = inb(i)
3212	cape1(idcum(i)) = cape(i)
3213	epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
3214	END DO
3215
3216	DO k = 1, nl
3217	DO i = 1, ncum
3218	vprecip1(idcum(i), k) = vprecip(i, k)
3219	evap1(idcum(i), k) = evap(i, k) !<<< RomP
3220	sig1(idcum(i), k) = sig(i, k)
3221	w01(idcum(i), k) = w0(i, k)
3222	ft1(idcum(i), k) = ft(i, k)
3223	fq1(idcum(i), k) = fq(i, k)
3224	fu1(idcum(i), k) = fu(i, k)
3225	fv1(idcum(i), k) = fv(i, k)
3226	ma1(idcum(i), k) = ma(i, k)
3227	upwd1(idcum(i), k) = upwd(i, k)
3228	dnwd1(idcum(i), k) = dnwd(i, k)
3229	dnwd01(idcum(i), k) = dnwd0(i, k)
3230	qcondc1(idcum(i), k) = qcondc(i, k)
3231	da1(idcum(i), k) = da(i, k)
3232	mp1(idcum(i), k) = mp(i, k)
3233	! RomP >>>
3234	ep1(idcum(i), k) = ep(i, k)
3235	d1a1(idcum(i), k) = d1a(i, k)
3236	dam1(idcum(i), k) = dam(i, k)
3237	clw1(idcum(i), k) = clw(i, k)
3238	eplamm1(idcum(i), k) = eplamm(i, k)
3239	wdtraina1(idcum(i), k) = wdtraina(i, k)
3240	wdtrainm1(idcum(i), k) = wdtrainm(i, k)
3241	! RomP <<<
3242	END DO
3243	END DO
3244
3245	DO i = 1, ncum
3246	sig1(idcum(i), nd) = sig(i, nd)
3247	END DO
3248
3249
3250	! do 2100 j=1,ntra
3251	! do 2110 k=1,nd ! oct3
3252	! do 2120 i=1,ncum
3253	! ftra1(idcum(i),k,j)=ftra(i,k,j)
3254	! 2120 continue
3255	! 2110 continue
3256	! 2100 continue
3257	DO j = 1, nd
3258	DO k = 1, nd
3259	DO i = 1, ncum
3260	sij1(idcum(i), k, j) = sij(i, k, j)
3261	phi1(idcum(i), k, j) = phi(i, k, j)
3262	phi21(idcum(i), k, j) = phi2(i, k, j)
3263	elij1(idcum(i), k, j) = elij(i, k, j)
3264	epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
3265	END DO
3266	END DO
3267	END DO
3268
3269	RETURN
3270	END SUBROUTINE cv30_uncompress
3271
3272	subroutine cv30_epmax_fn_cape(nloc,ncum,nd &
3273	,cape,ep,hp,icb,inb,clw,nk,t,h,lv &
3274	,epmax_diag)
3275	implicit none
3276
3277	! On fait varier epmax en fn de la cape
3278	! Il faut donc recalculer ep, et hp qui a déjà été calculé et
3279	! qui en dépend
3280	! Toutes les autres variables fn de ep sont calculées plus bas.
3281
3282	#include "cvthermo.h"
3283	#include "cv30param.h"
3284	#include "conema3.h"
3285
3286	! inputs:
3287	integer ncum, nd, nloc
3288	integer icb(nloc), inb(nloc)
3289	real cape(nloc)
3290	real clw(nloc,nd),lv(nloc,nd),t(nloc,nd),h(nloc,nd)
3291	integer nk(nloc)
3292	! inouts:
3293	real ep(nloc,nd)
3294	real hp(nloc,nd)
3295	! outputs ou local
3296	real epmax_diag(nloc)
3297	! locals
3298	integer i,k
3299	real hp_bak(nloc,nd)
3300
3301	! on recalcule ep et hp
3302
3303	if (coef_epmax_cape.gt.1e-12) then
3304	do i=1,ncum
3305	epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
3306	do k=1,nl
3307	ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
3308	ep(i,k)=amax1(ep(i,k),0.0)
3309	ep(i,k)=amin1(ep(i,k),epmax_diag(i))
3310	enddo
3311	enddo
3312
3313	! On recalcule hp:
3314	do k=1,nl
3315	do i=1,ncum
3316	hp_bak(i,k)=hp(i,k)
3317	enddo
3318	enddo
3319	do k=1,nlp
3320	do i=1,ncum
3321	hp(i,k)=h(i,k)
3322	enddo
3323	enddo
3324	do k=minorig+1,nl
3325	do i=1,ncum
3326	if((k.ge.icb(i)).and.(k.le.inb(i)))then
3327	hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
3328	endif
3329	enddo
3330	enddo !do k=minorig+1,n
3331	! write(,) 'cv30_routines 6218: hp(1,20)=',hp(1,20)
3332	do i=1,ncum
3333	do k=1,nl
3334	if (abs(hp_bak(i,k)-hp(i,k)).gt.0.01) then
3335	write(,) 'i,k=',i,k
3336	write(,) 'coef_epmax_cape=',coef_epmax_cape
3337	write(,) 'epmax_diag(i)=',epmax_diag(i)
3338	write(,) 'ep(i,k)=',ep(i,k)
3339	write(,) 'hp(i,k)=',hp(i,k)
3340	write(,) 'hp_bak(i,k)=',hp_bak(i,k)
3341	write(,) 'h(i,k)=',h(i,k)
3342	write(,) 'nk(i)=',nk(i)
3343	write(,) 'h(i,nk(i))=',h(i,nk(i))
3344	write(,) 'lv(i,k)=',lv(i,k)
3345	write(,) 't(i,k)=',t(i,k)
3346	write(,) 'clw(i,k)=',clw(i,k)
3347	write(,) 'cpd,cpv=',cpd,cpv
3348	stop
3349	endif
3350	enddo !do k=1,nl
3351	enddo !do i=1,ncum
3352	endif !if (coef_epmax_cape.gt.1e-12) then
3353
3354	return
3355	end subroutine cv30_epmax_fn_cape
3356
3357

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

Download in other formats:

Original Format