Context Navigation

cv3_routines.f90 @ 5622

Last change on this file since 5622 was 5622, checked in by evignon, 7 weeks ago
commission pour la conservation de l'eau dans la convection profonde. 1/ integre les corrections apportees par JYG dans cv3routines. 2/ par defaut, conserve dx*masse quand on appelle la convection tous les n pas de temps Cette commission conduit à une perte de convergence
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.9 KB

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

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

Download in other formats:

Original Format