Context Navigation

cv3_routines.f90 @ 5706

Last change on this file since 5706 was 5706, checked in by yann meurdesoif, 6 months ago

Convection GPU porting : suppress potential dependency between columns, may change results (cv3_undilute1)

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

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

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

Download in other formats:

Original Format