Context Navigation

cv3_routines.F @ 4359

Last change on this file since 4359 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: 105.7 KB

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

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

Download in other formats:

Original Format