Context Navigation

cv3_routines.F90 @ 2489

Last change on this file since 2489 was 2481, checked in by fhourdin, 8 years ago
Introduction d'une dependance epmax=f(Cape) sur proposition de Camille Risi
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: 140.6 KB

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

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

Download in other formats:

Original Format