Context Navigation

cv3_routines.F @ 624

Last change on this file since 624 was 619, checked in by lmdzadmin, 20 years ago
Rajout convection Kerry Emanuel pour traceurs- MAF+JYG
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 93.5 KB

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

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

Download in other formats:

Original Format