Context Navigation

cv3_routines.f90 @ 5695

Last change on this file since 5695 was 5695, checked in by yann meurdesoif, 5 weeks ago
Convection GPU porting : replace saved variables by PARAMETER when possible YM
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: 182.8 KB

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

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

Download in other formats:

Original Format