Context Navigation

cv30_routines.F @ 3357

Last change on this file since 3357 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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

File size: 93.8 KB

Line
1	!
2	! $Id: cv30_routines.F 1299 2010-01-20 14:27:21Z oboucher $
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)
1834	implicit none
1835
1836
1837	#include "cvthermo.h"
1838	#include "cv30param.h"
1839	#include "cvflag.h"
1840
1841	c inputs:
1842	integer ncum, nd, na, ntra, nloc
1843	integer icb(nloc), inb(nloc)
1844	real delt, plcl(nloc)
1845	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
1846	real u(nloc,nd), v(nloc,nd)
1847	real tra(nloc,nd,ntra)
1848	real p(nloc,nd), ph(nloc,nd+1)
1849	real th(nloc,na), gz(nloc,na)
1850	real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na)
1851	real cpn(nloc,na), tv(nloc,na)
1852	real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
1853
1854	c outputs:
1855	real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
1856	real water(nloc,na), evap(nloc,na), wt(nloc,na)
1857	real trap(nloc,na,ntra)
1858	real b(nloc,na)
1859
1860	c local variables
1861	integer i,j,k,il,num1
1862	real tinv, delti
1863	real awat, afac, afac1, afac2, bfac
1864	real pr1, pr2, sigt, b6, c6, revap, tevap, delth
1865	real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1866	real ampmax
1867	real lvcp(nloc,na)
1868	real wdtrain(nloc)
1869	logical lwork(nloc)
1870
1871
1872	c------------------------------------------------------
1873
1874	delti = 1./delt
1875	tinv=1./3.
1876
1877	mp(:,:)=0.
1878
1879	do i=1,nl
1880	do il=1,ncum
1881	mp(il,i)=0.0
1882	rp(il,i)=rr(il,i)
1883	up(il,i)=u(il,i)
1884	vp(il,i)=v(il,i)
1885	wt(il,i)=0.001
1886	water(il,i)=0.0
1887	evap(il,i)=0.0
1888	b(il,i)=0.0
1889	lvcp(il,i)=lv(il,i)/cpn(il,i)
1890	enddo
1891	enddo
1892
1893	c do k=1,ntra
1894	c do i=1,nd
1895	c do il=1,ncum
1896	c trap(il,i,k)=tra(il,i,k)
1897	c enddo
1898	c enddo
1899	c enddo
1900
1901	c
1902	c * check whether ep(inb)=0, if so, skip precipitating *
1903	c * downdraft calculation *
1904	c
1905
1906	do il=1,ncum
1907	lwork(il)=.TRUE.
1908	if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
1909	enddo
1910
1911	call zilch(wdtrain,ncum)
1912
1913	DO 400 i=nl+1,1,-1
1914
1915	num1=0
1916	do il=1,ncum
1917	if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1
1918	enddo
1919	if (num1.le.0) goto 400
1920
1921	c
1922	c * integrate liquid water equation to find condensed water *
1923	c * and condensed water flux *
1924	c
1925
1926	c
1927	c * begin downdraft loop *
1928	c
1929
1930	c
1931	c * calculate detrained precipitation *
1932	c
1933	do il=1,ncum
1934	if (i.le.inb(il) .and. lwork(il)) then
1935	if (cvflag_grav) then
1936	wdtrain(il)=gravep(il,i)m(il,i)*clw(il,i)
1937	else
1938	wdtrain(il)=10.0ep(il,i)m(il,i)*clw(il,i)
1939	endif
1940	endif
1941	enddo
1942
1943	if(i.gt.1)then
1944	do 320 j=1,i-1
1945	do il=1,ncum
1946	if (i.le.inb(il) .and. lwork(il)) then
1947	awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
1948	awat=amax1(awat,0.0)
1949	if (cvflag_grav) then
1950	wdtrain(il)=wdtrain(il)+gravawatment(il,j,i)
1951	else
1952	wdtrain(il)=wdtrain(il)+10.0awatment(il,j,i)
1953	endif
1954	endif
1955	enddo
1956	320 continue
1957	endif
1958
1959	c
1960	c * find rain water and evaporation using provisional *
1961	c * estimates of rp(i)and rp(i-1) *
1962	c
1963
1964	do 999 il=1,ncum
1965
1966	if (i.le.inb(il) .and. lwork(il)) then
1967
1968	wt(il,i)=45.0
1969
1970	if(i.lt.inb(il))then
1971	rp(il,i)=rp(il,i+1)
1972	: +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
1973	rp(il,i)=0.5*(rp(il,i)+rr(il,i))
1974	endif
1975	rp(il,i)=amax1(rp(il,i),0.0)
1976	rp(il,i)=amin1(rp(il,i),rs(il,i))
1977	rp(il,inb(il))=rr(il,inb(il))
1978
1979	if(i.eq.1)then
1980	afac=p(il,1)(rs(il,1)-rp(il,1))/(1.0e4+2000.0p(il,1)*rs(il,1))
1981	else
1982	rp(il,i-1)=rp(il,i)
1983	: +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i)
1984	rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
1985	rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
1986	rp(il,i-1)=amax1(rp(il,i-1),0.0)
1987	afac1=p(il,i)(rs(il,i)-rp(il,i))/(1.0e4+2000.0p(il,i)*rs(il,i))
1988	afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
1989	: /(1.0e4+2000.0p(il,i-1)rs(il,i-1))
1990	afac=0.5*(afac1+afac2)
1991	endif
1992	if(i.eq.inb(il))afac=0.0
1993	afac=amax1(afac,0.0)
1994	bfac=1./(sigd*wt(il,i))
1995	c
1996	cjyg1
1997	ccc sigt=1.0
1998	ccc if(i.ge.icb)sigt=sigp(i)
1999	c prise en compte de la variation progressive de sigt dans
2000	c les couches icb et icb-1:
2001	c pour plcl<ph(i+1), pr1=0 & pr2=1
2002	c pour plcl>ph(i), pr1=1 & pr2=0
2003	c pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2004	c sur le nuage, et pr2 est la proportion sous la base du
2005	c nuage.
2006	pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2007	pr1=max(0.,min(1.,pr1))
2008	pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2009	pr2=max(0.,min(1.,pr2))
2010	sigt=sigp(il,i)*pr1+pr2
2011	cjyg2
2012	c
2013	b6=bfac50.sigd(ph(il,i)-ph(il,i+1))sigt*afac
2014	c6=water(il,i+1)+bfac*wdtrain(il)
2015	: -50.sigdbfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2016	if(c6.gt.0.0)then
2017	revap=0.5(-b6+sqrt(b6b6+4.*c6))
2018	evap(il,i)=sigtafacrevap
2019	water(il,i)=revap*revap
2020	else
2021	evap(il,i)=-evap(il,i+1)
2022	: +0.02(wdtrain(il)+sigdwt(il,i)*water(il,i+1))
2023	: /(sigd*(ph(il,i)-ph(il,i+1)))
2024	end if
2025	c
2026	c * calculate precipitating downdraft mass flux under *
2027	c * hydrostatic approximation *
2028	c
2029	if (i.ne.1) then
2030
2031	tevap=amax1(0.0,evap(il,i))
2032	delth=amax1(0.001,(th(il,i)-th(il,i-1)))
2033	if (cvflag_grav) then
2034	mp(il,i)=100.ginvlvcp(il,i)sigdtevap
2035	: *(p(il,i-1)-p(il,i))/delth
2036	else
2037	mp(il,i)=10.lvcp(il,i)sigdtevap(p(il,i-1)-p(il,i))/delth
2038	endif
2039	c
2040	c * if hydrostatic assumption fails, *
2041	c * solve cubic difference equation for downdraft theta *
2042	c * and mass flux from two simultaneous differential eqns *
2043	c
2044	amfac=sigdsigd70.0ph(il,i)(p(il,i-1)-p(il,i))
2045	: (th(il,i)-th(il,i-1))/(tv(il,i)th(il,i))
2046	amp2=abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2047	if(amp2.gt.(0.1*amfac))then
2048	xf=100.0sigdsigdsigd(ph(il,i)-ph(il,i+1))
2049	tf=b(il,i)-5.0(th(il,i)-th(il,i-1))t(il,i)
2050	: /(lvcp(il,i)sigdth(il,i))
2051	af=xftf+mp(il,i+1)mp(il,i+1)*tinv
2052	bf=2.(tinvmp(il,i+1))*3+tinvmp(il,i+1)xftf
2053	: +50.(p(il,i-1)-p(il,i))xf*tevap
2054	fac2=1.0
2055	if(bf.lt.0.0)fac2=-1.0
2056	bf=abs(bf)
2057	ur=0.25bfbf-afafaftinvtinv*tinv
2058	if(ur.ge.0.0)then
2059	sru=sqrt(ur)
2060	fac=1.0
2061	if((0.5*bf-sru).lt.0.0)fac=-1.0
2062	mp(il,i)=mp(il,i+1)tinv+(0.5bf+sru)**tinv
2063	: +fac(abs(0.5bf-sru))**tinv
2064	else
2065	d=atan(2.*sqrt(-ur)/(bf+1.0e-28))
2066	if(fac2.lt.0.0)d=3.14159-d
2067	mp(il,i)=mp(il,i+1)tinv+2.sqrt(aftinv)cos(d*tinv)
2068	endif
2069	mp(il,i)=amax1(0.0,mp(il,i))
2070
2071	if (cvflag_grav) then
2072	Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2073	C il faut diviser par (mp(il,i)sigdgrav) et non par (mp(il,i)+sigd*0.1).
2074	C Et il faut bien revoir les facteurs 100.
2075	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2076	2 /(mp(il,i)+sigd*0.1)
2077	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
2078	else
2079	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
2080	2 /(mp(il,i)+sigd*0.1)
2081	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
2082	endif
2083	b(il,i-1)=amax1(b(il,i-1),0.0)
2084	endif
2085	c
2086	c * limit magnitude of mp(i) to meet cfl condition *
2087	c
2088	ampmax=2.0(ph(il,i)-ph(il,i+1))delti
2089	amp2=2.0(ph(il,i-1)-ph(il,i))delti
2090	ampmax=amin1(ampmax,amp2)
2091	mp(il,i)=amin1(mp(il,i),ampmax)
2092	c
2093	c * force mp to decrease linearly to zero *
2094	c * between cloud base and the surface *
2095	c
2096	if(p(il,i).gt.p(il,icb(il)))then
2097	mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2098	endif
2099
2100	360 continue
2101	endif ! i.eq.1
2102	c
2103	c * find mixing ratio of precipitating downdraft *
2104	c
2105
2106	if (i.ne.inb(il)) then
2107
2108	rp(il,i)=rr(il,i)
2109
2110	if(mp(il,i).gt.mp(il,i+1))then
2111
2112	if (cvflag_grav) then
2113	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2114	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
2115	: *(evap(il,i+1)+evap(il,i))
2116	else
2117	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
2118	: +5.sigd(ph(il,i)-ph(il,i+1))
2119	: *(evap(il,i+1)+evap(il,i))
2120	endif
2121	rp(il,i)=rp(il,i)/mp(il,i)
2122	up(il,i)=up(il,i+1)mp(il,i+1)+u(il,i)(mp(il,i)-mp(il,i+1))
2123	up(il,i)=up(il,i)/mp(il,i)
2124	vp(il,i)=vp(il,i+1)mp(il,i+1)+v(il,i)(mp(il,i)-mp(il,i+1))
2125	vp(il,i)=vp(il,i)/mp(il,i)
2126
2127	c do j=1,ntra
2128	c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2129	ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2130	c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2131	c trap(il,i,j)=trap(il,i,j)/mp(il,i)
2132	c end do
2133
2134	else
2135
2136	if(mp(il,i+1).gt.1.0e-16)then
2137	if (cvflag_grav) then
2138	rp(il,i)=rp(il,i+1)
2139	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
2140	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2141	else
2142	rp(il,i)=rp(il,i+1)
2143	: +5.sigd(ph(il,i)-ph(il,i+1))
2144	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2145	endif
2146	up(il,i)=up(il,i+1)
2147	vp(il,i)=vp(il,i+1)
2148
2149	c do j=1,ntra
2150	c trap(il,i,j)=trap(il,i+1,j)
2151	c end do
2152
2153	endif
2154	endif
2155	rp(il,i)=amin1(rp(il,i),rs(il,i))
2156	rp(il,i)=amax1(rp(il,i),0.0)
2157
2158	endif
2159	endif
2160	999 continue
2161
2162	400 continue
2163
2164	return
2165	end
2166
2167	SUBROUTINE cv30_yield(nloc,ncum,nd,na,ntra
2168	: ,icb,inb,delt
2169	: ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th
2170	: ,ep,clw,m,tp,mp,rp,up,vp,trap
2171	: ,wt,water,evap,b
2172	: ,ment,qent,uent,vent,nent,elij,traent,sig
2173	: ,tv,tvp
2174	: ,iflag,precip,VPrecip,ft,fr,fu,fv,ftra
2175	: ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
2176	implicit none
2177
2178	#include "cvthermo.h"
2179	#include "cv30param.h"
2180	#include "cvflag.h"
2181	#include "conema3.h"
2182
2183	c inputs:
2184	integer ncum,nd,na,ntra,nloc
2185	integer icb(nloc), inb(nloc)
2186	real delt
2187	real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
2188	real tra(nloc,nd,ntra), sig(nloc,nd)
2189	real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
2190	real th(nloc,na), p(nloc,nd), tp(nloc,na)
2191	real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
2192	real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
2193	real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
2194	real water(nloc,na), evap(nloc,na), b(nloc,na)
2195	real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
2196	cym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2197	real vent(nloc,na,na), elij(nloc,na,na)
2198	integer nent(nloc,na)
2199	real traent(nloc,na,na,ntra)
2200	real tv(nloc,nd), tvp(nloc,nd)
2201
2202	c input/output:
2203	integer iflag(nloc)
2204
2205	c outputs:
2206	real precip(nloc)
2207	real VPrecip(nloc,nd+1)
2208	real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
2209	real ftra(nloc,nd,ntra)
2210	real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
2211	real dnwd0(nloc,nd), mike(nloc,nd)
2212	real tls(nloc,nd), tps(nloc,nd)
2213	real qcondc(nloc,nd) ! cld
2214	real wd(nloc) ! gust
2215
2216	c local variables:
2217	integer i,k,il,n,j,num1
2218	real rat, awat, delti
2219	real ax, bx, cx, dx, ex
2220	real cpinv, rdcp, dpinv
2221	real lvcp(nloc,na), mke(nloc,na)
2222	real am(nloc), work(nloc), ad(nloc), amp1(nloc)
2223	c!! real up1(nloc), dn1(nloc)
2224	real up1(nloc,nd,nd), dn1(nloc,nd,nd)
2225	real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2226	real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
2227	real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
2228
2229
2230	c-------------------------------------------------------------
2231
2232	c initialization:
2233
2234	delti = 1.0/delt
2235
2236	do il=1,ncum
2237	precip(il)=0.0
2238	wd(il)=0.0 ! gust
2239	VPrecip(il,nd+1)=0.
2240	enddo
2241
2242	do i=1,nd
2243	do il=1,ncum
2244	VPrecip(il,i)=0.0
2245	ft(il,i)=0.0
2246	fr(il,i)=0.0
2247	fu(il,i)=0.0
2248	fv(il,i)=0.0
2249	qcondc(il,i)=0.0 ! cld
2250	qcond(il,i)=0.0 ! cld
2251	nqcond(il,i)=0.0 ! cld
2252	enddo
2253	enddo
2254
2255	c do j=1,ntra
2256	c do i=1,nd
2257	c do il=1,ncum
2258	c ftra(il,i,j)=0.0
2259	c enddo
2260	c enddo
2261	c enddo
2262
2263	do i=1,nl
2264	do il=1,ncum
2265	lvcp(il,i)=lv(il,i)/cpn(il,i)
2266	enddo
2267	enddo
2268
2269
2270	c
2271	c * calculate surface precipitation in mm/day *
2272	c
2273	do il=1,ncum
2274	if(ep(il,inb(il)).ge.0.0001)then
2275	if (cvflag_grav) then
2276	precip(il)=wt(il,1)sigdwater(il,1)86400.1000./(rowl*grav)
2277	else
2278	precip(il)=wt(il,1)sigdwater(il,1)*8640.
2279	endif
2280	endif
2281	enddo
2282
2283	C *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
2284	C
2285	c MAF rajout pour lessivage
2286	do k=1,nl
2287	do il=1,ncum
2288	if (k.le.inb(il)) then
2289	if (cvflag_grav) then
2290	VPrecip(il,k) = wt(il,k)sigdwater(il,k)/grav
2291	else
2292	VPrecip(il,k) = wt(il,k)sigdwater(il,k)/10.
2293	endif
2294	endif
2295	end do
2296	end do
2297	C
2298	c
2299	c * Calculate downdraft velocity scale *
2300	c * NE PAS UTILISER POUR L'INSTANT *
2301	c
2302	c! do il=1,ncum
2303	c! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
2304	c! : /(sigd*p(il,icb(il)))
2305	c! enddo
2306
2307	c
2308	c * calculate tendencies of lowest level potential temperature *
2309	c * and mixing ratio *
2310	c
2311	do il=1,ncum
2312	work(il)=1.0/(ph(il,1)-ph(il,2))
2313	am(il)=0.0
2314	enddo
2315
2316	do k=2,nl
2317	do il=1,ncum
2318	if (k.le.inb(il)) then
2319	am(il)=am(il)+m(il,k)
2320	endif
2321	enddo
2322	enddo
2323
2324	do il=1,ncum
2325
2326	c convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
2327	if (cvflag_grav) then
2328	if((0.01gravwork(il)*am(il)).ge.delti)iflag(il)=1!consist vect
2329	ft(il,1)=0.01gravwork(il)am(il)(t(il,2)-t(il,1)
2330	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2331	else
2332	if((0.1work(il)am(il)).ge.delti)iflag(il)=1 !consistency vect
2333	ft(il,1)=0.1work(il)am(il)*(t(il,2)-t(il,1)
2334	: +(gz(il,2)-gz(il,1))/cpn(il,1))
2335	endif
2336
2337	ft(il,1)=ft(il,1)-0.5lvcp(il,1)sigd*(evap(il,1)+evap(il,2))
2338
2339	if (cvflag_grav) then
2340	ft(il,1)=ft(il,1)-0.009gravsigd*mp(il,2)
2341	: t(il,1)b(il,1)*work(il)
2342	else
2343	ft(il,1)=ft(il,1)-0.09sigdmp(il,2)t(il,1)b(il,1)*work(il)
2344	endif
2345
2346	ft(il,1)=ft(il,1)+0.01sigdwt(il,1)(cl-cpd)water(il,2)*(t(il,2)
2347	:-t(il,1))*work(il)/cpn(il,1)
2348
2349	if (cvflag_grav) then
2350	Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
2351	c (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
2352	fr(il,1)=0.01gravmp(il,2)(rp(il,2)-rr(il,1))work(il)
2353	: +sigd0.5(evap(il,1)+evap(il,2))
2354	c+tard : +sigd*evap(il,1)
2355
2356	fr(il,1)=fr(il,1)+0.01gravam(il)(rr(il,2)-rr(il,1))work(il)
2357
2358	fu(il,1)=fu(il,1)+0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1))
2359	: +am(il)*(u(il,2)-u(il,1)))
2360	fv(il,1)=fv(il,1)+0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1))
2361	: +am(il)*(v(il,2)-v(il,1)))
2362	else ! cvflag_grav
2363	fr(il,1)=0.1mp(il,2)(rp(il,2)-rr(il,1))*work(il)
2364	: +sigd0.5(evap(il,1)+evap(il,2))
2365	fr(il,1)=fr(il,1)+0.1am(il)(rr(il,2)-rr(il,1))*work(il)
2366	fu(il,1)=fu(il,1)+0.1work(il)(mp(il,2)*(up(il,2)-u(il,1))
2367	: +am(il)*(u(il,2)-u(il,1)))
2368	fv(il,1)=fv(il,1)+0.1work(il)(mp(il,2)*(vp(il,2)-v(il,1))
2369	: +am(il)*(v(il,2)-v(il,1)))
2370	endif ! cvflag_grav
2371
2372	enddo ! il
2373
2374	c do j=1,ntra
2375	c do il=1,ncum
2376	c if (cvflag_grav) then
2377	c ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
2378	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2379	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2380	c else
2381	c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2382	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
2383	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2384	c endif
2385	c enddo
2386	c enddo
2387
2388	do j=2,nl
2389	do il=1,ncum
2390	if (j.le.inb(il)) then
2391	if (cvflag_grav) then
2392	fr(il,1)=fr(il,1)
2393	: +0.01gravwork(il)ment(il,j,1)(qent(il,j,1)-rr(il,1))
2394	fu(il,1)=fu(il,1)
2395	: +0.01gravwork(il)ment(il,j,1)(uent(il,j,1)-u(il,1))
2396	fv(il,1)=fv(il,1)
2397	: +0.01gravwork(il)ment(il,j,1)(vent(il,j,1)-v(il,1))
2398	else ! cvflag_grav
2399	fr(il,1)=fr(il,1)
2400	: +0.1work(il)ment(il,j,1)*(qent(il,j,1)-rr(il,1))
2401	fu(il,1)=fu(il,1)
2402	: +0.1work(il)ment(il,j,1)*(uent(il,j,1)-u(il,1))
2403	fv(il,1)=fv(il,1)
2404	: +0.1work(il)ment(il,j,1)*(vent(il,j,1)-v(il,1))
2405	endif ! cvflag_grav
2406	endif ! j
2407	enddo
2408	enddo
2409
2410	c do k=1,ntra
2411	c do j=2,nl
2412	c do il=1,ncum
2413	c if (j.le.inb(il)) then
2414
2415	c if (cvflag_grav) then
2416	c ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
2417	c : *(traent(il,j,1,k)-tra(il,1,k))
2418	c else
2419	c ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
2420	c : *(traent(il,j,1,k)-tra(il,1,k))
2421	c endif
2422
2423	c endif
2424	c enddo
2425	c enddo
2426	c enddo
2427
2428	c
2429	c * calculate tendencies of potential temperature and mixing ratio *
2430	c * at levels above the lowest level *
2431	c
2432	c * first find the net saturated updraft and downdraft mass fluxes *
2433	c * through each level *
2434	c
2435
2436	do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
2437
2438	num1=0
2439	do il=1,ncum
2440	if(i.le.inb(il))num1=num1+1
2441	enddo
2442	if(num1.le.0)go to 500
2443
2444	call zilch(amp1,ncum)
2445	call zilch(ad,ncum)
2446
2447	do 440 k=i+1,nl+1
2448	do 441 il=1,ncum
2449	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
2450	amp1(il)=amp1(il)+m(il,k)
2451	endif
2452	441 continue
2453	440 continue
2454
2455	do 450 k=1,i
2456	do 451 j=i+1,nl+1
2457	do 452 il=1,ncum
2458	if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
2459	amp1(il)=amp1(il)+ment(il,k,j)
2460	endif
2461	452 continue
2462	451 continue
2463	450 continue
2464
2465	do 470 k=1,i-1
2466	do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
2467	do 472 il=1,ncum
2468	if (i.le.inb(il) .and. j.le.inb(il)) then
2469	ad(il)=ad(il)+ment(il,j,k)
2470	endif
2471	472 continue
2472	471 continue
2473	470 continue
2474
2475	do 1350 il=1,ncum
2476	if (i.le.inb(il)) then
2477	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2478	cpinv=1.0/cpn(il,i)
2479
2480	c convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
2481	if (cvflag_grav) then
2482	if((0.01gravdpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
2483	else
2484	if((0.1dpinvamp1(il)).ge.delti)iflag(il)=1 ! vecto
2485	endif
2486
2487	if (cvflag_grav) then
2488	ft(il,i)=0.01gravdpinv(amp1(il)(t(il,i+1)-t(il,i)
2489	: +(gz(il,i+1)-gz(il,i))*cpinv)
2490	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2491	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
2492	rat=cpn(il,i-1)*cpinv
2493	ft(il,i)=ft(il,i)-0.009gravsigd(mp(il,i+1)t(il,i)*b(il,i)
2494	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2495	ft(il,i)=ft(il,i)+0.01gravdpinvment(il,i,i)(hp(il,i)-h(il,i)
2496	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2497	else ! cvflag_grav
2498	ft(il,i)=0.1dpinv(amp1(il)*(t(il,i+1)-t(il,i)
2499	: +(gz(il,i+1)-gz(il,i))*cpinv)
2500	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
2501	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
2502	rat=cpn(il,i-1)*cpinv
2503	ft(il,i)=ft(il,i)-0.09sigd(mp(il,i+1)t(il,i)b(il,i)
2504	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
2505	ft(il,i)=ft(il,i)+0.1dpinvment(il,i,i)*(hp(il,i)-h(il,i)
2506	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
2507	endif ! cvflag_grav
2508
2509
2510	ft(il,i)=ft(il,i)+0.01sigdwt(il,i)(cl-cpd)water(il,i+1)
2511	: (t(il,i+1)-t(il,i))dpinv*cpinv
2512
2513	if (cvflag_grav) then
2514	fr(il,i)=0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i))
2515	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2516	fu(il,i)=fu(il,i)+0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i))
2517	: -ad(il)*(u(il,i)-u(il,i-1)))
2518	fv(il,i)=fv(il,i)+0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i))
2519	: -ad(il)*(v(il,i)-v(il,i-1)))
2520	else ! cvflag_grav
2521	fr(il,i)=0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il,i))
2522	: -ad(il)*(rr(il,i)-rr(il,i-1)))
2523	fu(il,i)=fu(il,i)+0.1dpinv(amp1(il)*(u(il,i+1)-u(il,i))
2524	: -ad(il)*(u(il,i)-u(il,i-1)))
2525	fv(il,i)=fv(il,i)+0.1dpinv(amp1(il)*(v(il,i+1)-v(il,i))
2526	: -ad(il)*(v(il,i)-v(il,i-1)))
2527	endif ! cvflag_grav
2528
2529	endif ! i
2530	1350 continue
2531
2532	c do k=1,ntra
2533	c do il=1,ncum
2534	c if (i.le.inb(il)) then
2535	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2536	c cpinv=1.0/cpn(il,i)
2537	c if (cvflag_grav) then
2538	c ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
2539	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2540	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2541	c else
2542	c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2543	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
2544	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2545	c endif
2546	c endif
2547	c enddo
2548	c enddo
2549
2550	do 480 k=1,i-1
2551	do 1370 il=1,ncum
2552	if (i.le.inb(il)) then
2553	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2554	cpinv=1.0/cpn(il,i)
2555
2556	awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
2557	awat=amax1(awat,0.0)
2558
2559	if (cvflag_grav) then
2560	fr(il,i)=fr(il,i)
2561	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-awat-rr(il,i))
2562	fu(il,i)=fu(il,i)
2563	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2564	fv(il,i)=fv(il,i)
2565	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2566	else ! cvflag_grav
2567	fr(il,i)=fr(il,i)
2568	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
2569	fu(il,i)=fu(il,i)
2570	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2571	fv(il,i)=fv(il,i)
2572	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2573	endif ! cvflag_grav
2574
2575	c (saturated updrafts resulting from mixing) ! cld
2576	qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
2577	nqcond(il,i)=nqcond(il,i)+1. ! cld
2578	endif ! i
2579	1370 continue
2580	480 continue
2581
2582	c do j=1,ntra
2583	c do k=1,i-1
2584	c do il=1,ncum
2585	c if (i.le.inb(il)) then
2586	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2587	c cpinv=1.0/cpn(il,i)
2588	c if (cvflag_grav) then
2589	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2590	c : *(traent(il,k,i,j)-tra(il,i,j))
2591	c else
2592	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2593	c : *(traent(il,k,i,j)-tra(il,i,j))
2594	c endif
2595	c endif
2596	c enddo
2597	c enddo
2598	c enddo
2599
2600	do 490 k=i,nl+1
2601	do 1380 il=1,ncum
2602	if (i.le.inb(il) .and. k.le.inb(il)) then
2603	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2604	cpinv=1.0/cpn(il,i)
2605
2606	if (cvflag_grav) then
2607	fr(il,i)=fr(il,i)
2608	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-rr(il,i))
2609	fu(il,i)=fu(il,i)
2610	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
2611	fv(il,i)=fv(il,i)
2612	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
2613	else ! cvflag_grav
2614	fr(il,i)=fr(il,i)
2615	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-rr(il,i))
2616	fu(il,i)=fu(il,i)
2617	: +0.1dpinvment(il,k,i)*(uent(il,k,i)-u(il,i))
2618	fv(il,i)=fv(il,i)
2619	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
2620	endif ! cvflag_grav
2621	endif ! i and k
2622	1380 continue
2623	490 continue
2624
2625	c do j=1,ntra
2626	c do k=i,nl+1
2627	c do il=1,ncum
2628	c if (i.le.inb(il) .and. k.le.inb(il)) then
2629	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2630	c cpinv=1.0/cpn(il,i)
2631	c if (cvflag_grav) then
2632	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
2633	c : *(traent(il,k,i,j)-tra(il,i,j))
2634	c else
2635	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
2636	c : *(traent(il,k,i,j)-tra(il,i,j))
2637	c endif
2638	c endif ! i and k
2639	c enddo
2640	c enddo
2641	c enddo
2642
2643	do 1400 il=1,ncum
2644	if (i.le.inb(il)) then
2645	dpinv=1.0/(ph(il,i)-ph(il,i+1))
2646	cpinv=1.0/cpn(il,i)
2647
2648	if (cvflag_grav) then
2649	c sb: on ne fait pas encore la correction permettant de mieux
2650	c conserver l'eau:
2651	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
2652	: +0.01grav(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
2653	: (rp(il,i)-rr(il,i-1)))dpinv
2654
2655	fu(il,i)=fu(il,i)+0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i))
2656	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2657	fv(il,i)=fv(il,i)+0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i))
2658	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2659	else ! cvflag_grav
2660	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
2661	: +0.1(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)
2662	: (rp(il,i)-rr(il,i-1)))dpinv
2663	fu(il,i)=fu(il,i)+0.1(mp(il,i+1)(up(il,i+1)-u(il,i))
2664	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
2665	fv(il,i)=fv(il,i)+0.1(mp(il,i+1)(vp(il,i+1)-v(il,i))
2666	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
2667	endif ! cvflag_grav
2668
2669	endif ! i
2670	1400 continue
2671
2672	c sb: interface with the cloud parameterization: ! cld
2673
2674	do k=i+1,nl
2675	do il=1,ncum
2676	if (k.le.inb(il) .and. i.le.inb(il)) then ! cld
2677	C (saturated downdrafts resulting from mixing) ! cld
2678	qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
2679	nqcond(il,i)=nqcond(il,i)+1. ! cld
2680	endif ! cld
2681	enddo ! cld
2682	enddo ! cld
2683
2684	C (particular case: no detraining level is found) ! cld
2685	do il=1,ncum ! cld
2686	if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld
2687	qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
2688	nqcond(il,i)=nqcond(il,i)+1. ! cld
2689	endif ! cld
2690	enddo ! cld
2691
2692	do il=1,ncum ! cld
2693	if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld
2694	qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
2695	endif ! cld
2696	enddo
2697
2698	c do j=1,ntra
2699	c do il=1,ncum
2700	c if (i.le.inb(il)) then
2701	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2702	c cpinv=1.0/cpn(il,i)
2703
2704	c if (cvflag_grav) then
2705	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
2706	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2707	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2708	c else
2709	c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2710	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
2711	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2712	c endif
2713	c endif ! i
2714	c enddo
2715	c enddo
2716
2717	500 continue
2718
2719
2720	c * move the detrainment at level inb down to level inb-1 *
2721	c * in such a way as to preserve the vertically *
2722	c * integrated enthalpy and water tendencies *
2723	c
2724	do 503 il=1,ncum
2725
2726	ax=0.1ment(il,inb(il),inb(il))(hp(il,inb(il))-h(il,inb(il))
2727	: +t(il,inb(il))*(cpv-cpd)
2728	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
2729	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2730	ft(il,inb(il))=ft(il,inb(il))-ax
2731	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
2732	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
2733	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
2734
2735	bx=0.1ment(il,inb(il),inb(il))(qent(il,inb(il),inb(il))
2736	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2737	fr(il,inb(il))=fr(il,inb(il))-bx
2738	fr(il,inb(il)-1)=fr(il,inb(il)-1)
2739	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
2740	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2741
2742	cx=0.1ment(il,inb(il),inb(il))(uent(il,inb(il),inb(il))
2743	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2744	fu(il,inb(il))=fu(il,inb(il))-cx
2745	fu(il,inb(il)-1)=fu(il,inb(il)-1)
2746	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
2747	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2748
2749	dx=0.1ment(il,inb(il),inb(il))(vent(il,inb(il),inb(il))
2750	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2751	fv(il,inb(il))=fv(il,inb(il))-dx
2752	fv(il,inb(il)-1)=fv(il,inb(il)-1)
2753	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
2754	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
2755
2756	503 continue
2757
2758	c do j=1,ntra
2759	c do il=1,ncum
2760	c ex=0.1*ment(il,inb(il),inb(il))
2761	c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2762	c : /(ph(il,inb(il))-ph(il,inb(il)+1))
2763	c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2764	c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2765	c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2766	c : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2767	c enddo
2768	c enddo
2769
2770	c
2771	c * homoginize tendencies below cloud base *
2772	c
2773	c
2774	do il=1,ncum
2775	asum(il)=0.0
2776	bsum(il)=0.0
2777	csum(il)=0.0
2778	dsum(il)=0.0
2779	enddo
2780
2781	do i=1,nl
2782	do il=1,ncum
2783	if (i.le.(icb(il)-1)) then
2784	asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
2785	bsum(il)=bsum(il)+fr(il,i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))
2786	: *(ph(il,i)-ph(il,i+1))
2787	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
2788	: *(ph(il,i)-ph(il,i+1))
2789	dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
2790	endif
2791	enddo
2792	enddo
2793
2794	c!!! do 700 i=1,icb(il)-1
2795	do i=1,nl
2796	do il=1,ncum
2797	if (i.le.(icb(il)-1)) then
2798	ft(il,i)=asum(il)t(il,i)/(th(il,i)dsum(il))
2799	fr(il,i)=bsum(il)/csum(il)
2800	endif
2801	enddo
2802	enddo
2803
2804	c
2805	c * reset counter and return *
2806	c
2807	do il=1,ncum
2808	sig(il,nd)=2.0
2809	enddo
2810
2811
2812	do i=1,nd
2813	do il=1,ncum
2814	upwd(il,i)=0.0
2815	dnwd(il,i)=0.0
2816	enddo
2817	enddo
2818
2819	do i=1,nl
2820	do il=1,ncum
2821	dnwd0(il,i)=-mp(il,i)
2822	enddo
2823	enddo
2824	do i=nl+1,nd
2825	do il=1,ncum
2826	dnwd0(il,i)=0.
2827	enddo
2828	enddo
2829
2830
2831	do i=1,nl
2832	do il=1,ncum
2833	if (i.ge.icb(il) .and. i.le.inb(il)) then
2834	upwd(il,i)=0.0
2835	dnwd(il,i)=0.0
2836	endif
2837	enddo
2838	enddo
2839
2840	do i=1,nl
2841	do k=1,nl
2842	do il=1,ncum
2843	up1(il,k,i)=0.0
2844	dn1(il,k,i)=0.0
2845	enddo
2846	enddo
2847	enddo
2848
2849	do i=1,nl
2850	do k=i,nl
2851	do n=1,i-1
2852	do il=1,ncum
2853	if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2854	up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
2855	dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
2856	endif
2857	enddo
2858	enddo
2859	enddo
2860	enddo
2861
2862	do i=2,nl
2863	do k=i,nl
2864	do il=1,ncum
2865	ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2866	if (i.le.inb(il).and.k.le.inb(il)) then
2867	upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
2868	dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
2869	endif
2870	enddo
2871	enddo
2872	enddo
2873
2874
2875	c!!! DO il=1,ncum
2876	c!!! do i=icb(il),inb(il)
2877	c!!!
2878	c!!! upwd(il,i)=0.0
2879	c!!! dnwd(il,i)=0.0
2880	c!!! do k=i,inb(il)
2881	c!!! up1=0.0
2882	c!!! dn1=0.0
2883	c!!! do n=1,i-1
2884	c!!! up1=up1+ment(il,n,k)
2885	c!!! dn1=dn1-ment(il,k,n)
2886	c!!! enddo
2887	c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
2888	c!!! dnwd(il,i)=dnwd(il,i)+dn1
2889	c!!! enddo
2890	c!!! enddo
2891	c!!!
2892	c!!! ENDDO
2893
2894	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2895	c determination de la variation de flux ascendant entre
2896	c deux niveau non dilue mike
2897	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2898
2899	do i=1,nl
2900	do il=1,ncum
2901	mike(il,i)=m(il,i)
2902	enddo
2903	enddo
2904
2905	do i=nl+1,nd
2906	do il=1,ncum
2907	mike(il,i)=0.
2908	enddo
2909	enddo
2910
2911	do i=1,nd
2912	do il=1,ncum
2913	ma(il,i)=0
2914	enddo
2915	enddo
2916
2917	do i=1,nl
2918	do j=i,nl
2919	do il=1,ncum
2920	ma(il,i)=ma(il,i)+m(il,j)
2921	enddo
2922	enddo
2923	enddo
2924
2925	do i=nl+1,nd
2926	do il=1,ncum
2927	ma(il,i)=0.
2928	enddo
2929	enddo
2930
2931	do i=1,nl
2932	do il=1,ncum
2933	if (i.le.(icb(il)-1)) then
2934	ma(il,i)=0
2935	endif
2936	enddo
2937	enddo
2938
2939	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2940	c icb represente de niveau ou se trouve la
2941	c base du nuage , et inb le top du nuage
2942	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2943
2944	do i=1,nd
2945	do il=1,ncum
2946	mke(il,i)=upwd(il,i)+dnwd(il,i)
2947	enddo
2948	enddo
2949
2950	do i=1,nd
2951	DO 999 il=1,ncum
2952	rdcp=(rrd(1.-rr(il,i))-rr(il,i)rrv)
2953	: /(cpd(1.-rr(il,i))+rr(il,i)cpv)
2954	tls(il,i)=t(il,i)(1000.0/p(il,i))*rdcp
2955	tps(il,i)=tp(il,i)
2956	999 CONTINUE
2957	enddo
2958
2959	c
2960	c * diagnose the in-cloud mixing ratio * ! cld
2961	c * of condensed water * ! cld
2962	c ! cld
2963
2964	do i=1,nd ! cld
2965	do il=1,ncum ! cld
2966	mac(il,i)=0.0 ! cld
2967	wa(il,i)=0.0 ! cld
2968	siga(il,i)=0.0 ! cld
2969	sax(il,i)=0.0 ! cld
2970	enddo ! cld
2971	enddo ! cld
2972
2973	do i=minorig, nl ! cld
2974	do k=i+1,nl+1 ! cld
2975	do il=1,ncum ! cld
2976	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld
2977	mac(il,i)=mac(il,i)+m(il,k) ! cld
2978	endif ! cld
2979	enddo ! cld
2980	enddo ! cld
2981	enddo ! cld
2982
2983	do i=1,nl ! cld
2984	do j=1,i ! cld
2985	do il=1,ncum ! cld
2986	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
2987	: .and. j.ge.icb(il) ) then ! cld
2988	sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
2989	: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
2990	endif ! cld
2991	enddo ! cld
2992	enddo ! cld
2993	enddo ! cld
2994
2995	do i=1,nl ! cld
2996	do il=1,ncum ! cld
2997	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
2998	: .and. sax(il,i).gt.0.0 ) then ! cld
2999	wa(il,i)=sqrt(2.*sax(il,i)) ! cld
3000	endif ! cld
3001	enddo ! cld
3002	enddo ! cld
3003
3004	do i=1,nl ! cld
3005	do il=1,ncum ! cld
3006	if (wa(il,i).gt.0.0) ! cld
3007	: siga(il,i)=mac(il,i)/wa(il,i) ! cld
3008	: rrdtvp(il,i)/p(il,i)/100./delta ! cld
3009	siga(il,i) = min(siga(il,i),1.0) ! cld
3010	cIM cf. FH
3011	if (iflag_clw.eq.0) then
3012	qcondc(il,i)=siga(il,i)clw(il,i)(1.-ep(il,i)) ! cld
3013	: + (1.-siga(il,i))*qcond(il,i) ! cld
3014	else if (iflag_clw.eq.1) then
3015	qcondc(il,i)=qcond(il,i) ! cld
3016	endif
3017
3018	enddo ! cld
3019	enddo ! cld
3020
3021	return
3022	end
3023
3024	SUBROUTINE cv30_tracer(nloc,len,ncum,nd,na,
3025	& ment,sij,da,phi)
3026	implicit none
3027	c inputs:
3028	integer ncum, nd, na, nloc,len
3029	real ment(nloc,na,na),sij(nloc,na,na)
3030	c ouputs:
3031	real da(nloc,na),phi(nloc,na,na)
3032	c local variables:
3033	integer i,j,k
3034	c
3035	da(:,:)=0.
3036	c
3037	do j=1,na
3038	do k=1,na
3039	do i=1,ncum
3040	da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
3041	phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
3042	c print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j)
3043	end do
3044	end do
3045	end do
3046
3047	return
3048	end
3049
3050
3051	SUBROUTINE cv30_uncompress(nloc,len,ncum,nd,ntra,idcum
3052	: ,iflag
3053	: ,precip,VPrecip,sig,w0
3054	: ,ft,fq,fu,fv,ftra
3055	: ,inb
3056	: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
3057	: ,da,phi,mp
3058	: ,iflag1
3059	: ,precip1,VPrecip1,sig1,w01
3060	: ,ft1,fq1,fu1,fv1,ftra1
3061	: ,inb1
3062	: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
3063	: ,da1,phi1,mp1)
3064	implicit none
3065
3066	#include "cv30param.h"
3067
3068	c inputs:
3069	integer len, ncum, nd, ntra, nloc
3070	integer idcum(nloc)
3071	integer iflag(nloc)
3072	integer inb(nloc)
3073	real precip(nloc)
3074	real VPrecip(nloc,nd+1)
3075	real sig(nloc,nd), w0(nloc,nd)
3076	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
3077	real ftra(nloc,nd,ntra)
3078	real Ma(nloc,nd)
3079	real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
3080	real qcondc(nloc,nd)
3081	real wd(nloc),cape(nloc)
3082	real da(nloc,nd),phi(nloc,nd,nd),mp(nloc,nd)
3083
3084	c outputs:
3085	integer iflag1(len)
3086	integer inb1(len)
3087	real precip1(len)
3088	real VPrecip1(len,nd+1)
3089	real sig1(len,nd), w01(len,nd)
3090	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
3091	real ftra1(len,nd,ntra)
3092	real Ma1(len,nd)
3093	real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
3094	real qcondc1(nloc,nd)
3095	real wd1(nloc),cape1(nloc)
3096	real da1(nloc,nd),phi1(nloc,nd,nd),mp1(nloc,nd)
3097
3098	c local variables:
3099	integer i,k,j
3100
3101	do 2000 i=1,ncum
3102	precip1(idcum(i))=precip(i)
3103	iflag1(idcum(i))=iflag(i)
3104	wd1(idcum(i))=wd(i)
3105	inb1(idcum(i))=inb(i)
3106	cape1(idcum(i))=cape(i)
3107	2000 continue
3108
3109	do 2020 k=1,nl
3110	do 2010 i=1,ncum
3111	VPrecip1(idcum(i),k)=VPrecip(i,k)
3112	sig1(idcum(i),k)=sig(i,k)
3113	w01(idcum(i),k)=w0(i,k)
3114	ft1(idcum(i),k)=ft(i,k)
3115	fq1(idcum(i),k)=fq(i,k)
3116	fu1(idcum(i),k)=fu(i,k)
3117	fv1(idcum(i),k)=fv(i,k)
3118	Ma1(idcum(i),k)=Ma(i,k)
3119	upwd1(idcum(i),k)=upwd(i,k)
3120	dnwd1(idcum(i),k)=dnwd(i,k)
3121	dnwd01(idcum(i),k)=dnwd0(i,k)
3122	qcondc1(idcum(i),k)=qcondc(i,k)
3123	da1(idcum(i),k)=da(i,k)
3124	mp1(idcum(i),k)=mp(i,k)
3125	2010 continue
3126	2020 continue
3127
3128	do 2200 i=1,ncum
3129	sig1(idcum(i),nd)=sig(i,nd)
3130	2200 continue
3131
3132
3133	c do 2100 j=1,ntra
3134	c do 2110 k=1,nd ! oct3
3135	c do 2120 i=1,ncum
3136	c ftra1(idcum(i),k,j)=ftra(i,k,j)
3137	c 2120 continue
3138	c 2110 continue
3139	c 2100 continue
3140	do j=1,nd
3141	do k=1,nd
3142	do i=1,ncum
3143	phi1(idcum(i),k,j)=phi(i,k,j)
3144	end do
3145	end do
3146	end do
3147
3148	return
3149	end
3150

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

Download in other formats:

Original Format