Context Navigation

cv3_routines.F @ 1842

Last change on this file since 1842 was 1774, checked in by Laurent Fairhead, 11 years ago

Initialisation de variables dans la convection

JYG

Missing initialisations in convection

JYG

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

File size: 114.2 KB

Line
1	!
2	! $Id: cv3_routines.F 1774 2013-06-18 12:02:30Z idelkadi $
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	!AC! do 121 j=1,ntra
882	!AC!ccccc do 111 k=1,nl+1
883	!AC! do 111 k=1,nd
884	!AC! nn=0
885	!AC! do 101 i=1,len
886	!AC! if(iflag1(i).eq.0)then
887	!AC! nn=nn+1
888	!AC! tra(nn,k,j)=tra1(i,k,j)
889	!AC! endif
890	!AC! 101 continue
891	!AC! 111 continue
892	!AC! 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	!AC! do k=1,ntra
1636	!AC! do j=1,nd ! instead nlp
1637	!AC! do i=1,nd ! instead nlp
1638	!AC! do il=1,ncum
1639	!AC! traent(il,i,j,k)=tra(il,j,k)
1640	!AC! enddo
1641	!AC! enddo
1642	!AC! enddo
1643	!AC! 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	!AC! do k=1,ntra
1700	!AC! do j=minorig,nl
1701	!AC! do il=1,ncum
1702	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1703	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1704	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1705	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1706	!AC! endif
1707	!AC! enddo
1708	!AC! enddo
1709	!AC! 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	!AC! do j=1,ntra
1733	!AC! do i=minorig+1,nl
1734	!AC! do il=1,ncum
1735	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1736	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
1737	!AC! endif
1738	!AC! enddo
1739	!AC! enddo
1740	!AC! 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	!AC! do j=1,ntra
1907	!AC! do il=1,ncum
1908	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1909	!AC! : .and. csum(il,i).lt.m(il,i) ) then
1910	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
1911	!AC! endif
1912	!AC! enddo
1913	!AC! 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	o ,wdtrainA,wdtrainM) ! RomP
1954	implicit none
1955
1956
1957	#include "cvthermo.h"
1958	#include "cv3param.h"
1959	#include "cvflag.h"
1960
1961	c inputs:
1962	integer ncum, nd, na, ntra, nloc
1963	integer icb(nloc), inb(nloc)
1964	real delt, plcl(nloc)
1965	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd),gz(nloc,na)
1966	real u(nloc,nd), v(nloc,nd)
1967	real tra(nloc,nd,ntra)
1968	real p(nloc,nd), ph(nloc,nd+1)
1969	real ep(nloc,na), sigp(nloc,na), clw(nloc,na)
1970	real th(nloc,na),tv(nloc,na),lv(nloc,na),cpn(nloc,na)
1971	real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
1972	real coef_clos(nloc)
1973	c
1974	c input/output
1975	integer iflag(nloc)
1976	c
1977	c outputs:
1978	real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
1979	real water(nloc,na), evap(nloc,na), wt(nloc,na)
1980	real trap(nloc,na,ntra)
1981	real b(nloc,na), sigd(nloc)
1982	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
1983	! lascendance adiabatique et des flux melanges Pa et Pm.
1984	! Distinction des wdtrain
1985	! Pa = wdtrainA Pm = wdtrainM
1986	real wdtrainA(nloc,na), wdtrainM(nloc,na)
1987
1988	c local variables
1989	integer i,j,k,il,num1,ndp1
1990	real tinv, delti
1991	real awat, afac, afac1, afac2, bfac
1992	real pr1, pr2, sigt, b6, c6, revap, delth
1993	real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1994	real ampmax
1995	real tevap(nloc)
1996	real lvcp(nloc,na)
1997	real h(nloc,na),hm(nloc,na)
1998	real wdtrain(nloc)
1999	logical lwork(nloc),mplus(nloc)
2000
2001
2002	c------------------------------------------------------
2003
2004	delti = 1./delt
2005	tinv=1./3.
2006
2007	mp(:,:)=0.
2008
2009	do i=1,nl
2010	do il=1,ncum
2011	mp(il,i)=0.0
2012	rp(il,i)=rr(il,i)
2013	up(il,i)=u(il,i)
2014	vp(il,i)=v(il,i)
2015	wt(il,i)=0.001
2016	water(il,i)=0.0
2017	evap(il,i)=0.0
2018	b(il,i)=0.0
2019	lvcp(il,i)=lv(il,i)/cpn(il,i)
2020	enddo
2021	enddo
2022	!AC! do k=1,ntra
2023	!AC! do i=1,nd
2024	!AC! do il=1,ncum
2025	!AC! trap(il,i,k)=tra(il,i,k)
2026	!AC! enddo
2027	!AC! enddo
2028	!AC! enddo
2029	!! RomP >>>
2030	do i=1,nd
2031	do il=1,ncum
2032	wdtrainA(il,i)=0.0
2033	wdtrainM(il,i)=0.0
2034	enddo
2035	enddo
2036	!! RomP <<<
2037	c
2038	c * check whether ep(inb)=0, if so, skip precipitating *
2039	c * downdraft calculation *
2040	c
2041
2042	do il=1,ncum
2043	!! lwork(il)=.TRUE.
2044	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
2045	lwork(il)= ep(il,inb(il)) .ge. 0.0001
2046	enddo
2047
2048	c *** Set the fractionnal area sigd of precipitating downdraughts
2049	do il = 1,ncum
2050	sigd(il) = sigdz*coef_clos(il)
2051	enddo
2052
2053
2054	c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2055	c
2056	c * begin downdraft loop *
2057	c
2058	c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2059	c
2060	DO 400 i=nl+1,1,-1
2061
2062	num1=0
2063	do il=1,ncum
2064	if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1
2065	enddo
2066	if (num1.le.0) goto 400
2067
2068	call zilch(wdtrain,ncum)
2069
2070	c
2071	c * integrate liquid water equation to find condensed water *
2072	c * and condensed water flux *
2073	c
2074	c
2075	c * calculate detrained precipitation *
2076	c
2077	do il=1,ncum
2078	if (i.le.inb(il) .and. lwork(il)) then
2079	if (cvflag_grav) then
2080	wdtrain(il)=gravep(il,i)m(il,i)*clw(il,i)
2081	wdtrainA(il,i) = wdtrain(il)/grav ! Pa RomP
2082	else
2083	wdtrain(il)=10.0ep(il,i)m(il,i)*clw(il,i)
2084	wdtrainA(il,i) = wdtrain(il)/10. ! Pa RomP
2085	endif
2086	endif
2087	enddo
2088
2089	if(i.gt.1)then
2090	do 320 j=1,i-1
2091	do il=1,ncum
2092	if (i.le.inb(il) .and. lwork(il)) then
2093	awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
2094	awat=max(awat,0.0)
2095	if (cvflag_grav) then
2096	wdtrain(il)=wdtrain(il)+gravawatment(il,j,i)
2097	wdtrainM(il,i) = wdtrain(il)/grav-wdtrainA(il,i) ! Pm RomP
2098	else
2099	wdtrain(il)=wdtrain(il)+10.0awatment(il,j,i)
2100	wdtrainM(il,i) = wdtrain(il)/10.-wdtrainA(il,i) ! Pm RomP
2101	endif
2102	endif
2103	enddo
2104	320 continue
2105	endif
2106
2107	c
2108	c * find rain water and evaporation using provisional *
2109	c * estimates of rp(i)and rp(i-1) *
2110	c
2111
2112	do 995 il=1,ncum
2113	if (i.le.inb(il) .and. lwork(il)) then
2114
2115	wt(il,i)=45.0
2116
2117	if(i.lt.inb(il))then
2118	rp(il,i)=rp(il,i+1)
2119	: +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
2120	rp(il,i)=0.5*(rp(il,i)+rr(il,i))
2121	endif
2122	rp(il,i)=max(rp(il,i),0.0)
2123	rp(il,i)=amin1(rp(il,i),rs(il,i))
2124	rp(il,inb(il))=rr(il,inb(il))
2125
2126	if(i.eq.1)then
2127	afac=p(il,1)(rs(il,1)-rp(il,1))/(1.0e4+2000.0p(il,1)*rs(il,1))
2128	else
2129	rp(il,i-1)=rp(il,i)
2130	: +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i)
2131	rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
2132	rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
2133	rp(il,i-1)=max(rp(il,i-1),0.0)
2134	afac1=p(il,i)(rs(il,i)-rp(il,i))/(1.0e4+2000.0p(il,i)*rs(il,i))
2135	afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
2136	: /(1.0e4+2000.0p(il,i-1)rs(il,i-1))
2137	afac=0.5*(afac1+afac2)
2138	endif
2139	if(i.eq.inb(il))afac=0.0
2140	afac=max(afac,0.0)
2141	bfac=1./(sigd(il)*wt(il,i))
2142	c
2143	cjyg1
2144	ccc sigt=1.0
2145	ccc if(i.ge.icb)sigt=sigp(i)
2146	c prise en compte de la variation progressive de sigt dans
2147	c les couches icb et icb-1:
2148	c pour plcl<ph(i+1), pr1=0 & pr2=1
2149	c pour plcl>ph(i), pr1=1 & pr2=0
2150	c pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2151	c sur le nuage, et pr2 est la proportion sous la base du
2152	c nuage.
2153	pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2154	pr1=max(0.,min(1.,pr1))
2155	pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2156	pr2=max(0.,min(1.,pr2))
2157	sigt=sigp(il,i)*pr1+pr2
2158	cjyg2
2159	c
2160	cjyg----
2161	c b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2162	c c6 = water(il,i+1) + wdtrain(il)*bfac
2163	c revap=0.5(-b6+sqrt(b6b6+4.*c6))
2164	c evap(il,i)=sigtafacrevap
2165	c water(il,i)=revap*revap
2166	cc print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ',
2167	cc $ i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
2168	cc---end jyg---
2169	c
2170	c--------retour à la formulation originale d'Emanuel.
2171	b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2172	c6=water(il,i+1)+bfac*wdtrain(il)
2173	: -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2174	if(c6.gt.0.0)then
2175	revap=0.5(-b6+sqrt(b6b6+4.*c6))
2176	cjyg Dans sa formulation originale, Emanuel calcule l'evaporation par:
2177	cc evap(il,i)=sigtafacrevap
2178	c ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
2179	c Ici,l'evaporation evap est simplement calculee par l'equation de
2180	c conservation.
2181	water(il,i)=revap*revap
2182	else
2183	cjyg---- Correction : si c6 <= 0, water(il,i)=0.
2184	water(il,i) = 0.
2185	endif
2186	cJYG/IM : ci-dessous formulation originale de KE
2187	c evap(il,i)=-evap(il,i+1)
2188	c : +(wdtrain(il)+sigd(il)wt(il,i)water(il,i+1))
2189	c : /(sigd(il)(ph(il,i)-ph(il,i+1))50.)
2190	c
2191	cJYG/IM : ci-dessous modification formulation originale de KE
2192	c pour eliminer oscillations verticales de pluie se produisant
2193	c lorsqu'il y a evaporation totale de la pluie
2194	c
2195	c evap(il,i)= +(wdtrain(il)+sigd(il)wt(il,i)water(il,i+1)) !itlmd(jyg)
2196	c : /(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2197	c end if !itlmd(jyg)
2198	cjyg--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
2199	c moins [tt ce qui sort de la couche i]
2200	evap(il,i)=
2201	: (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))
2202	: /(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2203	c
2204	endif !(i.le.inb(il) .and. lwork(il))
2205	995 Continue
2206	c----------------------------------------------------------------
2207	c
2208	ccc
2209	c * calculate precipitating downdraft mass flux under *
2210	c * hydrostatic approximation *
2211	c
2212	Do 996 il = 1,ncum
2213	if (i.le.inb(il) .and. lwork(il) .and. i.ne.1) then
2214	c
2215	tevap(il)=max(0.0,evap(il,i))
2216	delth=max(0.001,(th(il,i)-th(il,i-1)))
2217	if (cvflag_grav) then
2218	mp(il,i)=100.ginvlvcp(il,i)sigd(il)tevap(il)
2219	: *(p(il,i-1)-p(il,i))/delth
2220	else
2221	mp(il,i)=10.lvcp(il,i)sigd(il)*tevap(il)
2222	: *(p(il,i-1)-p(il,i))/delth
2223	endif
2224	c
2225	endif !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2226	996 Continue
2227	c----------------------------------------------------------------
2228	c
2229	c * if hydrostatic assumption fails, *
2230	c * solve cubic difference equation for downdraft theta *
2231	c * and mass flux from two simultaneous differential eqns *
2232	c
2233	Do 997 il = 1,ncum
2234	if (i.le.inb(il) .and. lwork(il) .and. i.ne.1) then
2235	c
2236	amfac=sigd(il)sigd(il)70.0ph(il,i)(p(il,i-1)-p(il,i))
2237	: (th(il,i)-th(il,i-1))/(tv(il,i)th(il,i))
2238	amp2=abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2239	c
2240	if(amp2.gt.(0.1*amfac))then
2241	xf=100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
2242	tf=b(il,i)-5.0(th(il,i)-th(il,i-1))t(il,i)
2243	: /(lvcp(il,i)sigd(il)th(il,i))
2244	af=xftf+mp(il,i+1)mp(il,i+1)*tinv
2245	bf=2.(tinvmp(il,i+1))*3+tinvmp(il,i+1)xftf
2246	: +50.(p(il,i-1)-p(il,i))xf*tevap(il)
2247	fac2=1.0
2248	if(bf.lt.0.0)fac2=-1.0
2249	bf=abs(bf)
2250	ur=0.25bfbf-afafaftinvtinv*tinv
2251	if(ur.ge.0.0)then
2252	sru=sqrt(ur)
2253	fac=1.0
2254	if((0.5*bf-sru).lt.0.0)fac=-1.0
2255	mp(il,i)=mp(il,i+1)tinv+(0.5bf+sru)**tinv
2256	: +fac(abs(0.5bf-sru))**tinv
2257	else
2258	d=atan(2.*sqrt(-ur)/(bf+1.0e-28))
2259	if(fac2.lt.0.0)d=3.14159-d
2260	mp(il,i)=mp(il,i+1)tinv+2.sqrt(aftinv)cos(d*tinv)
2261	endif
2262	mp(il,i)=max(0.0,mp(il,i))
2263
2264	if (cvflag_grav) then
2265	Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2266	C il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
2267	C Et il faut bien revoir les facteurs 100.
2268	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap(il)
2269	2 /(mp(il,i)+sigd(il)*0.1)
2270	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)
2271	: sigd(il)th(il,i))
2272	else
2273	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap(il)
2274	2 /(mp(il,i)+sigd(il)*0.1)
2275	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)
2276	: sigd(il)th(il,i))
2277	endif
2278	b(il,i-1)=max(b(il,i-1),0.0)
2279	c
2280	endif !(amp2.gt.(0.1*amfac))
2281	c
2282	c * limit magnitude of mp(i) to meet cfl condition *
2283	c
2284	ampmax=2.0(ph(il,i)-ph(il,i+1))delti
2285	amp2=2.0(ph(il,i-1)-ph(il,i))delti
2286	ampmax=min(ampmax,amp2)
2287	mp(il,i)=min(mp(il,i),ampmax)
2288	c
2289	c * force mp to decrease linearly to zero *
2290	c * between cloud base and the surface *
2291	c
2292	c
2293	cc if(p(il,i).gt.p(il,icb(il)))then
2294	cc mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2295	cc endif
2296	if(ph(il,i) .gt. 0.9*plcl(il)) then
2297	mp(il,i) = mp(il,i)*(ph(il,1)-ph(il,i))/
2298	$ (ph(il,1)-0.9*plcl(il))
2299	endif
2300
2301	endif ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2302	997 Continue
2303	c----------------------------------------------------------------
2304	c
2305	c * find mixing ratio of precipitating downdraft *
2306	c
2307	Do il = 1,ncum
2308	if (i.lt.inb(il) .and. lwork(il)) then
2309	mplus(il) = mp(il,i).gt.mp(il,i+1)
2310	endif ! (i.lt.inb(il) .and. lwork(il))
2311	enddo
2312	c
2313	Do 999 il = 1,ncum
2314	if (i.lt.inb(il) .and. lwork(il)) then
2315	c
2316	rp(il,i)=rr(il,i)
2317
2318	if(mplus(il))then
2319
2320	if (cvflag_grav) then
2321	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2322	: +100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))
2323	: *(evap(il,i+1)+evap(il,i))
2324	else
2325	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2326	: +5.sigd(il)(ph(il,i)-ph(il,i+1))
2327	: *(evap(il,i+1)+evap(il,i))
2328	endif
2329	rp(il,i)=rp(il,i)/mp(il,i)
2330	up(il,i)=up(il,i+1)mp(il,i+1)+u(il,i)(mp(il,i)-mp(il,i+1))
2331	up(il,i)=up(il,i)/mp(il,i)
2332	vp(il,i)=vp(il,i+1)mp(il,i+1)+v(il,i)(mp(il,i)-mp(il,i+1))
2333	vp(il,i)=vp(il,i)/mp(il,i)
2334
2335	else ! if (mplus(il))
2336
2337	if(mp(il,i+1).gt.1.0e-16)then
2338	if (cvflag_grav) then
2339	rp(il,i)=rp(il,i+1)
2340	: +100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))
2341	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2342	else
2343	rp(il,i)=rp(il,i+1)
2344	: +5.sigd(il)(ph(il,i)-ph(il,i+1))
2345	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2346	endif
2347	up(il,i)=up(il,i+1)
2348	vp(il,i)=vp(il,i+1)
2349	endif ! (mp(il,i+1).gt.1.0e-16)
2350	endif ! (mplus(il)) else if (.not.mplus(il))
2351	c
2352	rp(il,i)=amin1(rp(il,i),rs(il,i))
2353	rp(il,i)=max(rp(il,i),0.0)
2354
2355	endif ! (i.lt.inb(il) .and. lwork(il))
2356	999 continue
2357	c----------------------------------------------------------------
2358	c
2359	c * find tracer concentrations in precipitating downdraft *
2360	c
2361	!AC! do j=1,ntra
2362	!AC! do il = 1,ncum
2363	!AC! if (i.lt.inb(il) .and. lwork(il)) then
2364	!AC!c
2365	!AC! if(mplus(il))then
2366	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2367	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2368	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2369	!AC! else ! if (mplus(il))
2370	!AC! if(mp(il,i+1).gt.1.0e-16)then
2371	!AC! trap(il,i,j)=trap(il,i+1,j)
2372	!AC! endif
2373	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
2374	!AC!c
2375	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
2376	!AC! enddo
2377	!AC! end do
2378
2379	400 continue
2380	c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2381	c
2382	c * end of downdraft loop *
2383	c
2384	c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2385	c
2386
2387	return
2388	end
2389
2390	SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
2391	: ,icb,inb,delt
2392	: ,t,rr,t_wake,rr_wake,s_wake,u,v,tra
2393	: ,gz,p,ph,h,hp,lv,cpn,th,th_wake
2394	: ,ep,clw,m,tp,mp,rp,up,vp,trap
2395	: ,wt,water,evap,b,sigd
2396	: ,ment,qent,hent,iflag_mix,uent,vent
2397	: ,nent,elij,traent,sig
2398	: ,tv,tvp,wghti
2399	: ,iflag,precip,Vprecip,ft,fr,fu,fv,ftra
2400	: ,cbmf,upwd,dnwd,dnwd0,ma,mip
2401	: ,tls,tps,qcondc,wd
2402	: ,ftd,fqd)
2403
2404	implicit none
2405
2406	#include "cvthermo.h"
2407	#include "cv3param.h"
2408	#include "cvflag.h"
2409	#include "conema3.h"
2410
2411	c inputs:
2412	c print*,'cv3_yield apres include'
2413	integer iflag_mix
2414	integer ncum,nd,na,ntra,nloc
2415	integer icb(nloc), inb(nloc)
2416	real delt
2417	real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
2418	real t_wake(nloc,nd), rr_wake(nloc,nd)
2419	real s_wake(nloc)
2420	real tra(nloc,nd,ntra), sig(nloc,nd)
2421	real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
2422	real th(nloc,na), p(nloc,nd), tp(nloc,na)
2423	real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
2424	real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
2425	real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
2426	real water(nloc,na), evap(nloc,na), b(nloc,na), sigd(nloc)
2427	real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
2428	real hent(nloc,na,na)
2429	cIM bug real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2430	real vent(nloc,na,na), elij(nloc,na,na)
2431	integer nent(nloc,nd)
2432	real traent(nloc,na,na,ntra)
2433	real tv(nloc,nd), tvp(nloc,nd), wghti(nloc,nd)
2434	c print*,'cv3_yield declarations 1'
2435	c input/output:
2436	integer iflag(nloc)
2437
2438	c outputs:
2439	real precip(nloc)
2440	real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
2441	real ftd(nloc,nd), fqd(nloc,nd)
2442	real ftra(nloc,nd,ntra)
2443	real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
2444	real dnwd0(nloc,nd), mip(nloc,nd)
2445	real Vprecip(nloc,nd+1)
2446	real tls(nloc,nd), tps(nloc,nd)
2447	real qcondc(nloc,nd) ! cld
2448	real wd(nloc) ! gust
2449	real cbmf(nloc)
2450	c print*,'cv3_yield declarations 2'
2451	c local variables:
2452	integer i,k,il,n,j,num1
2453	real rat, delti
2454	real ax, bx, cx, dx, ex
2455	real cpinv, rdcp, dpinv
2456	real awat(nloc)
2457	real lvcp(nloc,na), mke(nloc,na)
2458	real am(nloc), work(nloc), ad(nloc), amp1(nloc)
2459	c!! real up1(nloc), dn1(nloc)
2460	real up1(nloc,nd,nd), dn1(nloc,nd,nd)
2461	real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2462	real esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
2463	real th_wake(nloc,nd)
2464	real alpha_qpos(nloc),alpha_qpos1(nloc)
2465	real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
2466	real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
2467
2468	c print*,'cv3_yield declarations 3'
2469	c-------------------------------------------------------------
2470
2471	c initialization:
2472
2473	delti = 1.0/delt
2474	c print*,'cv3_yield initialisation delt', delt
2475	cprecip,Vprecip,ft,fr,fu,fv,ftra
2476	c : ,cbmf,upwd,dnwd,dnwd0,ma,mip
2477	c : ,tls,tps,qcondc,wd
2478	c : ,ftd,fqd )
2479	do il=1,ncum
2480	precip(il)=0.0
2481	Vprecip(il,nd+1)=0.0
2482	wd(il)=0.0 ! gust
2483	enddo
2484
2485	do i=1,nd
2486	do il=1,ncum
2487	Vprecip(il,i)=0.0
2488	ft(il,i)=0.0
2489	fr(il,i)=0.0
2490	fu(il,i)=0.0
2491	fv(il,i)=0.0
2492	upwd(il,i)=0.0
2493	dnwd(il,i)=0.0
2494	dnwd0(il,i)=0.0
2495	mip(il,i)=0.0
2496	ftd(il,i)=0.0
2497	fqd(il,i)=0.0
2498	qcondc(il,i)=0.0 ! cld
2499	qcond(il,i)=0.0 ! cld
2500	nqcond(il,i)=0.0 ! cld
2501	enddo
2502	enddo
2503	c print*,'cv3_yield initialisation 2'
2504	!AC! do j=1,ntra
2505	!AC! do i=1,nd
2506	!AC! do il=1,ncum
2507	!AC! ftra(il,i,j)=0.0
2508	!AC! enddo
2509	!AC! enddo
2510	!AC! enddo
2511	c print*,'cv3_yield initialisation 3'
2512	do i=1,nl
2513	do il=1,ncum
2514	lvcp(il,i)=lv(il,i)/cpn(il,i)
2515	enddo
2516	enddo
2517
2518
2519	c
2520	c * calculate surface precipitation in mm/day *
2521	c
2522	do il=1,ncum
2523	if(ep(il,inb(il)).ge.0.0001 .and. iflag(il) .le. 1)then
2524	if (cvflag_grav) then
2525	precip(il)=wt(il,1)sigd(il)water(il,1)86400.1000.
2526	: /(rowl*grav)
2527	else
2528	precip(il)=wt(il,1)sigd(il)water(il,1)*8640.
2529	endif
2530	endif
2531	enddo
2532	c print*,'cv3_yield apres calcul precip'
2533
2534	C
2535	C === calculate vertical profile of precipitation in kg/m2/s ===
2536	C
2537	do i = 1,nl
2538	do il=1,ncum
2539	if(ep(il,inb(il)).ge.0.0001 .and. i.le.inb(il)
2540	: .and. iflag(il) .le. 1)then
2541	if (cvflag_grav) then
2542	VPrecip(il,i) = wt(il,i)sigd(il)water(il,i)/grav
2543	else
2544	VPrecip(il,i) = wt(il,i)sigd(il)water(il,i)/10.
2545	endif
2546	endif
2547	enddo
2548	enddo
2549	C
2550	c
2551	c * Calculate downdraft velocity scale *
2552	c * NE PAS UTILISER POUR L'INSTANT *
2553	c
2554	c! do il=1,ncum
2555	c! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
2556	c! : /(sigd(il)*p(il,icb(il)))
2557	c! enddo
2558
2559	c
2560	c * calculate tendencies of lowest level potential temperature *
2561	c * and mixing ratio *
2562	c
2563	do il=1,ncum
2564	work(il)=1.0/(ph(il,1)-ph(il,2))
2565	cbmf(il)=0.0
2566	enddo
2567
2568	do k=2,nl
2569	do il=1,ncum
2570	if (k.ge.icb(il)) then
2571	cbmf(il)=cbmf(il)+m(il,k)
2572	endif
2573	enddo
2574	enddo
2575
2576	c print*,'cv3_yield avant ft'
2577	c AM is the part of cbmf taken from the first level
2578	do il=1,ncum
2579	am(il)=cbmf(il)*wghti(il,1)
2580	enddo
2581	c
2582	do il=1,ncum
2583	if (iflag(il) .le. 1) then
2584	c convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
2585	cjyg Correction pour conserver l'eau
2586	ccc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
2587	ft(il,1)=-lvcp(il,1)sigd(il)evap(il,1) !precip
2588
2589	if (cvflag_grav) then
2590	ft(il,1)=ft(il,1)-0.009gravsigd(il)*mp(il,2)
2591	: t_wake(il,1)b(il,1)*work(il)
2592	else
2593	ft(il,1)=ft(il,1)-0.09sigd(il)mp(il,2)
2594	: t_wake(il,1)b(il,1)*work(il)
2595	endif
2596
2597	ft(il,1)=ft(il,1)+0.01sigd(il)wt(il,1)(cl-cpd)water(il,2)
2598	: (t_wake(il,2)-t_wake(il,1))work(il)/cpn(il,1)
2599
2600	ftd(il,1) = ft(il,1) ! fin precip
2601
2602	if (cvflag_grav) then !sature
2603	if((0.01gravwork(il)*am(il)).ge.delti)iflag(il)=1!consist vect
2604	ft(il,1)=ft(il,1)+0.01gravwork(il)am(il)(t(il,2)-t(il,1)
2605	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2606	else
2607	if((0.1work(il)am(il)).ge.delti)iflag(il)=1 !consistency vect
2608	ft(il,1)=ft(il,1)+0.1work(il)am(il)*(t(il,2)-t(il,1)
2609	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2610	endif
2611	endif ! iflag
2612	enddo
2613
2614
2615	do j=2,nl
2616	IF (iflag_mix .gt. 0) then
2617	do il=1,ncum
2618	c FH WARNING a modifier :
2619	cpinv=0.
2620	c cpinv=1.0/cpn(il,1)
2621	if (j.le.inb(il) .and. iflag(il) .le. 1) then
2622	if (cvflag_grav) then
2623	ft(il,1)=ft(il,1)
2624	: +0.01gravwork(il)ment(il,j,1)(hent(il,j,1)-h(il,1)
2625	: +t(il,1)(cpv-cpd)(rr(il,1)-Qent(il,j,1)))*cpinv
2626	else
2627	ft(il,1)=ft(il,1)
2628	: +0.1work(il)ment(il,j,1)*(hent(il,j,1)-h(il,1)
2629	: +t(il,1)(cpv-cpd)(rr(il,1)-Qent(il,j,1)))*cpinv
2630	endif ! cvflag_grav
2631	endif ! j
2632	enddo
2633	ENDIF
2634	enddo
2635	! fin sature
2636
2637
2638	do il=1,ncum
2639	if (iflag(il) .le. 1) then
2640	if (cvflag_grav) then
2641	Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
2642	fr(il,1)=0.01gravmp(il,2)(rp(il,2)-rr_wake(il,1))work(il)
2643	: +sigd(il)*evap(il,1)
2644	ccc : +sigd(il)0.5(evap(il,1)+evap(il,2))
2645
2646	fqd(il,1)=fr(il,1) !precip
2647
2648	fr(il,1)=fr(il,1)+0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
2649
2650	fu(il,1)=fu(il,1)+0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1))
2651	: +am(il)*(u(il,2)-u(il,1)))
2652	fv(il,1)=fv(il,1)+0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1))
2653	: +am(il)*(v(il,2)-v(il,1)))
2654	else ! cvflag_grav
2655	fr(il,1)=0.1mp(il,2)(rp(il,2)-rr_wake(il,1))*work(il)
2656	: +sigd(il)*evap(il,1)
2657	ccc : +sigd(il)0.5(evap(il,1)+evap(il,2))
2658	fqd(il,1)=fr(il,1) !precip
2659	fr(il,1)=fr(il,1)+0.1am(il)(rr(il,2)-rr(il,1))*work(il)
2660	fu(il,1)=fu(il,1)+0.1work(il)(mp(il,2)*(up(il,2)-u(il,1))
2661	: +am(il)*(u(il,2)-u(il,1)))
2662	fv(il,1)=fv(il,1)+0.1work(il)(mp(il,2)*(vp(il,2)-v(il,1))
2663	: +am(il)*(v(il,2)-v(il,1)))
2664	endif ! cvflag_grav
2665	endif ! iflag
2666	enddo ! il
2667
2668
2669	!AC! do j=1,ntra
2670	!AC! do il=1,ncum
2671	!AC! if (iflag(il) .le. 1) then
2672	!AC! if (cvflag_grav) then
2673	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
2674	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2675	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2676	!AC! else
2677	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2678	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2679	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2680	!AC! endif
2681	!AC! endif ! iflag
2682	!AC! enddo
2683	!AC! enddo
2684
2685	do j=2,nl
2686	do il=1,ncum
2687	if (j.le.inb(il) .and. iflag(il) .le. 1) then
2688	if (cvflag_grav) then
2689	fr(il,1)=fr(il,1)
2690	: +0.01gravwork(il)ment(il,j,1)(qent(il,j,1)-rr(il,1))
2691	fu(il,1)=fu(il,1)
2692	: +0.01gravwork(il)ment(il,j,1)(uent(il,j,1)-u(il,1))
2693	fv(il,1)=fv(il,1)
2694	: +0.01gravwork(il)ment(il,j,1)(vent(il,j,1)-v(il,1))
2695	else ! cvflag_grav
2696	fr(il,1)=fr(il,1)
2697	: +0.1work(il)ment(il,j,1)*(qent(il,j,1)-rr(il,1))
2698	fu(il,1)=fu(il,1)
2699	: +0.1work(il)ment(il,j,1)*(uent(il,j,1)-u(il,1))
2700	fv(il,1)=fv(il,1)
2701	: +0.1work(il)ment(il,j,1)*(vent(il,j,1)-v(il,1)) ! fin sature
2702	endif ! cvflag_grav
2703	endif ! j
2704	enddo
2705	enddo
2706
2707	!AC! do k=1,ntra
2708	!AC! do j=2,nl
2709	!AC! do il=1,ncum
2710	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
2711	!AC!
2712	!AC! if (cvflag_grav) then
2713	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
2714	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
2715	!AC! else
2716	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
2717	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
2718	!AC! endif
2719	!AC!
2720	!AC! endif
2721	!AC! enddo
2722	!AC! enddo
2723	!AC! enddo
2724	c print*,'cv3_yield apres ft'
2725	c
2726	c * calculate tendencies of potential temperature and mixing ratio *
2727	c * at levels above the lowest level *
2728	c
2729	c * first find the net saturated updraft and downdraft mass fluxes *
2730	c * through each level *
2731	c
2732
2733	do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
2734
2735	num1=0
2736	do il=1,ncum
2737	if(i.le.inb(il) .and. iflag(il) .le. 1)num1=num1+1
2738	enddo
2739	if(num1.le.0)go to 500
2740
2741	call zilch(amp1,ncum)
2742	call zilch(ad,ncum)
2743
2744	do 440 k=1,nl+1
2745	do 441 il=1,ncum
2746	if(i.ge.icb(il)) then
2747	if(k.ge.i+1.and. k.le.(inb(il)+1)) then
2748	amp1(il)=amp1(il)+m(il,k)
2749	endif
2750	else
2751	c AMP1 is the part of cbmf taken from layers I and lower
2752	if(k.le.i) then
2753	amp1(il)=amp1(il)+cbmf(il)*wghti(il,k)
2754	endif
2755	endif
2756	441 continue
2757	440 continue
2758
2759	do 450 k=1,i
2760	do 451 j=i+1,nl+1
2761	do 452 il=1,ncum
2762	if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
2763	amp1(il)=amp1(il)+ment(il,k,j)
2764	endif
2765	452 continue
2766	451 continue
2767	450 continue
2768
2769	do 470 k=1,i-1
2770	do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
2771	do 472 il=1,ncum
2772	if (i.le.inb(il) .and. j.le.inb(il)) then
2773	ad(il)=ad(il)+ment(il,j,k)
2774	endif
2775	472 continue
2776	471 continue
2777	470 continue
2778
2779	do 1350 il=1,ncum
2780	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2781	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2782	cpinv=1.0/cpn(il,i)
2783
2784	c convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
2785	if (cvflag_grav) then
2786	if((0.01gravdpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
2787	else
2788	if((0.1dpinvamp1(il)).ge.delti)iflag(il)=1 ! vecto
2789	endif
2790
2791	! precip
2792	ccc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
2793	ft(il,i)= -sigd(il)lvcp(il,i)evap(il,i)
2794	rat=cpn(il,i-1)*cpinv
2795	c
2796	if (cvflag_grav) then
2797	ft(il,i)=ft(il,i)-0.009gravsigd(il)
2798	: (mp(il,i+1)t_wake(il,i)*b(il,i)
2799	: -mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
2800	ft(il,i)=ft(il,i)+0.01sigd(il)wt(il,i)(cl-cpd)water(il,i+1)
2801	: (t_wake(il,i+1)-t_wake(il,i))dpinv*cpinv
2802	ftd(il,i)=ft(il,i)
2803	! fin precip
2804	c
2805	! sature
2806	ft(il,i)=ft(il,i)+0.01gravdpinv(amp1(il)(t(il,i+1)-t(il,i)
2807	: +(gz(il,i+1)-gz(il,i))*cpinv)
2808	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2809
2810	c
2811	IF (iflag_mix .eq. 0) then
2812	ft(il,i)=ft(il,i)+0.01gravdpinvment(il,i,i)(hp(il,i)-h(il,i)
2813	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2814	ENDIF
2815	c
2816	else ! cvflag_grav
2817	ft(il,i)=ft(il,i)-0.09*sigd(il)
2818	: (mp(il,i+1)t_wake(il,i)*b(il,i)
2819	: -mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
2820	ft(il,i)=ft(il,i)+0.01sigd(il)wt(il,i)(cl-cpd)water(il,i+1)
2821	: (t_wake(il,i+1)-t_wake(il,i))dpinv*cpinv
2822	ftd(il,i)=ft(il,i)
2823	! fin precip
2824	c
2825	! sature
2826	ft(il,i)=ft(il,i)+0.1dpinv(amp1(il)*(t(il,i+1)-t(il,i)
2827	: +(gz(il,i+1)-gz(il,i))*cpinv)
2828	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2829
2830	c
2831	IF (iflag_mix .eq. 0) then
2832	ft(il,i)=ft(il,i)+0.1dpinvment(il,i,i)*(hp(il,i)-h(il,i)
2833	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2834	ENDIF
2835	endif ! cvflag_grav
2836
2837
2838	if (cvflag_grav) then
2839	c sb: on ne fait pas encore la correction permettant de mieux
2840	c conserver l'eau:
2841	c jyg: correction permettant de mieux conserver l'eau:
2842	ccc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
2843	fr(il,i)=sigd(il)*evap(il,i)
2844	: +0.01grav(mp(il,i+1)*(rp(il,i+1)-rr_wake(il,i))
2845	: -mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
2846	fqd(il,i)=fr(il,i) ! precip
2847
2848	fu(il,i)=0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i))
2849	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2850	fv(il,i)=0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i))
2851	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2852	else ! cvflag_grav
2853	ccc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
2854	fr(il,i)=sigd(il)*evap(il,i)
2855	: +0.1(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i))
2856	: -mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
2857	fqd(il,i)=fr(il,i) ! precip
2858
2859	fu(il,i)=0.1(mp(il,i+1)(up(il,i+1)-u(il,i))
2860	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2861	fv(il,i)=0.1(mp(il,i+1)(vp(il,i+1)-v(il,i))
2862	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2863	endif ! cvflag_grav
2864
2865
2866	if (cvflag_grav) then
2867	fr(il,i)=fr(il,i)+0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i))
2868	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2869	fu(il,i)=fu(il,i)+0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i))
2870	: -ad(il)*(u(il,i)-u(il,i-1)))
2871	fv(il,i)=fv(il,i)+0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i))
2872	: -ad(il)*(v(il,i)-v(il,i-1)))
2873	else ! cvflag_grav
2874	fr(il,i)=fr(il,i)+0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il,i))
2875	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2876	fu(il,i)=fu(il,i)+0.1dpinv(amp1(il)*(u(il,i+1)-u(il,i))
2877	: -ad(il)*(u(il,i)-u(il,i-1)))
2878	fv(il,i)=fv(il,i)+0.1dpinv(amp1(il)*(v(il,i+1)-v(il,i))
2879	: -ad(il)*(v(il,i)-v(il,i-1)))
2880	endif ! cvflag_grav
2881
2882	endif ! i
2883	1350 continue
2884
2885	!AC! do k=1,ntra
2886	!AC! do il=1,ncum
2887	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
2888	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2889	!AC! cpinv=1.0/cpn(il,i)
2890	!AC! if (cvflag_grav) then
2891	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
2892	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2893	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2894	!AC! else
2895	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2896	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2897	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2898	!AC! endif
2899	!AC! endif
2900	!AC! enddo
2901	!AC! enddo
2902
2903	do 480 k=1,i-1
2904	c
2905	do il = 1,ncum
2906	awat(il)=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
2907	awat(il)=max(awat(il),0.0)
2908	enddo
2909	c
2910	IF (iflag_mix .ne. 0) then
2911	do il=1,ncum
2912	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2913	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2914	cpinv=1.0/cpn(il,i)
2915	if (cvflag_grav) then
2916	ft(il,i)=ft(il,i)
2917	: +0.01gravdpinvment(il,k,i)(hent(il,k,i)-h(il,i)
2918	: +t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-Qent(il,k,i)))*cpinv
2919
2920	c
2921	c
2922	else
2923	ft(il,i)=ft(il,i)
2924	: +0.1dpinvment(il,k,i)*(hent(il,k,i)-h(il,i)
2925	: +t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-Qent(il,k,i)))*cpinv
2926	endif !cvflag_grav
2927	endif ! i
2928	enddo
2929	ENDIF
2930	c
2931	do 1370 il=1,ncum
2932	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2933	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2934	cpinv=1.0/cpn(il,i)
2935	if (cvflag_grav) then
2936	fr(il,i)=fr(il,i)
2937	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-awat(il)-rr(il,i))
2938	fu(il,i)=fu(il,i)
2939	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2940	fv(il,i)=fv(il,i)
2941	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2942	else ! cvflag_grav
2943	fr(il,i)=fr(il,i)
2944	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-awat(il)-rr(il,i))
2945	fu(il,i)=fu(il,i)
2946	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2947	fv(il,i)=fv(il,i)
2948	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2949	endif ! cvflag_grav
2950
2951	c (saturated updrafts resulting from mixing) ! cld
2952	qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat(il)) ! cld
2953	nqcond(il,i)=nqcond(il,i)+1. ! cld
2954	endif ! i
2955	1370 continue
2956	480 continue
2957
2958	!AC! do j=1,ntra
2959	!AC! do k=1,i-1
2960	!AC! do il=1,ncum
2961	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
2962	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2963	!AC! cpinv=1.0/cpn(il,i)
2964	!AC! if (cvflag_grav) then
2965	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2966	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
2967	!AC! else
2968	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2969	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
2970	!AC! endif
2971	!AC! endif
2972	!AC! enddo
2973	!AC! enddo
2974	!AC! enddo
2975
2976	do 490 k=i,nl+1
2977	c
2978	IF (iflag_mix .ne. 0) then
2979	do il=1,ncum
2980	if (i.le.inb(il) .and. k.le.inb(il)
2981	$ .and. iflag(il) .le. 1) then
2982	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2983	cpinv=1.0/cpn(il,i)
2984	if (cvflag_grav) then
2985	ft(il,i)=ft(il,i)
2986	: +0.01gravdpinvment(il,k,i)(hent(il,k,i)-h(il,i)
2987	: +t(il,i)(cpv-cpd)(rr(il,i)-Qent(il,k,i)))*cpinv
2988	c
2989	c
2990	else
2991	ft(il,i)=ft(il,i)
2992	: +0.1dpinvment(il,k,i)*(hent(il,k,i)-h(il,i)
2993	: +t(il,i)(cpv-cpd)(rr(il,i)-Qent(il,k,i)))*cpinv
2994	endif !cvflag_grav
2995	endif ! i
2996	enddo
2997	ENDIF
2998	c
2999	do 1380 il=1,ncum
3000	if (i.le.inb(il) .and. k.le.inb(il)
3001	$ .and. iflag(il) .le. 1) then
3002	dpinv=1.0/(ph(il,i)-ph(il,i+1))
3003	cpinv=1.0/cpn(il,i)
3004
3005	if (cvflag_grav) then
3006	fr(il,i)=fr(il,i)
3007	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-rr(il,i))
3008	fu(il,i)=fu(il,i)
3009	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
3010	fv(il,i)=fv(il,i)
3011	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
3012	else ! cvflag_grav
3013	fr(il,i)=fr(il,i)
3014	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-rr(il,i))
3015	fu(il,i)=fu(il,i)
3016	: +0.1dpinvment(il,k,i)*(uent(il,k,i)-u(il,i))
3017	fv(il,i)=fv(il,i)
3018	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
3019	endif ! cvflag_grav
3020	endif ! i and k
3021	1380 continue
3022	490 continue
3023
3024	!AC! do j=1,ntra
3025	!AC! do k=i,nl+1
3026	!AC! do il=1,ncum
3027	!AC! if (i.le.inb(il) .and. k.le.inb(il)
3028	!AC! $ .and. iflag(il) .le. 1) then
3029	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3030	!AC! cpinv=1.0/cpn(il,i)
3031	!AC! if (cvflag_grav) then
3032	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3033	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3034	!AC! else
3035	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3036	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3037	!AC! endif
3038	!AC! endif ! i and k
3039	!AC! enddo
3040	!AC! enddo
3041	!AC! enddo
3042
3043	c sb: interface with the cloud parameterization: ! cld
3044
3045	do k=i+1,nl
3046	do il=1,ncum
3047	if (k.le.inb(il) .and. i.le.inb(il)
3048	$ .and. iflag(il) .le. 1) then ! cld
3049	C (saturated downdrafts resulting from mixing) ! cld
3050	qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
3051	nqcond(il,i)=nqcond(il,i)+1. ! cld
3052	endif ! cld
3053	enddo ! cld
3054	enddo ! cld
3055
3056	C (particular case: no detraining level is found) ! cld
3057	do il=1,ncum ! cld
3058	if (i.le.inb(il) .and. nent(il,i).eq.0
3059	$ .and. iflag(il) .le. 1) then ! cld
3060	qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
3061	nqcond(il,i)=nqcond(il,i)+1. ! cld
3062	endif ! cld
3063	enddo ! cld
3064
3065	do il=1,ncum ! cld
3066	if (i.le.inb(il) .and. nqcond(il,i).ne.0
3067	$ .and. iflag(il) .le. 1) then ! cld
3068	qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
3069	endif ! cld
3070	enddo
3071
3072	!AC! do j=1,ntra
3073	!AC! do il=1,ncum
3074	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3075	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3076	!AC! cpinv=1.0/cpn(il,i)
3077	!AC!
3078	!AC! if (cvflag_grav) then
3079	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
3080	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3081	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3082	!AC! else
3083	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
3084	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3085	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3086	!AC! endif
3087	!AC! endif ! i
3088	!AC! enddo
3089	!AC! enddo
3090
3091
3092	500 continue
3093
3094
3095	c * move the detrainment at level inb down to level inb-1 *
3096	c * in such a way as to preserve the vertically *
3097	c * integrated enthalpy and water tendencies *
3098	c
3099	c Correction bug le 18-03-09
3100	do 503 il=1,ncum
3101	IF (iflag(il) .le. 1) THEN
3102	if (cvflag_grav) then
3103	ax=0.01gravment(il,inb(il),inb(il))*(hp(il,inb(il))
3104	: -h(il,inb(il))+t(il,inb(il))*(cpv-cpd)
3105	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
3106	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
3107	ft(il,inb(il))=ft(il,inb(il))-ax
3108	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
3109	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
3110	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
3111
3112	bx=0.01gravment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il))
3113	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
3114	fr(il,inb(il))=fr(il,inb(il))-bx
3115	fr(il,inb(il)-1)=fr(il,inb(il)-1)
3116	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
3117	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
3118
3119	cx=0.01gravment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il))
3120	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
3121	fu(il,inb(il))=fu(il,inb(il))-cx
3122	fu(il,inb(il)-1)=fu(il,inb(il)-1)
3123	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
3124	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
3125
3126	dx=0.01gravment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il))
3127	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
3128	fv(il,inb(il))=fv(il,inb(il))-dx
3129	fv(il,inb(il)-1)=fv(il,inb(il)-1)
3130	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
3131	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
3132	else
3133	ax=0.1ment(il,inb(il),inb(il))(hp(il,inb(il))
3134	: -h(il,inb(il))+t(il,inb(il))*(cpv-cpd)
3135	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
3136	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
3137	ft(il,inb(il))=ft(il,inb(il))-ax
3138	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
3139	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
3140	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
3141
3142	bx=0.1ment(il,inb(il),inb(il))(qent(il,inb(il),inb(il))
3143	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
3144	fr(il,inb(il))=fr(il,inb(il))-bx
3145	fr(il,inb(il)-1)=fr(il,inb(il)-1)
3146	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
3147	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
3148
3149	cx=0.1ment(il,inb(il),inb(il))(uent(il,inb(il),inb(il))
3150	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
3151	fu(il,inb(il))=fu(il,inb(il))-cx
3152	fu(il,inb(il)-1)=fu(il,inb(il)-1)
3153	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
3154	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
3155
3156	dx=0.1ment(il,inb(il),inb(il))(vent(il,inb(il),inb(il))
3157	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
3158	fv(il,inb(il))=fv(il,inb(il))-dx
3159	fv(il,inb(il)-1)=fv(il,inb(il)-1)
3160	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
3161	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
3162	endif
3163	ENDIF !iflag
3164	503 continue
3165
3166	!AC! do j=1,ntra
3167	!AC! do il=1,ncum
3168	!AC! IF (iflag(il) .le. 1) THEN
3169	!AC! IF (cvflag_grav) then
3170	!AC! ex=0.01gravment(il,inb(il),inb(il))
3171	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3172	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3173	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3174	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3175	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3176	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3177	!AC! else
3178	!AC! ex=0.1*ment(il,inb(il),inb(il))
3179	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3180	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3181	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3182	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3183	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3184	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3185	!AC! ENDIF !cvflag grav
3186	!AC! ENDIF !iflag
3187	!AC! enddo
3188	!AC! enddo
3189
3190	c
3191	c * homogenize tendencies below cloud base *
3192	c
3193	c
3194	do il=1,ncum
3195	asum(il)=0.0
3196	bsum(il)=0.0
3197	csum(il)=0.0
3198	dsum(il)=0.0
3199	esum(il)=0.0
3200	fsum(il)=0.0
3201	gsum(il)=0.0
3202	hsum(il)=0.0
3203	enddo
3204	c
3205	c do i=1,nl
3206	c do il=1,ncum
3207	c th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
3208	c enddo
3209	c enddo
3210	c
3211	do i=1,nl
3212	do il=1,ncum
3213	if (i.le.(icb(il)-1) .and. iflag(il) .le. 1) then
3214	cjyg Saturated part : use T profile
3215	asum(il)=asum(il)+(ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
3216	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))
3217	: (lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))
3218	: *(ph(il,i)-ph(il,i+1))
3219	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
3220	: *(ph(il,i)-ph(il,i+1))
3221	dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
3222	cjyg Unsaturated part : use T_wake profile
3223	esum(il)=esum(il)+ftd(il,i)*(ph(il,i)-ph(il,i+1))
3224	fsum(il)=fsum(il)+fqd(il,i)
3225	: (lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))
3226	: *(ph(il,i)-ph(il,i+1))
3227	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))
3228	: *(ph(il,i)-ph(il,i+1))
3229	hsum(il)=hsum(il)+t_wake(il,i)
3230	; *(ph(il,i)-ph(il,i+1))/th_wake(il,i)
3231	endif
3232	enddo
3233	enddo
3234
3235	c!!! do 700 i=1,icb(il)-1
3236	do i=1,nl
3237	do il=1,ncum
3238	if (i.le.(icb(il)-1) .and. iflag(il) .le. 1) then
3239	ftd(il,i)=esum(il)t_wake(il,i)/(th_wake(il,i)hsum(il))
3240	fqd(il,i)=fsum(il)/gsum(il)
3241	ft(il,i)=ftd(il,i)+asum(il)t(il,i)/(th(il,i)dsum(il))
3242	fr(il,i)=fqd(il,i)+bsum(il)/csum(il)
3243	endif
3244	enddo
3245	enddo
3246
3247	c
3248	c *** Check that moisture stays positive. If not, scale tendencies
3249	c in order to ensure moisture positivity
3250	DO il = 1,ncum
3251	alpha_qpos(il)=1.
3252	IF (iflag(il) .le. 1) THEN
3253	if (fr(il,1) .le. 0.) then
3254	alpha_qpos(il) = max(alpha_qpos(il) ,
3255	: (-delt*fr(il,1))/
3256	: (s_wake(il)rr_wake(il,1)+(1.-s_wake(il))rr(il,1)))
3257	end if
3258	ENDIF
3259	ENDDO
3260	DO i = 2,nl
3261	DO il = 1,ncum
3262	IF (iflag(il) .le. 1) THEN
3263	IF (fr(il,i) .le. 0.) THEN
3264	alpha_qpos1(il)=max(1. , (-delt*fr(il,i))/
3265	: (s_wake(il)rr_wake(il,i)+(1.-s_wake(il))rr(il,i)))
3266	IF (alpha_qpos1(il) .ge. alpha_qpos(il))
3267	: alpha_qpos(il)=alpha_qpos1(il)
3268	ENDIF
3269	ENDIF
3270	ENDDO
3271	ENDDO
3272	DO il = 1,ncum
3273	IF (iflag(il) .le. 1 .and. alpha_qpos(il) .gt. 1.001) THEN
3274	alpha_qpos(il) = alpha_qpos(il)*1.1
3275	ENDIF
3276	ENDDO
3277	DO il = 1,ncum
3278	IF (iflag(il) .le. 1) THEN
3279	sigd(il) = sigd(il)/alpha_qpos(il)
3280	precip(il) = precip(il)/alpha_qpos(il)
3281	ENDIF
3282	ENDDO
3283	DO i = 1,nl
3284	DO il = 1,ncum
3285	IF (iflag(il) .le. 1) THEN
3286	fr(il,i) = fr(il,i)/alpha_qpos(il)
3287	ft(il,i) = ft(il,i)/alpha_qpos(il)
3288	fqd(il,i) = fqd(il,i)/alpha_qpos(il)
3289	ftd(il,i) = ftd(il,i)/alpha_qpos(il)
3290	fu(il,i) = fu(il,i)/alpha_qpos(il)
3291	fv(il,i) = fv(il,i)/alpha_qpos(il)
3292	m(il,i) = m(il,i)/alpha_qpos(il)
3293	mp(il,i) = mp(il,i)/alpha_qpos(il)
3294	Vprecip(il,i) = Vprecip(il,i)/alpha_qpos(il)
3295	ENDIF
3296	ENDDO
3297	ENDDO
3298	DO i = 1,nl
3299	DO j = 1,nl
3300	DO il = 1,ncum
3301	IF (iflag(il) .le. 1) THEN
3302	ment(il,i,j) = ment(il,i,j)/alpha_qpos(il)
3303	ENDIF
3304	ENDDO
3305	ENDDO
3306	ENDDO
3307
3308	!AC! DO j = 1,ntra
3309	!AC! DO i = 1,nl
3310	!AC! DO il = 1,ncum
3311	!AC! IF (iflag(il) .le. 1) THEN
3312	!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
3313	!AC! ENDIF
3314	!AC! ENDDO
3315	!AC! ENDDO
3316	!AC! ENDDO
3317
3318	c
3319	c * reset counter and return *
3320	c
3321	do il=1,ncum
3322	sig(il,nd)=2.0
3323	enddo
3324
3325
3326	do i=1,nd
3327	do il=1,ncum
3328	upwd(il,i)=0.0
3329	dnwd(il,i)=0.0
3330	enddo
3331	enddo
3332
3333	do i=1,nl
3334	do il=1,ncum
3335	dnwd0(il,i)=-mp(il,i)
3336	enddo
3337	enddo
3338	do i=nl+1,nd
3339	do il=1,ncum
3340	dnwd0(il,i)=0.
3341	enddo
3342	enddo
3343
3344
3345	do i=1,nl
3346	do il=1,ncum
3347	if (i.ge.icb(il) .and. i.le.inb(il)) then
3348	upwd(il,i)=0.0
3349	dnwd(il,i)=0.0
3350	endif
3351	enddo
3352	enddo
3353
3354	do i=1,nl
3355	do k=1,nl
3356	do il=1,ncum
3357	up1(il,k,i)=0.0
3358	dn1(il,k,i)=0.0
3359	enddo
3360	enddo
3361	enddo
3362
3363	do i=1,nl
3364	do k=i,nl
3365	do n=1,i-1
3366	do il=1,ncum
3367	if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
3368	up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
3369	dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
3370	endif
3371	enddo
3372	enddo
3373	enddo
3374	enddo
3375
3376	do i=1,nl
3377	do k=1,nl
3378	do il=1,ncum
3379	if(i.ge.icb(il)) then
3380	if(k.ge.i.and. k.le.(inb(il))) then
3381	upwd(il,i)=upwd(il,i)+m(il,k)
3382	endif
3383	else
3384	if(k.lt.i) then
3385	upwd(il,i)=upwd(il,i)+cbmf(il)*wghti(il,k)
3386	endif
3387	endif
3388	cc print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
3389	end do
3390	end do
3391	end do
3392
3393	do i=2,nl
3394	do k=i,nl
3395	do il=1,ncum
3396	ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
3397	if (i.le.inb(il).and.k.le.inb(il)) then
3398	upwd(il,i)=upwd(il,i)+up1(il,k,i)
3399	dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
3400	endif
3401	cc print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
3402	enddo
3403	enddo
3404	enddo
3405
3406
3407	c!!! DO il=1,ncum
3408	c!!! do i=icb(il),inb(il)
3409	c!!!
3410	c!!! upwd(il,i)=0.0
3411	c!!! dnwd(il,i)=0.0
3412	c!!! do k=i,inb(il)
3413	c!!! up1=0.0
3414	c!!! dn1=0.0
3415	c!!! do n=1,i-1
3416	c!!! up1=up1+ment(il,n,k)
3417	c!!! dn1=dn1-ment(il,k,n)
3418	c!!! enddo
3419	c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
3420	c!!! dnwd(il,i)=dnwd(il,i)+dn1
3421	c!!! enddo
3422	c!!! enddo
3423	c!!!
3424	c!!! ENDDO
3425
3426	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3427	c determination de la variation de flux ascendant entre
3428	c deux niveau non dilue mip
3429	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3430
3431	do i=1,nl
3432	do il=1,ncum
3433	mip(il,i)=m(il,i)
3434	enddo
3435	enddo
3436
3437	do i=nl+1,nd
3438	do il=1,ncum
3439	mip(il,i)=0.
3440	enddo
3441	enddo
3442
3443	do i=1,nd
3444	do il=1,ncum
3445	ma(il,i)=0
3446	enddo
3447	enddo
3448
3449	do i=1,nl
3450	do j=i,nl
3451	do il=1,ncum
3452	ma(il,i)=ma(il,i)+m(il,j)
3453	enddo
3454	enddo
3455	enddo
3456
3457	do i=nl+1,nd
3458	do il=1,ncum
3459	ma(il,i)=0.
3460	enddo
3461	enddo
3462
3463	do i=1,nl
3464	do il=1,ncum
3465	if (i.le.(icb(il)-1)) then
3466	ma(il,i)=0
3467	endif
3468	enddo
3469	enddo
3470
3471	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3472	c icb represente de niveau ou se trouve la
3473	c base du nuage , et inb le top du nuage
3474	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3475
3476	do i=1,nd
3477	do il=1,ncum
3478	mke(il,i)=upwd(il,i)+dnwd(il,i)
3479	enddo
3480	enddo
3481
3482	do i=1,nd
3483	DO 999 il=1,ncum
3484	rdcp=(rrd(1.-rr(il,i))-rr(il,i)rrv)
3485	: /(cpd(1.-rr(il,i))+rr(il,i)cpv)
3486	tls(il,i)=t(il,i)(1000.0/p(il,i))*rdcp
3487	tps(il,i)=tp(il,i)
3488	999 CONTINUE
3489	enddo
3490
3491	c
3492	c * diagnose the in-cloud mixing ratio * ! cld
3493	c * of condensed water * ! cld
3494	c ! cld
3495
3496	do i=1,nd ! cld
3497	do il=1,ncum ! cld
3498	mac(il,i)=0.0 ! cld
3499	wa(il,i)=0.0 ! cld
3500	siga(il,i)=0.0 ! cld
3501	sax(il,i)=0.0 ! cld
3502	enddo ! cld
3503	enddo ! cld
3504
3505	do i=minorig, nl ! cld
3506	do k=i+1,nl+1 ! cld
3507	do il=1,ncum ! cld
3508	if (i.le.inb(il) .and. k.le.(inb(il)+1)
3509	$ .and. iflag(il) .le. 1) then ! cld
3510	mac(il,i)=mac(il,i)+m(il,k) ! cld
3511	endif ! cld
3512	enddo ! cld
3513	enddo ! cld
3514	enddo ! cld
3515
3516	do i=1,nl ! cld
3517	do j=1,i ! cld
3518	do il=1,ncum ! cld
3519	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
3520	: .and. j.ge.icb(il)
3521	: .and. iflag(il) .le. 1 ) then ! cld
3522	sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
3523	: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
3524	endif ! cld
3525	enddo ! cld
3526	enddo ! cld
3527	enddo ! cld
3528
3529	do i=1,nl ! cld
3530	do il=1,ncum ! cld
3531	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
3532	: .and. sax(il,i).gt.0.0
3533	: .and. iflag(il) .le. 1 ) then ! cld
3534	wa(il,i)=sqrt(2.*sax(il,i)) ! cld
3535	endif ! cld
3536	enddo ! cld
3537	enddo ! cld
3538
3539	do i=1,nl ! cld
3540	do il=1,ncum ! cld
3541	if (wa(il,i).gt.0.0 .and. iflag(il) .le. 1) ! cld
3542	: siga(il,i)=mac(il,i)/wa(il,i) ! cld
3543	: rrdtvp(il,i)/p(il,i)/100./delta ! cld
3544	siga(il,i) = min(siga(il,i),1.0) ! cld
3545	cIM cf. FH
3546	if (iflag_clw.eq.0) then
3547	qcondc(il,i)=siga(il,i)clw(il,i)(1.-ep(il,i)) ! cld
3548	: + (1.-siga(il,i))*qcond(il,i) ! cld
3549	else if (iflag_clw.eq.1) then
3550	qcondc(il,i)=qcond(il,i) ! cld
3551	endif
3552
3553	enddo ! cld
3554	enddo
3555	c print*,'cv3_yield fin'
3556	! cld
3557	return
3558	end
3559
3560	!AC! et !RomP >>>
3561	SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
3562	& ment,sigij,da,phi,phi2,d1a,dam,
3563	& ep,Vprecip,elij,clw,epmlmMm,eplaMm,
3564	& icb,inb)
3565	implicit none
3566
3567	#include "cv3param.h"
3568
3569	c inputs:
3570	integer ncum, nd, na, nloc,len
3571	real ment(nloc,na,na),sigij(nloc,na,na)
3572	real clw(nloc,nd),elij(nloc,na,na)
3573	real ep(nloc,na)
3574	integer icb(nloc),inb(nloc)
3575	real VPrecip(nloc,nd+1)
3576	c ouputs:
3577	real da(nloc,na),phi(nloc,na,na)
3578	real phi2(nloc,na,na)
3579	real d1a(nloc,na),dam(nloc,na)
3580	real epmlmMm(nloc,na,na),eplaMm(nloc,na)
3581	! variables pour tracer dans precip de l'AA et des mel
3582	c local variables:
3583	integer i,j,k
3584	real epm(nloc,na,na)
3585	c
3586	! variables d'Emanuel : du second indice au troisieme
3587	! ---> tab(i,k,j) -> de l origine k a l arrivee j
3588	! ment, sigij, elij
3589	! variables personnelles : du troisieme au second indice
3590	! ---> tab(i,j,k) -> de k a j
3591	! phi, phi2
3592	!
3593	! initialisations
3594	c
3595	da(:,:)=0.
3596	d1a(:,:)=0.
3597	dam(:,:)=0.
3598	epm(:,:,:)=0.
3599	eplaMm(:,:)=0.
3600	epmlmMm(:,:,:)=0.
3601	phi(:,:,:)=0.
3602	phi2(:,:,:)=0.
3603	c
3604	! fraction deau condensee dans les melanges convertie en precip : epm
3605	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3606	do j=1,na
3607	do k=1,na
3608	do i=1,ncum
3609	if(k.ge.icb(i).and.k.le.inb(i).and.
3610	!!jyg & j.ge.k.and.j.le.inb(i)) then
3611	!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3612	& j.gt.k.and.j.le.inb(i)) then
3613	epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/
3614	& max(elij(i,k,j),1.e-16)
3615	!!
3616	epm(i,j,k)=max(epm(i,j,k),0.0)
3617	endif
3618	end do
3619	end do
3620	end do
3621
3622	!
3623	do j=1,na
3624	do k=1,na
3625	do i=1,ncum
3626	if(k.ge.icb(i).and.k.le.inb(i)) then
3627	eplaMm(i,j)=eplaMm(i,j) + ep(i,j)*clw(i,j)
3628	& ment(i,j,k)(1.-sigij(i,j,k))
3629	endif
3630	end do
3631	end do
3632	end do
3633	!
3634	do j=1,na
3635	do k=1,j-1
3636	do i=1,ncum
3637	if(k.ge.icb(i).and.k.le.inb(i).and.
3638	& j.le.inb(i)) then
3639	epmlmMm(i,j,k)=epm(i,j,k)elij(i,k,j)ment(i,k,j)
3640	endif
3641	end do
3642	end do
3643	end do
3644
3645	! matrices pour calculer la tendance des concentrations dans cvltr.F90
3646	do j=1,na
3647	do k=1,na
3648	do i=1,ncum
3649	da(i,j)=da(i,j)+(1.-sigij(i,k,j))*ment(i,k,j)
3650	phi(i,j,k)=sigij(i,k,j)*ment(i,k,j)
3651	d1a(i,j)=d1a(i,j)+ment(i,k,j)*ep(i,k)
3652	& *(1.-sigij(i,k,j))
3653	if(k.le.j) then
3654	dam(i,j)=dam(i,j)+ment(i,k,j)
3655	& epm(i,k,j)(1.-ep(i,k))*(1.-sigij(i,k,j))
3656
3657	phi2(i,j,k)=phi(i,j,k)*epm(i,j,k)
3658	endif
3659	end do
3660	end do
3661	end do
3662
3663	return
3664	end
3665	!AC! et !RomP <<<
3666
3667	SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
3668	: ,iflag
3669	: ,precip,sig,w0
3670	: ,ft,fq,fu,fv,ftra
3671	: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
3672	: ,iflag1
3673	: ,precip1,sig1,w01
3674	: ,ft1,fq1,fu1,fv1,ftra1
3675	: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
3676	: )
3677	implicit none
3678
3679	#include "cv3param.h"
3680
3681	c inputs:
3682	integer len, ncum, nd, ntra, nloc
3683	integer idcum(nloc)
3684	integer iflag(nloc)
3685	real precip(nloc)
3686	real sig(nloc,nd), w0(nloc,nd)
3687	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
3688	real ftra(nloc,nd,ntra)
3689	real Ma(nloc,nd)
3690	real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
3691	real qcondc(nloc,nd)
3692	real wd(nloc),cape(nloc)
3693
3694	c outputs:
3695	integer iflag1(len)
3696	real precip1(len)
3697	real sig1(len,nd), w01(len,nd)
3698	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
3699	real ftra1(len,nd,ntra)
3700	real Ma1(len,nd)
3701	real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
3702	real qcondc1(nloc,nd)
3703	real wd1(nloc),cape1(nloc)
3704
3705	c local variables:
3706	integer i,k,j
3707
3708	do 2000 i=1,ncum
3709	precip1(idcum(i))=precip(i)
3710	iflag1(idcum(i))=iflag(i)
3711	wd1(idcum(i))=wd(i)
3712	cape1(idcum(i))=cape(i)
3713	2000 continue
3714
3715	do 2020 k=1,nl
3716	do 2010 i=1,ncum
3717	sig1(idcum(i),k)=sig(i,k)
3718	w01(idcum(i),k)=w0(i,k)
3719	ft1(idcum(i),k)=ft(i,k)
3720	fq1(idcum(i),k)=fq(i,k)
3721	fu1(idcum(i),k)=fu(i,k)
3722	fv1(idcum(i),k)=fv(i,k)
3723	Ma1(idcum(i),k)=Ma(i,k)
3724	upwd1(idcum(i),k)=upwd(i,k)
3725	dnwd1(idcum(i),k)=dnwd(i,k)
3726	dnwd01(idcum(i),k)=dnwd0(i,k)
3727	qcondc1(idcum(i),k)=qcondc(i,k)
3728	2010 continue
3729	2020 continue
3730
3731	do 2200 i=1,ncum
3732	sig1(idcum(i),nd)=sig(i,nd)
3733	2200 continue
3734
3735
3736	!AC! do 2100 j=1,ntra
3737	!AC!c oct3 do 2110 k=1,nl
3738	!AC! do 2110 k=1,nd ! oct3
3739	!AC! do 2120 i=1,ncum
3740	!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
3741	!AC! 2120 continue
3742	!AC! 2110 continue
3743	!AC! 2100 continue
3744	return
3745	end
3746

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

Download in other formats:

Original Format