Context Navigation

cv3_routines.F @ 1042

Last change on this file since 1042 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 93.6 KB

Line
1	!
2	! $Header$
3	!
4	c
5	c
6	SUBROUTINE cv3_param(nd,delt)
7	implicit none
8
9	c------------------------------------------------------------
10	c Set parameters for convectL for iflag_con = 3
11	c------------------------------------------------------------
12
13	C
14	C * PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *
15	C * PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *
16	C * PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *
17	C * EFFICIENCY IS ASSUMED TO BE UNITY *
18	C * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
19	C * SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
20	C * OF CLOUD *
21	C
22	C [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
23	C * ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *
24	C * APPROACH TO QUASI-EQUILIBRIUM *
25	C * (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *
26	C * (BETA MUST BE LESS THAN OR EQUAL TO 1) *
27	C
28	C * DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *
29	C * APPROACH TO QUASI-EQUILIBRIUM *
30	C * IT MUST BE LESS THAN 0 *
31
32	#include "cvparam3.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 cv3_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 "cvparam3.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 cv3_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 "cvparam3.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 cv3_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 "cvparam3.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 cv3_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 "cvparam3.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 cv3_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 "cvparam3.h"
751
752	c inputs:
753	integer len,ncum,nd,ntra,nloc
754	integer iflag1(len),nk1(len),icb1(len),icbs1(len)
755	real plcl1(len),tnk1(len),qnk1(len),gznk1(len)
756	real pbase1(len),buoybase1(len)
757	real t1(len,nd),q1(len,nd),qs1(len,nd),u1(len,nd),v1(len,nd)
758	real gz1(len,nd),h1(len,nd),lv1(len,nd),cpn1(len,nd)
759	real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd)
760	real tvp1(len,nd),clw1(len,nd)
761	real th1(len,nd)
762	real sig1(len,nd), w01(len,nd)
763	real tra1(len,nd,ntra)
764
765	c outputs:
766	c en fait, on a nloc=len pour l'instant (cf cv_driver)
767	integer iflag(nloc),nk(nloc),icb(nloc),icbs(nloc)
768	real plcl(nloc),tnk(nloc),qnk(nloc),gznk(nloc)
769	real pbase(nloc),buoybase(nloc)
770	real t(nloc,nd),q(nloc,nd),qs(nloc,nd),u(nloc,nd),v(nloc,nd)
771	real gz(nloc,nd),h(nloc,nd),lv(nloc,nd),cpn(nloc,nd)
772	real p(nloc,nd),ph(nloc,nd+1),tv(nloc,nd),tp(nloc,nd)
773	real tvp(nloc,nd),clw(nloc,nd)
774	real th(nloc,nd)
775	real sig(nloc,nd), w0(nloc,nd)
776	real tra(nloc,nd,ntra)
777
778	c local variables:
779	integer i,k,nn,j
780
781
782	do 110 k=1,nl+1
783	nn=0
784	do 100 i=1,len
785	if(iflag1(i).eq.0)then
786	nn=nn+1
787	sig(nn,k)=sig1(i,k)
788	w0(nn,k)=w01(i,k)
789	t(nn,k)=t1(i,k)
790	q(nn,k)=q1(i,k)
791	qs(nn,k)=qs1(i,k)
792	u(nn,k)=u1(i,k)
793	v(nn,k)=v1(i,k)
794	gz(nn,k)=gz1(i,k)
795	h(nn,k)=h1(i,k)
796	lv(nn,k)=lv1(i,k)
797	cpn(nn,k)=cpn1(i,k)
798	p(nn,k)=p1(i,k)
799	ph(nn,k)=ph1(i,k)
800	tv(nn,k)=tv1(i,k)
801	tp(nn,k)=tp1(i,k)
802	tvp(nn,k)=tvp1(i,k)
803	clw(nn,k)=clw1(i,k)
804	th(nn,k)=th1(i,k)
805	endif
806	100 continue
807	110 continue
808
809	c do 121 j=1,ntra
810	c do 111 k=1,nd
811	c nn=0
812	c do 101 i=1,len
813	c if(iflag1(i).eq.0)then
814	c nn=nn+1
815	c tra(nn,k,j)=tra1(i,k,j)
816	c endif
817	c 101 continue
818	c 111 continue
819	c 121 continue
820
821	if (nn.ne.ncum) then
822	print*,'strange! nn not equal to ncum: ',nn,ncum
823	stop
824	endif
825
826	nn=0
827	do 150 i=1,len
828	if(iflag1(i).eq.0)then
829	nn=nn+1
830	pbase(nn)=pbase1(i)
831	buoybase(nn)=buoybase1(i)
832	plcl(nn)=plcl1(i)
833	tnk(nn)=tnk1(i)
834	qnk(nn)=qnk1(i)
835	gznk(nn)=gznk1(i)
836	nk(nn)=nk1(i)
837	icb(nn)=icb1(i)
838	icbs(nn)=icbs1(i)
839	iflag(nn)=iflag1(i)
840	endif
841	150 continue
842
843	return
844	end
845
846	SUBROUTINE cv3_undilute2(nloc,ncum,nd,icb,icbs,nk
847	: ,tnk,qnk,gznk,t,q,qs,gz
848	: ,p,h,tv,lv,pbase,buoybase,plcl
849	o ,inb,tp,tvp,clw,hp,ep,sigp,buoy)
850	implicit none
851
852	C---------------------------------------------------------------------
853	C Purpose:
854	C FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
855	C &
856	C COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
857	C FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
858	C &
859	C FIND THE LEVEL OF NEUTRAL BUOYANCY
860	C
861	C Main differences convect3/convect4:
862	C - icbs (input) is the first level above LCL (may differ from icb)
863	C - many minor differences in the iterations
864	C - condensed water not removed from tvp in convect3
865	C - vertical profile of buoyancy computed here (use of buoybase)
866	C - the determination of inb is different
867	C - no inb1, only inb in output
868	C---------------------------------------------------------------------
869
870	#include "cvthermo.h"
871	#include "cvparam3.h"
872	#include "conema3.h"
873
874	c inputs:
875	integer ncum, nd, nloc
876	integer icb(nloc), icbs(nloc), nk(nloc)
877	real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
878	real p(nloc,nd)
879	real tnk(nloc), qnk(nloc), gznk(nloc)
880	real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
881	real pbase(nloc), buoybase(nloc), plcl(nloc)
882
883	c outputs:
884	integer inb(nloc)
885	real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)
886	real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)
887	real buoy(nloc,nd)
888
889	c local variables:
890	integer i, k
891	real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
892	real by, defrac, pden
893	real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
894	logical lcape(nloc)
895
896	!=====================================================================
897	! --- SOME INITIALIZATIONS
898	!=====================================================================
899
900	do 170 k=1,nl
901	do 160 i=1,ncum
902	ep(i,k)=0.0
903	sigp(i,k)=spfac
904	160 continue
905	170 continue
906
907	!=====================================================================
908	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
909	!=====================================================================
910	c
911	c --- The procedure is to solve the equation.
912	c cptp+Lqp+phi=cptnk+Lqnk+gznk.
913	c
914	c * Calculate certain parcel quantities, including static energy *
915	c
916	c
917	do 240 i=1,ncum
918	ah0(i)=(cpd(1.-qnk(i))+clqnk(i))*tnk(i)
919	cdebug & +qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
920	& +qnk(i)(lv0-clmcpv(tnk(i)-273.15))+gznk(i)
921	240 continue
922	c
923	c
924	c * Find lifted parcel quantities above cloud base *
925	c
926	c
927	do 300 k=minorig+1,nl
928	do 290 i=1,ncum
929	c ori if(k.ge.(icb(i)+1))then
930	if(k.ge.(icbs(i)+1))then ! convect3
931	tg=t(i,k)
932	qg=qs(i,k)
933	cdebug alv=lv0-clmcpv*(t(i,k)-t0)
934	alv=lv0-clmcpv*(t(i,k)-273.15)
935	c
936	c First iteration.
937	c
938	c ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
939	s=cpd(1.-qnk(i))+clqnk(i) ! convect3
940	: +alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
941	s=1./s
942	c ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
943	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gz(i,k) ! convect3
944	tg=tg+s*(ah0(i)-ahg)
945	c ori tg=max(tg,35.0)
946	cdebug tc=tg-t0
947	tc=tg-273.15
948	denom=243.5+tc
949	denom=MAX(denom,1.0) ! convect3
950	c ori if(tc.ge.0.0)then
951	es=6.112exp(17.67tc/denom)
952	c ori else
953	c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
954	c ori endif
955	qg=epses/(p(i,k)-es(1.-eps))
956	c
957	c Second iteration.
958	c
959	c ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
960	c ori s=1./s
961	c ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
962	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gz(i,k) ! convect3
963	tg=tg+s*(ah0(i)-ahg)
964	c ori tg=max(tg,35.0)
965	cdebug tc=tg-t0
966	tc=tg-273.15
967	denom=243.5+tc
968	denom=MAX(denom,1.0) ! convect3
969	c ori if(tc.ge.0.0)then
970	es=6.112exp(17.67tc/denom)
971	c ori else
972	c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
973	c ori endif
974	qg=epses/(p(i,k)-es(1.-eps))
975	c
976	cdebug alv=lv0-clmcpv*(t(i,k)-t0)
977	alv=lv0-clmcpv*(t(i,k)-273.15)
978	c print*,'cpd dans convect2 ',cpd
979	c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
980	c print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
981
982	c ori c approximation here:
983	c ori tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
984
985	c convect3: no approximation:
986	tp(i,k)=(ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
987
988	clw(i,k)=qnk(i)-qg
989	clw(i,k)=max(0.0,clw(i,k))
990	rg=qg/(1.-qnk(i))
991	c ori tvp(i,k)=tp(i,k)(1.+rgepsi)
992	c convect3: (qg utilise au lieu du vrai mixing ratio rg):
993	tvp(i,k)=tp(i,k)*(1.+qg/eps-qnk(i)) ! whole thing
994	endif
995	290 continue
996	300 continue
997	c
998	!=====================================================================
999	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
1000	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
1001	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
1002	!=====================================================================
1003	c
1004	c ori do 320 k=minorig+1,nl
1005	do 320 k=1,nl ! convect3
1006	do 310 i=1,ncum
1007	pden=ptcrit-pbcrit
1008	ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax
1009	ep(i,k)=amax1(ep(i,k),0.0)
1010	ep(i,k)=amin1(ep(i,k),epmax)
1011	sigp(i,k)=spfac
1012	c ori if(k.ge.(nk(i)+1))then
1013	c ori tca=tp(i,k)-t0
1014	c ori if(tca.ge.0.0)then
1015	c ori elacrit=elcrit
1016	c ori else
1017	c ori elacrit=elcrit*(1.0-tca/tlcrit)
1018	c ori endif
1019	c ori elacrit=max(elacrit,0.0)
1020	c ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
1021	c ori ep(i,k)=max(ep(i,k),0.0 )
1022	c ori ep(i,k)=min(ep(i,k),1.0 )
1023	c ori sigp(i,k)=sigs
1024	c ori endif
1025	310 continue
1026	320 continue
1027	c
1028	!=====================================================================
1029	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1030	! --- VIRTUAL TEMPERATURE
1031	!=====================================================================
1032	c
1033	c dans convect3, tvp est calcule en une seule fois, et sans retirer
1034	c l'eau condensee (~> reversible CAPE)
1035	c
1036	c ori do 340 k=minorig+1,nl
1037	c ori do 330 i=1,ncum
1038	c ori if(k.ge.(icb(i)+1))then
1039	c ori tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
1040	c oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
1041	c oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1042	c ori endif
1043	c ori 330 continue
1044	c ori 340 continue
1045
1046	c ori do 350 i=1,ncum
1047	c ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1048	c ori 350 continue
1049
1050	do 350 i=1,ncum ! convect3
1051	tp(i,nlp)=tp(i,nl) ! convect3
1052	350 continue ! convect3
1053	c
1054	c=====================================================================
1055	c --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1056	c=====================================================================
1057
1058	c-- this is for convect3 only:
1059
1060	c first estimate of buoyancy:
1061
1062	do 500 i=1,ncum
1063	do 501 k=1,nl
1064	buoy(i,k)=tvp(i,k)-tv(i,k)
1065	501 continue
1066	500 continue
1067
1068	c set buoyancy=buoybase for all levels below base
1069	c for safety, set buoy(icb)=buoybase
1070
1071	do 505 i=1,ncum
1072	do 506 k=1,nl
1073	if((k.ge.icb(i)).and.(k.le.nl).and.(p(i,k).ge.pbase(i)))then
1074	buoy(i,k)=buoybase(i)
1075	endif
1076	506 continue
1077	cIM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i)
1078	buoy(i,icb(i))=buoybase(i)
1079	505 continue
1080
1081	c-- end convect3
1082
1083	c=====================================================================
1084	c --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1085	c --- LEVEL OF NEUTRAL BUOYANCY
1086	c=====================================================================
1087	c
1088	c-- this is for convect3 only:
1089
1090	do 510 i=1,ncum
1091	inb(i)=nl-1
1092	510 continue
1093
1094	do 530 i=1,ncum
1095	do 535 k=1,nl-1
1096	if ((k.ge.icb(i)).and.(buoy(i,k).lt.dtovsh)) then
1097	inb(i)=MIN(inb(i),k)
1098	endif
1099	535 continue
1100	530 continue
1101
1102	c-- end convect3
1103
1104	c ori do 510 i=1,ncum
1105	c ori cape(i)=0.0
1106	c ori capem(i)=0.0
1107	c ori inb(i)=icb(i)+1
1108	c ori inb1(i)=inb(i)
1109	c ori 510 continue
1110	c
1111	c Originial Code
1112	c
1113	c do 530 k=minorig+1,nl-1
1114	c do 520 i=1,ncum
1115	c if(k.ge.(icb(i)+1))then
1116	c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1117	c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1118	c cape(i)=cape(i)+by
1119	c if(by.ge.0.0)inb1(i)=k+1
1120	c if(cape(i).gt.0.0)then
1121	c inb(i)=k+1
1122	c capem(i)=cape(i)
1123	c endif
1124	c endif
1125	c520 continue
1126	c530 continue
1127	c do 540 i=1,ncum
1128	c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1129	c cape(i)=capem(i)+byp
1130	c defrac=capem(i)-cape(i)
1131	c defrac=max(defrac,0.001)
1132	c frac(i)=-cape(i)/defrac
1133	c frac(i)=min(frac(i),1.0)
1134	c frac(i)=max(frac(i),0.0)
1135	c540 continue
1136	c
1137	c K Emanuel fix
1138	c
1139	c call zilch(byp,ncum)
1140	c do 530 k=minorig+1,nl-1
1141	c do 520 i=1,ncum
1142	c if(k.ge.(icb(i)+1))then
1143	c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1144	c cape(i)=cape(i)+by
1145	c if(by.ge.0.0)inb1(i)=k+1
1146	c if(cape(i).gt.0.0)then
1147	c inb(i)=k+1
1148	c capem(i)=cape(i)
1149	c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1150	c endif
1151	c endif
1152	c520 continue
1153	c530 continue
1154	c do 540 i=1,ncum
1155	c inb(i)=max(inb(i),inb1(i))
1156	c cape(i)=capem(i)+byp(i)
1157	c defrac=capem(i)-cape(i)
1158	c defrac=max(defrac,0.001)
1159	c frac(i)=-cape(i)/defrac
1160	c frac(i)=min(frac(i),1.0)
1161	c frac(i)=max(frac(i),0.0)
1162	c540 continue
1163	c
1164	c J Teixeira fix
1165	c
1166	c ori call zilch(byp,ncum)
1167	c ori do 515 i=1,ncum
1168	c ori lcape(i)=.true.
1169	c ori 515 continue
1170	c ori do 530 k=minorig+1,nl-1
1171	c ori do 520 i=1,ncum
1172	c ori if(cape(i).lt.0.0)lcape(i)=.false.
1173	c ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1174	c ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1175	c ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1176	c ori cape(i)=cape(i)+by
1177	c ori if(by.ge.0.0)inb1(i)=k+1
1178	c ori if(cape(i).gt.0.0)then
1179	c ori inb(i)=k+1
1180	c ori capem(i)=cape(i)
1181	c ori endif
1182	c ori endif
1183	c ori 520 continue
1184	c ori 530 continue
1185	c ori do 540 i=1,ncum
1186	c ori cape(i)=capem(i)+byp(i)
1187	c ori defrac=capem(i)-cape(i)
1188	c ori defrac=max(defrac,0.001)
1189	c ori frac(i)=-cape(i)/defrac
1190	c ori frac(i)=min(frac(i),1.0)
1191	c ori frac(i)=max(frac(i),0.0)
1192	c ori 540 continue
1193	c
1194	c=====================================================================
1195	c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1196	c=====================================================================
1197	c
1198	cym do i=1,ncum*nlp
1199	cym hp(i,1)=h(i,1)
1200	cym enddo
1201
1202	do k=1,nlp
1203	do i=1,ncum
1204	hp(i,k)=h(i,k)
1205	enddo
1206	enddo
1207
1208	do 600 k=minorig+1,nl
1209	do 590 i=1,ncum
1210	if((k.ge.icb(i)).and.(k.le.inb(i)))then
1211	hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
1212	endif
1213	590 continue
1214	600 continue
1215
1216	return
1217	end
1218
1219	SUBROUTINE cv3_closure(nloc,ncum,nd,icb,inb
1220	: ,pbase,p,ph,tv,buoy
1221	o ,sig,w0,cape,m)
1222	implicit none
1223
1224	!===================================================================
1225	! --- CLOSURE OF CONVECT3
1226	!
1227	! vectorization: S. Bony
1228	!===================================================================
1229
1230	#include "cvthermo.h"
1231	#include "cvparam3.h"
1232
1233	c input:
1234	integer ncum, nd, nloc
1235	integer icb(nloc), inb(nloc)
1236	real pbase(nloc)
1237	real p(nloc,nd), ph(nloc,nd+1)
1238	real tv(nloc,nd), buoy(nloc,nd)
1239
1240	c input/output:
1241	real sig(nloc,nd), w0(nloc,nd)
1242
1243	c output:
1244	real cape(nloc)
1245	real m(nloc,nd)
1246
1247	c local variables:
1248	integer i, j, k, icbmax
1249	real deltap, fac, w, amu
1250	real dtmin(nloc,nd), sigold(nloc,nd)
1251
1252
1253	c -------------------------------------------------------
1254	c -- Initialization
1255	c -------------------------------------------------------
1256
1257	do k=1,nl
1258	do i=1,ncum
1259	m(i,k)=0.0
1260	enddo
1261	enddo
1262
1263	c -------------------------------------------------------
1264	c -- Reset sig(i) and w0(i) for i>inb and i<icb
1265	c -------------------------------------------------------
1266
1267	c update sig and w0 above LNB:
1268
1269	do 100 k=1,nl-1
1270	do 110 i=1,ncum
1271	if ((inb(i).lt.(nl-1)).and.(k.ge.(inb(i)+1)))then
1272	sig(i,k)=beta*sig(i,k)
1273	: +2.alphabuoy(i,inb(i))*ABS(buoy(i,inb(i)))
1274	sig(i,k)=AMAX1(sig(i,k),0.0)
1275	w0(i,k)=beta*w0(i,k)
1276	endif
1277	110 continue
1278	100 continue
1279
1280	c compute icbmax:
1281
1282	icbmax=2
1283	do 200 i=1,ncum
1284	icbmax=MAX(icbmax,icb(i))
1285	200 continue
1286
1287	c update sig and w0 below cloud base:
1288
1289	do 300 k=1,icbmax
1290	do 310 i=1,ncum
1291	if (k.le.icb(i))then
1292	sig(i,k)=betasig(i,k)-2.alphabuoy(i,icb(i))buoy(i,icb(i))
1293	sig(i,k)=amax1(sig(i,k),0.0)
1294	w0(i,k)=beta*w0(i,k)
1295	endif
1296	310 continue
1297	300 continue
1298
1299	c! if(inb.lt.(nl-1))then
1300	c! do 85 i=inb+1,nl-1
1301	c! sig(i)=betasig(i)+2.alphabuoy(inb)
1302	c! 1 abs(buoy(inb))
1303	c! sig(i)=amax1(sig(i),0.0)
1304	c! w0(i)=beta*w0(i)
1305	c! 85 continue
1306	c! end if
1307
1308	c! do 87 i=1,icb
1309	c! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
1310	c! sig(i)=amax1(sig(i),0.0)
1311	c! w0(i)=beta*w0(i)
1312	c! 87 continue
1313
1314	c -------------------------------------------------------------
1315	c -- Reset fractional areas of updrafts and w0 at initial time
1316	c -- and after 10 time steps of no convection
1317	c -------------------------------------------------------------
1318
1319	do 400 k=1,nl-1
1320	do 410 i=1,ncum
1321	if (sig(i,nd).lt.1.5.or.sig(i,nd).gt.12.0)then
1322	sig(i,k)=0.0
1323	w0(i,k)=0.0
1324	endif
1325	410 continue
1326	400 continue
1327
1328	c -------------------------------------------------------------
1329	c -- Calculate convective available potential energy (cape),
1330	c -- vertical velocity (w), fractional area covered by
1331	c -- undilute updraft (sig), and updraft mass flux (m)
1332	c -------------------------------------------------------------
1333
1334	do 500 i=1,ncum
1335	cape(i)=0.0
1336	500 continue
1337
1338	c compute dtmin (minimum buoyancy between ICB and given level k):
1339
1340	do i=1,ncum
1341	do k=1,nl
1342	dtmin(i,k)=100.0
1343	enddo
1344	enddo
1345
1346	do 550 i=1,ncum
1347	do 560 k=1,nl
1348	do 570 j=minorig,nl
1349	if ( (k.ge.(icb(i)+1)).and.(k.le.inb(i)).and.
1350	: (j.ge.icb(i)).and.(j.le.(k-1)) )then
1351	dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1352	endif
1353	570 continue
1354	560 continue
1355	550 continue
1356
1357	c the interval on which cape is computed starts at pbase :
1358
1359	do 600 k=1,nl
1360	do 610 i=1,ncum
1361
1362	if ((k.ge.(icb(i)+1)).and.(k.le.inb(i))) then
1363
1364	deltap = MIN(pbase(i),ph(i,k-1))-MIN(pbase(i),ph(i,k))
1365	cape(i)=cape(i)+rrdbuoy(i,k-1)deltap/p(i,k-1)
1366	cape(i)=AMAX1(0.0,cape(i))
1367	sigold(i,k)=sig(i,k)
1368
1369	c dtmin(i,k)=100.0
1370	c do 97 j=icb(i),k-1 ! mauvaise vectorisation
1371	c dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1372	c 97 continue
1373
1374	sig(i,k)=betasig(i,k)+alphadtmin(i,k)*ABS(dtmin(i,k))
1375	sig(i,k)=amax1(sig(i,k),0.0)
1376	sig(i,k)=amin1(sig(i,k),0.01)
1377	fac=AMIN1(((dtcrit-dtmin(i,k))/dtcrit),1.0)
1378	w=(1.-beta)facSQRT(cape(i))+beta*w0(i,k)
1379	amu=0.5(sig(i,k)+sigold(i,k))w
1380	m(i,k)=amu0.007p(i,k)*(ph(i,k)-ph(i,k+1))/tv(i,k)
1381	w0(i,k)=w
1382	endif
1383
1384	610 continue
1385	600 continue
1386
1387	do 700 i=1,ncum
1388	w0(i,icb(i))=0.5*w0(i,icb(i)+1)
1389	m(i,icb(i))=0.5*m(i,icb(i)+1)
1390	: *(ph(i,icb(i))-ph(i,icb(i)+1))
1391	: /(ph(i,icb(i)+1)-ph(i,icb(i)+2))
1392	sig(i,icb(i))=sig(i,icb(i)+1)
1393	sig(i,icb(i)-1)=sig(i,icb(i))
1394	700 continue
1395
1396
1397	c! cape=0.0
1398	c! do 98 i=icb+1,inb
1399	c! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1400	c! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
1401	c! dcape=rrdbuoy(i-1)deltap/p(i-1)
1402	c! dlnp=deltap/p(i-1)
1403	c! cape=amax1(0.0,cape)
1404	c! sigold=sig(i)
1405
1406	c! dtmin=100.0
1407	c! do 97 j=icb,i-1
1408	c! dtmin=amin1(dtmin,buoy(j))
1409	c! 97 continue
1410
1411	c! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
1412	c! sig(i)=amax1(sig(i),0.0)
1413	c! sig(i)=amin1(sig(i),0.01)
1414	c! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1415	c! w=(1.-beta)facsqrt(cape)+beta*w0(i)
1416	c! amu=0.5(sig(i)+sigold)w
1417	c! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
1418	c! w0(i)=w
1419	c! 98 continue
1420	c! w0(icb)=0.5*w0(icb+1)
1421	c! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1422	c! sig(icb)=sig(icb+1)
1423	c! sig(icb-1)=sig(icb)
1424
1425	return
1426	end
1427
1428	SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb
1429	: ,ph,t,rr,rs,u,v,tra,h,lv,qnk
1430	: ,hp,tv,tvp,ep,clw,m,sig
1431	: ,ment,qent,uent,vent,sij,elij,ments,qents,traent)
1432	implicit none
1433
1434	!---------------------------------------------------------------------
1435	! a faire:
1436	! - changer rr(il,1) -> qnk(il)
1437	! - vectorisation de la partie normalisation des flux (do 789...)
1438	!---------------------------------------------------------------------
1439
1440	#include "cvthermo.h"
1441	#include "cvparam3.h"
1442
1443	c inputs:
1444	integer ncum, nd, na, ntra, nloc
1445	integer icb(nloc), inb(nloc), nk(nloc)
1446	real sig(nloc,nd)
1447	real qnk(nloc)
1448	real ph(nloc,nd+1)
1449	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
1450	real u(nloc,nd), v(nloc,nd)
1451	real tra(nloc,nd,ntra) ! input of convect3
1452	real lv(nloc,na), h(nloc,na), hp(nloc,na)
1453	real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
1454	real m(nloc,na) ! input of convect3
1455
1456	c outputs:
1457	real ment(nloc,na,na), qent(nloc,na,na)
1458	real uent(nloc,na,na), vent(nloc,na,na)
1459	real sij(nloc,na,na), elij(nloc,na,na)
1460	real traent(nloc,nd,nd,ntra)
1461	real ments(nloc,nd,nd), qents(nloc,nd,nd)
1462	real sigij(nloc,nd,nd)
1463
1464	c local variables:
1465	integer i, j, k, il, im, jm
1466	integer num1, num2
1467	integer nent(nloc,na)
1468	real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1469	real alt, smid, sjmin, sjmax, delp, delm
1470	real asij(nloc), smax(nloc), scrit(nloc)
1471	real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd)
1472	real wgh
1473	real zm(nloc,na)
1474	logical lwork(nloc)
1475
1476	c=====================================================================
1477	c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1478	c=====================================================================
1479
1480	c ori do 360 i=1,ncum*nlp
1481	do 361 j=1,nl
1482	do 360 i=1,ncum
1483	nent(i,j)=0
1484	c in convect3, m is computed in cv3_closure
1485	c ori m(i,1)=0.0
1486	360 continue
1487	361 continue
1488
1489	c ori do 400 k=1,nlp
1490	c ori do 390 j=1,nlp
1491	do 400 j=1,nl
1492	do 390 k=1,nl
1493	do 385 i=1,ncum
1494	qent(i,k,j)=rr(i,j)
1495	uent(i,k,j)=u(i,j)
1496	vent(i,k,j)=v(i,j)
1497	elij(i,k,j)=0.0
1498	cym ment(i,k,j)=0.0
1499	cym sij(i,k,j)=0.0
1500	385 continue
1501	390 continue
1502	400 continue
1503
1504	cym
1505	ment(1:ncum,1:nd,1:nd)=0.0
1506	sij(1:ncum,1:nd,1:nd)=0.0
1507
1508	c do k=1,ntra
1509	c do j=1,nd ! instead nlp
1510	c do i=1,nd ! instead nlp
1511	c do il=1,ncum
1512	c traent(il,i,j,k)=tra(il,j,k)
1513	c enddo
1514	c enddo
1515	c enddo
1516	c enddo
1517	zm(:,:)=0.
1518
1519	c=====================================================================
1520	c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1521	c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1522	c --- FRACTION (sij)
1523	c=====================================================================
1524
1525	do 750 i=minorig+1, nl
1526
1527	do 710 j=minorig,nl
1528	do 700 il=1,ncum
1529	if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1530	: (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1531
1532	rti=rr(il,1)-ep(il,i)*clw(il,i)
1533	bf2=1.+lv(il,j)lv(il,j)rs(il,j)/(rrvt(il,j)t(il,j)*cpd)
1534	anum=h(il,j)-hp(il,i)+(cpv-cpd)t(il,j)(rti-rr(il,j))
1535	denom=h(il,i)-hp(il,i)+(cpd-cpv)(rr(il,i)-rti)t(il,j)
1536	dei=denom
1537	if(abs(dei).lt.0.01)dei=0.01
1538	sij(il,i,j)=anum/dei
1539	sij(il,i,i)=1.0
1540	altem=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti-rs(il,j)
1541	altem=altem/bf2
1542	cwat=clw(il,j)*(1.-ep(il,j))
1543	stemp=sij(il,i,j)
1544	if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat)
1545	: .and.j.gt.i)then
1546	anum=anum-lv(il,j)(rti-rs(il,j)-cwatbf2)
1547	denom=denom+lv(il,j)*(rr(il,i)-rti)
1548	if(abs(denom).lt.0.01)denom=0.01
1549	sij(il,i,j)=anum/denom
1550	altem=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti-rs(il,j)
1551	altem=altem-(bf2-1.)*cwat
1552	end if
1553	if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then
1554	qent(il,i,j)=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti
1555	uent(il,i,j)=sij(il,i,j)u(il,i)+(1.-sij(il,i,j))u(il,nk(il))
1556	vent(il,i,j)=sij(il,i,j)v(il,i)+(1.-sij(il,i,j))v(il,nk(il))
1557	c!!! do k=1,ntra
1558	c!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1559	c!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1560	c!!! end do
1561	elij(il,i,j)=altem
1562	elij(il,i,j)=amax1(0.0,elij(il,i,j))
1563	ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
1564	nent(il,i)=nent(il,i)+1
1565	end if
1566	sij(il,i,j)=amax1(0.0,sij(il,i,j))
1567	sij(il,i,j)=amin1(1.0,sij(il,i,j))
1568	endif ! new
1569	700 continue
1570	710 continue
1571
1572	c do k=1,ntra
1573	c do j=minorig,nl
1574	c do il=1,ncum
1575	c if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1576	c : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1577	c traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1578	c : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1579	c endif
1580	c enddo
1581	c enddo
1582	c enddo
1583
1584	c
1585	c * if no air can entrain at level i assume that updraft detrains *
1586	c * at that level and calculate detrained air flux and properties *
1587	c
1588
1589	c@ do 170 i=icb(il),inb(il)
1590
1591	do 740 il=1,ncum
1592	if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then
1593	c@ if(nent(il,i).eq.0)then
1594	ment(il,i,i)=m(il,i)
1595	qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)
1596	uent(il,i,i)=u(il,nk(il))
1597	vent(il,i,i)=v(il,nk(il))
1598	elij(il,i,i)=clw(il,i)
1599	cMAF sij(il,i,i)=1.0
1600	sij(il,i,i)=0.0
1601	end if
1602	740 continue
1603	750 continue
1604
1605	c do j=1,ntra
1606	c do i=minorig+1,nl
1607	c do il=1,ncum
1608	c if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1609	c traent(il,i,i,j)=tra(il,nk(il),j)
1610	c endif
1611	c enddo
1612	c enddo
1613	c enddo
1614
1615	do 100 j=minorig,nl
1616	do 101 i=minorig,nl
1617	do 102 il=1,ncum
1618	if ((j.ge.(icb(il)-1)).and.(j.le.inb(il))
1619	: .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
1620	sigij(il,i,j)=sij(il,i,j)
1621	endif
1622	102 continue
1623	101 continue
1624	100 continue
1625	c@ enddo
1626
1627	c@170 continue
1628
1629	c=====================================================================
1630	c --- NORMALIZE ENTRAINED AIR MASS FLUXES
1631	c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
1632	c=====================================================================
1633
1634	cym call zilch(asum,ncum*nd)
1635	cym call zilch(bsum,ncum*nd)
1636	cym call zilch(csum,ncum*nd)
1637	call zilch(asum,nloc*nd)
1638	call zilch(csum,nloc*nd)
1639	call zilch(csum,nloc*nd)
1640
1641	do il=1,ncum
1642	lwork(il) = .FALSE.
1643	enddo
1644
1645	DO 789 i=minorig+1,nl
1646
1647	num1=0
1648	do il=1,ncum
1649	if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
1650	enddo
1651	if (num1.le.0) goto 789
1652
1653
1654	do 781 il=1,ncum
1655	if ( i.ge.icb(il) .and. i.le.inb(il) ) then
1656	lwork(il)=(nent(il,i).ne.0)
1657	qp=rr(il,1)-ep(il,i)*clw(il,i)
1658	anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
1659	: +(cpv-cpd)t(il,i)(qp-rr(il,i))
1660	denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp)
1661	: +(cpd-cpv)t(il,i)(rr(il,i)-qp)
1662	if(abs(denom).lt.0.01)denom=0.01
1663	scrit(il)=anum/denom
1664	alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
1665	if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0
1666	smax(il)=0.0
1667	asij(il)=0.0
1668	endif
1669	781 continue
1670
1671	do 175 j=nl,minorig,-1
1672
1673	num2=0
1674	do il=1,ncum
1675	if ( i.ge.icb(il) .and. i.le.inb(il) .and.
1676	: j.ge.(icb(il)-1) .and. j.le.inb(il)
1677	: .and. lwork(il) ) num2=num2+1
1678	enddo
1679	if (num2.le.0) goto 175
1680
1681	do 782 il=1,ncum
1682	if ( i.ge.icb(il) .and. i.le.inb(il) .and.
1683	: j.ge.(icb(il)-1) .and. j.le.inb(il)
1684	: .and. lwork(il) ) then
1685
1686	if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then
1687	wgh=1.0
1688	if(j.gt.i)then
1689	sjmax=amax1(sij(il,i,j+1),smax(il))
1690	sjmax=amin1(sjmax,scrit(il))
1691	smax(il)=amax1(sij(il,i,j),smax(il))
1692	sjmin=amax1(sij(il,i,j-1),smax(il))
1693	sjmin=amin1(sjmin,scrit(il))
1694	if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
1695	smid=amin1(sij(il,i,j),scrit(il))
1696	else
1697	sjmax=amax1(sij(il,i,j+1),scrit(il))
1698	smid=amax1(sij(il,i,j),scrit(il))
1699	sjmin=0.0
1700	if(j.gt.1)sjmin=sij(il,i,j-1)
1701	sjmin=amax1(sjmin,scrit(il))
1702	endif
1703	delp=abs(sjmax-smid)
1704	delm=abs(sjmin-smid)
1705	asij(il)=asij(il)+wgh*(delp+delm)
1706	ment(il,i,j)=ment(il,i,j)(delp+delm)wgh
1707	endif
1708	endif
1709	782 continue
1710
1711	175 continue
1712
1713	do il=1,ncum
1714	if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
1715	asij(il)=amax1(1.0e-16,asij(il))
1716	asij(il)=1.0/asij(il)
1717	asum(il,i)=0.0
1718	bsum(il,i)=0.0
1719	csum(il,i)=0.0
1720	endif
1721	enddo
1722
1723	do 180 j=minorig,nl
1724	do il=1,ncum
1725	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1726	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1727	ment(il,i,j)=ment(il,i,j)*asij(il)
1728	endif
1729	enddo
1730	180 continue
1731
1732	do 190 j=minorig,nl
1733	do il=1,ncum
1734	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1735	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1736	asum(il,i)=asum(il,i)+ment(il,i,j)
1737	ment(il,i,j)=ment(il,i,j)*sig(il,j)
1738	bsum(il,i)=bsum(il,i)+ment(il,i,j)
1739	endif
1740	enddo
1741	190 continue
1742
1743	do il=1,ncum
1744	if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
1745	bsum(il,i)=amax1(bsum(il,i),1.0e-16)
1746	bsum(il,i)=1.0/bsum(il,i)
1747	endif
1748	enddo
1749
1750	do 195 j=minorig,nl
1751	do il=1,ncum
1752	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1753	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1754	ment(il,i,j)=ment(il,i,j)asum(il,i)bsum(il,i)
1755	endif
1756	enddo
1757	195 continue
1758
1759	do 197 j=minorig,nl
1760	do il=1,ncum
1761	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1762	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1763	csum(il,i)=csum(il,i)+ment(il,i,j)
1764	endif
1765	enddo
1766	197 continue
1767
1768	do il=1,ncum
1769	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1770	: .and. csum(il,i).lt.m(il,i) ) then
1771	nent(il,i)=0
1772	ment(il,i,i)=m(il,i)
1773	qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)
1774	uent(il,i,i)=u(il,nk(il))
1775	vent(il,i,i)=v(il,nk(il))
1776	elij(il,i,i)=clw(il,i)
1777	cMAF sij(il,i,i)=1.0
1778	sij(il,i,i)=0.0
1779	endif
1780	enddo ! il
1781
1782	c do j=1,ntra
1783	c do il=1,ncum
1784	c if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1785	c : .and. csum(il,i).lt.m(il,i) ) then
1786	c traent(il,i,i,j)=tra(il,nk(il),j)
1787	c endif
1788	c enddo
1789	c enddo
1790	789 continue
1791	c
1792	c MAF: renormalisation de MENT
1793	do jm=1,nd
1794	do im=1,nd
1795	do il=1,ncum
1796	zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)
1797	end do
1798	end do
1799	end do
1800	c
1801	do jm=1,nd
1802	do im=1,nd
1803	do il=1,ncum
1804	if(zm(il,im).ne.0.) then
1805	ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)
1806	endif
1807	end do
1808	end do
1809	end do
1810	c
1811	do jm=1,nd
1812	do im=1,nd
1813	do 999 il=1,ncum
1814	qents(il,im,jm)=qent(il,im,jm)
1815	ments(il,im,jm)=ment(il,im,jm)
1816	999 continue
1817	enddo
1818	enddo
1819
1820	return
1821	end
1822
1823
1824	SUBROUTINE cv3_unsat(nloc,ncum,nd,na,ntra,icb,inb
1825	: ,t,rr,rs,gz,u,v,tra,p,ph
1826	: ,th,tv,lv,cpn,ep,sigp,clw
1827	: ,m,ment,elij,delt,plcl
1828	: ,mp,rp,up,vp,trap,wt,water,evap,b)
1829	implicit none
1830
1831
1832	#include "cvthermo.h"
1833	#include "cvparam3.h"
1834	#include "cvflag.h"
1835
1836	c inputs:
1837	integer ncum, nd, na, ntra, nloc
1838	integer icb(nloc), inb(nloc)
1839	real delt, plcl(nloc)
1840	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
1841	real u(nloc,nd), v(nloc,nd)
1842	real tra(nloc,nd,ntra)
1843	real p(nloc,nd), ph(nloc,nd+1)
1844	real th(nloc,na), gz(nloc,na)
1845	real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na)
1846	real cpn(nloc,na), tv(nloc,na)
1847	real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
1848
1849	c outputs:
1850	real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
1851	real water(nloc,na), evap(nloc,na), wt(nloc,na)
1852	real trap(nloc,na,ntra)
1853	real b(nloc,na)
1854
1855	c local variables
1856	integer i,j,k,il,num1
1857	real tinv, delti
1858	real awat, afac, afac1, afac2, bfac
1859	real pr1, pr2, sigt, b6, c6, revap, tevap, delth
1860	real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1861	real ampmax
1862	real lvcp(nloc,na)
1863	real wdtrain(nloc)
1864	logical lwork(nloc)
1865
1866
1867	c------------------------------------------------------
1868
1869	delti = 1./delt
1870	tinv=1./3.
1871
1872	mp(:,:)=0.
1873
1874	do i=1,nl
1875	do il=1,ncum
1876	mp(il,i)=0.0
1877	rp(il,i)=rr(il,i)
1878	up(il,i)=u(il,i)
1879	vp(il,i)=v(il,i)
1880	wt(il,i)=0.001
1881	water(il,i)=0.0
1882	evap(il,i)=0.0
1883	b(il,i)=0.0
1884	lvcp(il,i)=lv(il,i)/cpn(il,i)
1885	enddo
1886	enddo
1887
1888	c do k=1,ntra
1889	c do i=1,nd
1890	c do il=1,ncum
1891	c trap(il,i,k)=tra(il,i,k)
1892	c enddo
1893	c enddo
1894	c enddo
1895
1896	c
1897	c * check whether ep(inb)=0, if so, skip precipitating *
1898	c * downdraft calculation *
1899	c
1900
1901	do il=1,ncum
1902	lwork(il)=.TRUE.
1903	if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
1904	enddo
1905
1906	call zilch(wdtrain,ncum)
1907
1908	DO 400 i=nl+1,1,-1
1909
1910	num1=0
1911	do il=1,ncum
1912	if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1
1913	enddo
1914	if (num1.le.0) goto 400
1915
1916	c
1917	c * integrate liquid water equation to find condensed water *
1918	c * and condensed water flux *
1919	c
1920
1921	c
1922	c * begin downdraft loop *
1923	c
1924
1925	c
1926	c * calculate detrained precipitation *
1927	c
1928	do il=1,ncum
1929	if (i.le.inb(il) .and. lwork(il)) then
1930	if (cvflag_grav) then
1931	wdtrain(il)=gravep(il,i)m(il,i)*clw(il,i)
1932	else
1933	wdtrain(il)=10.0ep(il,i)m(il,i)*clw(il,i)
1934	endif
1935	endif
1936	enddo
1937
1938	if(i.gt.1)then
1939	do 320 j=1,i-1
1940	do il=1,ncum
1941	if (i.le.inb(il) .and. lwork(il)) then
1942	awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
1943	awat=amax1(awat,0.0)
1944	if (cvflag_grav) then
1945	wdtrain(il)=wdtrain(il)+gravawatment(il,j,i)
1946	else
1947	wdtrain(il)=wdtrain(il)+10.0awatment(il,j,i)
1948	endif
1949	endif
1950	enddo
1951	320 continue
1952	endif
1953
1954	c
1955	c * find rain water and evaporation using provisional *
1956	c * estimates of rp(i)and rp(i-1) *
1957	c
1958
1959	do 999 il=1,ncum
1960
1961	if (i.le.inb(il) .and. lwork(il)) then
1962
1963	wt(il,i)=45.0
1964
1965	if(i.lt.inb(il))then
1966	rp(il,i)=rp(il,i+1)
1967	: +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
1968	rp(il,i)=0.5*(rp(il,i)+rr(il,i))
1969	endif
1970	rp(il,i)=amax1(rp(il,i),0.0)
1971	rp(il,i)=amin1(rp(il,i),rs(il,i))
1972	rp(il,inb(il))=rr(il,inb(il))
1973
1974	if(i.eq.1)then
1975	afac=p(il,1)(rs(il,1)-rp(il,1))/(1.0e4+2000.0p(il,1)*rs(il,1))
1976	else
1977	rp(il,i-1)=rp(il,i)
1978	: +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i)
1979	rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
1980	rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
1981	rp(il,i-1)=amax1(rp(il,i-1),0.0)
1982	afac1=p(il,i)(rs(il,i)-rp(il,i))/(1.0e4+2000.0p(il,i)*rs(il,i))
1983	afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
1984	: /(1.0e4+2000.0p(il,i-1)rs(il,i-1))
1985	afac=0.5*(afac1+afac2)
1986	endif
1987	if(i.eq.inb(il))afac=0.0
1988	afac=amax1(afac,0.0)
1989	bfac=1./(sigd*wt(il,i))
1990	c
1991	cjyg1
1992	ccc sigt=1.0
1993	ccc if(i.ge.icb)sigt=sigp(i)
1994	c prise en compte de la variation progressive de sigt dans
1995	c les couches icb et icb-1:
1996	c pour plcl<ph(i+1), pr1=0 & pr2=1
1997	c pour plcl>ph(i), pr1=1 & pr2=0
1998	c pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
1999	c sur le nuage, et pr2 est la proportion sous la base du
2000	c nuage.
2001	pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2002	pr1=max(0.,min(1.,pr1))
2003	pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2004	pr2=max(0.,min(1.,pr2))
2005	sigt=sigp(il,i)*pr1+pr2
2006	cjyg2
2007	c
2008	b6=bfac50.sigd(ph(il,i)-ph(il,i+1))sigt*afac
2009	c6=water(il,i+1)+bfac*wdtrain(il)
2010	: -50.sigdbfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2011	if(c6.gt.0.0)then
2012	revap=0.5(-b6+sqrt(b6b6+4.*c6))
2013	evap(il,i)=sigtafacrevap
2014	water(il,i)=revap*revap
2015	else
2016	evap(il,i)=-evap(il,i+1)
2017	: +0.02(wdtrain(il)+sigdwt(il,i)*water(il,i+1))
2018	: /(sigd*(ph(il,i)-ph(il,i+1)))
2019	end if
2020	c
2021	c * calculate precipitating downdraft mass flux under *
2022	c * hydrostatic approximation *
2023	c
2024	if (i.ne.1) then
2025
2026	tevap=amax1(0.0,evap(il,i))
2027	delth=amax1(0.001,(th(il,i)-th(il,i-1)))
2028	if (cvflag_grav) then
2029	mp(il,i)=100.ginvlvcp(il,i)sigdtevap
2030	: *(p(il,i-1)-p(il,i))/delth
2031	else
2032	mp(il,i)=10.lvcp(il,i)sigdtevap(p(il,i-1)-p(il,i))/delth
2033	endif
2034	c
2035	c * if hydrostatic assumption fails, *
2036	c * solve cubic difference equation for downdraft theta *
2037	c * and mass flux from two simultaneous differential eqns *
2038	c
2039	amfac=sigdsigd70.0ph(il,i)(p(il,i-1)-p(il,i))
2040	: (th(il,i)-th(il,i-1))/(tv(il,i)th(il,i))
2041	amp2=abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2042	if(amp2.gt.(0.1*amfac))then
2043	xf=100.0sigdsigdsigd(ph(il,i)-ph(il,i+1))
2044	tf=b(il,i)-5.0(th(il,i)-th(il,i-1))t(il,i)
2045	: /(lvcp(il,i)sigdth(il,i))
2046	af=xftf+mp(il,i+1)mp(il,i+1)*tinv
2047	bf=2.(tinvmp(il,i+1))*3+tinvmp(il,i+1)xftf
2048	: +50.(p(il,i-1)-p(il,i))xf*tevap
2049	fac2=1.0
2050	if(bf.lt.0.0)fac2=-1.0
2051	bf=abs(bf)
2052	ur=0.25bfbf-afafaftinvtinv*tinv
2053	if(ur.ge.0.0)then
2054	sru=sqrt(ur)
2055	fac=1.0
2056	if((0.5*bf-sru).lt.0.0)fac=-1.0
2057	mp(il,i)=mp(il,i+1)tinv+(0.5bf+sru)**tinv
2058	: +fac(abs(0.5bf-sru))**tinv
2059	else
2060	d=atan(2.*sqrt(-ur)/(bf+1.0e-28))
2061	if(fac2.lt.0.0)d=3.14159-d
2062	mp(il,i)=mp(il,i+1)tinv+2.sqrt(aftinv)cos(d*tinv)
2063	endif
2064	mp(il,i)=amax1(0.0,mp(il,i))
2065
2066	if (cvflag_grav) then
2067	Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2068	C il faut diviser par (mp(il,i)sigdgrav) et non par (mp(il,i)+sigd*0.1).
2069	C Et il faut bien revoir les facteurs 100.
2070	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2071	2 /(mp(il,i)+sigd*0.1)
2072	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
2073	else
2074	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2075	2 /(mp(il,i)+sigd*0.1)
2076	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
2077	endif
2078	b(il,i-1)=amax1(b(il,i-1),0.0)
2079	endif
2080	c
2081	c * limit magnitude of mp(i) to meet cfl condition *
2082	c
2083	ampmax=2.0(ph(il,i)-ph(il,i+1))delti
2084	amp2=2.0(ph(il,i-1)-ph(il,i))delti
2085	ampmax=amin1(ampmax,amp2)
2086	mp(il,i)=amin1(mp(il,i),ampmax)
2087	c
2088	c * force mp to decrease linearly to zero *
2089	c * between cloud base and the surface *
2090	c
2091	if(p(il,i).gt.p(il,icb(il)))then
2092	mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2093	endif
2094
2095	360 continue
2096	endif ! i.eq.1
2097	c
2098	c * find mixing ratio of precipitating downdraft *
2099	c
2100
2101	if (i.ne.inb(il)) then
2102
2103	rp(il,i)=rr(il,i)
2104
2105	if(mp(il,i).gt.mp(il,i+1))then
2106
2107	if (cvflag_grav) then
2108	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2109	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
2110	: *(evap(il,i+1)+evap(il,i))
2111	else
2112	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2113	: +5.sigd(ph(il,i)-ph(il,i+1))
2114	: *(evap(il,i+1)+evap(il,i))
2115	endif
2116	rp(il,i)=rp(il,i)/mp(il,i)
2117	up(il,i)=up(il,i+1)mp(il,i+1)+u(il,i)(mp(il,i)-mp(il,i+1))
2118	up(il,i)=up(il,i)/mp(il,i)
2119	vp(il,i)=vp(il,i+1)mp(il,i+1)+v(il,i)(mp(il,i)-mp(il,i+1))
2120	vp(il,i)=vp(il,i)/mp(il,i)
2121
2122	c do j=1,ntra
2123	c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2124	ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2125	c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2126	c trap(il,i,j)=trap(il,i,j)/mp(il,i)
2127	c end do
2128
2129	else
2130
2131	if(mp(il,i+1).gt.1.0e-16)then
2132	if (cvflag_grav) then
2133	rp(il,i)=rp(il,i+1)
2134	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
2135	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2136	else
2137	rp(il,i)=rp(il,i+1)
2138	: +5.sigd(ph(il,i)-ph(il,i+1))
2139	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2140	endif
2141	up(il,i)=up(il,i+1)
2142	vp(il,i)=vp(il,i+1)
2143
2144	c do j=1,ntra
2145	c trap(il,i,j)=trap(il,i+1,j)
2146	c end do
2147
2148	endif
2149	endif
2150	rp(il,i)=amin1(rp(il,i),rs(il,i))
2151	rp(il,i)=amax1(rp(il,i),0.0)
2152
2153	endif
2154	endif
2155	999 continue
2156
2157	400 continue
2158
2159	return
2160	end
2161
2162	SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
2163	: ,icb,inb,delt
2164	: ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th
2165	: ,ep,clw,m,tp,mp,rp,up,vp,trap
2166	: ,wt,water,evap,b
2167	: ,ment,qent,uent,vent,nent,elij,traent,sig
2168	: ,tv,tvp
2169	: ,iflag,precip,VPrecip,ft,fr,fu,fv,ftra
2170	: ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
2171	implicit none
2172
2173	#include "cvthermo.h"
2174	#include "cvparam3.h"
2175	#include "cvflag.h"
2176	#include "conema3.h"
2177
2178	c inputs:
2179	integer ncum,nd,na,ntra,nloc
2180	integer icb(nloc), inb(nloc)
2181	real delt
2182	real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
2183	real tra(nloc,nd,ntra), sig(nloc,nd)
2184	real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
2185	real th(nloc,na), p(nloc,nd), tp(nloc,na)
2186	real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
2187	real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
2188	real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
2189	real water(nloc,na), evap(nloc,na), b(nloc,na)
2190	real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
2191	cym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2192	real vent(nloc,na,na), elij(nloc,na,na)
2193	integer nent(nloc,na)
2194	real traent(nloc,na,na,ntra)
2195	real tv(nloc,nd), tvp(nloc,nd)
2196
2197	c input/output:
2198	integer iflag(nloc)
2199
2200	c outputs:
2201	real precip(nloc)
2202	real VPrecip(nloc,nd+1)
2203	real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
2204	real ftra(nloc,nd,ntra)
2205	real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
2206	real dnwd0(nloc,nd), mike(nloc,nd)
2207	real tls(nloc,nd), tps(nloc,nd)
2208	real qcondc(nloc,nd) ! cld
2209	real wd(nloc) ! gust
2210
2211	c local variables:
2212	integer i,k,il,n,j,num1
2213	real rat, awat, delti
2214	real ax, bx, cx, dx, ex
2215	real cpinv, rdcp, dpinv
2216	real lvcp(nloc,na), mke(nloc,na)
2217	real am(nloc), work(nloc), ad(nloc), amp1(nloc)
2218	c!! real up1(nloc), dn1(nloc)
2219	real up1(nloc,nd,nd), dn1(nloc,nd,nd)
2220	real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2221	real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
2222	real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
2223
2224
2225	c-------------------------------------------------------------
2226
2227	c initialization:
2228
2229	delti = 1.0/delt
2230
2231	do il=1,ncum
2232	precip(il)=0.0
2233	wd(il)=0.0 ! gust
2234	VPrecip(il,nd+1)=0.
2235	enddo
2236
2237	do i=1,nd
2238	do il=1,ncum
2239	VPrecip(il,i)=0.0
2240	ft(il,i)=0.0
2241	fr(il,i)=0.0
2242	fu(il,i)=0.0
2243	fv(il,i)=0.0
2244	qcondc(il,i)=0.0 ! cld
2245	qcond(il,i)=0.0 ! cld
2246	nqcond(il,i)=0.0 ! cld
2247	enddo
2248	enddo
2249
2250	c do j=1,ntra
2251	c do i=1,nd
2252	c do il=1,ncum
2253	c ftra(il,i,j)=0.0
2254	c enddo
2255	c enddo
2256	c enddo
2257
2258	do i=1,nl
2259	do il=1,ncum
2260	lvcp(il,i)=lv(il,i)/cpn(il,i)
2261	enddo
2262	enddo
2263
2264
2265	c
2266	c * calculate surface precipitation in mm/day *
2267	c
2268	do il=1,ncum
2269	if(ep(il,inb(il)).ge.0.0001)then
2270	if (cvflag_grav) then
2271	precip(il)=wt(il,1)sigdwater(il,1)86400.1000./(rowl*grav)
2272	else
2273	precip(il)=wt(il,1)sigdwater(il,1)*8640.
2274	endif
2275	endif
2276	enddo
2277
2278	C *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
2279	C
2280	c MAF rajout pour lessivage
2281	do k=1,nl
2282	do il=1,ncum
2283	if (k.le.inb(il)) then
2284	if (cvflag_grav) then
2285	VPrecip(il,k) = wt(il,k)sigdwater(il,k)/grav
2286	else
2287	VPrecip(il,k) = wt(il,k)sigdwater(il,k)/10.
2288	endif
2289	endif
2290	end do
2291	end do
2292	C
2293	c
2294	c * Calculate downdraft velocity scale *
2295	c * NE PAS UTILISER POUR L'INSTANT *
2296	c
2297	c! do il=1,ncum
2298	c! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
2299	c! : /(sigd*p(il,icb(il)))
2300	c! enddo
2301
2302	c
2303	c * calculate tendencies of lowest level potential temperature *
2304	c * and mixing ratio *
2305	c
2306	do il=1,ncum
2307	work(il)=1.0/(ph(il,1)-ph(il,2))
2308	am(il)=0.0
2309	enddo
2310
2311	do k=2,nl
2312	do il=1,ncum
2313	if (k.le.inb(il)) then
2314	am(il)=am(il)+m(il,k)
2315	endif
2316	enddo
2317	enddo
2318
2319	do il=1,ncum
2320
2321	c convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
2322	if (cvflag_grav) then
2323	if((0.01gravwork(il)*am(il)).ge.delti)iflag(il)=1!consist vect
2324	ft(il,1)=0.01gravwork(il)am(il)(t(il,2)-t(il,1)
2325	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2326	else
2327	if((0.1work(il)am(il)).ge.delti)iflag(il)=1 !consistency vect
2328	ft(il,1)=0.1work(il)am(il)*(t(il,2)-t(il,1)
2329	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2330	endif
2331
2332	ft(il,1)=ft(il,1)-0.5lvcp(il,1)sigd*(evap(il,1)+evap(il,2))
2333
2334	if (cvflag_grav) then
2335	ft(il,1)=ft(il,1)-0.009gravsigd*mp(il,2)
2336	: t(il,1)b(il,1)*work(il)
2337	else
2338	ft(il,1)=ft(il,1)-0.09sigdmp(il,2)t(il,1)b(il,1)*work(il)
2339	endif
2340
2341	ft(il,1)=ft(il,1)+0.01sigdwt(il,1)(cl-cpd)water(il,2)*(t(il,2)
2342	:-t(il,1))*work(il)/cpn(il,1)
2343
2344	if (cvflag_grav) then
2345	Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
2346	c (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
2347	fr(il,1)=0.01gravmp(il,2)(rp(il,2)-rr(il,1))work(il)
2348	: +sigd0.5(evap(il,1)+evap(il,2))
2349	c+tard : +sigd*evap(il,1)
2350
2351	fr(il,1)=fr(il,1)+0.01gravam(il)(rr(il,2)-rr(il,1))work(il)
2352
2353	fu(il,1)=fu(il,1)+0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1))
2354	: +am(il)*(u(il,2)-u(il,1)))
2355	fv(il,1)=fv(il,1)+0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1))
2356	: +am(il)*(v(il,2)-v(il,1)))
2357	else ! cvflag_grav
2358	fr(il,1)=0.1mp(il,2)(rp(il,2)-rr(il,1))*work(il)
2359	: +sigd0.5(evap(il,1)+evap(il,2))
2360	fr(il,1)=fr(il,1)+0.1am(il)(rr(il,2)-rr(il,1))*work(il)
2361	fu(il,1)=fu(il,1)+0.1work(il)(mp(il,2)*(up(il,2)-u(il,1))
2362	: +am(il)*(u(il,2)-u(il,1)))
2363	fv(il,1)=fv(il,1)+0.1work(il)(mp(il,2)*(vp(il,2)-v(il,1))
2364	: +am(il)*(v(il,2)-v(il,1)))
2365	endif ! cvflag_grav
2366
2367	enddo ! il
2368
2369	c do j=1,ntra
2370	c do il=1,ncum
2371	c if (cvflag_grav) then
2372	c ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
2373	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2374	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2375	c else
2376	c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2377	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2378	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2379	c endif
2380	c enddo
2381	c enddo
2382
2383	do j=2,nl
2384	do il=1,ncum
2385	if (j.le.inb(il)) then
2386	if (cvflag_grav) then
2387	fr(il,1)=fr(il,1)
2388	: +0.01gravwork(il)ment(il,j,1)(qent(il,j,1)-rr(il,1))
2389	fu(il,1)=fu(il,1)
2390	: +0.01gravwork(il)ment(il,j,1)(uent(il,j,1)-u(il,1))
2391	fv(il,1)=fv(il,1)
2392	: +0.01gravwork(il)ment(il,j,1)(vent(il,j,1)-v(il,1))
2393	else ! cvflag_grav
2394	fr(il,1)=fr(il,1)
2395	: +0.1work(il)ment(il,j,1)*(qent(il,j,1)-rr(il,1))
2396	fu(il,1)=fu(il,1)
2397	: +0.1work(il)ment(il,j,1)*(uent(il,j,1)-u(il,1))
2398	fv(il,1)=fv(il,1)
2399	: +0.1work(il)ment(il,j,1)*(vent(il,j,1)-v(il,1))
2400	endif ! cvflag_grav
2401	endif ! j
2402	enddo
2403	enddo
2404
2405	c do k=1,ntra
2406	c do j=2,nl
2407	c do il=1,ncum
2408	c if (j.le.inb(il)) then
2409
2410	c if (cvflag_grav) then
2411	c ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
2412	c : *(traent(il,j,1,k)-tra(il,1,k))
2413	c else
2414	c ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
2415	c : *(traent(il,j,1,k)-tra(il,1,k))
2416	c endif
2417
2418	c endif
2419	c enddo
2420	c enddo
2421	c enddo
2422
2423	c
2424	c * calculate tendencies of potential temperature and mixing ratio *
2425	c * at levels above the lowest level *
2426	c
2427	c * first find the net saturated updraft and downdraft mass fluxes *
2428	c * through each level *
2429	c
2430
2431	do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
2432
2433	num1=0
2434	do il=1,ncum
2435	if(i.le.inb(il))num1=num1+1
2436	enddo
2437	if(num1.le.0)go to 500
2438
2439	call zilch(amp1,ncum)
2440	call zilch(ad,ncum)
2441
2442	do 440 k=i+1,nl+1
2443	do 441 il=1,ncum
2444	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
2445	amp1(il)=amp1(il)+m(il,k)
2446	endif
2447	441 continue
2448	440 continue
2449
2450	do 450 k=1,i
2451	do 451 j=i+1,nl+1
2452	do 452 il=1,ncum
2453	if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
2454	amp1(il)=amp1(il)+ment(il,k,j)
2455	endif
2456	452 continue
2457	451 continue
2458	450 continue
2459
2460	do 470 k=1,i-1
2461	do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
2462	do 472 il=1,ncum
2463	if (i.le.inb(il) .and. j.le.inb(il)) then
2464	ad(il)=ad(il)+ment(il,j,k)
2465	endif
2466	472 continue
2467	471 continue
2468	470 continue
2469
2470	do 1350 il=1,ncum
2471	if (i.le.inb(il)) then
2472	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2473	cpinv=1.0/cpn(il,i)
2474
2475	c convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
2476	if (cvflag_grav) then
2477	if((0.01gravdpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
2478	else
2479	if((0.1dpinvamp1(il)).ge.delti)iflag(il)=1 ! vecto
2480	endif
2481
2482	if (cvflag_grav) then
2483	ft(il,i)=0.01gravdpinv(amp1(il)(t(il,i+1)-t(il,i)
2484	: +(gz(il,i+1)-gz(il,i))*cpinv)
2485	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2486	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
2487	rat=cpn(il,i-1)*cpinv
2488	ft(il,i)=ft(il,i)-0.009gravsigd(mp(il,i+1)t(il,i)*b(il,i)
2489	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2490	ft(il,i)=ft(il,i)+0.01gravdpinvment(il,i,i)(hp(il,i)-h(il,i)
2491	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2492	else ! cvflag_grav
2493	ft(il,i)=0.1dpinv(amp1(il)*(t(il,i+1)-t(il,i)
2494	: +(gz(il,i+1)-gz(il,i))*cpinv)
2495	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2496	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
2497	rat=cpn(il,i-1)*cpinv
2498	ft(il,i)=ft(il,i)-0.09sigd(mp(il,i+1)t(il,i)b(il,i)
2499	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2500	ft(il,i)=ft(il,i)+0.1dpinvment(il,i,i)*(hp(il,i)-h(il,i)
2501	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2502	endif ! cvflag_grav
2503
2504
2505	ft(il,i)=ft(il,i)+0.01sigdwt(il,i)(cl-cpd)water(il,i+1)
2506	: (t(il,i+1)-t(il,i))dpinv*cpinv
2507
2508	if (cvflag_grav) then
2509	fr(il,i)=0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i))
2510	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2511	fu(il,i)=fu(il,i)+0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i))
2512	: -ad(il)*(u(il,i)-u(il,i-1)))
2513	fv(il,i)=fv(il,i)+0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i))
2514	: -ad(il)*(v(il,i)-v(il,i-1)))
2515	else ! cvflag_grav
2516	fr(il,i)=0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il,i))
2517	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2518	fu(il,i)=fu(il,i)+0.1dpinv(amp1(il)*(u(il,i+1)-u(il,i))
2519	: -ad(il)*(u(il,i)-u(il,i-1)))
2520	fv(il,i)=fv(il,i)+0.1dpinv(amp1(il)*(v(il,i+1)-v(il,i))
2521	: -ad(il)*(v(il,i)-v(il,i-1)))
2522	endif ! cvflag_grav
2523
2524	endif ! i
2525	1350 continue
2526
2527	c do k=1,ntra
2528	c do il=1,ncum
2529	c if (i.le.inb(il)) then
2530	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2531	c cpinv=1.0/cpn(il,i)
2532	c if (cvflag_grav) then
2533	c ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
2534	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2535	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2536	c else
2537	c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2538	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2539	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2540	c endif
2541	c endif
2542	c enddo
2543	c enddo
2544
2545	do 480 k=1,i-1
2546	do 1370 il=1,ncum
2547	if (i.le.inb(il)) then
2548	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2549	cpinv=1.0/cpn(il,i)
2550
2551	awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
2552	awat=amax1(awat,0.0)
2553
2554	if (cvflag_grav) then
2555	fr(il,i)=fr(il,i)
2556	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-awat-rr(il,i))
2557	fu(il,i)=fu(il,i)
2558	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2559	fv(il,i)=fv(il,i)
2560	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2561	else ! cvflag_grav
2562	fr(il,i)=fr(il,i)
2563	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
2564	fu(il,i)=fu(il,i)
2565	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2566	fv(il,i)=fv(il,i)
2567	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2568	endif ! cvflag_grav
2569
2570	c (saturated updrafts resulting from mixing) ! cld
2571	qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
2572	nqcond(il,i)=nqcond(il,i)+1. ! cld
2573	endif ! i
2574	1370 continue
2575	480 continue
2576
2577	c do j=1,ntra
2578	c do k=1,i-1
2579	c do il=1,ncum
2580	c if (i.le.inb(il)) then
2581	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2582	c cpinv=1.0/cpn(il,i)
2583	c if (cvflag_grav) then
2584	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2585	c : *(traent(il,k,i,j)-tra(il,i,j))
2586	c else
2587	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2588	c : *(traent(il,k,i,j)-tra(il,i,j))
2589	c endif
2590	c endif
2591	c enddo
2592	c enddo
2593	c enddo
2594
2595	do 490 k=i,nl+1
2596	do 1380 il=1,ncum
2597	if (i.le.inb(il) .and. k.le.inb(il)) then
2598	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2599	cpinv=1.0/cpn(il,i)
2600
2601	if (cvflag_grav) then
2602	fr(il,i)=fr(il,i)
2603	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-rr(il,i))
2604	fu(il,i)=fu(il,i)
2605	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2606	fv(il,i)=fv(il,i)
2607	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2608	else ! cvflag_grav
2609	fr(il,i)=fr(il,i)
2610	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-rr(il,i))
2611	fu(il,i)=fu(il,i)
2612	: +0.1dpinvment(il,k,i)*(uent(il,k,i)-u(il,i))
2613	fv(il,i)=fv(il,i)
2614	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2615	endif ! cvflag_grav
2616	endif ! i and k
2617	1380 continue
2618	490 continue
2619
2620	c do j=1,ntra
2621	c do k=i,nl+1
2622	c do il=1,ncum
2623	c if (i.le.inb(il) .and. k.le.inb(il)) then
2624	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2625	c cpinv=1.0/cpn(il,i)
2626	c if (cvflag_grav) then
2627	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2628	c : *(traent(il,k,i,j)-tra(il,i,j))
2629	c else
2630	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2631	c : *(traent(il,k,i,j)-tra(il,i,j))
2632	c endif
2633	c endif ! i and k
2634	c enddo
2635	c enddo
2636	c enddo
2637
2638	do 1400 il=1,ncum
2639	if (i.le.inb(il)) then
2640	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2641	cpinv=1.0/cpn(il,i)
2642
2643	if (cvflag_grav) then
2644	c sb: on ne fait pas encore la correction permettant de mieux
2645	c conserver l'eau:
2646	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
2647	: +0.01grav(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
2648	: (rp(il,i)-rr(il,i-1)))dpinv
2649
2650	fu(il,i)=fu(il,i)+0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i))
2651	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2652	fv(il,i)=fv(il,i)+0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i))
2653	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2654	else ! cvflag_grav
2655	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
2656	: +0.1(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)
2657	: (rp(il,i)-rr(il,i-1)))dpinv
2658	fu(il,i)=fu(il,i)+0.1(mp(il,i+1)(up(il,i+1)-u(il,i))
2659	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2660	fv(il,i)=fv(il,i)+0.1(mp(il,i+1)(vp(il,i+1)-v(il,i))
2661	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2662	endif ! cvflag_grav
2663
2664	endif ! i
2665	1400 continue
2666
2667	c sb: interface with the cloud parameterization: ! cld
2668
2669	do k=i+1,nl
2670	do il=1,ncum
2671	if (k.le.inb(il) .and. i.le.inb(il)) then ! cld
2672	C (saturated downdrafts resulting from mixing) ! cld
2673	qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
2674	nqcond(il,i)=nqcond(il,i)+1. ! cld
2675	endif ! cld
2676	enddo ! cld
2677	enddo ! cld
2678
2679	C (particular case: no detraining level is found) ! cld
2680	do il=1,ncum ! cld
2681	if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld
2682	qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
2683	nqcond(il,i)=nqcond(il,i)+1. ! cld
2684	endif ! cld
2685	enddo ! cld
2686
2687	do il=1,ncum ! cld
2688	if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld
2689	qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
2690	endif ! cld
2691	enddo
2692
2693	c do j=1,ntra
2694	c do il=1,ncum
2695	c if (i.le.inb(il)) then
2696	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2697	c cpinv=1.0/cpn(il,i)
2698
2699	c if (cvflag_grav) then
2700	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
2701	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2702	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2703	c else
2704	c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2705	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2706	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2707	c endif
2708	c endif ! i
2709	c enddo
2710	c enddo
2711
2712	500 continue
2713
2714
2715	c * move the detrainment at level inb down to level inb-1 *
2716	c * in such a way as to preserve the vertically *
2717	c * integrated enthalpy and water tendencies *
2718	c
2719	do 503 il=1,ncum
2720
2721	ax=0.1ment(il,inb(il),inb(il))(hp(il,inb(il))-h(il,inb(il))
2722	: +t(il,inb(il))*(cpv-cpd)
2723	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
2724	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2725	ft(il,inb(il))=ft(il,inb(il))-ax
2726	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
2727	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
2728	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
2729
2730	bx=0.1ment(il,inb(il),inb(il))(qent(il,inb(il),inb(il))
2731	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2732	fr(il,inb(il))=fr(il,inb(il))-bx
2733	fr(il,inb(il)-1)=fr(il,inb(il)-1)
2734	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
2735	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2736
2737	cx=0.1ment(il,inb(il),inb(il))(uent(il,inb(il),inb(il))
2738	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2739	fu(il,inb(il))=fu(il,inb(il))-cx
2740	fu(il,inb(il)-1)=fu(il,inb(il)-1)
2741	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
2742	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2743
2744	dx=0.1ment(il,inb(il),inb(il))(vent(il,inb(il),inb(il))
2745	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2746	fv(il,inb(il))=fv(il,inb(il))-dx
2747	fv(il,inb(il)-1)=fv(il,inb(il)-1)
2748	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
2749	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2750
2751	503 continue
2752
2753	c do j=1,ntra
2754	c do il=1,ncum
2755	c ex=0.1*ment(il,inb(il),inb(il))
2756	c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2757	c : /(ph(il,inb(il))-ph(il,inb(il)+1))
2758	c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2759	c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2760	c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2761	c : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2762	c enddo
2763	c enddo
2764
2765	c
2766	c * homoginize tendencies below cloud base *
2767	c
2768	c
2769	do il=1,ncum
2770	asum(il)=0.0
2771	bsum(il)=0.0
2772	csum(il)=0.0
2773	dsum(il)=0.0
2774	enddo
2775
2776	do i=1,nl
2777	do il=1,ncum
2778	if (i.le.(icb(il)-1)) then
2779	asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
2780	bsum(il)=bsum(il)+fr(il,i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))
2781	: *(ph(il,i)-ph(il,i+1))
2782	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
2783	: *(ph(il,i)-ph(il,i+1))
2784	dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
2785	endif
2786	enddo
2787	enddo
2788
2789	c!!! do 700 i=1,icb(il)-1
2790	do i=1,nl
2791	do il=1,ncum
2792	if (i.le.(icb(il)-1)) then
2793	ft(il,i)=asum(il)t(il,i)/(th(il,i)dsum(il))
2794	fr(il,i)=bsum(il)/csum(il)
2795	endif
2796	enddo
2797	enddo
2798
2799	c
2800	c * reset counter and return *
2801	c
2802	do il=1,ncum
2803	sig(il,nd)=2.0
2804	enddo
2805
2806
2807	do i=1,nd
2808	do il=1,ncum
2809	upwd(il,i)=0.0
2810	dnwd(il,i)=0.0
2811	enddo
2812	enddo
2813
2814	do i=1,nl
2815	do il=1,ncum
2816	dnwd0(il,i)=-mp(il,i)
2817	enddo
2818	enddo
2819	do i=nl+1,nd
2820	do il=1,ncum
2821	dnwd0(il,i)=0.
2822	enddo
2823	enddo
2824
2825
2826	do i=1,nl
2827	do il=1,ncum
2828	if (i.ge.icb(il) .and. i.le.inb(il)) then
2829	upwd(il,i)=0.0
2830	dnwd(il,i)=0.0
2831	endif
2832	enddo
2833	enddo
2834
2835	do i=1,nl
2836	do k=1,nl
2837	do il=1,ncum
2838	up1(il,k,i)=0.0
2839	dn1(il,k,i)=0.0
2840	enddo
2841	enddo
2842	enddo
2843
2844	do i=1,nl
2845	do k=i,nl
2846	do n=1,i-1
2847	do il=1,ncum
2848	if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2849	up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
2850	dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
2851	endif
2852	enddo
2853	enddo
2854	enddo
2855	enddo
2856
2857	do i=2,nl
2858	do k=i,nl
2859	do il=1,ncum
2860	ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2861	if (i.le.inb(il).and.k.le.inb(il)) then
2862	upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
2863	dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
2864	endif
2865	enddo
2866	enddo
2867	enddo
2868
2869
2870	c!!! DO il=1,ncum
2871	c!!! do i=icb(il),inb(il)
2872	c!!!
2873	c!!! upwd(il,i)=0.0
2874	c!!! dnwd(il,i)=0.0
2875	c!!! do k=i,inb(il)
2876	c!!! up1=0.0
2877	c!!! dn1=0.0
2878	c!!! do n=1,i-1
2879	c!!! up1=up1+ment(il,n,k)
2880	c!!! dn1=dn1-ment(il,k,n)
2881	c!!! enddo
2882	c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
2883	c!!! dnwd(il,i)=dnwd(il,i)+dn1
2884	c!!! enddo
2885	c!!! enddo
2886	c!!!
2887	c!!! ENDDO
2888
2889	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2890	c determination de la variation de flux ascendant entre
2891	c deux niveau non dilue mike
2892	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2893
2894	do i=1,nl
2895	do il=1,ncum
2896	mike(il,i)=m(il,i)
2897	enddo
2898	enddo
2899
2900	do i=nl+1,nd
2901	do il=1,ncum
2902	mike(il,i)=0.
2903	enddo
2904	enddo
2905
2906	do i=1,nd
2907	do il=1,ncum
2908	ma(il,i)=0
2909	enddo
2910	enddo
2911
2912	do i=1,nl
2913	do j=i,nl
2914	do il=1,ncum
2915	ma(il,i)=ma(il,i)+m(il,j)
2916	enddo
2917	enddo
2918	enddo
2919
2920	do i=nl+1,nd
2921	do il=1,ncum
2922	ma(il,i)=0.
2923	enddo
2924	enddo
2925
2926	do i=1,nl
2927	do il=1,ncum
2928	if (i.le.(icb(il)-1)) then
2929	ma(il,i)=0
2930	endif
2931	enddo
2932	enddo
2933
2934	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2935	c icb represente de niveau ou se trouve la
2936	c base du nuage , et inb le top du nuage
2937	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2938
2939	do i=1,nd
2940	do il=1,ncum
2941	mke(il,i)=upwd(il,i)+dnwd(il,i)
2942	enddo
2943	enddo
2944
2945	do i=1,nd
2946	DO 999 il=1,ncum
2947	rdcp=(rrd(1.-rr(il,i))-rr(il,i)rrv)
2948	: /(cpd(1.-rr(il,i))+rr(il,i)cpv)
2949	tls(il,i)=t(il,i)(1000.0/p(il,i))*rdcp
2950	tps(il,i)=tp(il,i)
2951	999 CONTINUE
2952	enddo
2953
2954	c
2955	c * diagnose the in-cloud mixing ratio * ! cld
2956	c * of condensed water * ! cld
2957	c ! cld
2958
2959	do i=1,nd ! cld
2960	do il=1,ncum ! cld
2961	mac(il,i)=0.0 ! cld
2962	wa(il,i)=0.0 ! cld
2963	siga(il,i)=0.0 ! cld
2964	sax(il,i)=0.0 ! cld
2965	enddo ! cld
2966	enddo ! cld
2967
2968	do i=minorig, nl ! cld
2969	do k=i+1,nl+1 ! cld
2970	do il=1,ncum ! cld
2971	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld
2972	mac(il,i)=mac(il,i)+m(il,k) ! cld
2973	endif ! cld
2974	enddo ! cld
2975	enddo ! cld
2976	enddo ! cld
2977
2978	do i=1,nl ! cld
2979	do j=1,i ! cld
2980	do il=1,ncum ! cld
2981	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
2982	: .and. j.ge.icb(il) ) then ! cld
2983	sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
2984	: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
2985	endif ! cld
2986	enddo ! cld
2987	enddo ! cld
2988	enddo ! cld
2989
2990	do i=1,nl ! cld
2991	do il=1,ncum ! cld
2992	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
2993	: .and. sax(il,i).gt.0.0 ) then ! cld
2994	wa(il,i)=sqrt(2.*sax(il,i)) ! cld
2995	endif ! cld
2996	enddo ! cld
2997	enddo ! cld
2998
2999	do i=1,nl ! cld
3000	do il=1,ncum ! cld
3001	if (wa(il,i).gt.0.0) ! cld
3002	: siga(il,i)=mac(il,i)/wa(il,i) ! cld
3003	: rrdtvp(il,i)/p(il,i)/100./delta ! cld
3004	siga(il,i) = min(siga(il,i),1.0) ! cld
3005	cIM cf. FH
3006	if (iflag_clw.eq.0) then
3007	qcondc(il,i)=siga(il,i)clw(il,i)(1.-ep(il,i)) ! cld
3008	: + (1.-siga(il,i))*qcond(il,i) ! cld
3009	else if (iflag_clw.eq.1) then
3010	qcondc(il,i)=qcond(il,i) ! cld
3011	endif
3012
3013	enddo ! cld
3014	enddo ! cld
3015
3016	return
3017	end
3018
3019	SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
3020	& ment,sij,da,phi)
3021	implicit none
3022	c inputs:
3023	integer ncum, nd, na, nloc,len
3024	real ment(nloc,na,na),sij(nloc,na,na)
3025	c ouputs:
3026	real da(nloc,na),phi(nloc,na,na)
3027	c local variables:
3028	integer i,j,k
3029	c
3030	da(:,:)=0.
3031	c
3032	do j=1,na
3033	do k=1,na
3034	do i=1,ncum
3035	da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
3036	phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
3037	c print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j)
3038	end do
3039	end do
3040	end do
3041
3042	return
3043	end
3044
3045
3046	SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
3047	: ,iflag
3048	: ,precip,VPrecip,sig,w0
3049	: ,ft,fq,fu,fv,ftra
3050	: ,inb
3051	: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
3052	: ,da,phi,mp
3053	: ,iflag1
3054	: ,precip1,VPrecip1,sig1,w01
3055	: ,ft1,fq1,fu1,fv1,ftra1
3056	: ,inb1
3057	: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
3058	: ,da1,phi1,mp1)
3059	implicit none
3060
3061	#include "cvparam3.h"
3062
3063	c inputs:
3064	integer len, ncum, nd, ntra, nloc
3065	integer idcum(nloc)
3066	integer iflag(nloc)
3067	integer inb(nloc)
3068	real precip(nloc)
3069	real VPrecip(nloc,nd+1)
3070	real sig(nloc,nd), w0(nloc,nd)
3071	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
3072	real ftra(nloc,nd,ntra)
3073	real Ma(nloc,nd)
3074	real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
3075	real qcondc(nloc,nd)
3076	real wd(nloc),cape(nloc)
3077	real da(nloc,nd),phi(nloc,nd,nd),mp(nloc,nd)
3078
3079	c outputs:
3080	integer iflag1(len)
3081	integer inb1(len)
3082	real precip1(len)
3083	real VPrecip1(len,nd+1)
3084	real sig1(len,nd), w01(len,nd)
3085	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
3086	real ftra1(len,nd,ntra)
3087	real Ma1(len,nd)
3088	real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
3089	real qcondc1(nloc,nd)
3090	real wd1(nloc),cape1(nloc)
3091	real da1(nloc,nd),phi1(nloc,nd,nd),mp1(nloc,nd)
3092
3093	c local variables:
3094	integer i,k,j
3095
3096	do 2000 i=1,ncum
3097	precip1(idcum(i))=precip(i)
3098	iflag1(idcum(i))=iflag(i)
3099	wd1(idcum(i))=wd(i)
3100	inb1(idcum(i))=inb(i)
3101	cape1(idcum(i))=cape(i)
3102	2000 continue
3103
3104	do 2020 k=1,nl
3105	do 2010 i=1,ncum
3106	VPrecip1(idcum(i),k)=VPrecip(i,k)
3107	sig1(idcum(i),k)=sig(i,k)
3108	w01(idcum(i),k)=w0(i,k)
3109	ft1(idcum(i),k)=ft(i,k)
3110	fq1(idcum(i),k)=fq(i,k)
3111	fu1(idcum(i),k)=fu(i,k)
3112	fv1(idcum(i),k)=fv(i,k)
3113	Ma1(idcum(i),k)=Ma(i,k)
3114	upwd1(idcum(i),k)=upwd(i,k)
3115	dnwd1(idcum(i),k)=dnwd(i,k)
3116	dnwd01(idcum(i),k)=dnwd0(i,k)
3117	qcondc1(idcum(i),k)=qcondc(i,k)
3118	da1(idcum(i),k)=da(i,k)
3119	mp1(idcum(i),k)=mp(i,k)
3120	2010 continue
3121	2020 continue
3122
3123	do 2200 i=1,ncum
3124	sig1(idcum(i),nd)=sig(i,nd)
3125	2200 continue
3126
3127
3128	c do 2100 j=1,ntra
3129	c do 2110 k=1,nd ! oct3
3130	c do 2120 i=1,ncum
3131	c ftra1(idcum(i),k,j)=ftra(i,k,j)
3132	c 2120 continue
3133	c 2110 continue
3134	c 2100 continue
3135	do j=1,nd
3136	do k=1,nd
3137	do i=1,ncum
3138	phi1(idcum(i),k,j)=phi(i,k,j)
3139	end do
3140	end do
3141	end do
3142
3143	return
3144	end
3145

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

Download in other formats:

Original Format