Context Navigation

cv3_routines.F @ 1850

Last change on this file since 1850 was 1849, checked in by Ehouarn Millour, 11 years ago
Including the thermodynamics of ice in the convection scheme (iactive only if iflag_ice_thermo==1). CR+JYG
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 126.3 KB

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

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

Download in other formats:

Original Format