Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

cv3_routines.F @ 4501

Last change on this file since 4501 was 416, checked in by lmdzadmin, 22 years ago
Inclusion initiale
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 90.5 KB

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

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

Download in other formats:

Original Format