Context Navigation

cv3_routines.F @ 996

Last change on this file since 996 was 993, checked in by Laurent Fairhead, 16 years ago
Oubli de declaration de variables sur la derniere correction LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 101.3 KB

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

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

Download in other formats:

Original Format