Context Navigation

cv3_routines.F90 @ 3809

Last change on this file since 3809 was 3809, checked in by ymipsl, 10 years ago
Add LMDZ in aquaplanet configuration YM
File size: 126.0 KB

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

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

Download in other formats:

Original Format