Context Navigation

cv3_routines.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 106.7 KB

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

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

Download in other formats:

Original Format