Context Navigation

cv3_routines.F @ 1501

Last change on this file since 1501 was 1501, checked in by idelkadi, 13 years ago
Rajout de la nouvelle formulation de l'efficacite des precipitations fonction de la temperature et contenu en eau liquide
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 109.4 KB

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

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

Download in other formats:

Original Format