Context Navigation

cv3_routines.F @ 1494

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

Optimization of routines from the new physics

Optimisation de routines de la nouvelle physique

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

File size: 108.7 KB

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

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

Download in other formats:

Original Format