Context Navigation

cv3_routines.F90 @ 2562

Last change on this file since 2562 was 2520, checked in by lguez, 8 years ago
Bug fix: icb1 was not made >= 2.
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: 147.4 KB

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

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

Download in other formats:

Original Format