Context Navigation

cv3_routines.F @ 1183

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

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

Download in other formats:

Original Format