Context Navigation

cv3_routines.F90 @ 1939

Last change on this file since 1939 was 265, checked in by lfita, 10 years ago
Adding control error values on convection, since with this it has ran!!!
File size: 136.1 KB

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

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

Download in other formats:

Original Format