Context Navigation

cv30_routines.F @ 5412

Last change on this file since 5412 was 1744, checked in by Ehouarn Millour, 12 years ago
Fix to avoid dividing by zero. JYG
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 98.4 KB

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

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

Download in other formats:

Original Format