Context Navigation

cv3_routines.f90 @ 5555

Last change on this file since 5555 was 5544, checked in by jyg, 7 days ago
Bug fix in cv3_routines.f90 (in cv3_unsat and cv3_yield) concerning water conservation. And small improvement in the coding of ice fraction in convective precipitation.
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.0 KB

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

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

Download in other formats:

Original Format