Context Navigation

cv3_routines.F @ 1160

Last change on this file since 1160 was 704, checked in by Laurent Fairhead, 18 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
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	mp(:,:)=0.
1872
1873	do i=1,nl
1874	do il=1,ncum
1875	mp(il,i)=0.0
1876	rp(il,i)=rr(il,i)
1877	up(il,i)=u(il,i)
1878	vp(il,i)=v(il,i)
1879	wt(il,i)=0.001
1880	water(il,i)=0.0
1881	evap(il,i)=0.0
1882	b(il,i)=0.0
1883	lvcp(il,i)=lv(il,i)/cpn(il,i)
1884	enddo
1885	enddo
1886
1887	c do k=1,ntra
1888	c do i=1,nd
1889	c do il=1,ncum
1890	c trap(il,i,k)=tra(il,i,k)
1891	c enddo
1892	c enddo
1893	c enddo
1894
1895	c
1896	c * check whether ep(inb)=0, if so, skip precipitating *
1897	c * downdraft calculation *
1898	c
1899
1900	do il=1,ncum
1901	lwork(il)=.TRUE.
1902	if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
1903	enddo
1904
1905	call zilch(wdtrain,ncum)
1906
1907	DO 400 i=nl+1,1,-1
1908
1909	num1=0
1910	do il=1,ncum
1911	if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1
1912	enddo
1913	if (num1.le.0) goto 400
1914
1915	c
1916	c * integrate liquid water equation to find condensed water *
1917	c * and condensed water flux *
1918	c
1919
1920	c
1921	c * begin downdraft loop *
1922	c
1923
1924	c
1925	c * calculate detrained precipitation *
1926	c
1927	do il=1,ncum
1928	if (i.le.inb(il) .and. lwork(il)) then
1929	if (cvflag_grav) then
1930	wdtrain(il)=gravep(il,i)m(il,i)*clw(il,i)
1931	else
1932	wdtrain(il)=10.0ep(il,i)m(il,i)*clw(il,i)
1933	endif
1934	endif
1935	enddo
1936
1937	if(i.gt.1)then
1938	do 320 j=1,i-1
1939	do il=1,ncum
1940	if (i.le.inb(il) .and. lwork(il)) then
1941	awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
1942	awat=amax1(awat,0.0)
1943	if (cvflag_grav) then
1944	wdtrain(il)=wdtrain(il)+gravawatment(il,j,i)
1945	else
1946	wdtrain(il)=wdtrain(il)+10.0awatment(il,j,i)
1947	endif
1948	endif
1949	enddo
1950	320 continue
1951	endif
1952
1953	c
1954	c * find rain water and evaporation using provisional *
1955	c * estimates of rp(i)and rp(i-1) *
1956	c
1957
1958	do 999 il=1,ncum
1959
1960	if (i.le.inb(il) .and. lwork(il)) then
1961
1962	wt(il,i)=45.0
1963
1964	if(i.lt.inb(il))then
1965	rp(il,i)=rp(il,i+1)
1966	: +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
1967	rp(il,i)=0.5*(rp(il,i)+rr(il,i))
1968	endif
1969	rp(il,i)=amax1(rp(il,i),0.0)
1970	rp(il,i)=amin1(rp(il,i),rs(il,i))
1971	rp(il,inb(il))=rr(il,inb(il))
1972
1973	if(i.eq.1)then
1974	afac=p(il,1)(rs(il,1)-rp(il,1))/(1.0e4+2000.0p(il,1)*rs(il,1))
1975	else
1976	rp(il,i-1)=rp(il,i)
1977	: +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i)
1978	rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
1979	rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
1980	rp(il,i-1)=amax1(rp(il,i-1),0.0)
1981	afac1=p(il,i)(rs(il,i)-rp(il,i))/(1.0e4+2000.0p(il,i)*rs(il,i))
1982	afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
1983	: /(1.0e4+2000.0p(il,i-1)rs(il,i-1))
1984	afac=0.5*(afac1+afac2)
1985	endif
1986	if(i.eq.inb(il))afac=0.0
1987	afac=amax1(afac,0.0)
1988	bfac=1./(sigd*wt(il,i))
1989	c
1990	cjyg1
1991	ccc sigt=1.0
1992	ccc if(i.ge.icb)sigt=sigp(i)
1993	c prise en compte de la variation progressive de sigt dans
1994	c les couches icb et icb-1:
1995	c pour plcl<ph(i+1), pr1=0 & pr2=1
1996	c pour plcl>ph(i), pr1=1 & pr2=0
1997	c pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
1998	c sur le nuage, et pr2 est la proportion sous la base du
1999	c nuage.
2000	pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2001	pr1=max(0.,min(1.,pr1))
2002	pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2003	pr2=max(0.,min(1.,pr2))
2004	sigt=sigp(il,i)*pr1+pr2
2005	cjyg2
2006	c
2007	b6=bfac50.sigd(ph(il,i)-ph(il,i+1))sigt*afac
2008	c6=water(il,i+1)+bfac*wdtrain(il)
2009	: -50.sigdbfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2010	if(c6.gt.0.0)then
2011	revap=0.5(-b6+sqrt(b6b6+4.*c6))
2012	evap(il,i)=sigtafacrevap
2013	water(il,i)=revap*revap
2014	else
2015	evap(il,i)=-evap(il,i+1)
2016	: +0.02(wdtrain(il)+sigdwt(il,i)*water(il,i+1))
2017	: /(sigd*(ph(il,i)-ph(il,i+1)))
2018	end if
2019	c
2020	c * calculate precipitating downdraft mass flux under *
2021	c * hydrostatic approximation *
2022	c
2023	if (i.ne.1) then
2024
2025	tevap=amax1(0.0,evap(il,i))
2026	delth=amax1(0.001,(th(il,i)-th(il,i-1)))
2027	if (cvflag_grav) then
2028	mp(il,i)=100.ginvlvcp(il,i)sigdtevap
2029	: *(p(il,i-1)-p(il,i))/delth
2030	else
2031	mp(il,i)=10.lvcp(il,i)sigdtevap(p(il,i-1)-p(il,i))/delth
2032	endif
2033	c
2034	c * if hydrostatic assumption fails, *
2035	c * solve cubic difference equation for downdraft theta *
2036	c * and mass flux from two simultaneous differential eqns *
2037	c
2038	amfac=sigdsigd70.0ph(il,i)(p(il,i-1)-p(il,i))
2039	: (th(il,i)-th(il,i-1))/(tv(il,i)th(il,i))
2040	amp2=abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2041	if(amp2.gt.(0.1*amfac))then
2042	xf=100.0sigdsigdsigd(ph(il,i)-ph(il,i+1))
2043	tf=b(il,i)-5.0(th(il,i)-th(il,i-1))t(il,i)
2044	: /(lvcp(il,i)sigdth(il,i))
2045	af=xftf+mp(il,i+1)mp(il,i+1)*tinv
2046	bf=2.(tinvmp(il,i+1))*3+tinvmp(il,i+1)xftf
2047	: +50.(p(il,i-1)-p(il,i))xf*tevap
2048	fac2=1.0
2049	if(bf.lt.0.0)fac2=-1.0
2050	bf=abs(bf)
2051	ur=0.25bfbf-afafaftinvtinv*tinv
2052	if(ur.ge.0.0)then
2053	sru=sqrt(ur)
2054	fac=1.0
2055	if((0.5*bf-sru).lt.0.0)fac=-1.0
2056	mp(il,i)=mp(il,i+1)tinv+(0.5bf+sru)**tinv
2057	: +fac(abs(0.5bf-sru))**tinv
2058	else
2059	d=atan(2.*sqrt(-ur)/(bf+1.0e-28))
2060	if(fac2.lt.0.0)d=3.14159-d
2061	mp(il,i)=mp(il,i+1)tinv+2.sqrt(aftinv)cos(d*tinv)
2062	endif
2063	mp(il,i)=amax1(0.0,mp(il,i))
2064
2065	if (cvflag_grav) then
2066	Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2067	C il faut diviser par (mp(il,i)sigdgrav) et non par (mp(il,i)+sigd*0.1).
2068	C Et il faut bien revoir les facteurs 100.
2069	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2070	2 /(mp(il,i)+sigd*0.1)
2071	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
2072	else
2073	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2074	2 /(mp(il,i)+sigd*0.1)
2075	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
2076	endif
2077	b(il,i-1)=amax1(b(il,i-1),0.0)
2078	endif
2079	c
2080	c * limit magnitude of mp(i) to meet cfl condition *
2081	c
2082	ampmax=2.0(ph(il,i)-ph(il,i+1))delti
2083	amp2=2.0(ph(il,i-1)-ph(il,i))delti
2084	ampmax=amin1(ampmax,amp2)
2085	mp(il,i)=amin1(mp(il,i),ampmax)
2086	c
2087	c * force mp to decrease linearly to zero *
2088	c * between cloud base and the surface *
2089	c
2090	if(p(il,i).gt.p(il,icb(il)))then
2091	mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2092	endif
2093
2094	360 continue
2095	endif ! i.eq.1
2096	c
2097	c * find mixing ratio of precipitating downdraft *
2098	c
2099
2100	if (i.ne.inb(il)) then
2101
2102	rp(il,i)=rr(il,i)
2103
2104	if(mp(il,i).gt.mp(il,i+1))then
2105
2106	if (cvflag_grav) then
2107	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2108	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
2109	: *(evap(il,i+1)+evap(il,i))
2110	else
2111	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2112	: +5.sigd(ph(il,i)-ph(il,i+1))
2113	: *(evap(il,i+1)+evap(il,i))
2114	endif
2115	rp(il,i)=rp(il,i)/mp(il,i)
2116	up(il,i)=up(il,i+1)mp(il,i+1)+u(il,i)(mp(il,i)-mp(il,i+1))
2117	up(il,i)=up(il,i)/mp(il,i)
2118	vp(il,i)=vp(il,i+1)mp(il,i+1)+v(il,i)(mp(il,i)-mp(il,i+1))
2119	vp(il,i)=vp(il,i)/mp(il,i)
2120
2121	c do j=1,ntra
2122	c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2123	ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2124	c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2125	c trap(il,i,j)=trap(il,i,j)/mp(il,i)
2126	c end do
2127
2128	else
2129
2130	if(mp(il,i+1).gt.1.0e-16)then
2131	if (cvflag_grav) then
2132	rp(il,i)=rp(il,i+1)
2133	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
2134	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2135	else
2136	rp(il,i)=rp(il,i+1)
2137	: +5.sigd(ph(il,i)-ph(il,i+1))
2138	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2139	endif
2140	up(il,i)=up(il,i+1)
2141	vp(il,i)=vp(il,i+1)
2142
2143	c do j=1,ntra
2144	c trap(il,i,j)=trap(il,i+1,j)
2145	c end do
2146
2147	endif
2148	endif
2149	rp(il,i)=amin1(rp(il,i),rs(il,i))
2150	rp(il,i)=amax1(rp(il,i),0.0)
2151
2152	endif
2153	endif
2154	999 continue
2155
2156	400 continue
2157
2158	return
2159	end
2160
2161	SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
2162	: ,icb,inb,delt
2163	: ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th
2164	: ,ep,clw,m,tp,mp,rp,up,vp,trap
2165	: ,wt,water,evap,b
2166	: ,ment,qent,uent,vent,nent,elij,traent,sig
2167	: ,tv,tvp
2168	: ,iflag,precip,VPrecip,ft,fr,fu,fv,ftra
2169	: ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
2170	implicit none
2171
2172	#include "cvthermo.h"
2173	#include "cvparam3.h"
2174	#include "cvflag.h"
2175	#include "conema3.h"
2176
2177	c inputs:
2178	integer ncum,nd,na,ntra,nloc
2179	integer icb(nloc), inb(nloc)
2180	real delt
2181	real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
2182	real tra(nloc,nd,ntra), sig(nloc,nd)
2183	real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
2184	real th(nloc,na), p(nloc,nd), tp(nloc,na)
2185	real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
2186	real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
2187	real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
2188	real water(nloc,na), evap(nloc,na), b(nloc,na)
2189	real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
2190	cym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2191	real vent(nloc,na,na), elij(nloc,na,na)
2192	integer nent(nloc,na)
2193	real traent(nloc,na,na,ntra)
2194	real tv(nloc,nd), tvp(nloc,nd)
2195
2196	c input/output:
2197	integer iflag(nloc)
2198
2199	c outputs:
2200	real precip(nloc)
2201	real VPrecip(nloc,nd+1)
2202	real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
2203	real ftra(nloc,nd,ntra)
2204	real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
2205	real dnwd0(nloc,nd), mike(nloc,nd)
2206	real tls(nloc,nd), tps(nloc,nd)
2207	real qcondc(nloc,nd) ! cld
2208	real wd(nloc) ! gust
2209
2210	c local variables:
2211	integer i,k,il,n,j,num1
2212	real rat, awat, delti
2213	real ax, bx, cx, dx, ex
2214	real cpinv, rdcp, dpinv
2215	real lvcp(nloc,na), mke(nloc,na)
2216	real am(nloc), work(nloc), ad(nloc), amp1(nloc)
2217	c!! real up1(nloc), dn1(nloc)
2218	real up1(nloc,nd,nd), dn1(nloc,nd,nd)
2219	real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2220	real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
2221	real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
2222
2223
2224	c-------------------------------------------------------------
2225
2226	c initialization:
2227
2228	delti = 1.0/delt
2229
2230	do il=1,ncum
2231	precip(il)=0.0
2232	wd(il)=0.0 ! gust
2233	VPrecip(il,nd+1)=0.
2234	enddo
2235
2236	do i=1,nd
2237	do il=1,ncum
2238	VPrecip(il,i)=0.0
2239	ft(il,i)=0.0
2240	fr(il,i)=0.0
2241	fu(il,i)=0.0
2242	fv(il,i)=0.0
2243	qcondc(il,i)=0.0 ! cld
2244	qcond(il,i)=0.0 ! cld
2245	nqcond(il,i)=0.0 ! cld
2246	enddo
2247	enddo
2248
2249	c do j=1,ntra
2250	c do i=1,nd
2251	c do il=1,ncum
2252	c ftra(il,i,j)=0.0
2253	c enddo
2254	c enddo
2255	c enddo
2256
2257	do i=1,nl
2258	do il=1,ncum
2259	lvcp(il,i)=lv(il,i)/cpn(il,i)
2260	enddo
2261	enddo
2262
2263
2264	c
2265	c * calculate surface precipitation in mm/day *
2266	c
2267	do il=1,ncum
2268	if(ep(il,inb(il)).ge.0.0001)then
2269	if (cvflag_grav) then
2270	precip(il)=wt(il,1)sigdwater(il,1)86400.1000./(rowl*grav)
2271	else
2272	precip(il)=wt(il,1)sigdwater(il,1)*8640.
2273	endif
2274	endif
2275	enddo
2276
2277	C *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
2278	C
2279	c MAF rajout pour lessivage
2280	do k=1,nl
2281	do il=1,ncum
2282	if (k.le.inb(il)) then
2283	if (cvflag_grav) then
2284	VPrecip(il,k) = wt(il,k)sigdwater(il,k)/grav
2285	else
2286	VPrecip(il,k) = wt(il,k)sigdwater(il,k)/10.
2287	endif
2288	endif
2289	end do
2290	end do
2291	C
2292	c
2293	c * Calculate downdraft velocity scale *
2294	c * NE PAS UTILISER POUR L'INSTANT *
2295	c
2296	c! do il=1,ncum
2297	c! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
2298	c! : /(sigd*p(il,icb(il)))
2299	c! enddo
2300
2301	c
2302	c * calculate tendencies of lowest level potential temperature *
2303	c * and mixing ratio *
2304	c
2305	do il=1,ncum
2306	work(il)=1.0/(ph(il,1)-ph(il,2))
2307	am(il)=0.0
2308	enddo
2309
2310	do k=2,nl
2311	do il=1,ncum
2312	if (k.le.inb(il)) then
2313	am(il)=am(il)+m(il,k)
2314	endif
2315	enddo
2316	enddo
2317
2318	do il=1,ncum
2319
2320	c convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
2321	if (cvflag_grav) then
2322	if((0.01gravwork(il)*am(il)).ge.delti)iflag(il)=1!consist vect
2323	ft(il,1)=0.01gravwork(il)am(il)(t(il,2)-t(il,1)
2324	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2325	else
2326	if((0.1work(il)am(il)).ge.delti)iflag(il)=1 !consistency vect
2327	ft(il,1)=0.1work(il)am(il)*(t(il,2)-t(il,1)
2328	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2329	endif
2330
2331	ft(il,1)=ft(il,1)-0.5lvcp(il,1)sigd*(evap(il,1)+evap(il,2))
2332
2333	if (cvflag_grav) then
2334	ft(il,1)=ft(il,1)-0.009gravsigd*mp(il,2)
2335	: t(il,1)b(il,1)*work(il)
2336	else
2337	ft(il,1)=ft(il,1)-0.09sigdmp(il,2)t(il,1)b(il,1)*work(il)
2338	endif
2339
2340	ft(il,1)=ft(il,1)+0.01sigdwt(il,1)(cl-cpd)water(il,2)*(t(il,2)
2341	:-t(il,1))*work(il)/cpn(il,1)
2342
2343	if (cvflag_grav) then
2344	Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
2345	c (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
2346	fr(il,1)=0.01gravmp(il,2)(rp(il,2)-rr(il,1))work(il)
2347	: +sigd0.5(evap(il,1)+evap(il,2))
2348	c+tard : +sigd*evap(il,1)
2349
2350	fr(il,1)=fr(il,1)+0.01gravam(il)(rr(il,2)-rr(il,1))work(il)
2351
2352	fu(il,1)=fu(il,1)+0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1))
2353	: +am(il)*(u(il,2)-u(il,1)))
2354	fv(il,1)=fv(il,1)+0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1))
2355	: +am(il)*(v(il,2)-v(il,1)))
2356	else ! cvflag_grav
2357	fr(il,1)=0.1mp(il,2)(rp(il,2)-rr(il,1))*work(il)
2358	: +sigd0.5(evap(il,1)+evap(il,2))
2359	fr(il,1)=fr(il,1)+0.1am(il)(rr(il,2)-rr(il,1))*work(il)
2360	fu(il,1)=fu(il,1)+0.1work(il)(mp(il,2)*(up(il,2)-u(il,1))
2361	: +am(il)*(u(il,2)-u(il,1)))
2362	fv(il,1)=fv(il,1)+0.1work(il)(mp(il,2)*(vp(il,2)-v(il,1))
2363	: +am(il)*(v(il,2)-v(il,1)))
2364	endif ! cvflag_grav
2365
2366	enddo ! il
2367
2368	c do j=1,ntra
2369	c do il=1,ncum
2370	c if (cvflag_grav) then
2371	c ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
2372	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2373	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2374	c else
2375	c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2376	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2377	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2378	c endif
2379	c enddo
2380	c enddo
2381
2382	do j=2,nl
2383	do il=1,ncum
2384	if (j.le.inb(il)) then
2385	if (cvflag_grav) then
2386	fr(il,1)=fr(il,1)
2387	: +0.01gravwork(il)ment(il,j,1)(qent(il,j,1)-rr(il,1))
2388	fu(il,1)=fu(il,1)
2389	: +0.01gravwork(il)ment(il,j,1)(uent(il,j,1)-u(il,1))
2390	fv(il,1)=fv(il,1)
2391	: +0.01gravwork(il)ment(il,j,1)(vent(il,j,1)-v(il,1))
2392	else ! cvflag_grav
2393	fr(il,1)=fr(il,1)
2394	: +0.1work(il)ment(il,j,1)*(qent(il,j,1)-rr(il,1))
2395	fu(il,1)=fu(il,1)
2396	: +0.1work(il)ment(il,j,1)*(uent(il,j,1)-u(il,1))
2397	fv(il,1)=fv(il,1)
2398	: +0.1work(il)ment(il,j,1)*(vent(il,j,1)-v(il,1))
2399	endif ! cvflag_grav
2400	endif ! j
2401	enddo
2402	enddo
2403
2404	c do k=1,ntra
2405	c do j=2,nl
2406	c do il=1,ncum
2407	c if (j.le.inb(il)) then
2408
2409	c if (cvflag_grav) then
2410	c ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
2411	c : *(traent(il,j,1,k)-tra(il,1,k))
2412	c else
2413	c ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
2414	c : *(traent(il,j,1,k)-tra(il,1,k))
2415	c endif
2416
2417	c endif
2418	c enddo
2419	c enddo
2420	c enddo
2421
2422	c
2423	c * calculate tendencies of potential temperature and mixing ratio *
2424	c * at levels above the lowest level *
2425	c
2426	c * first find the net saturated updraft and downdraft mass fluxes *
2427	c * through each level *
2428	c
2429
2430	do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
2431
2432	num1=0
2433	do il=1,ncum
2434	if(i.le.inb(il))num1=num1+1
2435	enddo
2436	if(num1.le.0)go to 500
2437
2438	call zilch(amp1,ncum)
2439	call zilch(ad,ncum)
2440
2441	do 440 k=i+1,nl+1
2442	do 441 il=1,ncum
2443	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
2444	amp1(il)=amp1(il)+m(il,k)
2445	endif
2446	441 continue
2447	440 continue
2448
2449	do 450 k=1,i
2450	do 451 j=i+1,nl+1
2451	do 452 il=1,ncum
2452	if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
2453	amp1(il)=amp1(il)+ment(il,k,j)
2454	endif
2455	452 continue
2456	451 continue
2457	450 continue
2458
2459	do 470 k=1,i-1
2460	do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
2461	do 472 il=1,ncum
2462	if (i.le.inb(il) .and. j.le.inb(il)) then
2463	ad(il)=ad(il)+ment(il,j,k)
2464	endif
2465	472 continue
2466	471 continue
2467	470 continue
2468
2469	do 1350 il=1,ncum
2470	if (i.le.inb(il)) then
2471	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2472	cpinv=1.0/cpn(il,i)
2473
2474	c convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
2475	if (cvflag_grav) then
2476	if((0.01gravdpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
2477	else
2478	if((0.1dpinvamp1(il)).ge.delti)iflag(il)=1 ! vecto
2479	endif
2480
2481	if (cvflag_grav) then
2482	ft(il,i)=0.01gravdpinv(amp1(il)(t(il,i+1)-t(il,i)
2483	: +(gz(il,i+1)-gz(il,i))*cpinv)
2484	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2485	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
2486	rat=cpn(il,i-1)*cpinv
2487	ft(il,i)=ft(il,i)-0.009gravsigd(mp(il,i+1)t(il,i)*b(il,i)
2488	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2489	ft(il,i)=ft(il,i)+0.01gravdpinvment(il,i,i)(hp(il,i)-h(il,i)
2490	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2491	else ! cvflag_grav
2492	ft(il,i)=0.1dpinv(amp1(il)*(t(il,i+1)-t(il,i)
2493	: +(gz(il,i+1)-gz(il,i))*cpinv)
2494	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2495	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
2496	rat=cpn(il,i-1)*cpinv
2497	ft(il,i)=ft(il,i)-0.09sigd(mp(il,i+1)t(il,i)b(il,i)
2498	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2499	ft(il,i)=ft(il,i)+0.1dpinvment(il,i,i)*(hp(il,i)-h(il,i)
2500	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2501	endif ! cvflag_grav
2502
2503
2504	ft(il,i)=ft(il,i)+0.01sigdwt(il,i)(cl-cpd)water(il,i+1)
2505	: (t(il,i+1)-t(il,i))dpinv*cpinv
2506
2507	if (cvflag_grav) then
2508	fr(il,i)=0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i))
2509	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2510	fu(il,i)=fu(il,i)+0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i))
2511	: -ad(il)*(u(il,i)-u(il,i-1)))
2512	fv(il,i)=fv(il,i)+0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i))
2513	: -ad(il)*(v(il,i)-v(il,i-1)))
2514	else ! cvflag_grav
2515	fr(il,i)=0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il,i))
2516	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2517	fu(il,i)=fu(il,i)+0.1dpinv(amp1(il)*(u(il,i+1)-u(il,i))
2518	: -ad(il)*(u(il,i)-u(il,i-1)))
2519	fv(il,i)=fv(il,i)+0.1dpinv(amp1(il)*(v(il,i+1)-v(il,i))
2520	: -ad(il)*(v(il,i)-v(il,i-1)))
2521	endif ! cvflag_grav
2522
2523	endif ! i
2524	1350 continue
2525
2526	c do k=1,ntra
2527	c do il=1,ncum
2528	c if (i.le.inb(il)) then
2529	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2530	c cpinv=1.0/cpn(il,i)
2531	c if (cvflag_grav) then
2532	c ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
2533	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2534	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2535	c else
2536	c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2537	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2538	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2539	c endif
2540	c endif
2541	c enddo
2542	c enddo
2543
2544	do 480 k=1,i-1
2545	do 1370 il=1,ncum
2546	if (i.le.inb(il)) then
2547	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2548	cpinv=1.0/cpn(il,i)
2549
2550	awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
2551	awat=amax1(awat,0.0)
2552
2553	if (cvflag_grav) then
2554	fr(il,i)=fr(il,i)
2555	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-awat-rr(il,i))
2556	fu(il,i)=fu(il,i)
2557	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2558	fv(il,i)=fv(il,i)
2559	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2560	else ! cvflag_grav
2561	fr(il,i)=fr(il,i)
2562	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
2563	fu(il,i)=fu(il,i)
2564	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2565	fv(il,i)=fv(il,i)
2566	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2567	endif ! cvflag_grav
2568
2569	c (saturated updrafts resulting from mixing) ! cld
2570	qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
2571	nqcond(il,i)=nqcond(il,i)+1. ! cld
2572	endif ! i
2573	1370 continue
2574	480 continue
2575
2576	c do j=1,ntra
2577	c do k=1,i-1
2578	c do il=1,ncum
2579	c if (i.le.inb(il)) then
2580	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2581	c cpinv=1.0/cpn(il,i)
2582	c if (cvflag_grav) then
2583	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2584	c : *(traent(il,k,i,j)-tra(il,i,j))
2585	c else
2586	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2587	c : *(traent(il,k,i,j)-tra(il,i,j))
2588	c endif
2589	c endif
2590	c enddo
2591	c enddo
2592	c enddo
2593
2594	do 490 k=i,nl+1
2595	do 1380 il=1,ncum
2596	if (i.le.inb(il) .and. k.le.inb(il)) then
2597	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2598	cpinv=1.0/cpn(il,i)
2599
2600	if (cvflag_grav) then
2601	fr(il,i)=fr(il,i)
2602	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-rr(il,i))
2603	fu(il,i)=fu(il,i)
2604	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2605	fv(il,i)=fv(il,i)
2606	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2607	else ! cvflag_grav
2608	fr(il,i)=fr(il,i)
2609	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-rr(il,i))
2610	fu(il,i)=fu(il,i)
2611	: +0.1dpinvment(il,k,i)*(uent(il,k,i)-u(il,i))
2612	fv(il,i)=fv(il,i)
2613	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2614	endif ! cvflag_grav
2615	endif ! i and k
2616	1380 continue
2617	490 continue
2618
2619	c do j=1,ntra
2620	c do k=i,nl+1
2621	c do il=1,ncum
2622	c if (i.le.inb(il) .and. k.le.inb(il)) then
2623	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2624	c cpinv=1.0/cpn(il,i)
2625	c if (cvflag_grav) then
2626	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2627	c : *(traent(il,k,i,j)-tra(il,i,j))
2628	c else
2629	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2630	c : *(traent(il,k,i,j)-tra(il,i,j))
2631	c endif
2632	c endif ! i and k
2633	c enddo
2634	c enddo
2635	c enddo
2636
2637	do 1400 il=1,ncum
2638	if (i.le.inb(il)) then
2639	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2640	cpinv=1.0/cpn(il,i)
2641
2642	if (cvflag_grav) then
2643	c sb: on ne fait pas encore la correction permettant de mieux
2644	c conserver l'eau:
2645	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
2646	: +0.01grav(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
2647	: (rp(il,i)-rr(il,i-1)))dpinv
2648
2649	fu(il,i)=fu(il,i)+0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i))
2650	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2651	fv(il,i)=fv(il,i)+0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i))
2652	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2653	else ! cvflag_grav
2654	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
2655	: +0.1(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)
2656	: (rp(il,i)-rr(il,i-1)))dpinv
2657	fu(il,i)=fu(il,i)+0.1(mp(il,i+1)(up(il,i+1)-u(il,i))
2658	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2659	fv(il,i)=fv(il,i)+0.1(mp(il,i+1)(vp(il,i+1)-v(il,i))
2660	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2661	endif ! cvflag_grav
2662
2663	endif ! i
2664	1400 continue
2665
2666	c sb: interface with the cloud parameterization: ! cld
2667
2668	do k=i+1,nl
2669	do il=1,ncum
2670	if (k.le.inb(il) .and. i.le.inb(il)) then ! cld
2671	C (saturated downdrafts resulting from mixing) ! cld
2672	qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
2673	nqcond(il,i)=nqcond(il,i)+1. ! cld
2674	endif ! cld
2675	enddo ! cld
2676	enddo ! cld
2677
2678	C (particular case: no detraining level is found) ! cld
2679	do il=1,ncum ! cld
2680	if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld
2681	qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
2682	nqcond(il,i)=nqcond(il,i)+1. ! cld
2683	endif ! cld
2684	enddo ! cld
2685
2686	do il=1,ncum ! cld
2687	if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld
2688	qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
2689	endif ! cld
2690	enddo
2691
2692	c do j=1,ntra
2693	c do il=1,ncum
2694	c if (i.le.inb(il)) then
2695	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2696	c cpinv=1.0/cpn(il,i)
2697
2698	c if (cvflag_grav) then
2699	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
2700	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2701	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2702	c else
2703	c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2704	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2705	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2706	c endif
2707	c endif ! i
2708	c enddo
2709	c enddo
2710
2711	500 continue
2712
2713
2714	c * move the detrainment at level inb down to level inb-1 *
2715	c * in such a way as to preserve the vertically *
2716	c * integrated enthalpy and water tendencies *
2717	c
2718	do 503 il=1,ncum
2719
2720	ax=0.1ment(il,inb(il),inb(il))(hp(il,inb(il))-h(il,inb(il))
2721	: +t(il,inb(il))*(cpv-cpd)
2722	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
2723	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2724	ft(il,inb(il))=ft(il,inb(il))-ax
2725	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
2726	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
2727	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
2728
2729	bx=0.1ment(il,inb(il),inb(il))(qent(il,inb(il),inb(il))
2730	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2731	fr(il,inb(il))=fr(il,inb(il))-bx
2732	fr(il,inb(il)-1)=fr(il,inb(il)-1)
2733	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
2734	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2735
2736	cx=0.1ment(il,inb(il),inb(il))(uent(il,inb(il),inb(il))
2737	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2738	fu(il,inb(il))=fu(il,inb(il))-cx
2739	fu(il,inb(il)-1)=fu(il,inb(il)-1)
2740	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
2741	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2742
2743	dx=0.1ment(il,inb(il),inb(il))(vent(il,inb(il),inb(il))
2744	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2745	fv(il,inb(il))=fv(il,inb(il))-dx
2746	fv(il,inb(il)-1)=fv(il,inb(il)-1)
2747	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
2748	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2749
2750	503 continue
2751
2752	c do j=1,ntra
2753	c do il=1,ncum
2754	c ex=0.1*ment(il,inb(il),inb(il))
2755	c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2756	c : /(ph(il,inb(il))-ph(il,inb(il)+1))
2757	c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2758	c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2759	c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2760	c : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2761	c enddo
2762	c enddo
2763
2764	c
2765	c * homoginize tendencies below cloud base *
2766	c
2767	c
2768	do il=1,ncum
2769	asum(il)=0.0
2770	bsum(il)=0.0
2771	csum(il)=0.0
2772	dsum(il)=0.0
2773	enddo
2774
2775	do i=1,nl
2776	do il=1,ncum
2777	if (i.le.(icb(il)-1)) then
2778	asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
2779	bsum(il)=bsum(il)+fr(il,i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))
2780	: *(ph(il,i)-ph(il,i+1))
2781	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
2782	: *(ph(il,i)-ph(il,i+1))
2783	dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
2784	endif
2785	enddo
2786	enddo
2787
2788	c!!! do 700 i=1,icb(il)-1
2789	do i=1,nl
2790	do il=1,ncum
2791	if (i.le.(icb(il)-1)) then
2792	ft(il,i)=asum(il)t(il,i)/(th(il,i)dsum(il))
2793	fr(il,i)=bsum(il)/csum(il)
2794	endif
2795	enddo
2796	enddo
2797
2798	c
2799	c * reset counter and return *
2800	c
2801	do il=1,ncum
2802	sig(il,nd)=2.0
2803	enddo
2804
2805
2806	do i=1,nd
2807	do il=1,ncum
2808	upwd(il,i)=0.0
2809	dnwd(il,i)=0.0
2810	enddo
2811	enddo
2812
2813	do i=1,nl
2814	do il=1,ncum
2815	dnwd0(il,i)=-mp(il,i)
2816	enddo
2817	enddo
2818	do i=nl+1,nd
2819	do il=1,ncum
2820	dnwd0(il,i)=0.
2821	enddo
2822	enddo
2823
2824
2825	do i=1,nl
2826	do il=1,ncum
2827	if (i.ge.icb(il) .and. i.le.inb(il)) then
2828	upwd(il,i)=0.0
2829	dnwd(il,i)=0.0
2830	endif
2831	enddo
2832	enddo
2833
2834	do i=1,nl
2835	do k=1,nl
2836	do il=1,ncum
2837	up1(il,k,i)=0.0
2838	dn1(il,k,i)=0.0
2839	enddo
2840	enddo
2841	enddo
2842
2843	do i=1,nl
2844	do k=i,nl
2845	do n=1,i-1
2846	do il=1,ncum
2847	if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2848	up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
2849	dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
2850	endif
2851	enddo
2852	enddo
2853	enddo
2854	enddo
2855
2856	do i=2,nl
2857	do k=i,nl
2858	do il=1,ncum
2859	ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2860	if (i.le.inb(il).and.k.le.inb(il)) then
2861	upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
2862	dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
2863	endif
2864	enddo
2865	enddo
2866	enddo
2867
2868
2869	c!!! DO il=1,ncum
2870	c!!! do i=icb(il),inb(il)
2871	c!!!
2872	c!!! upwd(il,i)=0.0
2873	c!!! dnwd(il,i)=0.0
2874	c!!! do k=i,inb(il)
2875	c!!! up1=0.0
2876	c!!! dn1=0.0
2877	c!!! do n=1,i-1
2878	c!!! up1=up1+ment(il,n,k)
2879	c!!! dn1=dn1-ment(il,k,n)
2880	c!!! enddo
2881	c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
2882	c!!! dnwd(il,i)=dnwd(il,i)+dn1
2883	c!!! enddo
2884	c!!! enddo
2885	c!!!
2886	c!!! ENDDO
2887
2888	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2889	c determination de la variation de flux ascendant entre
2890	c deux niveau non dilue mike
2891	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2892
2893	do i=1,nl
2894	do il=1,ncum
2895	mike(il,i)=m(il,i)
2896	enddo
2897	enddo
2898
2899	do i=nl+1,nd
2900	do il=1,ncum
2901	mike(il,i)=0.
2902	enddo
2903	enddo
2904
2905	do i=1,nd
2906	do il=1,ncum
2907	ma(il,i)=0
2908	enddo
2909	enddo
2910
2911	do i=1,nl
2912	do j=i,nl
2913	do il=1,ncum
2914	ma(il,i)=ma(il,i)+m(il,j)
2915	enddo
2916	enddo
2917	enddo
2918
2919	do i=nl+1,nd
2920	do il=1,ncum
2921	ma(il,i)=0.
2922	enddo
2923	enddo
2924
2925	do i=1,nl
2926	do il=1,ncum
2927	if (i.le.(icb(il)-1)) then
2928	ma(il,i)=0
2929	endif
2930	enddo
2931	enddo
2932
2933	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2934	c icb represente de niveau ou se trouve la
2935	c base du nuage , et inb le top du nuage
2936	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2937
2938	do i=1,nd
2939	do il=1,ncum
2940	mke(il,i)=upwd(il,i)+dnwd(il,i)
2941	enddo
2942	enddo
2943
2944	do i=1,nd
2945	DO 999 il=1,ncum
2946	rdcp=(rrd(1.-rr(il,i))-rr(il,i)rrv)
2947	: /(cpd(1.-rr(il,i))+rr(il,i)cpv)
2948	tls(il,i)=t(il,i)(1000.0/p(il,i))*rdcp
2949	tps(il,i)=tp(il,i)
2950	999 CONTINUE
2951	enddo
2952
2953	c
2954	c * diagnose the in-cloud mixing ratio * ! cld
2955	c * of condensed water * ! cld
2956	c ! cld
2957
2958	do i=1,nd ! cld
2959	do il=1,ncum ! cld
2960	mac(il,i)=0.0 ! cld
2961	wa(il,i)=0.0 ! cld
2962	siga(il,i)=0.0 ! cld
2963	sax(il,i)=0.0 ! cld
2964	enddo ! cld
2965	enddo ! cld
2966
2967	do i=minorig, nl ! cld
2968	do k=i+1,nl+1 ! cld
2969	do il=1,ncum ! cld
2970	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld
2971	mac(il,i)=mac(il,i)+m(il,k) ! cld
2972	endif ! cld
2973	enddo ! cld
2974	enddo ! cld
2975	enddo ! cld
2976
2977	do i=1,nl ! cld
2978	do j=1,i ! cld
2979	do il=1,ncum ! cld
2980	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
2981	: .and. j.ge.icb(il) ) then ! cld
2982	sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
2983	: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
2984	endif ! cld
2985	enddo ! cld
2986	enddo ! cld
2987	enddo ! cld
2988
2989	do i=1,nl ! cld
2990	do il=1,ncum ! cld
2991	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
2992	: .and. sax(il,i).gt.0.0 ) then ! cld
2993	wa(il,i)=sqrt(2.*sax(il,i)) ! cld
2994	endif ! cld
2995	enddo ! cld
2996	enddo ! cld
2997
2998	do i=1,nl ! cld
2999	do il=1,ncum ! cld
3000	if (wa(il,i).gt.0.0) ! cld
3001	: siga(il,i)=mac(il,i)/wa(il,i) ! cld
3002	: rrdtvp(il,i)/p(il,i)/100./delta ! cld
3003	siga(il,i) = min(siga(il,i),1.0) ! cld
3004	cIM cf. FH
3005	if (iflag_clw.eq.0) then
3006	qcondc(il,i)=siga(il,i)clw(il,i)(1.-ep(il,i)) ! cld
3007	: + (1.-siga(il,i))*qcond(il,i) ! cld
3008	else if (iflag_clw.eq.1) then
3009	qcondc(il,i)=qcond(il,i) ! cld
3010	endif
3011
3012	enddo ! cld
3013	enddo ! cld
3014
3015	return
3016	end
3017
3018	SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
3019	& ment,sij,da,phi)
3020	implicit none
3021	c inputs:
3022	integer ncum, nd, na, nloc,len
3023	real ment(nloc,na,na),sij(nloc,na,na)
3024	c ouputs:
3025	real da(nloc,na),phi(nloc,na,na)
3026	c local variables:
3027	integer i,j,k
3028	c
3029	da(:,:)=0.
3030	c
3031	do j=1,na
3032	do k=1,na
3033	do i=1,ncum
3034	da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
3035	phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
3036	c print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j)
3037	end do
3038	end do
3039	end do
3040
3041	return
3042	end
3043
3044
3045	SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
3046	: ,iflag
3047	: ,precip,VPrecip,sig,w0
3048	: ,ft,fq,fu,fv,ftra
3049	: ,inb
3050	: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
3051	: ,da,phi,mp
3052	: ,iflag1
3053	: ,precip1,VPrecip1,sig1,w01
3054	: ,ft1,fq1,fu1,fv1,ftra1
3055	: ,inb1
3056	: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
3057	: ,da1,phi1,mp1)
3058	implicit none
3059
3060	#include "cvparam3.h"
3061
3062	c inputs:
3063	integer len, ncum, nd, ntra, nloc
3064	integer idcum(nloc)
3065	integer iflag(nloc)
3066	integer inb(nloc)
3067	real precip(nloc)
3068	real VPrecip(nloc,nd+1)
3069	real sig(nloc,nd), w0(nloc,nd)
3070	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
3071	real ftra(nloc,nd,ntra)
3072	real Ma(nloc,nd)
3073	real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
3074	real qcondc(nloc,nd)
3075	real wd(nloc),cape(nloc)
3076	real da(nloc,nd),phi(nloc,nd,nd),mp(nloc,nd)
3077
3078	c outputs:
3079	integer iflag1(len)
3080	integer inb1(len)
3081	real precip1(len)
3082	real VPrecip1(len,nd+1)
3083	real sig1(len,nd), w01(len,nd)
3084	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
3085	real ftra1(len,nd,ntra)
3086	real Ma1(len,nd)
3087	real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
3088	real qcondc1(nloc,nd)
3089	real wd1(nloc),cape1(nloc)
3090	real da1(nloc,nd),phi1(nloc,nd,nd),mp1(nloc,nd)
3091
3092	c local variables:
3093	integer i,k,j
3094
3095	do 2000 i=1,ncum
3096	precip1(idcum(i))=precip(i)
3097	iflag1(idcum(i))=iflag(i)
3098	wd1(idcum(i))=wd(i)
3099	inb1(idcum(i))=inb(i)
3100	cape1(idcum(i))=cape(i)
3101	2000 continue
3102
3103	do 2020 k=1,nl
3104	do 2010 i=1,ncum
3105	VPrecip1(idcum(i),k)=VPrecip(i,k)
3106	sig1(idcum(i),k)=sig(i,k)
3107	w01(idcum(i),k)=w0(i,k)
3108	ft1(idcum(i),k)=ft(i,k)
3109	fq1(idcum(i),k)=fq(i,k)
3110	fu1(idcum(i),k)=fu(i,k)
3111	fv1(idcum(i),k)=fv(i,k)
3112	Ma1(idcum(i),k)=Ma(i,k)
3113	upwd1(idcum(i),k)=upwd(i,k)
3114	dnwd1(idcum(i),k)=dnwd(i,k)
3115	dnwd01(idcum(i),k)=dnwd0(i,k)
3116	qcondc1(idcum(i),k)=qcondc(i,k)
3117	da1(idcum(i),k)=da(i,k)
3118	mp1(idcum(i),k)=mp(i,k)
3119	2010 continue
3120	2020 continue
3121
3122	do 2200 i=1,ncum
3123	sig1(idcum(i),nd)=sig(i,nd)
3124	2200 continue
3125
3126
3127	c do 2100 j=1,ntra
3128	c do 2110 k=1,nd ! oct3
3129	c do 2120 i=1,ncum
3130	c ftra1(idcum(i),k,j)=ftra(i,k,j)
3131	c 2120 continue
3132	c 2110 continue
3133	c 2100 continue
3134	do j=1,nd
3135	do k=1,nd
3136	do i=1,ncum
3137	phi1(idcum(i),k,j)=phi(i,k,j)
3138	end do
3139	end do
3140	end do
3141
3142	return
3143	end
3144

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

Download in other formats:

Original Format