Context Navigation

cv3_routines.F @ 5416

Last change on this file since 5416 was 1467, checked in by idelkadi, 14 years ago
Input file conv_param.data added to change convection parameters without changing them inside the code
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 107.1 KB

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

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

Download in other formats:

Original Format