Context Navigation

cv3_routines.F @ 558

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

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

Download in other formats:

Original Format