Context Navigation

cv3_routines.F @ 879

Last change on this file since 879 was 879, checked in by Laurent Fairhead, 16 years ago
Suite de la bascule vers une physique avec thermiques, nouvelle convection, poche froide ... LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 100.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,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
1528	c local variables:
1529	integer i, j, k, il, im, jm
1530	integer num1, num2
1531	integer nent(nloc,na)
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	b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2078	c6 = water(il,i+1) + wdtrain(il)*bfac
2079	revap=0.5(-b6+sqrt(b6b6+4.*c6))
2080	evap(il,i)=sigtafacrevap
2081	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	if (1.eq.0) then
2088	b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2089	c6=water(il,i+1)+bfac*wdtrain(il)
2090	: -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2091	if(c6.gt.0.0)then
2092	revap=0.5(-b6+sqrt(b6b6+4.*c6))
2093	evap(il,i)=sigtafacrevap
2094	water(il,i)=revap*revap
2095	else
2096	evap(il,i)=-evap(il,i+1)
2097	: +(wdtrain(il)+sigd(il)wt(il,i)water(il,i+1))
2098	: /(sigd(il)(ph(il,i)-ph(il,i+1))50.)
2099	end if
2100	endif
2101	ccc
2102	c * calculate precipitating downdraft mass flux under *
2103	c * hydrostatic approximation *
2104	c
2105	if (i.ne.1) then
2106
2107	tevap=amax1(0.0,evap(il,i))
2108	delth=amax1(0.001,(th(il,i)-th(il,i-1)))
2109	if (cvflag_grav) then
2110	mp(il,i)=100.ginvlvcp(il,i)sigd(il)tevap
2111	: *(p(il,i-1)-p(il,i))/delth
2112	else
2113	mp(il,i)=10.lvcp(il,i)sigd(il)*tevap
2114	: *(p(il,i-1)-p(il,i))/delth
2115	endif
2116	c
2117	c * if hydrostatic assumption fails, *
2118	c * solve cubic difference equation for downdraft theta *
2119	c * and mass flux from two simultaneous differential eqns *
2120	c
2121	amfac=sigd(il)sigd(il)70.0ph(il,i)(p(il,i-1)-p(il,i))
2122	: (th(il,i)-th(il,i-1))/(tv(il,i)th(il,i))
2123	amp2=abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2124	if(amp2.gt.(0.1*amfac))then
2125	xf=100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
2126	tf=b(il,i)-5.0(th(il,i)-th(il,i-1))t(il,i)
2127	: /(lvcp(il,i)sigd(il)th(il,i))
2128	af=xftf+mp(il,i+1)mp(il,i+1)*tinv
2129	bf=2.(tinvmp(il,i+1))*3+tinvmp(il,i+1)xftf
2130	: +50.(p(il,i-1)-p(il,i))xf*tevap
2131	fac2=1.0
2132	if(bf.lt.0.0)fac2=-1.0
2133	bf=abs(bf)
2134	ur=0.25bfbf-afafaftinvtinv*tinv
2135	if(ur.ge.0.0)then
2136	sru=sqrt(ur)
2137	fac=1.0
2138	if((0.5*bf-sru).lt.0.0)fac=-1.0
2139	mp(il,i)=mp(il,i+1)tinv+(0.5bf+sru)**tinv
2140	: +fac(abs(0.5bf-sru))**tinv
2141	else
2142	d=atan(2.*sqrt(-ur)/(bf+1.0e-28))
2143	if(fac2.lt.0.0)d=3.14159-d
2144	mp(il,i)=mp(il,i+1)tinv+2.sqrt(aftinv)cos(d*tinv)
2145	endif
2146	mp(il,i)=amax1(0.0,mp(il,i))
2147
2148	if (cvflag_grav) then
2149	Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2150	C il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
2151	C Et il faut bien revoir les facteurs 100.
2152	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2153	2 /(mp(il,i)+sigd(il)*0.1)
2154	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)
2155	: sigd(il)th(il,i))
2156	else
2157	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2158	2 /(mp(il,i)+sigd(il)*0.1)
2159	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)
2160	: sigd(il)th(il,i))
2161	endif
2162	b(il,i-1)=amax1(b(il,i-1),0.0)
2163	endif
2164	c
2165	c * limit magnitude of mp(i) to meet cfl condition *
2166	c
2167	ampmax=2.0(ph(il,i)-ph(il,i+1))delti
2168	amp2=2.0(ph(il,i-1)-ph(il,i))delti
2169	ampmax=amin1(ampmax,amp2)
2170	mp(il,i)=amin1(mp(il,i),ampmax)
2171	c
2172	c * force mp to decrease linearly to zero *
2173	c * between cloud base and the surface *
2174	c
2175	c
2176	cc if(p(il,i).gt.p(il,icb(il)))then
2177	cc mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2178	cc endif
2179	if(ph(il,i) .gt. 0.9*plcl(il)) then
2180	mp(il,i) = mp(il,i)*(ph(il,1)-ph(il,i))/
2181	$ (ph(il,1)-0.9*plcl(il))
2182	endif
2183
2184	360 continue
2185	endif ! i.eq.1
2186	c
2187	c * find mixing ratio of precipitating downdraft *
2188	c
2189
2190	if (i.ne.inb(il)) then
2191
2192	rp(il,i)=rr(il,i)
2193
2194	if(mp(il,i).gt.mp(il,i+1))then
2195
2196	if (cvflag_grav) then
2197	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2198	: +100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))
2199	: *(evap(il,i+1)+evap(il,i))
2200	else
2201	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2202	: +5.sigd(il)(ph(il,i)-ph(il,i+1))
2203	: *(evap(il,i+1)+evap(il,i))
2204	endif
2205	rp(il,i)=rp(il,i)/mp(il,i)
2206	up(il,i)=up(il,i+1)mp(il,i+1)+u(il,i)(mp(il,i)-mp(il,i+1))
2207	up(il,i)=up(il,i)/mp(il,i)
2208	vp(il,i)=vp(il,i+1)mp(il,i+1)+v(il,i)(mp(il,i)-mp(il,i+1))
2209	vp(il,i)=vp(il,i)/mp(il,i)
2210
2211	do j=1,ntra
2212	trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2213	: +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2214	trap(il,i,j)=trap(il,i,j)/mp(il,i)
2215	end do
2216
2217	else
2218
2219	if(mp(il,i+1).gt.1.0e-16)then
2220	if (cvflag_grav) then
2221	rp(il,i)=rp(il,i+1)
2222	: +100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))
2223	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2224	else
2225	rp(il,i)=rp(il,i+1)
2226	: +5.sigd(il)(ph(il,i)-ph(il,i+1))
2227	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2228	endif
2229	up(il,i)=up(il,i+1)
2230	vp(il,i)=vp(il,i+1)
2231
2232	do j=1,ntra
2233	trap(il,i,j)=trap(il,i+1,j)
2234	end do
2235
2236	endif
2237	endif
2238	rp(il,i)=amin1(rp(il,i),rs(il,i))
2239	rp(il,i)=amax1(rp(il,i),0.0)
2240
2241	endif
2242	endif
2243	999 continue
2244
2245	400 continue
2246
2247	return
2248	end
2249
2250	SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
2251	: ,icb,inb,delt
2252	: ,t,rr,t_wake,rr_wake,u,v,tra
2253	: ,gz,p,ph,h,hp,lv,cpn,th
2254	: ,ep,clw,m,tp,mp,rp,up,vp,trap
2255	: ,wt,water,evap,b,sigd
2256	: ,ment,qent,hent,iflag_mix,uent,vent
2257	: ,nent,elij,traent,sig
2258	: ,tv,tvp,wghti
2259	: ,iflag,precip,Vprecip,ft,fr,fu,fv,ftra
2260	: ,cbmf,upwd,dnwd,dnwd0,ma,mip
2261	: ,tls,tps,qcondc,wd
2262	: ,ftd,fqd)
2263
2264	implicit none
2265
2266	#include "cvthermo.h"
2267	#include "cv3param.h"
2268	#include "cvflag.h"
2269	#include "conema3.h"
2270
2271	c inputs:
2272	c print*,'cv3_yield apres include'
2273	integer iflag_mix
2274	integer ncum,nd,na,ntra,nloc
2275	integer icb(nloc), inb(nloc)
2276	real delt
2277	real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
2278	real t_wake(nloc,nd), rr_wake(nloc,nd)
2279	real tra(nloc,nd,ntra), sig(nloc,nd)
2280	real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
2281	real th(nloc,na), p(nloc,nd), tp(nloc,na)
2282	real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
2283	real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
2284	real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
2285	real water(nloc,na), evap(nloc,na), b(nloc,na), sigd(nloc)
2286	real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
2287	real hent(nloc,na,na)
2288	real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2289	real traent(nloc,na,na,ntra)
2290	real tv(nloc,nd), tvp(nloc,nd), wghti(nloc,nd)
2291	c print*,'cv3_yield declarations 1'
2292	c input/output:
2293	integer iflag(nloc)
2294
2295	c outputs:
2296	real precip(nloc)
2297	real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
2298	real ftd(nloc,nd), fqd(nloc,nd)
2299	real ftra(nloc,nd,ntra)
2300	real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
2301	real dnwd0(nloc,nd), mip(nloc,nd)
2302	real Vprecip(nloc,nd)
2303	real tls(nloc,nd), tps(nloc,nd)
2304	real qcondc(nloc,nd) ! cld
2305	real wd(nloc) ! gust
2306	real cbmf(nloc)
2307	c print*,'cv3_yield declarations 2'
2308	c local variables:
2309	integer i,k,il,n,j,num1
2310	real rat, delti
2311	real ax, bx, cx, dx, ex
2312	real cpinv, rdcp, dpinv
2313	real awat(nloc)
2314	real lvcp(nloc,na), mke(nloc,na)
2315	real am(nloc), work(nloc), ad(nloc), amp1(nloc)
2316	c!! real up1(nloc), dn1(nloc)
2317	real up1(nloc,nd,nd), dn1(nloc,nd,nd)
2318	real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
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	do j=2,nl
2469	IF (iflag_mix .gt. 0) then
2470	do il=1,ncum
2471	if (j.le.inb(il) .and. iflag(il) .le. 1) then
2472	if (cvflag_grav) then
2473	ft(il,1)=ft(il,1)
2474	: +0.01gravwork(il)ment(il,j,1)(hent(il,j,1)-h(il,1)
2475	: +t(il,1)(cpv-cpd)(rr(il,1)-Qent(il,j,1)))*cpinv
2476	else
2477	ft(il,1)=ft(il,1)
2478	: +0.1work(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	endif ! cvflag_grav
2481	endif ! j
2482	enddo
2483	ENDIF
2484	enddo
2485	! fin sature
2486
2487
2488	do il=1,ncum
2489	if (iflag(il) .le. 1) then
2490	if (cvflag_grav) then
2491	Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
2492	fr(il,1)=0.01gravmp(il,2)(rp(il,2)-rr_wake(il,1))work(il)
2493	: +sigd(il)*evap(il,1)
2494	ccc : +sigd(il)0.5(evap(il,1)+evap(il,2))
2495
2496	fqd(il,1)=fr(il,1) !precip
2497
2498	fr(il,1)=fr(il,1)+0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
2499
2500	fu(il,1)=fu(il,1)+0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1))
2501	: +am(il)*(u(il,2)-u(il,1)))
2502	fv(il,1)=fv(il,1)+0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1))
2503	: +am(il)*(v(il,2)-v(il,1)))
2504	else ! cvflag_grav
2505	fr(il,1)=0.1mp(il,2)(rp(il,2)-rr_wake(il,1))*work(il)
2506	: +sigd(il)*evap(il,1)
2507	ccc : +sigd(il)0.5(evap(il,1)+evap(il,2))
2508	fqd(il,1)=fr(il,1) !precip
2509	fr(il,1)=fr(il,1)+0.1am(il)(rr(il,2)-rr(il,1))*work(il)
2510	fu(il,1)=fu(il,1)+0.1work(il)(mp(il,2)*(up(il,2)-u(il,1))
2511	: +am(il)*(u(il,2)-u(il,1)))
2512	fv(il,1)=fv(il,1)+0.1work(il)(mp(il,2)*(vp(il,2)-v(il,1))
2513	: +am(il)*(v(il,2)-v(il,1)))
2514	endif ! cvflag_grav
2515	endif ! iflag
2516	enddo ! il
2517
2518
2519	do j=1,ntra
2520	do il=1,ncum
2521	if (iflag(il) .le. 1) then
2522	if (cvflag_grav) then
2523	ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
2524	: (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2525	: +am(il)*(tra(il,2,j)-tra(il,1,j)))
2526	else
2527	ftra(il,1,j)=ftra(il,1,j)+0.1*work(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	endif
2531	endif ! iflag
2532	enddo
2533	enddo
2534
2535	do j=2,nl
2536	do il=1,ncum
2537	if (j.le.inb(il) .and. iflag(il) .le. 1) then
2538	if (cvflag_grav) then
2539	fr(il,1)=fr(il,1)
2540	: +0.01gravwork(il)ment(il,j,1)(qent(il,j,1)-rr(il,1))
2541	fu(il,1)=fu(il,1)
2542	: +0.01gravwork(il)ment(il,j,1)(uent(il,j,1)-u(il,1))
2543	fv(il,1)=fv(il,1)
2544	: +0.01gravwork(il)ment(il,j,1)(vent(il,j,1)-v(il,1))
2545	else ! cvflag_grav
2546	fr(il,1)=fr(il,1)
2547	: +0.1work(il)ment(il,j,1)*(qent(il,j,1)-rr(il,1))
2548	fu(il,1)=fu(il,1)
2549	: +0.1work(il)ment(il,j,1)*(uent(il,j,1)-u(il,1))
2550	fv(il,1)=fv(il,1)
2551	: +0.1work(il)ment(il,j,1)*(vent(il,j,1)-v(il,1)) ! fin sature
2552	endif ! cvflag_grav
2553	endif ! j
2554	enddo
2555	enddo
2556
2557	do k=1,ntra
2558	do j=2,nl
2559	do il=1,ncum
2560	if (j.le.inb(il) .and. iflag(il) .le. 1) then
2561
2562	if (cvflag_grav) then
2563	ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
2564	: *(traent(il,j,1,k)-tra(il,1,k))
2565	else
2566	ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
2567	: *(traent(il,j,1,k)-tra(il,1,k))
2568	endif
2569
2570	endif
2571	enddo
2572	enddo
2573	enddo
2574	c print*,'cv3_yield apres ft'
2575	c
2576	c * calculate tendencies of potential temperature and mixing ratio *
2577	c * at levels above the lowest level *
2578	c
2579	c * first find the net saturated updraft and downdraft mass fluxes *
2580	c * through each level *
2581	c
2582
2583	do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
2584
2585	num1=0
2586	do il=1,ncum
2587	if(i.le.inb(il) .and. iflag(il) .le. 1)num1=num1+1
2588	enddo
2589	if(num1.le.0)go to 500
2590
2591	call zilch(amp1,ncum)
2592	call zilch(ad,ncum)
2593
2594	do 440 k=1,nl+1
2595	do 441 il=1,ncum
2596	if(i.ge.icb(il)) then
2597	if(k.ge.i+1.and. k.le.(inb(il)+1)) then
2598	amp1(il)=amp1(il)+m(il,k)
2599	endif
2600	else
2601	c AMP1 is the part of cbmf taken from layers I and lower
2602	if(k.le.i) then
2603	amp1(il)=amp1(il)+cbmf(il)*wghti(il,k)
2604	endif
2605	endif
2606	441 continue
2607	440 continue
2608
2609	do 450 k=1,i
2610	do 451 j=i+1,nl+1
2611	do 452 il=1,ncum
2612	if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
2613	amp1(il)=amp1(il)+ment(il,k,j)
2614	endif
2615	452 continue
2616	451 continue
2617	450 continue
2618
2619	do 470 k=1,i-1
2620	do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
2621	do 472 il=1,ncum
2622	if (i.le.inb(il) .and. j.le.inb(il)) then
2623	ad(il)=ad(il)+ment(il,j,k)
2624	endif
2625	472 continue
2626	471 continue
2627	470 continue
2628
2629	do 1350 il=1,ncum
2630	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2631	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2632	cpinv=1.0/cpn(il,i)
2633
2634	c convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
2635	if (cvflag_grav) then
2636	if((0.01gravdpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
2637	else
2638	if((0.1dpinvamp1(il)).ge.delti)iflag(il)=1 ! vecto
2639	endif
2640
2641	! precip
2642	ccc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
2643	ft(il,i)= -sigd(il)lvcp(il,i)evap(il,i)
2644	rat=cpn(il,i-1)*cpinv
2645	c
2646	if (cvflag_grav) then
2647	ft(il,i)=ft(il,i)-0.009gravsigd(il)
2648	: (mp(il,i+1)t_wake(il,i)*b(il,i)
2649	: -mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
2650	ft(il,i)=ft(il,i)+0.01sigd(il)wt(il,i)(cl-cpd)water(il,i+1)
2651	: (t_wake(il,i+1)-t_wake(il,i))dpinv*cpinv
2652	ftd(il,i)=ft(il,i)
2653	! fin precip
2654	c
2655	! sature
2656	ft(il,i)=ft(il,i)+0.01gravdpinv(amp1(il)(t(il,i+1)-t(il,i)
2657	: +(gz(il,i+1)-gz(il,i))*cpinv)
2658	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2659
2660	c
2661	IF (iflag_mix .eq. 0) then
2662	ft(il,i)=ft(il,i)+0.01gravdpinvment(il,i,i)(hp(il,i)-h(il,i)
2663	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2664	ENDIF
2665	c
2666	else ! cvflag_grav
2667	ft(il,i)=ft(il,i)-0.09*sigd(il)
2668	: (mp(il,i+1)t_wake(il,i)*b(il,i)
2669	: -mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
2670	ft(il,i)=ft(il,i)+0.01sigd(il)wt(il,i)(cl-cpd)water(il,i+1)
2671	: (t_wake(il,i+1)-t_wake(il,i))dpinv*cpinv
2672	ftd(il,i)=ft(il,i)
2673	! fin precip
2674	c
2675	! sature
2676	ft(il,i)=ft(il,i)+0.1dpinv(amp1(il)*(t(il,i+1)-t(il,i)
2677	: +(gz(il,i+1)-gz(il,i))*cpinv)
2678	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2679
2680	c
2681	IF (iflag_mix .eq. 0) then
2682	ft(il,i)=ft(il,i)+0.1dpinvment(il,i,i)*(hp(il,i)-h(il,i)
2683	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2684	ENDIF
2685	endif ! cvflag_grav
2686
2687
2688	if (cvflag_grav) then
2689	c sb: on ne fait pas encore la correction permettant de mieux
2690	c conserver l'eau:
2691	c jyg: correction permettant de mieux conserver l'eau:
2692	ccc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
2693	fr(il,i)=sigd(il)*evap(il,i)
2694	: +0.01grav(mp(il,i+1)*(rp(il,i+1)-rr_wake(il,i))
2695	: -mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
2696	fqd(il,i)=fr(il,i) ! precip
2697
2698	fu(il,i)=0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i))
2699	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2700	fv(il,i)=0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i))
2701	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2702	else ! cvflag_grav
2703	ccc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
2704	fr(il,i)=sigd(il)*evap(il,i)
2705	: +0.1(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i))
2706	: -mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
2707	fqd(il,i)=fr(il,i) ! precip
2708
2709	fu(il,i)=0.1(mp(il,i+1)(up(il,i+1)-u(il,i))
2710	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2711	fv(il,i)=0.1(mp(il,i+1)(vp(il,i+1)-v(il,i))
2712	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2713	endif ! cvflag_grav
2714
2715
2716	if (cvflag_grav) then
2717	fr(il,i)=fr(il,i)+0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i))
2718	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2719	fu(il,i)=fu(il,i)+0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i))
2720	: -ad(il)*(u(il,i)-u(il,i-1)))
2721	fv(il,i)=fv(il,i)+0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i))
2722	: -ad(il)*(v(il,i)-v(il,i-1)))
2723	else ! cvflag_grav
2724	fr(il,i)=fr(il,i)+0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il,i))
2725	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2726	fu(il,i)=fu(il,i)+0.1dpinv(amp1(il)*(u(il,i+1)-u(il,i))
2727	: -ad(il)*(u(il,i)-u(il,i-1)))
2728	fv(il,i)=fv(il,i)+0.1dpinv(amp1(il)*(v(il,i+1)-v(il,i))
2729	: -ad(il)*(v(il,i)-v(il,i-1)))
2730	endif ! cvflag_grav
2731
2732	endif ! i
2733	1350 continue
2734
2735	do k=1,ntra
2736	do il=1,ncum
2737	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2738	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2739	cpinv=1.0/cpn(il,i)
2740	if (cvflag_grav) then
2741	ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
2742	: (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2743	: -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2744	else
2745	ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
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	endif
2749	endif
2750	enddo
2751	enddo
2752
2753	do 480 k=1,i-1
2754	c
2755	do il = 1,ncum
2756	awat(il)=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
2757	awat(il)=amax1(awat(il),0.0)
2758	enddo
2759	c
2760	IF (iflag_mix .ne. 0) then
2761	do il=1,ncum
2762	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2763	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2764	cpinv=1.0/cpn(il,i)
2765	if (cvflag_grav) then
2766	ft(il,i)=ft(il,i)
2767	: +0.01gravdpinvment(il,k,i)(hent(il,k,i)-h(il,i)
2768	: +t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-Qent(il,k,i)))*cpinv
2769
2770	c
2771	c
2772	else
2773	ft(il,i)=ft(il,i)
2774	: +0.1dpinvment(il,k,i)*(hent(il,k,i)-h(il,i)
2775	: +t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-Qent(il,k,i)))*cpinv
2776	endif !cvflag_grav
2777	endif ! i
2778	enddo
2779	ENDIF
2780	c
2781	do 1370 il=1,ncum
2782	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2783	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2784	cpinv=1.0/cpn(il,i)
2785	if (cvflag_grav) then
2786	fr(il,i)=fr(il,i)
2787	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-awat(il)-rr(il,i))
2788	fu(il,i)=fu(il,i)
2789	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2790	fv(il,i)=fv(il,i)
2791	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2792	else ! cvflag_grav
2793	fr(il,i)=fr(il,i)
2794	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-awat(il)-rr(il,i))
2795	fu(il,i)=fu(il,i)
2796	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2797	fv(il,i)=fv(il,i)
2798	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2799	endif ! cvflag_grav
2800
2801	c (saturated updrafts resulting from mixing) ! cld
2802	qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat(il)) ! cld
2803	nqcond(il,i)=nqcond(il,i)+1. ! cld
2804	endif ! i
2805	1370 continue
2806	480 continue
2807
2808	do j=1,ntra
2809	do k=1,i-1
2810	do il=1,ncum
2811	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2812	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2813	cpinv=1.0/cpn(il,i)
2814	if (cvflag_grav) then
2815	ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2816	: *(traent(il,k,i,j)-tra(il,i,j))
2817	else
2818	ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2819	: *(traent(il,k,i,j)-tra(il,i,j))
2820	endif
2821	endif
2822	enddo
2823	enddo
2824	enddo
2825
2826	do 490 k=i,nl+1
2827	c
2828	IF (iflag_mix .ne. 0) then
2829	do il=1,ncum
2830	if (i.le.inb(il) .and. k.le.inb(il)
2831	$ .and. iflag(il) .le. 1) then
2832	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2833	cpinv=1.0/cpn(il,i)
2834	if (cvflag_grav) then
2835	ft(il,i)=ft(il,i)
2836	: +0.01gravdpinvment(il,k,i)(hent(il,k,i)-h(il,i)
2837	: +t(il,i)(cpv-cpd)(rr(il,i)-Qent(il,k,i)))*cpinv
2838	c
2839	c
2840	else
2841	ft(il,i)=ft(il,i)
2842	: +0.1dpinvment(il,k,i)*(hent(il,k,i)-h(il,i)
2843	: +t(il,i)(cpv-cpd)(rr(il,i)-Qent(il,k,i)))*cpinv
2844	endif !cvflag_grav
2845	endif ! i
2846	enddo
2847	ENDIF
2848	c
2849	do 1380 il=1,ncum
2850	if (i.le.inb(il) .and. k.le.inb(il)
2851	$ .and. iflag(il) .le. 1) then
2852	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2853	cpinv=1.0/cpn(il,i)
2854
2855	if (cvflag_grav) then
2856	fr(il,i)=fr(il,i)
2857	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-rr(il,i))
2858	fu(il,i)=fu(il,i)
2859	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2860	fv(il,i)=fv(il,i)
2861	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2862	else ! cvflag_grav
2863	fr(il,i)=fr(il,i)
2864	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-rr(il,i))
2865	fu(il,i)=fu(il,i)
2866	: +0.1dpinvment(il,k,i)*(uent(il,k,i)-u(il,i))
2867	fv(il,i)=fv(il,i)
2868	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2869	endif ! cvflag_grav
2870	endif ! i and k
2871	1380 continue
2872	490 continue
2873
2874	do j=1,ntra
2875	do k=i,nl+1
2876	do il=1,ncum
2877	if (i.le.inb(il) .and. k.le.inb(il)
2878	$ .and. iflag(il) .le. 1) then
2879	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2880	cpinv=1.0/cpn(il,i)
2881	if (cvflag_grav) then
2882	ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2883	: *(traent(il,k,i,j)-tra(il,i,j))
2884	else
2885	ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2886	: *(traent(il,k,i,j)-tra(il,i,j))
2887	endif
2888	endif ! i and k
2889	enddo
2890	enddo
2891	enddo
2892
2893	c sb: interface with the cloud parameterization: ! cld
2894
2895	do k=i+1,nl
2896	do il=1,ncum
2897	if (k.le.inb(il) .and. i.le.inb(il)
2898	$ .and. iflag(il) .le. 1) then ! cld
2899	C (saturated downdrafts resulting from mixing) ! cld
2900	qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
2901	nqcond(il,i)=nqcond(il,i)+1. ! cld
2902	endif ! cld
2903	enddo ! cld
2904	enddo ! cld
2905
2906	C (particular case: no detraining level is found) ! cld
2907	do il=1,ncum ! cld
2908	if (i.le.inb(il) .and. nent(il,i).eq.0
2909	$ .and. iflag(il) .le. 1) then ! cld
2910	qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
2911	nqcond(il,i)=nqcond(il,i)+1. ! cld
2912	endif ! cld
2913	enddo ! cld
2914
2915	do il=1,ncum ! cld
2916	if (i.le.inb(il) .and. nqcond(il,i).ne.0
2917	$ .and. iflag(il) .le. 1) then ! cld
2918	qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
2919	endif ! cld
2920	enddo
2921
2922	do j=1,ntra
2923	do il=1,ncum
2924	if (i.le.inb(il) .and. iflag(il) .le. 1) then
2925	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2926	cpinv=1.0/cpn(il,i)
2927
2928	if (cvflag_grav) then
2929	ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
2930	: (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2931	: -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
2932	else
2933	ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
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	endif
2937	endif ! i
2938	enddo
2939	enddo
2940
2941
2942	500 continue
2943
2944
2945	c * move the detrainment at level inb down to level inb-1 *
2946	c * in such a way as to preserve the vertically *
2947	c * integrated enthalpy and water tendencies *
2948	c
2949	do 503 il=1,ncum
2950	IF (iflag(il) .le. 1) THEN
2951	ax=0.1ment(il,inb(il),inb(il))(hp(il,inb(il))-h(il,inb(il))
2952	: +t(il,inb(il))*(cpv-cpd)
2953	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
2954	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2955	ft(il,inb(il))=ft(il,inb(il))-ax
2956	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
2957	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
2958	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
2959
2960	bx=0.1ment(il,inb(il),inb(il))(qent(il,inb(il),inb(il))
2961	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2962	fr(il,inb(il))=fr(il,inb(il))-bx
2963	fr(il,inb(il)-1)=fr(il,inb(il)-1)
2964	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
2965	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2966
2967	cx=0.1ment(il,inb(il),inb(il))(uent(il,inb(il),inb(il))
2968	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2969	fu(il,inb(il))=fu(il,inb(il))-cx
2970	fu(il,inb(il)-1)=fu(il,inb(il)-1)
2971	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
2972	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2973
2974	dx=0.1ment(il,inb(il),inb(il))(vent(il,inb(il),inb(il))
2975	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2976	fv(il,inb(il))=fv(il,inb(il))-dx
2977	fv(il,inb(il)-1)=fv(il,inb(il)-1)
2978	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
2979	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2980	ENDIF !iflag
2981	503 continue
2982
2983	do j=1,ntra
2984	do il=1,ncum
2985	IF (iflag(il) .le. 1) THEN
2986	ex=0.1*ment(il,inb(il),inb(il))
2987	: *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2988	: /(ph(i l,inb(il))-ph(il,inb(il)+1))
2989	ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2990	ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2991	: +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2992	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2993	ENDIF !iflag
2994	enddo
2995	enddo
2996
2997	c
2998	c * homogenize tendencies below cloud base *
2999	c
3000	c
3001	do il=1,ncum
3002	asum(il)=0.0
3003	bsum(il)=0.0
3004	csum(il)=0.0
3005	dsum(il)=0.0
3006	enddo
3007
3008	do i=1,nl
3009	do il=1,ncum
3010	if (i.le.(icb(il)-1) .and. iflag(il) .le. 1) then
3011	asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
3012	bsum(il)=bsum(il)+fr(il,i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))
3013	: *(ph(il,i)-ph(il,i+1))
3014	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
3015	: *(ph(il,i)-ph(il,i+1))
3016	dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
3017	endif
3018	enddo
3019	enddo
3020
3021	c!!! do 700 i=1,icb(il)-1
3022	do i=1,nl
3023	do il=1,ncum
3024	if (i.le.(icb(il)-1) .and. iflag(il) .le. 1) then
3025	ft(il,i)=asum(il)t(il,i)/(th(il,i)dsum(il))
3026	fr(il,i)=bsum(il)/csum(il)
3027	endif
3028	enddo
3029	enddo
3030
3031	c
3032	c * reset counter and return *
3033	c
3034	do il=1,ncum
3035	sig(il,nd)=2.0
3036	enddo
3037
3038
3039	do i=1,nd
3040	do il=1,ncum
3041	upwd(il,i)=0.0
3042	dnwd(il,i)=0.0
3043	enddo
3044	enddo
3045
3046	do i=1,nl
3047	do il=1,ncum
3048	dnwd0(il,i)=-mp(il,i)
3049	enddo
3050	enddo
3051	do i=nl+1,nd
3052	do il=1,ncum
3053	dnwd0(il,i)=0.
3054	enddo
3055	enddo
3056
3057
3058	do i=1,nl
3059	do il=1,ncum
3060	if (i.ge.icb(il) .and. i.le.inb(il)) then
3061	upwd(il,i)=0.0
3062	dnwd(il,i)=0.0
3063	endif
3064	enddo
3065	enddo
3066
3067	do i=1,nl
3068	do k=1,nl
3069	do il=1,ncum
3070	up1(il,k,i)=0.0
3071	dn1(il,k,i)=0.0
3072	enddo
3073	enddo
3074	enddo
3075
3076	do i=1,nl
3077	do k=i,nl
3078	do n=1,i-1
3079	do il=1,ncum
3080	if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
3081	up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
3082	dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
3083	endif
3084	enddo
3085	enddo
3086	enddo
3087	enddo
3088
3089	do i=1,nl
3090	do k=1,nl
3091	do il=1,ncum
3092	if(i.ge.icb(il)) then
3093	if(k.ge.i.and. k.le.(inb(il))) then
3094	upwd(il,i)=upwd(il,i)+m(il,k)
3095	endif
3096	else
3097	if(k.lt.i) then
3098	upwd(il,i)=upwd(il,i)+cbmf(il)*wghti(il,k)
3099	endif
3100	endif
3101	cc print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
3102	end do
3103	end do
3104	end do
3105
3106	do i=2,nl
3107	do k=i,nl
3108	do il=1,ncum
3109	ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
3110	if (i.le.inb(il).and.k.le.inb(il)) then
3111	upwd(il,i)=upwd(il,i)+up1(il,k,i)
3112	dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
3113	endif
3114	cc print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
3115	enddo
3116	enddo
3117	enddo
3118
3119
3120	c!!! DO il=1,ncum
3121	c!!! do i=icb(il),inb(il)
3122	c!!!
3123	c!!! upwd(il,i)=0.0
3124	c!!! dnwd(il,i)=0.0
3125	c!!! do k=i,inb(il)
3126	c!!! up1=0.0
3127	c!!! dn1=0.0
3128	c!!! do n=1,i-1
3129	c!!! up1=up1+ment(il,n,k)
3130	c!!! dn1=dn1-ment(il,k,n)
3131	c!!! enddo
3132	c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
3133	c!!! dnwd(il,i)=dnwd(il,i)+dn1
3134	c!!! enddo
3135	c!!! enddo
3136	c!!!
3137	c!!! ENDDO
3138
3139	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3140	c determination de la variation de flux ascendant entre
3141	c deux niveau non dilue mip
3142	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3143
3144	do i=1,nl
3145	do il=1,ncum
3146	mip(il,i)=m(il,i)
3147	enddo
3148	enddo
3149
3150	do i=nl+1,nd
3151	do il=1,ncum
3152	mip(il,i)=0.
3153	enddo
3154	enddo
3155
3156	do i=1,nd
3157	do il=1,ncum
3158	ma(il,i)=0
3159	enddo
3160	enddo
3161
3162	do i=1,nl
3163	do j=i,nl
3164	do il=1,ncum
3165	ma(il,i)=ma(il,i)+m(il,j)
3166	enddo
3167	enddo
3168	enddo
3169
3170	do i=nl+1,nd
3171	do il=1,ncum
3172	ma(il,i)=0.
3173	enddo
3174	enddo
3175
3176	do i=1,nl
3177	do il=1,ncum
3178	if (i.le.(icb(il)-1)) then
3179	ma(il,i)=0
3180	endif
3181	enddo
3182	enddo
3183
3184	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3185	c icb represente de niveau ou se trouve la
3186	c base du nuage , et inb le top du nuage
3187	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3188
3189	do i=1,nd
3190	do il=1,ncum
3191	mke(il,i)=upwd(il,i)+dnwd(il,i)
3192	enddo
3193	enddo
3194
3195	do i=1,nd
3196	DO 999 il=1,ncum
3197	rdcp=(rrd(1.-rr(il,i))-rr(il,i)rrv)
3198	: /(cpd(1.-rr(il,i))+rr(il,i)cpv)
3199	tls(il,i)=t(il,i)(1000.0/p(il,i))*rdcp
3200	tps(il,i)=tp(il,i)
3201	999 CONTINUE
3202	enddo
3203
3204	c
3205	c * diagnose the in-cloud mixing ratio * ! cld
3206	c * of condensed water * ! cld
3207	c ! cld
3208
3209	do i=1,nd ! cld
3210	do il=1,ncum ! cld
3211	mac(il,i)=0.0 ! cld
3212	wa(il,i)=0.0 ! cld
3213	siga(il,i)=0.0 ! cld
3214	sax(il,i)=0.0 ! cld
3215	enddo ! cld
3216	enddo ! cld
3217
3218	do i=minorig, nl ! cld
3219	do k=i+1,nl+1 ! cld
3220	do il=1,ncum ! cld
3221	if (i.le.inb(il) .and. k.le.(inb(il)+1)
3222	$ .and. iflag(il) .le. 1) then ! cld
3223	mac(il,i)=mac(il,i)+m(il,k) ! cld
3224	endif ! cld
3225	enddo ! cld
3226	enddo ! cld
3227	enddo ! cld
3228
3229	do i=1,nl ! cld
3230	do j=1,i ! cld
3231	do il=1,ncum ! cld
3232	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
3233	: .and. j.ge.icb(il)
3234	: .and. iflag(il) .le. 1 ) then ! cld
3235	sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
3236	: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
3237	endif ! cld
3238	enddo ! cld
3239	enddo ! cld
3240	enddo ! cld
3241
3242	do i=1,nl ! cld
3243	do il=1,ncum ! cld
3244	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
3245	: .and. sax(il,i).gt.0.0
3246	: .and. iflag(il) .le. 1 ) then ! cld
3247	wa(il,i)=sqrt(2.*sax(il,i)) ! cld
3248	endif ! cld
3249	enddo ! cld
3250	enddo ! cld
3251
3252	do i=1,nl ! cld
3253	do il=1,ncum ! cld
3254	if (wa(il,i).gt.0.0 .and. iflag(il) .le. 1) ! cld
3255	: siga(il,i)=mac(il,i)/wa(il,i) ! cld
3256	: rrdtvp(il,i)/p(il,i)/100./delta ! cld
3257	siga(il,i) = min(siga(il,i),1.0) ! cld
3258	cIM cf. FH
3259	if (iflag_clw.eq.0) then
3260	qcondc(il,i)=siga(il,i)clw(il,i)(1.-ep(il,i)) ! cld
3261	: + (1.-siga(il,i))*qcond(il,i) ! cld
3262	else if (iflag_clw.eq.1) then
3263	qcondc(il,i)=qcond(il,i) ! cld
3264	endif
3265
3266	enddo ! cld
3267	enddo
3268	c print*,'cv3_yield fin'
3269	! cld
3270	return
3271	end
3272
3273
3274	SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
3275	: ,iflag
3276	: ,precip,sig,w0
3277	: ,ft,fq,fu,fv,ftra
3278	: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
3279	: ,iflag1
3280	: ,precip1,sig1,w01
3281	: ,ft1,fq1,fu1,fv1,ftra1
3282	: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
3283	: )
3284	implicit none
3285
3286	#include "cv3param.h"
3287
3288	c inputs:
3289	integer len, ncum, nd, ntra, nloc
3290	integer idcum(nloc)
3291	integer iflag(nloc)
3292	real precip(nloc)
3293	real sig(nloc,nd), w0(nloc,nd)
3294	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
3295	real ftra(nloc,nd,ntra)
3296	real Ma(nloc,nd)
3297	real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
3298	real qcondc(nloc,nd)
3299	real wd(nloc),cape(nloc)
3300
3301	c outputs:
3302	integer iflag1(len)
3303	real precip1(len)
3304	real sig1(len,nd), w01(len,nd)
3305	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
3306	real ftra1(len,nd,ntra)
3307	real Ma1(len,nd)
3308	real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
3309	real qcondc1(nloc,nd)
3310	real wd1(nloc),cape1(nloc)
3311
3312	c local variables:
3313	integer i,k,j
3314
3315	do 2000 i=1,ncum
3316	precip1(idcum(i))=precip(i)
3317	iflag1(idcum(i))=iflag(i)
3318	wd1(idcum(i))=wd(i)
3319	cape1(idcum(i))=cape(i)
3320	2000 continue
3321
3322	do 2020 k=1,nl
3323	do 2010 i=1,ncum
3324	sig1(idcum(i),k)=sig(i,k)
3325	w01(idcum(i),k)=w0(i,k)
3326	ft1(idcum(i),k)=ft(i,k)
3327	fq1(idcum(i),k)=fq(i,k)
3328	fu1(idcum(i),k)=fu(i,k)
3329	fv1(idcum(i),k)=fv(i,k)
3330	Ma1(idcum(i),k)=Ma(i,k)
3331	upwd1(idcum(i),k)=upwd(i,k)
3332	dnwd1(idcum(i),k)=dnwd(i,k)
3333	dnwd01(idcum(i),k)=dnwd0(i,k)
3334	qcondc1(idcum(i),k)=qcondc(i,k)
3335	2010 continue
3336	2020 continue
3337
3338	do 2200 i=1,ncum
3339	sig1(idcum(i),nd)=sig(i,nd)
3340	2200 continue
3341
3342
3343	do 2100 j=1,ntra
3344	c oct3 do 2110 k=1,nl
3345	do 2110 k=1,nd ! oct3
3346	do 2120 i=1,ncum
3347	ftra1(idcum(i),k,j)=ftra(i,k,j)
3348	2120 continue
3349	2110 continue
3350	2100 continue
3351	return
3352	end
3353

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

Download in other formats:

Original Format