Context Navigation

cv3_routines.F @ 1483

Last change on this file since 1483 was 1334, checked in by idelkadi, 14 years ago
Rajout de champs de sorties Correction dans la partie convection (nouvelle physique)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 105.5 KB

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

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

Download in other formats:

Original Format