Context Navigation

cv3_routines.F90 @ 2257

Last change on this file since 2257 was 2257, checked in by jyg, 9 years ago
Improvement of three loops in cv3_undilute2 in cv3_routines.F90
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 128.8 KB

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

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

Download in other formats:

Original Format