Context Navigation

cv3_routines.F @ 973

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

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

Download in other formats:

Original Format