Context Navigation

cv3_routines.F @ 1479

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

Un peu de ménage sur les prints

A little cleanup on the prints

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

File size: 107.2 KB

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

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

Download in other formats:

Original Format