Context Navigation

cv3_routines.f90 @ 5623

Last change on this file since 5623 was 5618, checked in by aborella, 3 months ago
Merge with trunk testing r5597. We have convergence in prod and debug in NPv7.0.1c
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: 185.0 KB

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

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

Download in other formats:

Original Format