Context Navigation

cv3_routines.F @ 601

Last change on this file since 601 was 495, checked in by lmdzadmin, 20 years ago
Optimisation de differentes routines, IM, MAF, FH LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 92.1 KB

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

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

Download in other formats:

Original Format