Context Navigation

cv3_routines.F @ 1046

Last change on this file since 1046 was 1044, checked in by lmdzadmin, 16 years ago
Correction oscillations verticalles de la pluie lorsque celle-ci s'evapore en totalite (correction a la programmation originale de Emanuel) JYG/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 101.8 KB

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

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

Download in other formats:

Original Format