Context Navigation

cv3_routines.f90 @ 5447

Last change on this file since 5447 was 5346, checked in by fhourdin, 4 weeks ago

Debut de replaysation de la convection profonde.

Regroupement de cvparam, cv3param et cvthermo (récemment
passés de statut de .h à module, dans un unique module
lmdz_cv_ini.f90

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: 180.6 KB

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

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

Download in other formats:

Original Format