Context Navigation

cv3_routines.F @ 3287

Last change on this file since 3287 was 1554, checked in by Laurent Fairhead, 13 years ago

Changes in some default valeus of parameters

Modifications des valeurs par défaut de certains paramètres

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 109.7 KB

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

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

Download in other formats:

Original Format