Context Navigation

cv3_routines.F @ 524

Last change on this file since 524 was 524, checked in by lmdzadmin, 21 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 92.1 KB

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

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

Download in other formats:

Original Format