Context Navigation

cv3_routines.F @ 433

Last change on this file since 433 was 433, checked in by lmdzadmin, 21 years ago

Convergence avec la version de Ionela dec 2002

YOMCST.? : suppression RI0 (IM)
albedo.F : facteur 1.2 sur le nouveau calcul (IM)
clesphys.h : rajout de différentes ctes (concentration des gaz) (IM)
clmain.F : separation des flux LW, SW (JLD)

remplace qsurf par yqsol (IM)

conf_phys.F90 : rajout de différentes ctes (gaz + orbite) (IM)
convect3.F : DPINV+SIGD*0.5*(EVAP(1)+EVAP(2)) (SBL)
cv3_routines.F:
cvparam3.h : compatibilite avec conema3 TEMPORAIRE (FH)
phyetat0.F : lecture de co2_ppm et solaire pour tests de coherence
phyredem.F : co2_ppm et solaire passé en common
physiq.F : separation flux LW, SW

rajout diagnostiques (slp, w500)
suppression iflag_con = 4
clwcon0=qcondc (FH)
position dU "ENDIF ! ok_cvl"

radlwsw.F : passage des concentrations gaz dans un common (IM)

PEMIS(i) = 1.0 (JLD pour cohérence ORCHIDEE)

stdlevvar.F90 :
suphec.F : suppression init. des ctes orbitales (IM)

nouvelles E/S (ini_hist..., write_hist...)

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 90.7 KB

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

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

Download in other formats:

Original Format