Context Navigation

cv3_routines.f90 @ 5506

Last change on this file since 5506 was 5502, checked in by fhourdin, 35 hours ago
Flag moved to lmdz_cv_ini/cv3_param.
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: 181.4 KB

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

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

Download in other formats:

Original Format