Context Navigation

cv3_routines.F @ 486

Last change on this file since 486 was 486, checked in by lmdzadmin, 20 years ago
Phasage avec la version de Ionela IM/LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 91.3 KB

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

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

Download in other formats:

Original Format