Context Navigation

cv3_routines.F @ 593

Last change on this file since 593 was 559, checked in by lmdzadmin, 20 years ago
Initialisations diverses YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 92.3 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	do i=1,ncum*nlp
1198	hp(i,1)=h(i,1)
1199	enddo
1200
1201	do 600 k=minorig+1,nl
1202	do 590 i=1,ncum
1203	if((k.ge.icb(i)).and.(k.le.inb(i)))then
1204	hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
1205	endif
1206	590 continue
1207	600 continue
1208
1209	return
1210	end
1211
1212	SUBROUTINE cv3_closure(nloc,ncum,nd,icb,inb
1213	: ,pbase,p,ph,tv,buoy
1214	o ,sig,w0,cape,m)
1215	implicit none
1216
1217	!===================================================================
1218	! --- CLOSURE OF CONVECT3
1219	!
1220	! vectorization: S. Bony
1221	!===================================================================
1222
1223	#include "cvthermo.h"
1224	#include "cvparam3.h"
1225
1226	c input:
1227	integer ncum, nd, nloc
1228	integer icb(nloc), inb(nloc)
1229	real pbase(nloc)
1230	real p(nloc,nd), ph(nloc,nd+1)
1231	real tv(nloc,nd), buoy(nloc,nd)
1232
1233	c input/output:
1234	real sig(nloc,nd), w0(nloc,nd)
1235
1236	c output:
1237	real cape(nloc)
1238	real m(nloc,nd)
1239
1240	c local variables:
1241	integer i, j, k, icbmax
1242	real deltap, fac, w, amu
1243	real dtmin(nloc,nd), sigold(nloc,nd)
1244
1245
1246	c -------------------------------------------------------
1247	c -- Initialization
1248	c -------------------------------------------------------
1249
1250	do k=1,nl
1251	do i=1,ncum
1252	m(i,k)=0.0
1253	enddo
1254	enddo
1255
1256	c -------------------------------------------------------
1257	c -- Reset sig(i) and w0(i) for i>inb and i<icb
1258	c -------------------------------------------------------
1259
1260	c update sig and w0 above LNB:
1261
1262	do 100 k=1,nl-1
1263	do 110 i=1,ncum
1264	if ((inb(i).lt.(nl-1)).and.(k.ge.(inb(i)+1)))then
1265	sig(i,k)=beta*sig(i,k)
1266	: +2.alphabuoy(i,inb(i))*ABS(buoy(i,inb(i)))
1267	sig(i,k)=AMAX1(sig(i,k),0.0)
1268	w0(i,k)=beta*w0(i,k)
1269	endif
1270	110 continue
1271	100 continue
1272
1273	c compute icbmax:
1274
1275	icbmax=2
1276	do 200 i=1,ncum
1277	icbmax=MAX(icbmax,icb(i))
1278	200 continue
1279
1280	c update sig and w0 below cloud base:
1281
1282	do 300 k=1,icbmax
1283	do 310 i=1,ncum
1284	if (k.le.icb(i))then
1285	sig(i,k)=betasig(i,k)-2.alphabuoy(i,icb(i))buoy(i,icb(i))
1286	sig(i,k)=amax1(sig(i,k),0.0)
1287	w0(i,k)=beta*w0(i,k)
1288	endif
1289	310 continue
1290	300 continue
1291
1292	c! if(inb.lt.(nl-1))then
1293	c! do 85 i=inb+1,nl-1
1294	c! sig(i)=betasig(i)+2.alphabuoy(inb)
1295	c! 1 abs(buoy(inb))
1296	c! sig(i)=amax1(sig(i),0.0)
1297	c! w0(i)=beta*w0(i)
1298	c! 85 continue
1299	c! end if
1300
1301	c! do 87 i=1,icb
1302	c! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
1303	c! sig(i)=amax1(sig(i),0.0)
1304	c! w0(i)=beta*w0(i)
1305	c! 87 continue
1306
1307	c -------------------------------------------------------------
1308	c -- Reset fractional areas of updrafts and w0 at initial time
1309	c -- and after 10 time steps of no convection
1310	c -------------------------------------------------------------
1311
1312	do 400 k=1,nl-1
1313	do 410 i=1,ncum
1314	if (sig(i,nd).lt.1.5.or.sig(i,nd).gt.12.0)then
1315	sig(i,k)=0.0
1316	w0(i,k)=0.0
1317	endif
1318	410 continue
1319	400 continue
1320
1321	c -------------------------------------------------------------
1322	c -- Calculate convective available potential energy (cape),
1323	c -- vertical velocity (w), fractional area covered by
1324	c -- undilute updraft (sig), and updraft mass flux (m)
1325	c -------------------------------------------------------------
1326
1327	do 500 i=1,ncum
1328	cape(i)=0.0
1329	500 continue
1330
1331	c compute dtmin (minimum buoyancy between ICB and given level k):
1332
1333	do i=1,ncum
1334	do k=1,nl
1335	dtmin(i,k)=100.0
1336	enddo
1337	enddo
1338
1339	do 550 i=1,ncum
1340	do 560 k=1,nl
1341	do 570 j=minorig,nl
1342	if ( (k.ge.(icb(i)+1)).and.(k.le.inb(i)).and.
1343	: (j.ge.icb(i)).and.(j.le.(k-1)) )then
1344	dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1345	endif
1346	570 continue
1347	560 continue
1348	550 continue
1349
1350	c the interval on which cape is computed starts at pbase :
1351
1352	do 600 k=1,nl
1353	do 610 i=1,ncum
1354
1355	if ((k.ge.(icb(i)+1)).and.(k.le.inb(i))) then
1356
1357	deltap = MIN(pbase(i),ph(i,k-1))-MIN(pbase(i),ph(i,k))
1358	cape(i)=cape(i)+rrdbuoy(i,k-1)deltap/p(i,k-1)
1359	cape(i)=AMAX1(0.0,cape(i))
1360	sigold(i,k)=sig(i,k)
1361
1362	c dtmin(i,k)=100.0
1363	c do 97 j=icb(i),k-1 ! mauvaise vectorisation
1364	c dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1365	c 97 continue
1366
1367	sig(i,k)=betasig(i,k)+alphadtmin(i,k)*ABS(dtmin(i,k))
1368	sig(i,k)=amax1(sig(i,k),0.0)
1369	sig(i,k)=amin1(sig(i,k),0.01)
1370	fac=AMIN1(((dtcrit-dtmin(i,k))/dtcrit),1.0)
1371	w=(1.-beta)facSQRT(cape(i))+beta*w0(i,k)
1372	amu=0.5(sig(i,k)+sigold(i,k))w
1373	m(i,k)=amu0.007p(i,k)*(ph(i,k)-ph(i,k+1))/tv(i,k)
1374	w0(i,k)=w
1375	endif
1376
1377	610 continue
1378	600 continue
1379
1380	do 700 i=1,ncum
1381	w0(i,icb(i))=0.5*w0(i,icb(i)+1)
1382	m(i,icb(i))=0.5*m(i,icb(i)+1)
1383	: *(ph(i,icb(i))-ph(i,icb(i)+1))
1384	: /(ph(i,icb(i)+1)-ph(i,icb(i)+2))
1385	sig(i,icb(i))=sig(i,icb(i)+1)
1386	sig(i,icb(i)-1)=sig(i,icb(i))
1387	700 continue
1388
1389
1390	c! cape=0.0
1391	c! do 98 i=icb+1,inb
1392	c! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1393	c! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
1394	c! dcape=rrdbuoy(i-1)deltap/p(i-1)
1395	c! dlnp=deltap/p(i-1)
1396	c! cape=amax1(0.0,cape)
1397	c! sigold=sig(i)
1398
1399	c! dtmin=100.0
1400	c! do 97 j=icb,i-1
1401	c! dtmin=amin1(dtmin,buoy(j))
1402	c! 97 continue
1403
1404	c! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
1405	c! sig(i)=amax1(sig(i),0.0)
1406	c! sig(i)=amin1(sig(i),0.01)
1407	c! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1408	c! w=(1.-beta)facsqrt(cape)+beta*w0(i)
1409	c! amu=0.5(sig(i)+sigold)w
1410	c! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
1411	c! w0(i)=w
1412	c! 98 continue
1413	c! w0(icb)=0.5*w0(icb+1)
1414	c! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1415	c! sig(icb)=sig(icb+1)
1416	c! sig(icb-1)=sig(icb)
1417
1418	return
1419	end
1420
1421	SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb
1422	: ,ph,t,rr,rs,u,v,tra,h,lv,qnk
1423	: ,hp,tv,tvp,ep,clw,m,sig
1424	: ,ment,qent,uent,vent,sij,elij,ments,qents,traent)
1425	implicit none
1426
1427	!---------------------------------------------------------------------
1428	! a faire:
1429	! - changer rr(il,1) -> qnk(il)
1430	! - vectorisation de la partie normalisation des flux (do 789...)
1431	!---------------------------------------------------------------------
1432
1433	#include "cvthermo.h"
1434	#include "cvparam3.h"
1435
1436	c inputs:
1437	integer ncum, nd, na, ntra, nloc
1438	integer icb(nloc), inb(nloc), nk(nloc)
1439	real sig(nloc,nd)
1440	real qnk(nloc)
1441	real ph(nloc,nd+1)
1442	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
1443	real u(nloc,nd), v(nloc,nd)
1444	real tra(nloc,nd,ntra) ! input of convect3
1445	real lv(nloc,na), h(nloc,na), hp(nloc,na)
1446	real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
1447	real m(nloc,na) ! input of convect3
1448
1449	c outputs:
1450	real ment(nloc,na,na), qent(nloc,na,na)
1451	real uent(nloc,na,na), vent(nloc,na,na)
1452	real sij(nloc,na,na), elij(nloc,na,na)
1453	real traent(nloc,nd,nd,ntra)
1454	real ments(nloc,nd,nd), qents(nloc,nd,nd)
1455	real sigij(nloc,nd,nd)
1456
1457	c local variables:
1458	integer i, j, k, il, im, jm
1459	integer num1, num2
1460	integer nent(nloc,na)
1461	real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1462	real alt, smid, sjmin, sjmax, delp, delm
1463	real asij(nloc), smax(nloc), scrit(nloc)
1464	real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd)
1465	real wgh
1466	real zm(nloc,na)
1467	logical lwork(nloc)
1468
1469	c=====================================================================
1470	c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1471	c=====================================================================
1472
1473	c ori do 360 i=1,ncum*nlp
1474	do 361 j=1,nl
1475	do 360 i=1,ncum
1476	nent(i,j)=0
1477	c in convect3, m is computed in cv3_closure
1478	c ori m(i,1)=0.0
1479	360 continue
1480	361 continue
1481
1482	c ori do 400 k=1,nlp
1483	c ori do 390 j=1,nlp
1484	do 400 j=1,nl
1485	do 390 k=1,nl
1486	do 385 i=1,ncum
1487	qent(i,k,j)=rr(i,j)
1488	uent(i,k,j)=u(i,j)
1489	vent(i,k,j)=v(i,j)
1490	elij(i,k,j)=0.0
1491	cym ment(i,k,j)=0.0
1492	cym sij(i,k,j)=0.0
1493	385 continue
1494	390 continue
1495	400 continue
1496
1497	cym
1498	ment(1:ncum,1:nd,1:nd)=0.0
1499	sij(1:ncum,1:nd,1:nd)=0.0
1500
1501	c do k=1,ntra
1502	c do j=1,nd ! instead nlp
1503	c do i=1,nd ! instead nlp
1504	c do il=1,ncum
1505	c traent(il,i,j,k)=tra(il,j,k)
1506	c enddo
1507	c enddo
1508	c enddo
1509	c enddo
1510	zm(:,:)=0.
1511
1512	c=====================================================================
1513	c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1514	c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1515	c --- FRACTION (sij)
1516	c=====================================================================
1517
1518	do 750 i=minorig+1, nl
1519
1520	do 710 j=minorig,nl
1521	do 700 il=1,ncum
1522	if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1523	: (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1524
1525	rti=rr(il,1)-ep(il,i)*clw(il,i)
1526	bf2=1.+lv(il,j)lv(il,j)rs(il,j)/(rrvt(il,j)t(il,j)*cpd)
1527	anum=h(il,j)-hp(il,i)+(cpv-cpd)t(il,j)(rti-rr(il,j))
1528	denom=h(il,i)-hp(il,i)+(cpd-cpv)(rr(il,i)-rti)t(il,j)
1529	dei=denom
1530	if(abs(dei).lt.0.01)dei=0.01
1531	sij(il,i,j)=anum/dei
1532	sij(il,i,i)=1.0
1533	altem=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti-rs(il,j)
1534	altem=altem/bf2
1535	cwat=clw(il,j)*(1.-ep(il,j))
1536	stemp=sij(il,i,j)
1537	if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat)
1538	: .and.j.gt.i)then
1539	anum=anum-lv(il,j)(rti-rs(il,j)-cwatbf2)
1540	denom=denom+lv(il,j)*(rr(il,i)-rti)
1541	if(abs(denom).lt.0.01)denom=0.01
1542	sij(il,i,j)=anum/denom
1543	altem=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti-rs(il,j)
1544	altem=altem-(bf2-1.)*cwat
1545	end if
1546	if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then
1547	qent(il,i,j)=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti
1548	uent(il,i,j)=sij(il,i,j)u(il,i)+(1.-sij(il,i,j))u(il,nk(il))
1549	vent(il,i,j)=sij(il,i,j)v(il,i)+(1.-sij(il,i,j))v(il,nk(il))
1550	c!!! do k=1,ntra
1551	c!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1552	c!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1553	c!!! end do
1554	elij(il,i,j)=altem
1555	elij(il,i,j)=amax1(0.0,elij(il,i,j))
1556	ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
1557	nent(il,i)=nent(il,i)+1
1558	end if
1559	sij(il,i,j)=amax1(0.0,sij(il,i,j))
1560	sij(il,i,j)=amin1(1.0,sij(il,i,j))
1561	endif ! new
1562	700 continue
1563	710 continue
1564
1565	c do k=1,ntra
1566	c do j=minorig,nl
1567	c do il=1,ncum
1568	c if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1569	c : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1570	c traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1571	c : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1572	c endif
1573	c enddo
1574	c enddo
1575	c enddo
1576
1577	c
1578	c * if no air can entrain at level i assume that updraft detrains *
1579	c * at that level and calculate detrained air flux and properties *
1580	c
1581
1582	c@ do 170 i=icb(il),inb(il)
1583
1584	do 740 il=1,ncum
1585	if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then
1586	c@ if(nent(il,i).eq.0)then
1587	ment(il,i,i)=m(il,i)
1588	qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)
1589	uent(il,i,i)=u(il,nk(il))
1590	vent(il,i,i)=v(il,nk(il))
1591	elij(il,i,i)=clw(il,i)
1592	cMAF sij(il,i,i)=1.0
1593	sij(il,i,i)=0.0
1594	end if
1595	740 continue
1596	750 continue
1597
1598	c do j=1,ntra
1599	c do i=minorig+1,nl
1600	c do il=1,ncum
1601	c if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1602	c traent(il,i,i,j)=tra(il,nk(il),j)
1603	c endif
1604	c enddo
1605	c enddo
1606	c enddo
1607
1608	do 100 j=minorig,nl
1609	do 101 i=minorig,nl
1610	do 102 il=1,ncum
1611	if ((j.ge.(icb(il)-1)).and.(j.le.inb(il))
1612	: .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
1613	sigij(il,i,j)=sij(il,i,j)
1614	endif
1615	102 continue
1616	101 continue
1617	100 continue
1618	c@ enddo
1619
1620	c@170 continue
1621
1622	c=====================================================================
1623	c --- NORMALIZE ENTRAINED AIR MASS FLUXES
1624	c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
1625	c=====================================================================
1626
1627	cym call zilch(asum,ncum*nd)
1628	cym call zilch(bsum,ncum*nd)
1629	cym call zilch(csum,ncum*nd)
1630	call zilch(asum,nloc*nd)
1631	call zilch(csum,nloc*nd)
1632	call zilch(csum,nloc*nd)
1633
1634	do il=1,ncum
1635	lwork(il) = .FALSE.
1636	enddo
1637
1638	DO 789 i=minorig+1,nl
1639
1640	num1=0
1641	do il=1,ncum
1642	if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
1643	enddo
1644	if (num1.le.0) goto 789
1645
1646
1647	do 781 il=1,ncum
1648	if ( i.ge.icb(il) .and. i.le.inb(il) ) then
1649	lwork(il)=(nent(il,i).ne.0)
1650	qp=rr(il,1)-ep(il,i)*clw(il,i)
1651	anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
1652	: +(cpv-cpd)t(il,i)(qp-rr(il,i))
1653	denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp)
1654	: +(cpd-cpv)t(il,i)(rr(il,i)-qp)
1655	if(abs(denom).lt.0.01)denom=0.01
1656	scrit(il)=anum/denom
1657	alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
1658	if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0
1659	smax(il)=0.0
1660	asij(il)=0.0
1661	endif
1662	781 continue
1663
1664	do 175 j=nl,minorig,-1
1665
1666	num2=0
1667	do il=1,ncum
1668	if ( i.ge.icb(il) .and. i.le.inb(il) .and.
1669	: j.ge.(icb(il)-1) .and. j.le.inb(il)
1670	: .and. lwork(il) ) num2=num2+1
1671	enddo
1672	if (num2.le.0) goto 175
1673
1674	do 782 il=1,ncum
1675	if ( i.ge.icb(il) .and. i.le.inb(il) .and.
1676	: j.ge.(icb(il)-1) .and. j.le.inb(il)
1677	: .and. lwork(il) ) then
1678
1679	if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then
1680	wgh=1.0
1681	if(j.gt.i)then
1682	sjmax=amax1(sij(il,i,j+1),smax(il))
1683	sjmax=amin1(sjmax,scrit(il))
1684	smax(il)=amax1(sij(il,i,j),smax(il))
1685	sjmin=amax1(sij(il,i,j-1),smax(il))
1686	sjmin=amin1(sjmin,scrit(il))
1687	if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
1688	smid=amin1(sij(il,i,j),scrit(il))
1689	else
1690	sjmax=amax1(sij(il,i,j+1),scrit(il))
1691	smid=amax1(sij(il,i,j),scrit(il))
1692	sjmin=0.0
1693	if(j.gt.1)sjmin=sij(il,i,j-1)
1694	sjmin=amax1(sjmin,scrit(il))
1695	endif
1696	delp=abs(sjmax-smid)
1697	delm=abs(sjmin-smid)
1698	asij(il)=asij(il)+wgh*(delp+delm)
1699	ment(il,i,j)=ment(il,i,j)(delp+delm)wgh
1700	endif
1701	endif
1702	782 continue
1703
1704	175 continue
1705
1706	do il=1,ncum
1707	if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
1708	asij(il)=amax1(1.0e-16,asij(il))
1709	asij(il)=1.0/asij(il)
1710	asum(il,i)=0.0
1711	bsum(il,i)=0.0
1712	csum(il,i)=0.0
1713	endif
1714	enddo
1715
1716	do 180 j=minorig,nl
1717	do il=1,ncum
1718	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1719	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1720	ment(il,i,j)=ment(il,i,j)*asij(il)
1721	endif
1722	enddo
1723	180 continue
1724
1725	do 190 j=minorig,nl
1726	do il=1,ncum
1727	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1728	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1729	asum(il,i)=asum(il,i)+ment(il,i,j)
1730	ment(il,i,j)=ment(il,i,j)*sig(il,j)
1731	bsum(il,i)=bsum(il,i)+ment(il,i,j)
1732	endif
1733	enddo
1734	190 continue
1735
1736	do il=1,ncum
1737	if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
1738	bsum(il,i)=amax1(bsum(il,i),1.0e-16)
1739	bsum(il,i)=1.0/bsum(il,i)
1740	endif
1741	enddo
1742
1743	do 195 j=minorig,nl
1744	do il=1,ncum
1745	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1746	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1747	ment(il,i,j)=ment(il,i,j)asum(il,i)bsum(il,i)
1748	endif
1749	enddo
1750	195 continue
1751
1752	do 197 j=minorig,nl
1753	do il=1,ncum
1754	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1755	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1756	csum(il,i)=csum(il,i)+ment(il,i,j)
1757	endif
1758	enddo
1759	197 continue
1760
1761	do il=1,ncum
1762	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1763	: .and. csum(il,i).lt.m(il,i) ) then
1764	nent(il,i)=0
1765	ment(il,i,i)=m(il,i)
1766	qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)
1767	uent(il,i,i)=u(il,nk(il))
1768	vent(il,i,i)=v(il,nk(il))
1769	elij(il,i,i)=clw(il,i)
1770	cMAF sij(il,i,i)=1.0
1771	sij(il,i,i)=0.0
1772	endif
1773	enddo ! il
1774
1775	c do j=1,ntra
1776	c do il=1,ncum
1777	c if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1778	c : .and. csum(il,i).lt.m(il,i) ) then
1779	c traent(il,i,i,j)=tra(il,nk(il),j)
1780	c endif
1781	c enddo
1782	c enddo
1783	789 continue
1784	c
1785	c MAF: renormalisation de MENT
1786	do jm=1,nd
1787	do im=1,nd
1788	do il=1,ncum
1789	zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)
1790	end do
1791	end do
1792	end do
1793	c
1794	do jm=1,nd
1795	do im=1,nd
1796	do il=1,ncum
1797	if(zm(il,im).ne.0.) then
1798	ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)
1799	endif
1800	end do
1801	end do
1802	end do
1803	c
1804	do jm=1,nd
1805	do im=1,nd
1806	do 999 il=1,ncum
1807	qents(il,im,jm)=qent(il,im,jm)
1808	ments(il,im,jm)=ment(il,im,jm)
1809	999 continue
1810	enddo
1811	enddo
1812
1813	return
1814	end
1815
1816
1817	SUBROUTINE cv3_unsat(nloc,ncum,nd,na,ntra,icb,inb
1818	: ,t,rr,rs,gz,u,v,tra,p,ph
1819	: ,th,tv,lv,cpn,ep,sigp,clw
1820	: ,m,ment,elij,delt,plcl
1821	: ,mp,rp,up,vp,trap,wt,water,evap,b)
1822	implicit none
1823
1824
1825	#include "cvthermo.h"
1826	#include "cvparam3.h"
1827	#include "cvflag.h"
1828
1829	c inputs:
1830	integer ncum, nd, na, ntra, nloc
1831	integer icb(nloc), inb(nloc)
1832	real delt, plcl(nloc)
1833	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
1834	real u(nloc,nd), v(nloc,nd)
1835	real tra(nloc,nd,ntra)
1836	real p(nloc,nd), ph(nloc,nd+1)
1837	real th(nloc,na), gz(nloc,na)
1838	real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na)
1839	real cpn(nloc,na), tv(nloc,na)
1840	real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
1841
1842	c outputs:
1843	real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
1844	real water(nloc,na), evap(nloc,na), wt(nloc,na)
1845	real trap(nloc,na,ntra)
1846	real b(nloc,na)
1847
1848	c local variables
1849	integer i,j,k,il,num1
1850	real tinv, delti
1851	real awat, afac, afac1, afac2, bfac
1852	real pr1, pr2, sigt, b6, c6, revap, tevap, delth
1853	real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1854	real ampmax
1855	real lvcp(nloc,na)
1856	real wdtrain(nloc)
1857	logical lwork(nloc)
1858
1859
1860	c------------------------------------------------------
1861
1862	delti = 1./delt
1863	tinv=1./3.
1864
1865	do i=1,nl
1866	do il=1,ncum
1867	mp(il,i)=0.0
1868	rp(il,i)=rr(il,i)
1869	up(il,i)=u(il,i)
1870	vp(il,i)=v(il,i)
1871	wt(il,i)=0.001
1872	water(il,i)=0.0
1873	evap(il,i)=0.0
1874	b(il,i)=0.0
1875	lvcp(il,i)=lv(il,i)/cpn(il,i)
1876	enddo
1877	enddo
1878
1879	c do k=1,ntra
1880	c do i=1,nd
1881	c do il=1,ncum
1882	c trap(il,i,k)=tra(il,i,k)
1883	c enddo
1884	c enddo
1885	c enddo
1886
1887	c
1888	c * check whether ep(inb)=0, if so, skip precipitating *
1889	c * downdraft calculation *
1890	c
1891
1892	do il=1,ncum
1893	lwork(il)=.TRUE.
1894	if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
1895	enddo
1896
1897	call zilch(wdtrain,ncum)
1898
1899	DO 400 i=nl+1,1,-1
1900
1901	num1=0
1902	do il=1,ncum
1903	if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1
1904	enddo
1905	if (num1.le.0) goto 400
1906
1907	c
1908	c * integrate liquid water equation to find condensed water *
1909	c * and condensed water flux *
1910	c
1911
1912	c
1913	c * begin downdraft loop *
1914	c
1915
1916	c
1917	c * calculate detrained precipitation *
1918	c
1919	do il=1,ncum
1920	if (i.le.inb(il) .and. lwork(il)) then
1921	if (cvflag_grav) then
1922	wdtrain(il)=gravep(il,i)m(il,i)*clw(il,i)
1923	else
1924	wdtrain(il)=10.0ep(il,i)m(il,i)*clw(il,i)
1925	endif
1926	endif
1927	enddo
1928
1929	if(i.gt.1)then
1930	do 320 j=1,i-1
1931	do il=1,ncum
1932	if (i.le.inb(il) .and. lwork(il)) then
1933	awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
1934	awat=amax1(awat,0.0)
1935	if (cvflag_grav) then
1936	wdtrain(il)=wdtrain(il)+gravawatment(il,j,i)
1937	else
1938	wdtrain(il)=wdtrain(il)+10.0awatment(il,j,i)
1939	endif
1940	endif
1941	enddo
1942	320 continue
1943	endif
1944
1945	c
1946	c * find rain water and evaporation using provisional *
1947	c * estimates of rp(i)and rp(i-1) *
1948	c
1949
1950	do 999 il=1,ncum
1951
1952	if (i.le.inb(il) .and. lwork(il)) then
1953
1954	wt(il,i)=45.0
1955
1956	if(i.lt.inb(il))then
1957	rp(il,i)=rp(il,i+1)
1958	: +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
1959	rp(il,i)=0.5*(rp(il,i)+rr(il,i))
1960	endif
1961	rp(il,i)=amax1(rp(il,i),0.0)
1962	rp(il,i)=amin1(rp(il,i),rs(il,i))
1963	rp(il,inb(il))=rr(il,inb(il))
1964
1965	if(i.eq.1)then
1966	afac=p(il,1)(rs(il,1)-rp(il,1))/(1.0e4+2000.0p(il,1)*rs(il,1))
1967	else
1968	rp(il,i-1)=rp(il,i)
1969	: +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i)
1970	rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
1971	rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
1972	rp(il,i-1)=amax1(rp(il,i-1),0.0)
1973	afac1=p(il,i)(rs(il,i)-rp(il,i))/(1.0e4+2000.0p(il,i)*rs(il,i))
1974	afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
1975	: /(1.0e4+2000.0p(il,i-1)rs(il,i-1))
1976	afac=0.5*(afac1+afac2)
1977	endif
1978	if(i.eq.inb(il))afac=0.0
1979	afac=amax1(afac,0.0)
1980	bfac=1./(sigd*wt(il,i))
1981	c
1982	cjyg1
1983	ccc sigt=1.0
1984	ccc if(i.ge.icb)sigt=sigp(i)
1985	c prise en compte de la variation progressive de sigt dans
1986	c les couches icb et icb-1:
1987	c pour plcl<ph(i+1), pr1=0 & pr2=1
1988	c pour plcl>ph(i), pr1=1 & pr2=0
1989	c pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
1990	c sur le nuage, et pr2 est la proportion sous la base du
1991	c nuage.
1992	pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
1993	pr1=max(0.,min(1.,pr1))
1994	pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
1995	pr2=max(0.,min(1.,pr2))
1996	sigt=sigp(il,i)*pr1+pr2
1997	cjyg2
1998	c
1999	b6=bfac50.sigd(ph(il,i)-ph(il,i+1))sigt*afac
2000	c6=water(il,i+1)+bfac*wdtrain(il)
2001	: -50.sigdbfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2002	if(c6.gt.0.0)then
2003	revap=0.5(-b6+sqrt(b6b6+4.*c6))
2004	evap(il,i)=sigtafacrevap
2005	water(il,i)=revap*revap
2006	else
2007	evap(il,i)=-evap(il,i+1)
2008	: +0.02(wdtrain(il)+sigdwt(il,i)*water(il,i+1))
2009	: /(sigd*(ph(il,i)-ph(il,i+1)))
2010	end if
2011	c
2012	c * calculate precipitating downdraft mass flux under *
2013	c * hydrostatic approximation *
2014	c
2015	if (i.ne.1) then
2016
2017	tevap=amax1(0.0,evap(il,i))
2018	delth=amax1(0.001,(th(il,i)-th(il,i-1)))
2019	if (cvflag_grav) then
2020	mp(il,i)=100.ginvlvcp(il,i)sigdtevap
2021	: *(p(il,i-1)-p(il,i))/delth
2022	else
2023	mp(il,i)=10.lvcp(il,i)sigdtevap(p(il,i-1)-p(il,i))/delth
2024	endif
2025	c
2026	c * if hydrostatic assumption fails, *
2027	c * solve cubic difference equation for downdraft theta *
2028	c * and mass flux from two simultaneous differential eqns *
2029	c
2030	amfac=sigdsigd70.0ph(il,i)(p(il,i-1)-p(il,i))
2031	: (th(il,i)-th(il,i-1))/(tv(il,i)th(il,i))
2032	amp2=abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2033	if(amp2.gt.(0.1*amfac))then
2034	xf=100.0sigdsigdsigd(ph(il,i)-ph(il,i+1))
2035	tf=b(il,i)-5.0(th(il,i)-th(il,i-1))t(il,i)
2036	: /(lvcp(il,i)sigdth(il,i))
2037	af=xftf+mp(il,i+1)mp(il,i+1)*tinv
2038	bf=2.(tinvmp(il,i+1))*3+tinvmp(il,i+1)xftf
2039	: +50.(p(il,i-1)-p(il,i))xf*tevap
2040	fac2=1.0
2041	if(bf.lt.0.0)fac2=-1.0
2042	bf=abs(bf)
2043	ur=0.25bfbf-afafaftinvtinv*tinv
2044	if(ur.ge.0.0)then
2045	sru=sqrt(ur)
2046	fac=1.0
2047	if((0.5*bf-sru).lt.0.0)fac=-1.0
2048	mp(il,i)=mp(il,i+1)tinv+(0.5bf+sru)**tinv
2049	: +fac(abs(0.5bf-sru))**tinv
2050	else
2051	d=atan(2.*sqrt(-ur)/(bf+1.0e-28))
2052	if(fac2.lt.0.0)d=3.14159-d
2053	mp(il,i)=mp(il,i+1)tinv+2.sqrt(aftinv)cos(d*tinv)
2054	endif
2055	mp(il,i)=amax1(0.0,mp(il,i))
2056
2057	if (cvflag_grav) then
2058	Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2059	C il faut diviser par (mp(il,i)sigdgrav) et non par (mp(il,i)+sigd*0.1).
2060	C Et il faut bien revoir les facteurs 100.
2061	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2062	2 /(mp(il,i)+sigd*0.1)
2063	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
2064	else
2065	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2066	2 /(mp(il,i)+sigd*0.1)
2067	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
2068	endif
2069	b(il,i-1)=amax1(b(il,i-1),0.0)
2070	endif
2071	c
2072	c * limit magnitude of mp(i) to meet cfl condition *
2073	c
2074	ampmax=2.0(ph(il,i)-ph(il,i+1))delti
2075	amp2=2.0(ph(il,i-1)-ph(il,i))delti
2076	ampmax=amin1(ampmax,amp2)
2077	mp(il,i)=amin1(mp(il,i),ampmax)
2078	c
2079	c * force mp to decrease linearly to zero *
2080	c * between cloud base and the surface *
2081	c
2082	if(p(il,i).gt.p(il,icb(il)))then
2083	mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2084	endif
2085
2086	360 continue
2087	endif ! i.eq.1
2088	c
2089	c * find mixing ratio of precipitating downdraft *
2090	c
2091
2092	if (i.ne.inb(il)) then
2093
2094	rp(il,i)=rr(il,i)
2095
2096	if(mp(il,i).gt.mp(il,i+1))then
2097
2098	if (cvflag_grav) then
2099	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2100	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
2101	: *(evap(il,i+1)+evap(il,i))
2102	else
2103	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2104	: +5.sigd(ph(il,i)-ph(il,i+1))
2105	: *(evap(il,i+1)+evap(il,i))
2106	endif
2107	rp(il,i)=rp(il,i)/mp(il,i)
2108	up(il,i)=up(il,i+1)mp(il,i+1)+u(il,i)(mp(il,i)-mp(il,i+1))
2109	up(il,i)=up(il,i)/mp(il,i)
2110	vp(il,i)=vp(il,i+1)mp(il,i+1)+v(il,i)(mp(il,i)-mp(il,i+1))
2111	vp(il,i)=vp(il,i)/mp(il,i)
2112
2113	c do j=1,ntra
2114	c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2115	ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2116	c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2117	c trap(il,i,j)=trap(il,i,j)/mp(il,i)
2118	c end do
2119
2120	else
2121
2122	if(mp(il,i+1).gt.1.0e-16)then
2123	if (cvflag_grav) then
2124	rp(il,i)=rp(il,i+1)
2125	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
2126	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2127	else
2128	rp(il,i)=rp(il,i+1)
2129	: +5.sigd(ph(il,i)-ph(il,i+1))
2130	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2131	endif
2132	up(il,i)=up(il,i+1)
2133	vp(il,i)=vp(il,i+1)
2134
2135	c do j=1,ntra
2136	c trap(il,i,j)=trap(il,i+1,j)
2137	c end do
2138
2139	endif
2140	endif
2141	rp(il,i)=amin1(rp(il,i),rs(il,i))
2142	rp(il,i)=amax1(rp(il,i),0.0)
2143
2144	endif
2145	endif
2146	999 continue
2147
2148	400 continue
2149
2150	return
2151	end
2152
2153	SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
2154	: ,icb,inb,delt
2155	: ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th
2156	: ,ep,clw,m,tp,mp,rp,up,vp,trap
2157	: ,wt,water,evap,b
2158	: ,ment,qent,uent,vent,nent,elij,traent,sig
2159	: ,tv,tvp
2160	: ,iflag,precip,ft,fr,fu,fv,ftra
2161	: ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
2162	implicit none
2163
2164	#include "cvthermo.h"
2165	#include "cvparam3.h"
2166	#include "cvflag.h"
2167	#include "conema3.h"
2168
2169	c inputs:
2170	integer ncum,nd,na,ntra,nloc
2171	integer icb(nloc), inb(nloc)
2172	real delt
2173	real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
2174	real tra(nloc,nd,ntra), sig(nloc,nd)
2175	real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
2176	real th(nloc,na), p(nloc,nd), tp(nloc,na)
2177	real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
2178	real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
2179	real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
2180	real water(nloc,na), evap(nloc,na), b(nloc,na)
2181	real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
2182	real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2183	real traent(nloc,na,na,ntra)
2184	real tv(nloc,nd), tvp(nloc,nd)
2185
2186	c input/output:
2187	integer iflag(nloc)
2188
2189	c outputs:
2190	real precip(nloc)
2191	real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
2192	real ftra(nloc,nd,ntra)
2193	real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
2194	real dnwd0(nloc,nd), mike(nloc,nd)
2195	real tls(nloc,nd), tps(nloc,nd)
2196	real qcondc(nloc,nd) ! cld
2197	real wd(nloc) ! gust
2198
2199	c local variables:
2200	integer i,k,il,n,j,num1
2201	real rat, awat, delti
2202	real ax, bx, cx, dx, ex
2203	real cpinv, rdcp, dpinv
2204	real lvcp(nloc,na), mke(nloc,na)
2205	real am(nloc), work(nloc), ad(nloc), amp1(nloc)
2206	c!! real up1(nloc), dn1(nloc)
2207	real up1(nloc,nd,nd), dn1(nloc,nd,nd)
2208	real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2209	real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
2210	real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
2211
2212
2213	c-------------------------------------------------------------
2214
2215	c initialization:
2216
2217	delti = 1.0/delt
2218
2219	do il=1,ncum
2220	precip(il)=0.0
2221	wd(il)=0.0 ! gust
2222	enddo
2223
2224	do i=1,nd
2225	do il=1,ncum
2226	ft(il,i)=0.0
2227	fr(il,i)=0.0
2228	fu(il,i)=0.0
2229	fv(il,i)=0.0
2230	qcondc(il,i)=0.0 ! cld
2231	qcond(il,i)=0.0 ! cld
2232	nqcond(il,i)=0.0 ! cld
2233	enddo
2234	enddo
2235
2236	c do j=1,ntra
2237	c do i=1,nd
2238	c do il=1,ncum
2239	c ftra(il,i,j)=0.0
2240	c enddo
2241	c enddo
2242	c enddo
2243
2244	do i=1,nl
2245	do il=1,ncum
2246	lvcp(il,i)=lv(il,i)/cpn(il,i)
2247	enddo
2248	enddo
2249
2250
2251	c
2252	c * calculate surface precipitation in mm/day *
2253	c
2254	do il=1,ncum
2255	if(ep(il,inb(il)).ge.0.0001)then
2256	if (cvflag_grav) then
2257	precip(il)=wt(il,1)sigdwater(il,1)86400.1000./(rowl*grav)
2258	else
2259	precip(il)=wt(il,1)sigdwater(il,1)*8640.
2260	endif
2261	endif
2262	enddo
2263
2264	c
2265	c * Calculate downdraft velocity scale *
2266	c * NE PAS UTILISER POUR L'INSTANT *
2267	c
2268	c! do il=1,ncum
2269	c! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
2270	c! : /(sigd*p(il,icb(il)))
2271	c! enddo
2272
2273	c
2274	c * calculate tendencies of lowest level potential temperature *
2275	c * and mixing ratio *
2276	c
2277	do il=1,ncum
2278	work(il)=1.0/(ph(il,1)-ph(il,2))
2279	am(il)=0.0
2280	enddo
2281
2282	do k=2,nl
2283	do il=1,ncum
2284	if (k.le.inb(il)) then
2285	am(il)=am(il)+m(il,k)
2286	endif
2287	enddo
2288	enddo
2289
2290	do il=1,ncum
2291
2292	c convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
2293	if (cvflag_grav) then
2294	if((0.01gravwork(il)*am(il)).ge.delti)iflag(il)=1!consist vect
2295	ft(il,1)=0.01gravwork(il)am(il)(t(il,2)-t(il,1)
2296	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2297	else
2298	if((0.1work(il)am(il)).ge.delti)iflag(il)=1 !consistency vect
2299	ft(il,1)=0.1work(il)am(il)*(t(il,2)-t(il,1)
2300	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2301	endif
2302
2303	ft(il,1)=ft(il,1)-0.5lvcp(il,1)sigd*(evap(il,1)+evap(il,2))
2304
2305	if (cvflag_grav) then
2306	ft(il,1)=ft(il,1)-0.009gravsigd*mp(il,2)
2307	: t(il,1)b(il,1)*work(il)
2308	else
2309	ft(il,1)=ft(il,1)-0.09sigdmp(il,2)t(il,1)b(il,1)*work(il)
2310	endif
2311
2312	ft(il,1)=ft(il,1)+0.01sigdwt(il,1)(cl-cpd)water(il,2)*(t(il,2)
2313	:-t(il,1))*work(il)/cpn(il,1)
2314
2315	if (cvflag_grav) then
2316	Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
2317	c (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
2318	fr(il,1)=0.01gravmp(il,2)(rp(il,2)-rr(il,1))work(il)
2319	: +sigd0.5(evap(il,1)+evap(il,2))
2320	c+tard : +sigd*evap(il,1)
2321
2322	fr(il,1)=fr(il,1)+0.01gravam(il)(rr(il,2)-rr(il,1))work(il)
2323
2324	fu(il,1)=fu(il,1)+0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1))
2325	: +am(il)*(u(il,2)-u(il,1)))
2326	fv(il,1)=fv(il,1)+0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1))
2327	: +am(il)*(v(il,2)-v(il,1)))
2328	else ! cvflag_grav
2329	fr(il,1)=0.1mp(il,2)(rp(il,2)-rr(il,1))*work(il)
2330	: +sigd0.5(evap(il,1)+evap(il,2))
2331	fr(il,1)=fr(il,1)+0.1am(il)(rr(il,2)-rr(il,1))*work(il)
2332	fu(il,1)=fu(il,1)+0.1work(il)(mp(il,2)*(up(il,2)-u(il,1))
2333	: +am(il)*(u(il,2)-u(il,1)))
2334	fv(il,1)=fv(il,1)+0.1work(il)(mp(il,2)*(vp(il,2)-v(il,1))
2335	: +am(il)*(v(il,2)-v(il,1)))
2336	endif ! cvflag_grav
2337
2338	enddo ! il
2339
2340	c do j=1,ntra
2341	c do il=1,ncum
2342	c if (cvflag_grav) then
2343	c ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
2344	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2345	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2346	c else
2347	c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2348	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2349	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2350	c endif
2351	c enddo
2352	c enddo
2353
2354	do j=2,nl
2355	do il=1,ncum
2356	if (j.le.inb(il)) then
2357	if (cvflag_grav) then
2358	fr(il,1)=fr(il,1)
2359	: +0.01gravwork(il)ment(il,j,1)(qent(il,j,1)-rr(il,1))
2360	fu(il,1)=fu(il,1)
2361	: +0.01gravwork(il)ment(il,j,1)(uent(il,j,1)-u(il,1))
2362	fv(il,1)=fv(il,1)
2363	: +0.01gravwork(il)ment(il,j,1)(vent(il,j,1)-v(il,1))
2364	else ! cvflag_grav
2365	fr(il,1)=fr(il,1)
2366	: +0.1work(il)ment(il,j,1)*(qent(il,j,1)-rr(il,1))
2367	fu(il,1)=fu(il,1)
2368	: +0.1work(il)ment(il,j,1)*(uent(il,j,1)-u(il,1))
2369	fv(il,1)=fv(il,1)
2370	: +0.1work(il)ment(il,j,1)*(vent(il,j,1)-v(il,1))
2371	endif ! cvflag_grav
2372	endif ! j
2373	enddo
2374	enddo
2375
2376	c do k=1,ntra
2377	c do j=2,nl
2378	c do il=1,ncum
2379	c if (j.le.inb(il)) then
2380
2381	c if (cvflag_grav) then
2382	c ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
2383	c : *(traent(il,j,1,k)-tra(il,1,k))
2384	c else
2385	c ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
2386	c : *(traent(il,j,1,k)-tra(il,1,k))
2387	c endif
2388
2389	c endif
2390	c enddo
2391	c enddo
2392	c enddo
2393
2394	c
2395	c * calculate tendencies of potential temperature and mixing ratio *
2396	c * at levels above the lowest level *
2397	c
2398	c * first find the net saturated updraft and downdraft mass fluxes *
2399	c * through each level *
2400	c
2401
2402	do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
2403
2404	num1=0
2405	do il=1,ncum
2406	if(i.le.inb(il))num1=num1+1
2407	enddo
2408	if(num1.le.0)go to 500
2409
2410	call zilch(amp1,ncum)
2411	call zilch(ad,ncum)
2412
2413	do 440 k=i+1,nl+1
2414	do 441 il=1,ncum
2415	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
2416	amp1(il)=amp1(il)+m(il,k)
2417	endif
2418	441 continue
2419	440 continue
2420
2421	do 450 k=1,i
2422	do 451 j=i+1,nl+1
2423	do 452 il=1,ncum
2424	if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
2425	amp1(il)=amp1(il)+ment(il,k,j)
2426	endif
2427	452 continue
2428	451 continue
2429	450 continue
2430
2431	do 470 k=1,i-1
2432	do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
2433	do 472 il=1,ncum
2434	if (i.le.inb(il) .and. j.le.inb(il)) then
2435	ad(il)=ad(il)+ment(il,j,k)
2436	endif
2437	472 continue
2438	471 continue
2439	470 continue
2440
2441	do 1350 il=1,ncum
2442	if (i.le.inb(il)) then
2443	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2444	cpinv=1.0/cpn(il,i)
2445
2446	c convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
2447	if (cvflag_grav) then
2448	if((0.01gravdpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
2449	else
2450	if((0.1dpinvamp1(il)).ge.delti)iflag(il)=1 ! vecto
2451	endif
2452
2453	if (cvflag_grav) then
2454	ft(il,i)=0.01gravdpinv(amp1(il)(t(il,i+1)-t(il,i)
2455	: +(gz(il,i+1)-gz(il,i))*cpinv)
2456	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2457	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
2458	rat=cpn(il,i-1)*cpinv
2459	ft(il,i)=ft(il,i)-0.009gravsigd(mp(il,i+1)t(il,i)*b(il,i)
2460	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2461	ft(il,i)=ft(il,i)+0.01gravdpinvment(il,i,i)(hp(il,i)-h(il,i)
2462	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2463	else ! cvflag_grav
2464	ft(il,i)=0.1dpinv(amp1(il)*(t(il,i+1)-t(il,i)
2465	: +(gz(il,i+1)-gz(il,i))*cpinv)
2466	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2467	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
2468	rat=cpn(il,i-1)*cpinv
2469	ft(il,i)=ft(il,i)-0.09sigd(mp(il,i+1)t(il,i)b(il,i)
2470	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2471	ft(il,i)=ft(il,i)+0.1dpinvment(il,i,i)*(hp(il,i)-h(il,i)
2472	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2473	endif ! cvflag_grav
2474
2475
2476	ft(il,i)=ft(il,i)+0.01sigdwt(il,i)(cl-cpd)water(il,i+1)
2477	: (t(il,i+1)-t(il,i))dpinv*cpinv
2478
2479	if (cvflag_grav) then
2480	fr(il,i)=0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i))
2481	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2482	fu(il,i)=fu(il,i)+0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i))
2483	: -ad(il)*(u(il,i)-u(il,i-1)))
2484	fv(il,i)=fv(il,i)+0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i))
2485	: -ad(il)*(v(il,i)-v(il,i-1)))
2486	else ! cvflag_grav
2487	fr(il,i)=0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il,i))
2488	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2489	fu(il,i)=fu(il,i)+0.1dpinv(amp1(il)*(u(il,i+1)-u(il,i))
2490	: -ad(il)*(u(il,i)-u(il,i-1)))
2491	fv(il,i)=fv(il,i)+0.1dpinv(amp1(il)*(v(il,i+1)-v(il,i))
2492	: -ad(il)*(v(il,i)-v(il,i-1)))
2493	endif ! cvflag_grav
2494
2495	endif ! i
2496	1350 continue
2497
2498	c do k=1,ntra
2499	c do il=1,ncum
2500	c if (i.le.inb(il)) then
2501	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2502	c cpinv=1.0/cpn(il,i)
2503	c if (cvflag_grav) then
2504	c ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
2505	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2506	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2507	c else
2508	c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2509	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2510	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2511	c endif
2512	c endif
2513	c enddo
2514	c enddo
2515
2516	do 480 k=1,i-1
2517	do 1370 il=1,ncum
2518	if (i.le.inb(il)) then
2519	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2520	cpinv=1.0/cpn(il,i)
2521
2522	awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
2523	awat=amax1(awat,0.0)
2524
2525	if (cvflag_grav) then
2526	fr(il,i)=fr(il,i)
2527	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-awat-rr(il,i))
2528	fu(il,i)=fu(il,i)
2529	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2530	fv(il,i)=fv(il,i)
2531	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2532	else ! cvflag_grav
2533	fr(il,i)=fr(il,i)
2534	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
2535	fu(il,i)=fu(il,i)
2536	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2537	fv(il,i)=fv(il,i)
2538	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2539	endif ! cvflag_grav
2540
2541	c (saturated updrafts resulting from mixing) ! cld
2542	qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
2543	nqcond(il,i)=nqcond(il,i)+1. ! cld
2544	endif ! i
2545	1370 continue
2546	480 continue
2547
2548	c do j=1,ntra
2549	c do k=1,i-1
2550	c do il=1,ncum
2551	c if (i.le.inb(il)) then
2552	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2553	c cpinv=1.0/cpn(il,i)
2554	c if (cvflag_grav) then
2555	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2556	c : *(traent(il,k,i,j)-tra(il,i,j))
2557	c else
2558	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2559	c : *(traent(il,k,i,j)-tra(il,i,j))
2560	c endif
2561	c endif
2562	c enddo
2563	c enddo
2564	c enddo
2565
2566	do 490 k=i,nl+1
2567	do 1380 il=1,ncum
2568	if (i.le.inb(il) .and. k.le.inb(il)) then
2569	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2570	cpinv=1.0/cpn(il,i)
2571
2572	if (cvflag_grav) then
2573	fr(il,i)=fr(il,i)
2574	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-rr(il,i))
2575	fu(il,i)=fu(il,i)
2576	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2577	fv(il,i)=fv(il,i)
2578	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2579	else ! cvflag_grav
2580	fr(il,i)=fr(il,i)
2581	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-rr(il,i))
2582	fu(il,i)=fu(il,i)
2583	: +0.1dpinvment(il,k,i)*(uent(il,k,i)-u(il,i))
2584	fv(il,i)=fv(il,i)
2585	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2586	endif ! cvflag_grav
2587	endif ! i and k
2588	1380 continue
2589	490 continue
2590
2591	c do j=1,ntra
2592	c do k=i,nl+1
2593	c do il=1,ncum
2594	c if (i.le.inb(il) .and. k.le.inb(il)) then
2595	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2596	c cpinv=1.0/cpn(il,i)
2597	c if (cvflag_grav) then
2598	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2599	c : *(traent(il,k,i,j)-tra(il,i,j))
2600	c else
2601	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2602	c : *(traent(il,k,i,j)-tra(il,i,j))
2603	c endif
2604	c endif ! i and k
2605	c enddo
2606	c enddo
2607	c enddo
2608
2609	do 1400 il=1,ncum
2610	if (i.le.inb(il)) then
2611	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2612	cpinv=1.0/cpn(il,i)
2613
2614	if (cvflag_grav) then
2615	c sb: on ne fait pas encore la correction permettant de mieux
2616	c conserver l'eau:
2617	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
2618	: +0.01grav(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
2619	: (rp(il,i)-rr(il,i-1)))dpinv
2620
2621	fu(il,i)=fu(il,i)+0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i))
2622	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2623	fv(il,i)=fv(il,i)+0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i))
2624	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2625	else ! cvflag_grav
2626	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
2627	: +0.1(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)
2628	: (rp(il,i)-rr(il,i-1)))dpinv
2629	fu(il,i)=fu(il,i)+0.1(mp(il,i+1)(up(il,i+1)-u(il,i))
2630	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2631	fv(il,i)=fv(il,i)+0.1(mp(il,i+1)(vp(il,i+1)-v(il,i))
2632	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2633	endif ! cvflag_grav
2634
2635	endif ! i
2636	1400 continue
2637
2638	c sb: interface with the cloud parameterization: ! cld
2639
2640	do k=i+1,nl
2641	do il=1,ncum
2642	if (k.le.inb(il) .and. i.le.inb(il)) then ! cld
2643	C (saturated downdrafts resulting from mixing) ! cld
2644	qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
2645	nqcond(il,i)=nqcond(il,i)+1. ! cld
2646	endif ! cld
2647	enddo ! cld
2648	enddo ! cld
2649
2650	C (particular case: no detraining level is found) ! cld
2651	do il=1,ncum ! cld
2652	if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld
2653	qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
2654	nqcond(il,i)=nqcond(il,i)+1. ! cld
2655	endif ! cld
2656	enddo ! cld
2657
2658	do il=1,ncum ! cld
2659	if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld
2660	qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
2661	endif ! cld
2662	enddo
2663
2664	c do j=1,ntra
2665	c do il=1,ncum
2666	c if (i.le.inb(il)) then
2667	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2668	c cpinv=1.0/cpn(il,i)
2669
2670	c if (cvflag_grav) then
2671	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
2672	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2673	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2674	c else
2675	c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2676	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2677	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2678	c endif
2679	c endif ! i
2680	c enddo
2681	c enddo
2682
2683	500 continue
2684
2685
2686	c * move the detrainment at level inb down to level inb-1 *
2687	c * in such a way as to preserve the vertically *
2688	c * integrated enthalpy and water tendencies *
2689	c
2690	do 503 il=1,ncum
2691
2692	ax=0.1ment(il,inb(il),inb(il))(hp(il,inb(il))-h(il,inb(il))
2693	: +t(il,inb(il))*(cpv-cpd)
2694	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
2695	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2696	ft(il,inb(il))=ft(il,inb(il))-ax
2697	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
2698	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
2699	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
2700
2701	bx=0.1ment(il,inb(il),inb(il))(qent(il,inb(il),inb(il))
2702	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2703	fr(il,inb(il))=fr(il,inb(il))-bx
2704	fr(il,inb(il)-1)=fr(il,inb(il)-1)
2705	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
2706	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2707
2708	cx=0.1ment(il,inb(il),inb(il))(uent(il,inb(il),inb(il))
2709	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2710	fu(il,inb(il))=fu(il,inb(il))-cx
2711	fu(il,inb(il)-1)=fu(il,inb(il)-1)
2712	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
2713	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2714
2715	dx=0.1ment(il,inb(il),inb(il))(vent(il,inb(il),inb(il))
2716	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2717	fv(il,inb(il))=fv(il,inb(il))-dx
2718	fv(il,inb(il)-1)=fv(il,inb(il)-1)
2719	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
2720	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2721
2722	503 continue
2723
2724	c do j=1,ntra
2725	c do il=1,ncum
2726	c ex=0.1*ment(il,inb(il),inb(il))
2727	c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2728	c : /(ph(il,inb(il))-ph(il,inb(il)+1))
2729	c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2730	c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2731	c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2732	c : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2733	c enddo
2734	c enddo
2735
2736	c
2737	c * homoginize tendencies below cloud base *
2738	c
2739	c
2740	do il=1,ncum
2741	asum(il)=0.0
2742	bsum(il)=0.0
2743	csum(il)=0.0
2744	dsum(il)=0.0
2745	enddo
2746
2747	do i=1,nl
2748	do il=1,ncum
2749	if (i.le.(icb(il)-1)) then
2750	asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
2751	bsum(il)=bsum(il)+fr(il,i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))
2752	: *(ph(il,i)-ph(il,i+1))
2753	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
2754	: *(ph(il,i)-ph(il,i+1))
2755	dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
2756	endif
2757	enddo
2758	enddo
2759
2760	c!!! do 700 i=1,icb(il)-1
2761	do i=1,nl
2762	do il=1,ncum
2763	if (i.le.(icb(il)-1)) then
2764	ft(il,i)=asum(il)t(il,i)/(th(il,i)dsum(il))
2765	fr(il,i)=bsum(il)/csum(il)
2766	endif
2767	enddo
2768	enddo
2769
2770	c
2771	c * reset counter and return *
2772	c
2773	do il=1,ncum
2774	sig(il,nd)=2.0
2775	enddo
2776
2777
2778	do i=1,nd
2779	do il=1,ncum
2780	upwd(il,i)=0.0
2781	dnwd(il,i)=0.0
2782	enddo
2783	enddo
2784
2785	do i=1,nl
2786	do il=1,ncum
2787	dnwd0(il,i)=-mp(il,i)
2788	enddo
2789	enddo
2790	do i=nl+1,nd
2791	do il=1,ncum
2792	dnwd0(il,i)=0.
2793	enddo
2794	enddo
2795
2796
2797	do i=1,nl
2798	do il=1,ncum
2799	if (i.ge.icb(il) .and. i.le.inb(il)) then
2800	upwd(il,i)=0.0
2801	dnwd(il,i)=0.0
2802	endif
2803	enddo
2804	enddo
2805
2806	do i=1,nl
2807	do k=1,nl
2808	do il=1,ncum
2809	up1(il,k,i)=0.0
2810	dn1(il,k,i)=0.0
2811	enddo
2812	enddo
2813	enddo
2814
2815	do i=1,nl
2816	do k=i,nl
2817	do n=1,i-1
2818	do il=1,ncum
2819	if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2820	up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
2821	dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
2822	endif
2823	enddo
2824	enddo
2825	enddo
2826	enddo
2827
2828	do i=2,nl
2829	do k=i,nl
2830	do il=1,ncum
2831	ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2832	if (i.le.inb(il).and.k.le.inb(il)) then
2833	upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
2834	dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
2835	endif
2836	enddo
2837	enddo
2838	enddo
2839
2840
2841	c!!! DO il=1,ncum
2842	c!!! do i=icb(il),inb(il)
2843	c!!!
2844	c!!! upwd(il,i)=0.0
2845	c!!! dnwd(il,i)=0.0
2846	c!!! do k=i,inb(il)
2847	c!!! up1=0.0
2848	c!!! dn1=0.0
2849	c!!! do n=1,i-1
2850	c!!! up1=up1+ment(il,n,k)
2851	c!!! dn1=dn1-ment(il,k,n)
2852	c!!! enddo
2853	c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
2854	c!!! dnwd(il,i)=dnwd(il,i)+dn1
2855	c!!! enddo
2856	c!!! enddo
2857	c!!!
2858	c!!! ENDDO
2859
2860	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2861	c determination de la variation de flux ascendant entre
2862	c deux niveau non dilue mike
2863	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2864
2865	do i=1,nl
2866	do il=1,ncum
2867	mike(il,i)=m(il,i)
2868	enddo
2869	enddo
2870
2871	do i=nl+1,nd
2872	do il=1,ncum
2873	mike(il,i)=0.
2874	enddo
2875	enddo
2876
2877	do i=1,nd
2878	do il=1,ncum
2879	ma(il,i)=0
2880	enddo
2881	enddo
2882
2883	do i=1,nl
2884	do j=i,nl
2885	do il=1,ncum
2886	ma(il,i)=ma(il,i)+m(il,j)
2887	enddo
2888	enddo
2889	enddo
2890
2891	do i=nl+1,nd
2892	do il=1,ncum
2893	ma(il,i)=0.
2894	enddo
2895	enddo
2896
2897	do i=1,nl
2898	do il=1,ncum
2899	if (i.le.(icb(il)-1)) then
2900	ma(il,i)=0
2901	endif
2902	enddo
2903	enddo
2904
2905	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2906	c icb represente de niveau ou se trouve la
2907	c base du nuage , et inb le top du nuage
2908	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2909
2910	do i=1,nd
2911	do il=1,ncum
2912	mke(il,i)=upwd(il,i)+dnwd(il,i)
2913	enddo
2914	enddo
2915
2916	do i=1,nd
2917	DO 999 il=1,ncum
2918	rdcp=(rrd(1.-rr(il,i))-rr(il,i)rrv)
2919	: /(cpd(1.-rr(il,i))+rr(il,i)cpv)
2920	tls(il,i)=t(il,i)(1000.0/p(il,i))*rdcp
2921	tps(il,i)=tp(il,i)
2922	999 CONTINUE
2923	enddo
2924
2925	c
2926	c * diagnose the in-cloud mixing ratio * ! cld
2927	c * of condensed water * ! cld
2928	c ! cld
2929
2930	do i=1,nd ! cld
2931	do il=1,ncum ! cld
2932	mac(il,i)=0.0 ! cld
2933	wa(il,i)=0.0 ! cld
2934	siga(il,i)=0.0 ! cld
2935	sax(il,i)=0.0 ! cld
2936	enddo ! cld
2937	enddo ! cld
2938
2939	do i=minorig, nl ! cld
2940	do k=i+1,nl+1 ! cld
2941	do il=1,ncum ! cld
2942	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld
2943	mac(il,i)=mac(il,i)+m(il,k) ! cld
2944	endif ! cld
2945	enddo ! cld
2946	enddo ! cld
2947	enddo ! cld
2948
2949	do i=1,nl ! cld
2950	do j=1,i ! cld
2951	do il=1,ncum ! cld
2952	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
2953	: .and. j.ge.icb(il) ) then ! cld
2954	sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
2955	: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
2956	endif ! cld
2957	enddo ! cld
2958	enddo ! cld
2959	enddo ! cld
2960
2961	do i=1,nl ! cld
2962	do il=1,ncum ! cld
2963	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
2964	: .and. sax(il,i).gt.0.0 ) then ! cld
2965	wa(il,i)=sqrt(2.*sax(il,i)) ! cld
2966	endif ! cld
2967	enddo ! cld
2968	enddo ! cld
2969
2970	do i=1,nl ! cld
2971	do il=1,ncum ! cld
2972	if (wa(il,i).gt.0.0) ! cld
2973	: siga(il,i)=mac(il,i)/wa(il,i) ! cld
2974	: rrdtvp(il,i)/p(il,i)/100./delta ! cld
2975	siga(il,i) = min(siga(il,i),1.0) ! cld
2976	cIM cf. FH
2977	if (iflag_clw.eq.0) then
2978	qcondc(il,i)=siga(il,i)clw(il,i)(1.-ep(il,i)) ! cld
2979	: + (1.-siga(il,i))*qcond(il,i) ! cld
2980	else if (iflag_clw.eq.1) then
2981	qcondc(il,i)=qcond(il,i) ! cld
2982	endif
2983
2984	enddo ! cld
2985	enddo ! cld
2986
2987	return
2988	end
2989
2990	SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
2991	& ment,sij,da,phi)
2992	implicit none
2993	c inputs:
2994	integer ncum, nd, na, nloc,len
2995	real ment(nloc,na,na),sij(nloc,na,na)
2996	c ouputs:
2997	real da(nloc,na),phi(nloc,na,na)
2998	c local variables:
2999	integer i,j,k
3000	c
3001	da(:,:)=0.
3002	c
3003	do j=1,na
3004	do k=1,na
3005	do i=1,ncum
3006	da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
3007	phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
3008	c print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j)
3009	end do
3010	end do
3011	end do
3012
3013	return
3014	end
3015
3016
3017	SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
3018	: ,iflag
3019	: ,precip,sig,w0
3020	: ,ft,fq,fu,fv,ftra
3021	: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
3022	: ,iflag1
3023	: ,precip1,sig1,w01
3024	: ,ft1,fq1,fu1,fv1,ftra1
3025	: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
3026	: )
3027	implicit none
3028
3029	#include "cvparam3.h"
3030
3031	c inputs:
3032	integer len, ncum, nd, ntra, nloc
3033	integer idcum(nloc)
3034	integer iflag(nloc)
3035	real precip(nloc)
3036	real sig(nloc,nd), w0(nloc,nd)
3037	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
3038	real ftra(nloc,nd,ntra)
3039	real Ma(nloc,nd)
3040	real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
3041	real qcondc(nloc,nd)
3042	real wd(nloc),cape(nloc)
3043
3044	c outputs:
3045	integer iflag1(len)
3046	real precip1(len)
3047	real sig1(len,nd), w01(len,nd)
3048	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
3049	real ftra1(len,nd,ntra)
3050	real Ma1(len,nd)
3051	real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
3052	real qcondc1(nloc,nd)
3053	real wd1(nloc),cape1(nloc)
3054
3055	c local variables:
3056	integer i,k,j
3057
3058	do 2000 i=1,ncum
3059	precip1(idcum(i))=precip(i)
3060	iflag1(idcum(i))=iflag(i)
3061	wd1(idcum(i))=wd(i)
3062	cape1(idcum(i))=cape(i)
3063	2000 continue
3064
3065	do 2020 k=1,nl
3066	do 2010 i=1,ncum
3067	sig1(idcum(i),k)=sig(i,k)
3068	w01(idcum(i),k)=w0(i,k)
3069	ft1(idcum(i),k)=ft(i,k)
3070	fq1(idcum(i),k)=fq(i,k)
3071	fu1(idcum(i),k)=fu(i,k)
3072	fv1(idcum(i),k)=fv(i,k)
3073	Ma1(idcum(i),k)=Ma(i,k)
3074	upwd1(idcum(i),k)=upwd(i,k)
3075	dnwd1(idcum(i),k)=dnwd(i,k)
3076	dnwd01(idcum(i),k)=dnwd0(i,k)
3077	qcondc1(idcum(i),k)=qcondc(i,k)
3078	2010 continue
3079	2020 continue
3080
3081	do 2200 i=1,ncum
3082	sig1(idcum(i),nd)=sig(i,nd)
3083	2200 continue
3084
3085
3086	c do 2100 j=1,ntra
3087	c do 2110 k=1,nd ! oct3
3088	c do 2120 i=1,ncum
3089	c ftra1(idcum(i),k,j)=ftra(i,k,j)
3090	c 2120 continue
3091	c 2110 continue
3092	c 2100 continue
3093
3094	return
3095	end
3096

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

Download in other formats:

Original Format