Context Navigation

cv3_routines.F @ 597

Last change on this file since 597 was 597, checked in by Laurent Fairhead, 20 years ago
Ne passe pas en parallele YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 92.4 KB

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

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

Download in other formats:

Original Format