Context Navigation

cv3_routines.f90 @ 5747

Last change on this file since 5747 was 5707, checked in by yann meurdesoif, 3 weeks ago

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

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.6 KB

Line
1
2	! $Id: cv3_routines.f90 5707 2025-06-16 14:20:42Z dcugnet $
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	!ym icbmax = 2
2137	!ym DO i = 1, ncum
2138	!ym icbmax = max(icbmax, icb(i))
2139	!ym END DO
2140
2141	! update sig and w0 below cloud base:
2142
2143	!ym DO k = 1, icbmax
2144	DO k = 1, nd
2145	DO i = 1, ncum
2146	IF (k<=MAX(2,icb(i))) THEN
2147	IF (k<=icb(i)) THEN
2148	sig(i, k) = beta*sig(i, k) - &
2149	2.alphabuoy(i, icb(i))*buoy(i, icb(i))
2150	sig(i, k) = max(sig(i,k), 0.0)
2151	w0(i, k) = beta*w0(i, k)
2152	END IF
2153	ENDIF
2154	END DO
2155	END DO
2156
2157	!! if(inb.lt.(nl-1))then
2158	!! do 85 i=inb+1,nl-1
2159	!! sig(i)=betasig(i)+2.alphabuoy(inb)
2160	!! 1 abs(buoy(inb))
2161	!! sig(i)=max(sig(i),0.0)
2162	!! w0(i)=beta*w0(i)
2163	!! 85 continue
2164	!! end if
2165
2166	!! do 87 i=1,icb
2167	!! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
2168	!! sig(i)=max(sig(i),0.0)
2169	!! w0(i)=beta*w0(i)
2170	!! 87 continue
2171
2172	! -------------------------------------------------------------
2173	! -- Reset fractional areas of updrafts and w0 at initial time
2174	! -- and after 10 time steps of no convection
2175	! -------------------------------------------------------------
2176
2177	DO k = 1, nl - 1
2178	DO i = 1, ncum
2179	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
2180	sig(i, k) = 0.0
2181	w0(i, k) = 0.0
2182	END IF
2183	END DO
2184	END DO
2185
2186	! -------------------------------------------------------------
2187	! -- Calculate convective available potential energy (cape),
2188	! -- vertical velocity (w), fractional area covered by
2189	! -- undilute updraft (sig), and updraft mass flux (m)
2190	! -------------------------------------------------------------
2191
2192	DO i = 1, ncum
2193	cape(i) = 0.0
2194	END DO
2195
2196	! compute dtmin (minimum buoyancy between ICB and given level k):
2197
2198	DO i = 1, ncum
2199	DO k = 1, nl
2200	dtmin(i, k) = 100.0
2201	END DO
2202	END DO
2203
2204	DO i = 1, ncum
2205	DO k = 1, nl
2206	DO j = minorig, nl
2207	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
2208	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
2209	END IF
2210	END DO
2211	END DO
2212	END DO
2213
2214	! the interval on which cape is computed starts at pbase :
2215
2216	DO k = 1, nl
2217	DO i = 1, ncum
2218
2219	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
2220
2221	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
2222	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
2223	cape(i) = amax1(0.0, cape(i))
2224	sigold(i, k) = sig(i, k)
2225
2226	! dtmin(i,k)=100.0
2227	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
2228	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
2229	! 97 continue
2230
2231	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
2232	sig(i, k) = max(sig(i,k), 0.0)
2233	sig(i, k) = amin1(sig(i,k), 0.01)
2234	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
2235	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
2236	amu = 0.5(sig(i,k)+sigold(i,k))w
2237	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
2238	w0(i, k) = w
2239	END IF
2240
2241	END DO
2242	END DO
2243
2244	DO i = 1, ncum
2245	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
2246	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))
2247	sig(i, icb(i)) = sig(i, icb(i)+1)
2248	sig(i, icb(i)-1) = sig(i, icb(i))
2249	END DO
2250
2251	! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
2252	! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
2253	! ccc the final mass flux (cbmflast) is greater than the target mass flux
2254	! ccc (cbmf) ??).
2255	! cc
2256	! c do i = 1,ncum
2257	! c cbmflast(i) = 0.
2258	! c enddo
2259	! cc
2260	! c do k= 1,nl
2261	! c do i = 1,ncum
2262	! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
2263	! c cbmflast(i) = cbmflast(i)+M(i,k)
2264	! c ENDIF
2265	! c enddo
2266	! c enddo
2267	! cc
2268	! c do i = 1,ncum
2269	! c IF (cbmflast(i) .lt. 1.e-6) THEN
2270	! c iflag(i) = 3
2271	! c ENDIF
2272	! c enddo
2273	! cc
2274	! c do k= 1,nl
2275	! c do i = 1,ncum
2276	! c IF (iflag(i) .ge. 3) THEN
2277	! c M(i,k) = 0.
2278	! c sig(i,k) = 0.
2279	! c w0(i,k) = 0.
2280	! c ENDIF
2281	! c enddo
2282	! c enddo
2283	! cc
2284	!! cape=0.0
2285	!! do 98 i=icb+1,inb
2286	!! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
2287	!! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
2288	!! dcape=rrdbuoy(i-1)deltap/p(i-1)
2289	!! dlnp=deltap/p(i-1)
2290	!! cape=max(0.0,cape)
2291	!! sigold=sig(i)
2292
2293	!! dtmin=100.0
2294	!! do 97 j=icb,i-1
2295	!! dtmin=amin1(dtmin,buoy(j))
2296	!! 97 continue
2297
2298	!! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
2299	!! sig(i)=max(sig(i),0.0)
2300	!! sig(i)=amin1(sig(i),0.01)
2301	!! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
2302	!! w=(1.-beta)facsqrt(cape)+beta*w0(i)
2303	!! amu=0.5(sig(i)+sigold)w
2304	!! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
2305	!! w0(i)=w
2306	!! 98 continue
2307	!! w0(icb)=0.5*w0(icb+1)
2308	!! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
2309	!! sig(icb)=sig(icb+1)
2310	!! sig(icb-1)=sig(icb)
2311
2312	RETURN
2313	END SUBROUTINE cv3_closure
2314
2315	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
2316	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
2317	unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
2318	ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
2319	USE cvflag_mod_h
2320	USE lmdz_cv_ini, ONLY : cpd,cpv,minorig,nl,rrv,cpd,ginv,grav
2321	IMPLICIT NONE
2322
2323	! ---------------------------------------------------------------------
2324	! a faire:
2325	! - vectorisation de la partie normalisation des flux (do 789...)
2326	! ---------------------------------------------------------------------
2327
2328	!inputs:
2329	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2330	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
2331	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
2332	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
2333	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2334	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
2335	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
2336	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
2337	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, h, hp
2338	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf, frac
2339	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp, ep, clw
2340	REAL, DIMENSION (nloc, na), INTENT (IN) :: m ! input of convect3
2341
2342	!outputs:
2343	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: ment, qent
2344	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
2345	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: sij, elij
2346	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
2347	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: ments, qents
2348	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
2349
2350	!local variables:
2351	INTEGER i, j, k, il, im, jm
2352	INTEGER num1, num2
2353	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
2354	REAL alt, smid, sjmin, sjmax, delp, delm
2355	REAL asij(nloc), smax(nloc), scrit(nloc)
2356	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
2357	REAL sigij(nloc, nd, nd)
2358	REAL wgh
2359	REAL zm(nloc, na)
2360	LOGICAL lwork(nloc)
2361
2362	! =====================================================================
2363	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
2364	! =====================================================================
2365
2366	! ori do 360 i=1,ncum*nlp
2367	DO j = 1, nl
2368	DO il = 1, ncum
2369	nent(il, j) = 0
2370	! in convect3, m is computed in cv3_closure
2371	! ori m(i,1)=0.0
2372	END DO
2373	END DO
2374
2375	! ori do 400 k=1,nlp
2376	! ori do 390 j=1,nlp
2377	DO j = 1, nl
2378	DO k = 1, nl
2379	DO il = 1, ncum
2380	qent(il, k, j) = rr(il, j)
2381	uent(il, k, j) = u(il, j)
2382	vent(il, k, j) = v(il, j)
2383	elij(il, k, j) = 0.0
2384	!ym ment(i,k,j)=0.0
2385	!ym sij(i,k,j)=0.0
2386	END DO
2387	END DO
2388	END DO
2389
2390	!ym
2391	ment(1:ncum, 1:nd, 1:nd) = 0.0
2392	sij(1:ncum, 1:nd, 1:nd) = 0.0
2393
2394	!AC! do k=1,ntra
2395	!AC! do j=1,nd ! instead nlp
2396	!AC! do i=1,nd ! instead nlp
2397	!AC! do il=1,ncum
2398	!AC! traent(il,i,j,k)=tra(il,j,k)
2399	!AC! enddo
2400	!AC! enddo
2401	!AC! enddo
2402	!AC! enddo
2403	zm(:, :) = 0.
2404
2405	! =====================================================================
2406	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
2407	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
2408	! --- FRACTION (sij)
2409	! =====================================================================
2410
2411	DO i = minorig + 1, nl
2412
2413	DO j = minorig, nl
2414	DO il = 1, ncum
2415	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
2416
2417	rti = qnk(il) - ep(il, i)*clw(il, i)
2418	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
2419
2420
2421	IF (cvflag_ice) THEN
2422	! print*,cvflag_ice,'cvflag_ice dans do 700'
2423	IF (t(il,j)<=263.15) THEN
2424	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
2425	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
2426	END IF
2427	END IF
2428
2429	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
2430	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
2431	dei = denom
2432	IF (abs(dei)<0.01) dei = 0.01
2433	sij(il, i, j) = anum/dei
2434	sij(il, i, i) = 1.0
2435	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
2436	altem = altem/bf2
2437	cwat = clw(il, j)*(1.-ep(il,j))
2438	stemp = sij(il, i, j)
2439	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
2440
2441	IF (cvflag_ice) THEN
2442	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
2443	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
2444	ELSE
2445	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
2446	denom = denom + lv(il, j)*(rr(il,i)-rti)
2447	END IF
2448
2449	IF (abs(denom)<0.01) denom = 0.01
2450	sij(il, i, j) = anum/denom
2451	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
2452	altem = altem - (bf2-1.)*cwat
2453	END IF
2454	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
2455	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
2456	uent(il, i, j) = sij(il, i, j)u(il, i) + (1.-sij(il,i,j))unk(il)
2457	vent(il, i, j) = sij(il, i, j)v(il, i) + (1.-sij(il,i,j))vnk(il)
2458	!!!! do k=1,ntra
2459	!!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2460	!!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2461	!!!! end do
2462	elij(il, i, j) = altem
2463	elij(il, i, j) = max(0.0, elij(il,i,j))
2464	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
2465	nent(il, i) = nent(il, i) + 1
2466	END IF
2467	sij(il, i, j) = max(0.0, sij(il,i,j))
2468	sij(il, i, j) = amin1(1.0, sij(il,i,j))
2469	END IF ! new
2470	END DO
2471	END DO
2472
2473	!AC! do k=1,ntra
2474	!AC! do j=minorig,nl
2475	!AC! do il=1,ncum
2476	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
2477	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
2478	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2479	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2480	!AC! endif
2481	!AC! enddo
2482	!AC! enddo
2483	!AC! enddo
2484
2485
2486	! * if no air can entrain at level i assume that updraft detrains *
2487	! * at that level and calculate detrained air flux and properties *
2488
2489
2490	! @ do 170 i=icb(il),inb(il)
2491
2492	DO il = 1, ncum
2493	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
2494	! @ if(nent(il,i).eq.0)then
2495	ment(il, i, i) = m(il, i)
2496	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2497	uent(il, i, i) = unk(il)
2498	vent(il, i, i) = vnk(il)
2499	elij(il, i, i) = clw(il, i)
2500	! MAF sij(il,i,i)=1.0
2501	sij(il, i, i) = 0.0
2502	END IF
2503	END DO
2504	END DO
2505
2506	!AC! do j=1,ntra
2507	!AC! do i=minorig+1,nl
2508	!AC! do il=1,ncum
2509	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
2510	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2511	!AC! endif
2512	!AC! enddo
2513	!AC! enddo
2514	!AC! enddo
2515
2516	DO j = minorig, nl
2517	DO i = minorig, nl
2518	DO il = 1, ncum
2519	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
2520	sigij(il, i, j) = sij(il, i, j)
2521	END IF
2522	END DO
2523	END DO
2524	END DO
2525	! @ enddo
2526
2527	! @170 continue
2528
2529	! =====================================================================
2530	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
2531	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
2532	! =====================================================================
2533	asum(1:nloc,1:nd) = 0.
2534	csum(1:nloc,1:nd) = 0.
2535
2536	DO il = 1, ncum
2537	lwork(il) = .FALSE.
2538	END DO
2539
2540	DO i = minorig + 1, nl
2541
2542	num1 = 0
2543	DO il = 1, ncum
2544	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
2545	END DO
2546	!ym IF (num1<=0) GO TO 789
2547	IF (num1<=0) CYCLE
2548
2549	DO il = 1, ncum
2550	IF (i>=icb(il) .AND. i<=inb(il)) THEN
2551	lwork(il) = (nent(il,i)/=0)
2552	qp = qnk(il) - ep(il, i)*clw(il, i)
2553
2554	IF (cvflag_ice) THEN
2555
2556	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
2557	(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
2558	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
2559	(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
2560	ELSE
2561
2562	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
2563	(cpv-cpd)t(il, i)(qp-rr(il,i))
2564	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
2565	(cpd-cpv)t(il, i)(rr(il,i)-qp)
2566	END IF
2567
2568	IF (abs(denom)<0.01) denom = 0.01
2569	scrit(il) = anum/denom
2570	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
2571	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
2572	smax(il) = 0.0
2573	asij(il) = 0.0
2574	END IF
2575	END DO
2576
2577	DO j = nl, minorig, -1
2578
2579	num2 = 0
2580	DO il = 1, ncum
2581	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2582	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2583	lwork(il)) num2 = num2 + 1
2584	END DO
2585	!ym IF (num2<=0) GO TO 175
2586	IF (num2<=0) CYCLE
2587
2588	DO il = 1, ncum
2589	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2590	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2591	lwork(il)) THEN
2592
2593	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
2594	wgh = 1.0
2595	IF (j>i) THEN
2596	sjmax = max(sij(il,i,j+1), smax(il))
2597	sjmax = amin1(sjmax, scrit(il))
2598	smax(il) = max(sij(il,i,j), smax(il))
2599	sjmin = max(sij(il,i,j-1), smax(il))
2600	sjmin = amin1(sjmin, scrit(il))
2601	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
2602	smid = amin1(sij(il,i,j), scrit(il))
2603	ELSE
2604	sjmax = max(sij(il,i,j+1), scrit(il))
2605	smid = max(sij(il,i,j), scrit(il))
2606	sjmin = 0.0
2607	IF (j>1) sjmin = sij(il, i, j-1)
2608	sjmin = max(sjmin, scrit(il))
2609	END IF
2610	delp = abs(sjmax-smid)
2611	delm = abs(sjmin-smid)
2612	asij(il) = asij(il) + wgh*(delp+delm)
2613	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
2614	END IF
2615	END IF
2616	END DO
2617
2618	175 END DO
2619
2620	DO il = 1, ncum
2621	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2622	asij(il) = max(1.0E-16, asij(il))
2623	asij(il) = 1.0/asij(il)
2624	asum(il, i) = 0.0
2625	bsum(il, i) = 0.0
2626	csum(il, i) = 0.0
2627	END IF
2628	END DO
2629
2630	DO j = minorig, nl
2631	DO il = 1, ncum
2632	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2633	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2634	ment(il, i, j) = ment(il, i, j)*asij(il)
2635	END IF
2636	END DO
2637	END DO
2638
2639	DO j = minorig, nl
2640	DO il = 1, ncum
2641	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2642	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2643	asum(il, i) = asum(il, i) + ment(il, i, j)
2644	ment(il, i, j) = ment(il, i, j)*sig(il, j)
2645	bsum(il, i) = bsum(il, i) + ment(il, i, j)
2646	END IF
2647	END DO
2648	END DO
2649
2650	DO il = 1, ncum
2651	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2652	bsum(il, i) = max(bsum(il,i), 1.0E-16)
2653	bsum(il, i) = 1.0/bsum(il, i)
2654	END IF
2655	END DO
2656
2657	DO j = minorig, nl
2658	DO il = 1, ncum
2659	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2660	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2661	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
2662	END IF
2663	END DO
2664	END DO
2665
2666	DO j = minorig, nl
2667	DO il = 1, ncum
2668	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2669	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2670	csum(il, i) = csum(il, i) + ment(il, i, j)
2671	END IF
2672	END DO
2673	END DO
2674
2675	DO il = 1, ncum
2676	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2677	csum(il,i)<m(il,i)) THEN
2678	nent(il, i) = 0
2679	ment(il, i, i) = m(il, i)
2680	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2681	uent(il, i, i) = unk(il)
2682	vent(il, i, i) = vnk(il)
2683	elij(il, i, i) = clw(il, i)
2684	! MAF sij(il,i,i)=1.0
2685	sij(il, i, i) = 0.0
2686	END IF
2687	END DO ! il
2688
2689	!AC! do j=1,ntra
2690	!AC! do il=1,ncum
2691	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
2692	!AC! : .and. csum(il,i).lt.m(il,i) ) then
2693	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2694	!AC! endif
2695	!AC! enddo
2696	!AC! enddo
2697	789 END DO
2698
2699	! MAF: renormalisation de MENT
2700	zm(1:nloc,1:na) = 0.
2701
2702	DO jm = 1, nl
2703	DO im = 1, nl
2704	DO il = 1, ncum
2705	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
2706	END DO
2707	END DO
2708	END DO
2709
2710	DO jm = 1, nl
2711	DO im = 1, nl
2712	DO il = 1, ncum
2713	IF (zm(il,im)/=0.) THEN
2714	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
2715	END IF
2716	END DO
2717	END DO
2718	END DO
2719
2720	DO jm = 1, nl
2721	DO im = 1, nl
2722	DO il = 1, ncum
2723	qents(il, im, jm) = qent(il, im, jm)
2724	ments(il, im, jm) = ment(il, im, jm)
2725	END DO
2726	END DO
2727	END DO
2728
2729	RETURN
2730	END SUBROUTINE cv3_mixing
2731
2732	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
2733	t, rr, rs, gz, u, v, tra, p, ph, &
2734	th, tv, lv, lf, cpn, ep, sigp, clw, frac_s, qpreca, frac_a, qta , & !!jygprl
2735	m, ment, elij, delt, plcl, coef_clos, &
2736	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
2737	faci, b, sigd, &
2738	wdtrainA, wdtrainS, wdtrainM) ! RomP
2739	USE lmdz_cv_ini, ONLY : cpd,ginv,grav,nl,nlp,sigdz
2740	USE cvflag_mod_h
2741	USE print_control_mod, ONLY: prt_level, lunout
2742	USE nuage_params_mod_h
2743	IMPLICIT NONE
2744
2745	!inputs:
2746	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2747	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
2748	REAL, INTENT(IN) :: delt
2749	REAL, DIMENSION (nloc), INTENT (IN) :: plcl
2750	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
2751	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz
2752	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
2753	REAL, DIMENSION (nloc, nd, ntra), INTENT(IN) :: tra
2754	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
2755	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2756	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, sigp, clw !adiab ascent shedding
2757	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_s !ice fraction in adiab ascent shedding !!jygprl
2758	REAL, DIMENSION (nloc, na), INTENT (IN) :: qpreca !adiab ascent precip !!jygprl
2759	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_a !ice fraction in adiab ascent precip !!jygprl
2760	REAL, DIMENSION (nloc, na), INTENT (IN) :: qta !adiab ascent specific total water !!jygprl
2761	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tv, lv, cpn
2762	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
2763	REAL, DIMENSION (nloc, na), INTENT (IN) :: m
2764	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: ment, elij
2765	REAL, DIMENSION (nloc), INTENT (IN) :: coef_clos
2766
2767	!input/output
2768	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag(nloc)
2769
2770	!outputs:
2771	REAL, DIMENSION (nloc, na), INTENT (OUT) :: mp, rp, up, vp
2772	REAL, DIMENSION (nloc, na), INTENT (OUT) :: water, evap, wt
2773	REAL, DIMENSION (nloc, na), INTENT (OUT) :: ice, fondue
2774	REAL, DIMENSION (nloc, na), INTENT (OUT) :: faci ! ice fraction in precipitation
2775	REAL, DIMENSION (nloc, na, ntra), INTENT (OUT) :: trap
2776	REAL, DIMENSION (nloc, na), INTENT (OUT) :: b
2777	REAL, DIMENSION (nloc), INTENT (OUT) :: sigd
2778	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
2779	! de l ascendance adiabatique et des flux melanges Pa et Pm.
2780	! Distinction des wdtrain
2781	! Pa = wdtrainA Pm = wdtrainM
2782	REAL, DIMENSION (nloc, na), INTENT (OUT) :: wdtrainA, wdtrainS, wdtrainM
2783
2784	!local variables
2785	INTEGER i, j, k, il, num1, ndp1
2786	REAL smallestreal
2787	REAL tinv, delti, coef
2788	REAL awat, afac, afac1, afac2, bfac
2789	REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
2790	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
2791	REAL ampmax, thaw
2792	REAL tevap(nloc)
2793	REAL, DIMENSION (nloc, na) :: lvcp, lfcp
2794	REAL, DIMENSION (nloc, na) :: h, hm
2795	REAL, DIMENSION (nloc, na) :: ma
2796	REAL, DIMENSION (nloc, na) :: frac ! ice fraction in precipitation source
2797	REAL, DIMENSION (nloc, na) :: fraci ! provisionnal ice fraction in precipitation
2798	REAL, DIMENSION (nloc, na) :: prec
2799	REAL wdtrain(nloc)
2800	LOGICAL lwork(nloc), mplus(nloc)
2801
2802
2803	! ------------------------------------------------------
2804	IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1)
2805
2806	smallestreal=tiny(smallestreal)
2807
2808	! =============================
2809	! --- INITIALIZE OUTPUT ARRAYS
2810	! =============================
2811	! (loops up to nl+1)
2812	mp(:,:) = 0.
2813	rp(:,:) = 0.
2814	up(:,:) = 0.
2815	vp(:,:) = 0.
2816	water(:,:) = 0.
2817	evap(:,:) = 0.
2818	wt(:,:) = 0.
2819	ice(:,:) = 0.
2820	fondue(:,:) = 0.
2821	faci(:,:) = 0.
2822	b(:,:) = 0.
2823	sigd(:) = 0.
2824	!! RomP >>>
2825	wdtrainA(:,:) = 0.
2826	wdtrainS(:,:) = 0.
2827	wdtrainM(:,:) = 0.
2828	!! RomP <<<
2829
2830	DO i = 1, nlp
2831	DO il = 1, ncum
2832	rp(il, i) = rr(il, i)
2833	up(il, i) = u(il, i)
2834	vp(il, i) = v(il, i)
2835	wt(il, i) = 0.001
2836	END DO
2837	END DO
2838
2839	! *** Set the fractionnal area sigd of precipitating downdraughts
2840	DO il = 1, ncum
2841	sigd(il) = sigdz*coef_clos(il)
2842	END DO
2843
2844	! =====================================================================
2845	! --- INITIALIZE VARIOUS ARRAYS AND PARAMETERS USED IN THE COMPUTATIONS
2846	! =====================================================================
2847	! (loops up to nl+1)
2848
2849	delti = 1./delt
2850	tinv = 1./3.
2851
2852	DO i = 1, nlp
2853	DO il = 1, ncum
2854	frac(il, i) = 0.0
2855	fraci(il, i) = 0.0
2856	prec(il, i) = 0.0
2857	lvcp(il, i) = lv(il, i)/cpn(il, i)
2858	lfcp(il, i) = lf(il, i)/cpn(il, i)
2859	END DO
2860	END DO
2861
2862	!AC! do k=1,ntra
2863	!AC! do i=1,nd
2864	!AC! do il=1,ncum
2865	!AC! trap(il,i,k)=tra(il,i,k)
2866	!AC! enddo
2867	!AC! enddo
2868	!AC! enddo
2869
2870	! * check whether ep(inb)=0, if so, skip precipitating *
2871	! * downdraft calculation *
2872
2873
2874	DO il = 1, ncum
2875	!! lwork(il)=.TRUE.
2876	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
2877	!jyg<
2878	!! lwork(il) = ep(il, inb(il)) >= 0.0001
2879	lwork(il) = ep(il, inb(il)) >= 0.0001 .AND. iflag(il) <= 2
2880	END DO
2881
2882	!
2883	! Get adiabatic ascent mass flux
2884	!
2885	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2886	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
2887	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2888	!!! Warning : this option leads to water conservation violation
2889	!!! Expert only
2890	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2891	DO il = 1, ncum
2892	ma(il, nlp) = 0.
2893	ma(il, 1) = 0.
2894	END DO
2895
2896	DO i = nl, 2, -1
2897	DO il = 1, ncum
2898	ma(il, i) = ma(il, i+1)*(1.-qta(il,i))/(1.-qta(il,i-1)) + m(il, i)
2899	END DO
2900	END DO
2901	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2902	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
2903	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2904	DO il = 1, ncum
2905	ma(il, nlp) = 0.
2906	ma(il, 1) = 0.
2907	END DO
2908
2909	DO i = nl, 2, -1
2910	DO il = 1, ncum
2911	ma(il, i) = ma(il, i+1) + m(il, i)
2912	END DO
2913	END DO
2914
2915	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
2916	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2917
2918	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2919	!
2920	! * begin downdraft loop *
2921	!
2922	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2923
2924	DO i = nl + 1, 1, -1
2925
2926	num1 = 0
2927	DO il = 1, ncum
2928	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
2929	END DO
2930	!ym IF (num1<=0) GO TO 400
2931	IF (num1<=0) CYCLE
2932
2933	wdtrain(1:ncum) = 0.0
2934
2935
2936	! * integrate liquid water equation to find condensed water *
2937	! * and condensed water flux *
2938	!
2939	!
2940	! * calculate detrained precipitation *
2941
2942
2943	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2944	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2945	DO il = 1, ncum
2946	IF (i<=inb(il) .AND. lwork(il)) THEN
2947	wdtrainS(il, i) = ep(il, i)m(il, i)clw(il, i) ! jyg
2948	END IF
2949	END DO
2950
2951	IF (i>1) THEN
2952	DO j = 1, i - 1
2953	DO il = 1, ncum
2954	IF (i<=inb(il) .AND. lwork(il)) THEN
2955	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
2956	awat = max(awat, 0.0)
2957	wdtrainM(il, i) = wdtrainM(il, i) + awat*ment(il, j, i) ! jyg
2958	END IF
2959	END DO
2960	END DO
2961	END IF
2962
2963	IF (cvflag_prec_eject) THEN
2964	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2965	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
2966	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2967	!!! Warning : this option leads to water conservation violation
2968	!!! Expert only
2969	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2970	IF ( i > 1) THEN
2971	DO il = 1, ncum
2972	IF (i<=inb(il) .AND. lwork(il)) THEN
2973	wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))/(1. - qta(il, i-1)) ! Pa jygprl
2974	END IF
2975	END DO
2976	ENDIF ! ( i > 1)
2977	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2978	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
2979	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2980	IF ( i > 1) THEN
2981	DO il = 1, ncum
2982	IF (i<=inb(il) .AND. lwork(il)) THEN
2983	wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i)) ! Pa jygprl
2984	END IF
2985	END DO
2986	ENDIF ! ( i > 1)
2987
2988	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
2989	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2990	ENDIF ! (cvflag_prec_eject)
2991
2992	IF ( i > 1) THEN
2993	DO il = 1, ncum
2994	IF (i<=inb(il) .AND. lwork(il)) THEN
2995	wdtrain(il) = grav*(wdtrainS(il,i) + wdtrainM(il,i) + wdtrainA(il,i))
2996	END IF
2997	END DO
2998	ENDIF ! ( i > 1)
2999
3000	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3001	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3002
3003	! * find rain water and evaporation using provisional *
3004	! * estimates of rp(i)and rp(i-1) *
3005
3006
3007	IF (cvflag_ice) THEN !!jygprl
3008	IF (cvflag_prec_eject) THEN
3009	DO il = 1, ncum !!jygprl
3010	IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl
3011	frac(il, i) = (frac_a(il,i)wdtrainA(il,i)+frac_s(il,i)(wdtrainS(il,i)+wdtrainM(il,i))) / & !!jygprl
3012	max(wdtrainA(il,i)+wdtrainS(il,i)+wdtrainM(il,i),smallestreal) !!jygprl
3013	fraci(il, i) = frac(il, i) !!jygprl
3014	END IF !!jygprl
3015	END DO !!jygprl
3016	ELSE ! (cvflag_prec_eject)
3017	DO il = 1, ncum !!jygprl
3018	IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl
3019	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3020	IF (keepbug_ice_frac) THEN
3021	frac(il, i) = frac_s(il, i)
3022	! Ice fraction computed again here as a function of the temperature seen by unsaturated downdraughts
3023	! (i.e. the cold pool temperature) for compatibility with earlier versions.
3024	fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
3025	fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
3026	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3027	ELSE ! (keepbug_ice_frac)
3028	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3029	frac(il, i) = frac_s(il, i)
3030	fraci(il, i) = frac(il, i) !!jygprl
3031	ENDIF ! (keepbug_ice_frac)
3032	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3033	END IF !!jygprl
3034	END DO !!jygprl
3035	ENDIF ! (cvflag_prec_eject)
3036	END IF !!jygprl
3037
3038
3039	DO il = 1, ncum
3040	IF (i<=inb(il) .AND. lwork(il)) THEN
3041
3042	wt(il, i) = 45.0
3043
3044	IF (i<inb(il)) THEN
3045	rp(il, i) = rp(il, i+1) + &
3046	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
3047	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
3048	END IF
3049	rp(il, i) = max(rp(il,i), 0.0)
3050	rp(il, i) = amin1(rp(il,i), rs(il,i))
3051	rp(il, inb(il)) = rr(il, inb(il))
3052
3053	IF (i==1) THEN
3054	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
3055	IF (cvflag_ice) THEN
3056	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
3057	END IF
3058	ELSE
3059	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)
3060	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
3061	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
3062	rp(il, i-1) = max(rp(il,i-1), 0.0)
3063	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
3064	afac2 = p(il, i-1)(rs(il,i-1)-rp(il,i-1))/(1.0E4+2000.0p(il,i-1)*rs(il,i-1))
3065	afac = 0.5*(afac1+afac2)
3066	END IF
3067	IF (i==inb(il)) afac = 0.0
3068	afac = max(afac, 0.0)
3069	bfac = 1./(sigd(il)*wt(il,i))
3070
3071	!
3072	IF (prt_level >= 20) THEN
3073	Print*, 'cv3_unsat after provisional rp estimate: rp, afac, bfac ', &
3074	i, rp(1, i), afac,bfac
3075	ENDIF
3076	!
3077	!JYG1
3078	! cc sigt=1.0
3079	! cc if(i.ge.icb)sigt=sigp(i)
3080	! prise en compte de la variation progressive de sigt dans
3081	! les couches icb et icb-1:
3082	! pour plcl<ph(i+1), pr1=0 & pr2=1
3083	! pour plcl>ph(i), pr1=1 & pr2=0
3084	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
3085	! sur le nuage, et pr2 est la proportion sous la base du
3086	! nuage.
3087	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
3088	pr1 = max(0., min(1.,pr1))
3089	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
3090	pr2 = max(0., min(1.,pr2))
3091	sigt = sigp(il, i)*pr1 + pr2
3092	!JYG2
3093
3094	!JYG----
3095	! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3096	! c6 = water(il,i+1) + wdtrain(il)*bfac
3097	! c6 = prec(il,i+1) + wdtrain(il)*bfac
3098	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
3099	! evap(il,i)=sigtafacrevap
3100	! water(il,i)=revap*revap
3101	! prec(il,i)=revap*revap
3102	!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
3103	!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
3104	!!---end jyg---
3105
3106	! --------retour � la formulation originale d'Emanuel.
3107	IF (cvflag_ice) THEN
3108
3109	! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3110	! c6=prec(il,i+1)+bfac*wdtrain(il) &
3111	! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
3112	! if(c6.gt.0.0)then
3113	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
3114
3115	!JAM Attention: evap=sigt*E
3116	! Modification: evap devient l'�vaporation en milieu de couche
3117	! car n�cessaire dans cv3_yield
3118	! Du coup, il faut modifier pas mal d'�quations...
3119	! et l'expression de afac qui devient afac1
3120	! revap=sqrt((prec(i+1)+prec(i))/2)
3121
3122	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
3123	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
3124	! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
3125	! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
3126	! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
3127	IF (c6>b6*b6+1.E-20) THEN
3128	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
3129	ELSE
3130	revap = (-b6+sqrt(b6b6+4.c6))/2.
3131	END IF
3132	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
3133	! print*,prec(il,i),'neige'
3134
3135	!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
3136	! c evap(il,i)=sigtafacrevap
3137	! ce qui n'est pas correct. Dans cv_routines, la formulation a �t� modifiee.
3138	! Ici,l'evaporation evap est simplement calculee par l'equation de
3139	! conservation.
3140	! prec(il,i)=revap*revap
3141	! else
3142	!JYG---- Correction : si c6 <= 0, water(il,i)=0.
3143	! prec(il,i)=0.
3144	! endif
3145
3146	!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
3147	! moins [tt ce qui sort de la couche i]
3148	! print *, 'evap avec ice'
3149	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
3150	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
3151	!
3152	IF (prt_level >= 20) THEN
3153	Print*, 'cv3_unsat after evap computation: wdtrain, sigd, wt, prec(i+1),prec(i) ', &
3154	i, wdtrain(1), sigd(1), wt(1,i), prec(1,i+1),prec(1,i)
3155	ENDIF
3156	!
3157
3158	!jyg<
3159	d6 = prec(il,i)-prec(il,i+1)
3160
3161	!! d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
3162	!! e6 = bfac*wdtrain(il)
3163	!! f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
3164	!>jyg
3165	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
3166	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
3167	thaw = min(max(thaw,0.0), 1.0)
3168	!jyg<
3169	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
3170	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
3171	water(il, i) = min(prec(il,i), max(water(il,i), 0.))
3172	ice(il, i) = min(prec(il,i), max(ice(il,i), 0.))
3173
3174	!! water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
3175	!! water(il, i) = max(water(il,i), 0.)
3176	!! ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
3177	!! ice(il, i) = max(ice(il,i), 0.)
3178	!>jyg
3179	fondue(il, i) = ice(il, i)*thaw
3180	water(il, i) = water(il, i) + fondue(il, i)
3181	ice(il, i) = ice(il, i) - fondue(il, i)
3182
3183	!! IF (water(il,i)+ice(il,i)<1.E-30) THEN
3184	!! faci(il, i) = 0.
3185	!! ELSE
3186	!! faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
3187	!! END IF
3188
3189	faci(il,i) = ice(il, i)/max((water(il,i)+ice(il,i)), smallestreal)
3190
3191	! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
3192	! water(il,i)=max(water(il,i),0.)
3193	! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
3194	! ice(il,i)=max(ice(il,i),0.)
3195	! fondue(il,i)=ice(il,i)*thaw
3196	! water(il,i)=water(il,i)+fondue(il,i)
3197	! ice(il,i)=ice(il,i)-fondue(il,i)
3198
3199	! if((water(il,i)+ice(il,i)).lt.1.e-30)then
3200	! faci(il,i)=0.
3201	! else
3202	! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
3203	! endif
3204
3205	ELSE
3206	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3207	c6 = water(il, i+1) + bfac*wdtrain(il) - &
3208	50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
3209	IF (c6>0.0) THEN
3210	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
3211	water(il, i) = revap*revap
3212	ELSE
3213	water(il, i) = 0.
3214	END IF
3215	! print *, 'evap sans ice'
3216	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
3217	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
3218
3219	END IF
3220	END IF !(i.le.inb(il) .and. lwork(il))
3221	END DO
3222	! ----------------------------------------------------------------
3223
3224	! cc
3225	! * calculate precipitating downdraft mass flux under *
3226	! * hydrostatic approximation *
3227
3228	DO il = 1, ncum
3229	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
3230
3231	tevap(il) = max(0.0, evap(il,i))
3232	delth = max(0.001, (th(il,i)-th(il,i-1)))
3233	IF (cvflag_ice) THEN
3234	IF (cvflag_grav) THEN
3235	mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
3236	(p(il,i-1)-p(il,i))/delth + &
3237	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
3238	(p(il,i-1)-p(il,i))/delth + &
3239	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
3240	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
3241	ELSE
3242	mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
3243	(p(il,i-1)-p(il,i))/delth + &
3244	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
3245	(p(il,i-1)-p(il,i))/delth + &
3246	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
3247	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
3248
3249	END IF
3250	ELSE
3251	IF (cvflag_grav) THEN
3252	mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
3253	(p(il,i-1)-p(il,i))/delth
3254	ELSE
3255	mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
3256	(p(il,i-1)-p(il,i))/delth
3257	END IF
3258
3259	END IF
3260
3261	END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
3262	IF (prt_level .GE. 20) THEN
3263	PRINT *,'cv3_unsat, mp hydrostatic ', i, mp(il,i)
3264	ENDIF
3265	END DO
3266	! ----------------------------------------------------------------
3267
3268	! * if hydrostatic assumption fails, *
3269	! * solve cubic difference equation for downdraft theta *
3270	! * and mass flux from two simultaneous differential eqns *
3271
3272	DO il = 1, ncum
3273	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
3274
3275	amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
3276	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
3277	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
3278
3279	IF (amp2>(0.1*amfac)) THEN
3280	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
3281	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
3282	(lvcp(il,i)sigd(il)th(il,i))
3283	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
3284
3285	IF (cvflag_ice) THEN
3286	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
3287	50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
3288	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
3289	ELSE
3290
3291	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
3292	50.(p(il,i-1)-p(il,i))xf*tevap(il)
3293	END IF
3294
3295	fac2 = 1.0
3296	IF (bf<0.0) fac2 = -1.0
3297	bf = abs(bf)
3298	ur = 0.25bfbf - afafaftinvtinv*tinv
3299	IF (ur>=0.0) THEN
3300	sru = sqrt(ur)
3301	fac = 1.0
3302	IF ((0.5*bf-sru)<0.0) fac = -1.0
3303	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
3304	fac(abs(0.5bf-sru))**tinv
3305	ELSE
3306	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
3307	IF (fac2<0.0) d = 3.14159 - d
3308	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
3309	END IF
3310	mp(il, i) = max(0.0, mp(il,i))
3311	IF (prt_level .GE. 20) THEN
3312	PRINT *,'cv3_unsat, mp cubic ', i, mp(il,i)
3313	ENDIF
3314
3315	IF (cvflag_ice) THEN
3316	IF (cvflag_grav) THEN
3317	!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
3318	! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
3319	! Et il faut bien revoir les facteurs 100.
3320	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i)) &
3321	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
3322	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
3323	(ph(il,i)-ph(il,i+1))) / &
3324	(mp(il,i)+sigd(il)*0.1) - &
3325	10.0(th(il,i)-th(il,i-1))t(il, i) / &
3326	(lvcp(il,i)sigd(il)th(il,i))
3327	ELSE
3328	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
3329	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
3330	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
3331	(ph(il,i)-ph(il,i+1))) / &
3332	(mp(il,i)+sigd(il)*0.1) - &
3333	10.0(th(il,i)-th(il,i-1))t(il, i) / &
3334	(lvcp(il,i)sigd(il)th(il,i))
3335	END IF
3336	ELSE
3337	IF (cvflag_grav) THEN
3338	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
3339	(mp(il,i)+sigd(il)*0.1) - &
3340	10.0(th(il,i)-th(il,i-1))t(il, i) / &
3341	(lvcp(il,i)sigd(il)th(il,i))
3342	ELSE
3343	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
3344	(mp(il,i)+sigd(il)*0.1) - &
3345	10.0(th(il,i)-th(il,i-1))t(il, i) / &
3346	(lvcp(il,i)sigd(il)th(il,i))
3347	END IF
3348	END IF
3349	b(il, i-1) = max(b(il,i-1), 0.0)
3350
3351	END IF !(amp2.gt.(0.1*amfac))
3352
3353	!jyg< This part shifted 10 lines farther
3354	!!! * limit magnitude of mp(i) to meet cfl condition *
3355	!!
3356	!! ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
3357	!! amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
3358	!! ampmax = min(ampmax, amp2)
3359	!! mp(il, i) = min(mp(il,i), ampmax)
3360	!>jyg
3361
3362	! * force mp to decrease linearly to zero *
3363	! * between cloud base and the surface *
3364
3365
3366	! c if(p(il,i).gt.p(il,icb(il)))then
3367	! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
3368	! c endif
3369	IF (ph(il,i)>0.9*plcl(il)) THEN
3370	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
3371	END IF
3372
3373	!jyg< Shifted part
3374	! * limit magnitude of mp(i) to meet cfl condition *
3375
3376	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
3377	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
3378	ampmax = min(ampmax, amp2)
3379	mp(il, i) = min(mp(il,i), ampmax)
3380	!>jyg
3381
3382	END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
3383	END DO
3384	! ----------------------------------------------------------------
3385	!
3386	IF (prt_level >= 20) THEN
3387	Print*, 'cv3_unsat after mp computation: mp, b(i), b(i-1) ', &
3388	i, mp(1, i), b(1,i), b(1,max(i-1,1))
3389	ENDIF
3390	!
3391
3392	! * find mixing ratio of precipitating downdraft *
3393
3394	DO il = 1, ncum
3395	IF (i<inb(il) .AND. lwork(il)) THEN
3396	mplus(il) = mp(il, i) > mp(il, i+1)
3397	END IF ! (i.lt.inb(il) .and. lwork(il))
3398	END DO
3399
3400	DO il = 1, ncum
3401	IF (i<inb(il) .AND. lwork(il)) THEN
3402
3403	rp(il, i) = rr(il, i)
3404
3405	IF (mplus(il)) THEN
3406
3407	IF (cvflag_grav) THEN
3408	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
3409	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
3410	ELSE
3411	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
3412	5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
3413	END IF
3414	rp(il, i) = rp(il, i)/mp(il, i)
3415	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
3416	up(il, i) = up(il, i)/mp(il, i)
3417	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
3418	vp(il, i) = vp(il, i)/mp(il, i)
3419
3420	ELSE ! if (mplus(il))
3421
3422	IF (mp(il,i+1)>1.0E-16) THEN
3423	IF (cvflag_grav) THEN
3424	rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
3425	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
3426	ELSE
3427	rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
3428	(evap(il,i+1)+evap(il,i))/mp(il, i+1)
3429	END IF
3430	up(il, i) = up(il, i+1)
3431	vp(il, i) = vp(il, i+1)
3432	END IF ! (mp(il,i+1).gt.1.0e-16)
3433	END IF ! (mplus(il)) else if (.not.mplus(il))
3434
3435	rp(il, i) = amin1(rp(il,i), rs(il,i))
3436	rp(il, i) = max(rp(il,i), 0.0)
3437
3438	END IF ! (i.lt.inb(il) .and. lwork(il))
3439	END DO
3440	! ----------------------------------------------------------------
3441
3442	! * find tracer concentrations in precipitating downdraft *
3443
3444	!AC! do j=1,ntra
3445	!AC! do il = 1,ncum
3446	!AC! if (i.lt.inb(il) .and. lwork(il)) then
3447	!AC!c
3448	!AC! if(mplus(il))then
3449	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
3450	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
3451	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
3452	!AC! else ! if (mplus(il))
3453	!AC! if(mp(il,i+1).gt.1.0e-16)then
3454	!AC! trap(il,i,j)=trap(il,i+1,j)
3455	!AC! endif
3456	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
3457	!AC!c
3458	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
3459	!AC! enddo
3460	!AC! end do
3461
3462	400 END DO
3463	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3464
3465	! * end of downdraft loop *
3466
3467	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3468
3469	RETURN
3470
3471	END SUBROUTINE cv3_unsat
3472
3473	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
3474	icb, inb, delt, &
3475	t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
3476	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
3477	ep, clw, qpreca, m, tp, mp, rp, up, vp, trap, &
3478	wt, water, ice, evap, fondue, faci, b, sigd, &
3479	ment, qent, hent, iflag_mix, uent, vent, &
3480	nent, elij, traent, sig, &
3481	tv, tvp, wghti, &
3482	iflag, precip, Vprecip, Vprecipi, & ! jyg: Vprecipi
3483	ft, fr, fr_comp, fu, fv, ftra, & ! jyg
3484	cbmf, upwd, dnwd, dnwd0, ma, mip, &
3485	!! tls, tps, ! useless . jyg
3486	qcondc, wd, &
3487	ftd, fqd, qta, qtc, sigt, detrain, tau_cld_cv, coefw_cld_cv)
3488
3489	USE conema3_mod_h
3490	USE print_control_mod, ONLY: lunout, prt_level
3491	USE add_phys_tend_mod, only : fl_cor_ebil
3492	USE cvflag_mod_h
3493	USE lmdz_cv_ini, ONLY : grav,minorig,nl,nlp,rowl,rrd,nl,ci,cl,cpd,cpv
3494	USE lmdz_cv_ini, ONLY : restore_bug_cvdn
3495
3496	IMPLICIT NONE
3497
3498
3499	!inputs:
3500	INTEGER, INTENT (IN) :: iflag_mix
3501	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
3502	LOGICAL, INTENT (IN) :: ok_conserv_q
3503	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
3504	REAL, INTENT (IN) :: delt
3505	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, u, v
3506	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t_wake, rr_wake
3507	REAL, DIMENSION (nloc), INTENT (IN) :: s_wake
3508	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra
3509	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
3510	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
3511	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz, h, hp
3512	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tp
3513	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, cpn, ep, clw
3514	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
3515	REAL, DIMENSION (nloc, na), INTENT (IN) :: rp, up
3516	REAL, DIMENSION (nloc, na), INTENT (IN) :: vp
3517	REAL, DIMENSION (nloc, nd), INTENT (IN) :: wt
3518	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: trap
3519	REAL, DIMENSION (nloc, na), INTENT (IN) :: water, evap, b
3520	REAL, DIMENSION (nloc, na), INTENT (IN) :: fondue, faci, ice
3521	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: qent, uent
3522	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: hent
3523	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: vent, elij
3524	INTEGER, DIMENSION (nloc, nd), INTENT (IN) :: nent
3525	REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent
3526	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti
3527	REAL, DIMENSION (nloc, nd), INTENT (IN) :: qta
3528	REAL, DIMENSION (nloc, na),INTENT(IN) :: qpreca
3529	REAL, INTENT(IN) :: tau_cld_cv, coefw_cld_cv
3530	!
3531	!input/output:
3532	REAL, DIMENSION (nloc, na), INTENT (INOUT) :: m, mp
3533	REAL, DIMENSION (nloc, na, na), INTENT (INOUT) :: ment
3534	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
3535	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig
3536	REAL, DIMENSION (nloc), INTENT (INOUT) :: sigd
3537	!
3538	!outputs:
3539	REAL, DIMENSION (nloc), INTENT (OUT) :: precip
3540	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ft, fr, fu, fv , fr_comp
3541	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ftd, fqd
3542	REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT) :: ftra
3543	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: upwd, dnwd, ma
3544	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: dnwd0, mip
3545	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecip
3546	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecipi
3547	!! REAL tls(nloc, nd), tps(nloc, nd) ! useless . jyg
3548	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qcondc ! cld
3549	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qtc, sigt ! cld
3550	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: detrain ! Louis : pour le calcul de Klein du terme de variance qui détraine dans lenvironnement
3551	REAL, DIMENSION (nloc), INTENT (OUT) :: wd ! gust
3552	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf
3553	!
3554	!local variables:
3555	INTEGER :: i, k, il, n, j, num1
3556	REAL :: rat, delti
3557	REAL :: ax, bx, cx, dx, ex
3558	REAL :: cpinv, rdcp, dpinv
3559	REAL :: sigaq
3560	REAL, DIMENSION (nloc) :: awat
3561	REAL, DIMENSION (nloc, nd) :: lvcp, lfcp ! , mke ! unused . jyg
3562	REAL, DIMENSION (nloc) :: am, work, ad, amp1
3563	!! real up1(nloc), dn1(nloc)
3564	REAL, DIMENSION (nloc, nd, nd) :: up1, dn1
3565	!jyg<
3566	REAL, DIMENSION (nloc, nd) :: up_to, up_from
3567	REAL, DIMENSION (nloc, nd) :: dn_to, dn_from
3568	!>jyg
3569	REAL, DIMENSION (nloc) :: asum, bsum, csum, dsum
3570	REAL, DIMENSION (nloc) :: esum, fsum, gsum, hsum
3571	REAL, DIMENSION (nloc, nd) :: th_wake
3572	REAL, DIMENSION (nloc) :: alpha_qpos, alpha_qpos1
3573	REAL, DIMENSION (nloc, nd) :: qcond, nqcond, wa ! cld
3574	REAL, DIMENSION (nloc, nd) :: siga, sax, mac ! cld
3575	REAL, DIMENSION (nloc) :: sument
3576	REAL, DIMENSION (nloc, nd) :: sigment, qtment ! cld
3577	REAL, DIMENSION (nloc, nd, nd) :: qdet
3578	!! REAL sumdq !jyg
3579	!
3580	! -------------------------------------------------------------
3581
3582
3583	! initialization:
3584
3585	delti = 1.0/delt
3586	! print*,'cv3_yield initialisation delt', delt
3587
3588	DO il = 1, ncum
3589	precip(il) = 0.0
3590	wd(il) = 0.0 ! gust
3591	END DO
3592
3593	! Fluxes are on a staggered grid : loops extend up to nl+1
3594	DO i = 1, nlp
3595	DO il = 1, ncum
3596	Vprecip(il, i) = 0.0
3597	Vprecipi(il, i) = 0.0 ! jyg
3598	upwd(il, i) = 0.0
3599	dnwd(il, i) = 0.0
3600	dnwd0(il, i) = 0.0
3601	mip(il, i) = 0.0
3602	END DO
3603	END DO
3604	DO i = 1, nl
3605	DO il = 1, ncum
3606	ft(il, i) = 0.0
3607	fr(il, i) = 0.0
3608	fr_comp(il,i) = 0.0
3609	fu(il, i) = 0.0
3610	fv(il, i) = 0.0
3611	ftd(il, i) = 0.0
3612	fqd(il, i) = 0.0
3613	qcondc(il, i) = 0.0 ! cld
3614	qcond(il, i) = 0.0 ! cld
3615	qtc(il, i) = 0.0 ! cld
3616	qtment(il, i) = 0.0 ! cld
3617	sigment(il, i) = 0.0 ! cld
3618	sigt(il, i) = 0.0 ! cld
3619	qdet(il,i,:) = 0.0 ! cld
3620	detrain(il, i) = 0.0 ! cld
3621	nqcond(il, i) = 0.0 ! cld
3622	END DO
3623	END DO
3624	! print*,'cv3_yield initialisation 2'
3625	!AC! do j=1,ntra
3626	!AC! do i=1,nd
3627	!AC! do il=1,ncum
3628	!AC! ftra(il,i,j)=0.0
3629	!AC! enddo
3630	!AC! enddo
3631	!AC! enddo
3632	! print*,'cv3_yield initialisation 3'
3633	DO i = 1, nl
3634	DO il = 1, ncum
3635	lvcp(il, i) = lv(il, i)/cpn(il, i)
3636	lfcp(il, i) = lf(il, i)/cpn(il, i)
3637	END DO
3638	END DO
3639
3640
3641
3642	! * calculate surface precipitation in mm/day *
3643
3644	DO il = 1, ncum
3645	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
3646	IF (cvflag_ice) THEN
3647	precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
3648	86400.1000./(rowl*grav)
3649	ELSE
3650	precip(il) = wt(il, 1)sigd(il)water(il, 1) &
3651	86400.1000./(rowl*grav)
3652	END IF
3653	END IF
3654	END DO
3655	! print*,'cv3_yield apres calcul precip'
3656
3657
3658	! === calculate vertical profile of precipitation in kg/m2/s ===
3659
3660	DO i = 1, nl
3661	DO il = 1, ncum
3662	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
3663	IF (cvflag_ice) THEN
3664	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
3665	Vprecipi(il, i) = wt(il, i)sigd(il)ice(il,i)/grav ! jyg
3666	ELSE
3667	Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
3668	Vprecipi(il, i) = 0. ! jyg
3669	END IF
3670	END IF
3671	END DO
3672	END DO
3673
3674
3675	! * Calculate downdraft velocity scale *
3676	! * NE PAS UTILISER POUR L'INSTANT *
3677
3678	!! do il=1,ncum
3679	!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
3680	!! /(sigd(il)*p(il,icb(il)))
3681	!! enddo
3682
3683
3684	! * calculate tendencies of lowest level potential temperature *
3685	! * and mixing ratio *
3686
3687	DO il = 1, ncum
3688	work(il) = 1.0/(ph(il,1)-ph(il,2))
3689	cbmf(il) = 0.0
3690	END DO
3691
3692	! - Adiabatic ascent mass flux "ma" and cloud base mass flux "cbmf"
3693	!-----------------------------------------------------------------
3694	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3695	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
3696	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3697	!!! Warning : this option leads to water conservation violation
3698	!!! Expert only
3699	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3700	DO il = 1, ncum
3701	ma(il, nlp) = 0.
3702	ma(il, 1) = 0.
3703	END DO
3704	DO k = nl, 2, -1
3705	DO il = 1, ncum
3706	ma(il, k) = ma(il, k+1)*(1.-qta(il, k))/(1.-qta(il, k-1)) + m(il, k)
3707	cbmf(il) = max(cbmf(il), ma(il,k))
3708	END DO
3709	END DO
3710	DO k = 2,nl
3711	DO il = 1, ncum
3712	IF (k <icb(il)) THEN
3713	ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
3714	ENDIF
3715	END DO
3716	END DO
3717	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3718	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
3719	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3720	!! Line kept for compatibility with earlier versions
3721	DO k = 2, nl
3722	DO il = 1, ncum
3723	IF (k>=icb(il)) THEN
3724	cbmf(il) = cbmf(il) + m(il, k)
3725	END IF
3726	END DO
3727	END DO
3728
3729	DO il = 1, ncum
3730	ma(il, nlp) = 0.
3731	ma(il, 1) = 0.
3732	END DO
3733	DO k = nl, 2, -1
3734	DO il = 1, ncum
3735	ma(il, k) = ma(il, k+1) + m(il, k)
3736	END DO
3737	END DO
3738	DO k = 2,nl
3739	DO il = 1, ncum
3740	IF (k <icb(il)) THEN
3741	ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
3742	ENDIF
3743	END DO
3744	END DO
3745
3746	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
3747	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3748	!
3749	! print*,'cv3_yield avant ft'
3750	! am is the part of cbmf taken from the first level
3751	DO il = 1, ncum
3752	am(il) = cbmf(il)*wghti(il, 1)
3753	END DO
3754
3755	DO il = 1, ncum
3756	IF (iflag(il)<=1) THEN
3757	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
3758	!JYG Correction pour conserver l'eau
3759	! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
3760	IF (cvflag_ice) THEN
3761	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
3762	lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
3763	lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
3764	(100.*(ph(il,1)-ph(il,2))) !precip
3765	ELSE
3766	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
3767	END IF
3768
3769	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
3770
3771	IF (cvflag_ice) THEN
3772	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3773	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
3774	0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
3775	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3776	ELSE
3777	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3778	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3779	END IF
3780
3781	ftd(il, 1) = ft(il, 1) ! fin precip
3782
3783	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
3784	!jyg<
3785	IF (fl_cor_ebil >= 2) THEN
3786	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3787	((t(il,2)-t(il,1))*cpn(il,2)+gz(il,2)-gz(il,1))/cpn(il,1)
3788	ELSE
3789	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3790	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
3791	ENDIF
3792	!>jyg
3793	END IF ! iflag
3794	END DO
3795
3796
3797	DO j = 2, nl
3798	IF (iflag_mix>0) THEN
3799	DO il = 1, ncum
3800	! FH WARNING a modifier :
3801	cpinv = 0.
3802	! cpinv=1.0/cpn(il,1)
3803	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3804	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
3805	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
3806	END IF ! j
3807	END DO
3808	END IF
3809	END DO
3810	! fin sature
3811
3812
3813	DO il = 1, ncum
3814	IF (iflag(il)<=1) THEN
3815	!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
3816	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
3817	sigd(il)*evap(il, 1)
3818	!!! sigd(il)0.5(evap(il,1)+evap(il,2))
3819
3820	fqd(il, 1) = fr(il, 1) !precip
3821
3822	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
3823
3824	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
3825	am(il)*(u(il,2)-u(il,1)))
3826	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
3827	am(il)*(v(il,2)-v(il,1)))
3828	END IF ! iflag
3829	END DO ! il
3830
3831
3832	!AC! do j=1,ntra
3833	!AC! do il=1,ncum
3834	!AC! if (iflag(il) .le. 1) then
3835	!AC! if (cvflag_grav) then
3836	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
3837	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3838	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3839	!AC! else
3840	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
3841	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3842	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3843	!AC! endif
3844	!AC! endif ! iflag
3845	!AC! enddo
3846	!AC! enddo
3847
3848	DO j = 2, nl
3849	DO il = 1, ncum
3850	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3851	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3852	fr_comp(il,1) = fr_comp(il,1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3853	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
3854	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
3855	END IF ! j
3856	END DO
3857	END DO
3858
3859	!AC! do k=1,ntra
3860	!AC! do j=2,nl
3861	!AC! do il=1,ncum
3862	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
3863	!AC!
3864	!AC! if (cvflag_grav) then
3865	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
3866	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3867	!AC! else
3868	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
3869	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3870	!AC! endif
3871	!AC!
3872	!AC! endif
3873	!AC! enddo
3874	!AC! enddo
3875	!AC! enddo
3876	! print*,'cv3_yield apres ft'
3877
3878	!jyg<
3879	!-----------------------------------------------------------
3880	IF (ok_optim_yield) THEN !\|
3881	!-----------------------------------------------------------
3882	!
3883	!* *
3884	!* Compute convective mass fluxes upwd and dnwd *
3885
3886	!
3887	! =================================================
3888	! upward fluxes \|
3889	! ------------------------------------------------
3890	!
3891	upwd(:,:) = 0.
3892	up_to(:,:) = 0.
3893	up_from(:,:) = 0.
3894	!
3895	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3896	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
3897	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3898	!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
3899	!! is taken into account.
3900	!! WARNING : in the present version, taking into account the mass-flux decrease due to
3901	!! precipitation ejection leads to water conservation violation.
3902	!
3903	! - Upward mass flux of mixed draughts
3904	!---------------------------------------
3905	DO i = 2, nl
3906	DO j = 1, i-1
3907	DO il = 1, ncum
3908	IF (i<=inb(il)) THEN
3909	up_to(il,i) = up_to(il,i) + ment(il,j,i)
3910	ENDIF
3911	ENDDO
3912	ENDDO
3913	ENDDO
3914	!
3915	DO j = 3, nl
3916	DO i = 2, j-1
3917	DO il = 1, ncum
3918	IF (j<=inb(il)) THEN
3919	up_from(il,i) = up_from(il,i) + ment(il,i,j)
3920	ENDIF
3921	ENDDO
3922	ENDDO
3923	ENDDO
3924	!
3925	! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
3926	!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
3927	!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
3928	!
3929	DO i = 2, nlp
3930	DO il = 1, ncum
3931	IF (i<=inb(il)+1) THEN
3932	upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
3933	ENDIF
3934	ENDDO
3935	ENDDO
3936	!
3937	! - Total upward mass flux
3938	!---------------------------
3939	DO i = 2, nlp
3940	DO il = 1, ncum
3941	IF (i<=inb(il)+1) THEN
3942	upwd(il,i) = upwd(il,i) + ma(il,i)
3943	ENDIF
3944	ENDDO
3945	ENDDO
3946	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3947	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
3948	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3949	!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
3950	!! is not taken into account.
3951	!
3952	! - Upward mass flux
3953	!-------------------
3954	DO i = 2, nl
3955	DO il = 1, ncum
3956	IF (i<=inb(il)) THEN
3957	up_to(il,i) = m(il,i)
3958	ENDIF
3959	ENDDO
3960	DO j = 1, i-1
3961	DO il = 1, ncum
3962	IF (i<=inb(il)) THEN
3963	up_to(il,i) = up_to(il,i) + ment(il,j,i)
3964	ENDIF
3965	ENDDO
3966	ENDDO
3967	ENDDO
3968	!
3969	DO i = 1, nl
3970	DO il = 1, ncum
3971	IF (i<=inb(il)) THEN
3972	up_from(il,i) = cbmf(il)*wghti(il,i)
3973	ENDIF
3974	ENDDO
3975	ENDDO
3976	!
3977	DO j = 3, nl
3978	DO i = 2, j-1
3979	DO il = 1, ncum
3980	IF (j<=inb(il)) THEN
3981	up_from(il,i) = up_from(il,i) + ment(il,i,j)
3982	ENDIF
3983	ENDDO
3984	ENDDO
3985	ENDDO
3986	!
3987	! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
3988	!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
3989	!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
3990	!
3991	DO i = 2, nlp
3992	DO il = 1, ncum
3993	IF (i<=inb(il)+1) THEN
3994	upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
3995	ENDIF
3996	ENDDO
3997	ENDDO
3998
3999
4000	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
4001	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4002
4003	!
4004	! =================================================
4005	! downward fluxes \|
4006	! ------------------------------------------------
4007	dnwd(:,:) = 0.
4008	dn_to(:,:) = 0.
4009	dn_from(:,:) = 0.
4010	DO i = 1, nl
4011	DO j = i+1, nl
4012	DO il = 1, ncum
4013	IF (j<=inb(il)) THEN
4014	!! dn_to(il,i) = dn_to(il,i) + ment(il,j,i) !jyg,20220202
4015	dn_to(il,i) = dn_to(il,i) - ment(il,j,i)
4016	ENDIF
4017	ENDDO
4018	ENDDO
4019	ENDDO
4020	!
4021	DO j = 1, nl
4022	DO i = j+1, nl
4023	DO il = 1, ncum
4024	IF (i<=inb(il)) THEN
4025	!! dn_from(il,i) = dn_from(il,i) + ment(il,i,j) !jyg,20220202
4026	dn_from(il,i) = dn_from(il,i) - ment(il,i,j)
4027	ENDIF
4028	ENDDO
4029	ENDDO
4030	ENDDO
4031	!
4032	! The difference between dnwd(il,i) and dnwd(il,i+1) is due to downdrafts ending in layer
4033	!(i) (theses drafts cross interface (i+1) but not interface(i)) and to downdrafts
4034	!starting from layer (i) (theses drafts cross interface (i) but not interface(i+1)):
4035	!
4036	DO i = nl-1, 1, -1
4037	DO il = 1, ncum
4038	!! dnwd(il,i) = max(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i)) !jyg,20220202
4039	dnwd(il,i) = min(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i))
4040	ENDDO
4041	ENDDO
4042	! =================================================
4043	!
4044	!-----------------------------------------------------------
4045	ENDIF !(ok_optim_yield) !\|
4046	!-----------------------------------------------------------
4047	!>jyg
4048
4049	! * calculate tendencies of potential temperature and mixing ratio *
4050	! * at levels above the lowest level *
4051
4052	! * first find the net saturated updraft and downdraft mass fluxes *
4053	! * through each level *
4054
4055	!jyg<
4056	!! DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
4057	DO i = 2, nl
4058	!>jyg
4059
4060	num1 = 0
4061	DO il = 1, ncum
4062	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
4063	END DO
4064	!ym IF (num1<=0) GO TO 500
4065	IF (num1<=0) CYCLE
4066
4067	!
4068	!jyg<
4069	!-----------------------------------------------------------
4070	IF (ok_optim_yield) THEN !\|
4071	!-----------------------------------------------------------
4072
4073	! Restoring a bug that was found and corrected in svn release
4074	! 5544; which appears to have a much stronger impact than initially
4075	! thought
4076
4077	if ( restore_bug_cvdn ) then
4078	DO il = 1, ncum
4079	amp1(il) = upwd(il,i+1)
4080	ad(il) = dnwd(il,i)
4081	ENDDO
4082	else
4083	DO il = 1, ncum
4084	amp1(il) = upwd(il,i+1)
4085	ad(il) = - dnwd(il,i)
4086	ENDDO
4087	endif
4088	!-----------------------------------------------------------
4089	ELSE !(ok_optim_yield) !\|
4090	!-----------------------------------------------------------
4091	!>jyg
4092	DO il = 1,ncum
4093	amp1(il) = 0.
4094	ad(il) = 0.
4095	ENDDO
4096
4097	DO k = 1, nl + 1
4098	DO il = 1, ncum
4099	IF (i>=icb(il)) THEN
4100	IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
4101	amp1(il) = amp1(il) + m(il, k)
4102	END IF
4103	ELSE
4104	! AMP1 is the part of cbmf taken from layers I and lower
4105	IF (k<=i) THEN
4106	amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
4107	END IF
4108	END IF
4109	END DO
4110	END DO
4111
4112	DO j = i + 1, nl + 1
4113	DO k = 1, i
4114	!yor! reverted j and k loops
4115	DO il = 1, ncum
4116	!yor! IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN ! the second condition implies the first !
4117	IF (j<=(inb(il)+1)) THEN
4118	amp1(il) = amp1(il) + ment(il, k, j)
4119	END IF
4120	END DO
4121	END DO
4122	END DO
4123
4124	DO k = 1, i - 1
4125	!jyg<
4126	!! DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
4127	DO j = i, nl
4128	!>jyg
4129	DO il = 1, ncum
4130	!yor! IF (i<=inb(il) .AND. j<=inb(il)) THEN ! the second condition implies the 1st !
4131	IF (j<=inb(il)) THEN
4132	ad(il) = ad(il) + ment(il, j, k)
4133	END IF
4134	END DO
4135	END DO
4136	END DO
4137	!
4138	!-----------------------------------------------------------
4139	ENDIF !(ok_optim_yield) !\|
4140	!-----------------------------------------------------------
4141	!
4142	!! print *,'yield, i, amp1, ad', i, amp1(1), ad(1)
4143
4144	DO il = 1, ncum
4145	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4146	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4147	cpinv = 1.0/cpn(il, i)
4148
4149	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
4150	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
4151
4152	! precip
4153	! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
4154	IF (cvflag_ice) THEN
4155	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
4156	sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
4157	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
4158	ELSE
4159	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
4160	END IF
4161
4162	rat = cpn(il, i-1)*cpinv
4163
4164	ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
4165	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
4166	IF (cvflag_ice) THEN
4167	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
4168	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
4169	0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
4170	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
4171	ELSE
4172	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
4173	(t_wake(il,i+1)-t_wake(il,i))dpinv &
4174	cpinv
4175	END IF
4176
4177	ftd(il, i) = ft(il, i)
4178	! fin precip
4179
4180	! sature
4181	!jyg<
4182	IF (fl_cor_ebil >= 2) THEN
4183	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
4184	( amp1(il)( (t(il,i+1)-t(il,i))cpn(il,i+1) + gz(il,i+1)-gz(il,i))*cpinv - &
4185	ad(il)( (t(il,i)-t(il,i-1))cpn(il,i-1) + gz(il,i)-gz(il,i-1))*cpinv)
4186	ELSE
4187	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
4188	(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
4189	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
4190	ENDIF
4191	!>jyg
4192
4193
4194	IF (iflag_mix==0) THEN
4195	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
4196	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
4197	END IF
4198	!
4199	! sb: on ne fait pas encore la correction permettant de mieux
4200	! conserver l'eau:
4201	!JYG: correction permettant de mieux conserver l'eau:
4202	! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
4203	fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
4204	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
4205	fqd(il, i) = fr(il, i) ! precip
4206
4207	fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
4208	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
4209	fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
4210	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
4211
4212
4213	fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
4214	ad(il)*(rr(il,i)-rr(il,i-1)))
4215	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
4216	ad(il)*(u(il,i)-u(il,i-1)))
4217	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
4218	ad(il)*(v(il,i)-v(il,i-1)))
4219
4220	END IF ! i
4221	END DO
4222
4223	!AC! do k=1,ntra
4224	!AC! do il=1,ncum
4225	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4226	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4227	!AC! cpinv=1.0/cpn(il,i)
4228	!AC! if (cvflag_grav) then
4229	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
4230	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
4231	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
4232	!AC! else
4233	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
4234	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
4235	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
4236	!AC! endif
4237	!AC! endif
4238	!AC! enddo
4239	!AC! enddo
4240
4241	DO k = 1, i - 1
4242
4243	DO il = 1, ncum
4244	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
4245	awat(il) = max(awat(il), 0.0)
4246	END DO
4247
4248	IF (iflag_mix/=0) THEN
4249	DO il = 1, ncum
4250	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4251	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4252	cpinv = 1.0/cpn(il, i)
4253	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
4254	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
4255	!
4256	!
4257	END IF ! i
4258	END DO
4259	END IF
4260
4261	DO il = 1, ncum
4262	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4263	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4264	cpinv = 1.0/cpn(il, i)
4265	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
4266	(qent(il,k,i)-awat(il)-rr(il,i))
4267	fr_comp(il,i) = fr_comp(il,i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-awat(il)-rr(il,i))
4268	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
4269	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
4270
4271	! (saturated updrafts resulting from mixing) ! cld
4272	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
4273	qdet(il,k,i) = (qent(il,k,i)-awat(il)) ! cld Louis : specific humidity in detraining water
4274	qtment(il, i) = qtment(il, i) + qent(il,k,i) ! cld
4275	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4276	END IF ! i
4277	END DO
4278	END DO
4279
4280	!AC! do j=1,ntra
4281	!AC! do k=1,i-1
4282	!AC! do il=1,ncum
4283	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4284	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4285	!AC! cpinv=1.0/cpn(il,i)
4286	!AC! if (cvflag_grav) then
4287	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
4288	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4289	!AC! else
4290	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
4291	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4292	!AC! endif
4293	!AC! endif
4294	!AC! enddo
4295	!AC! enddo
4296	!AC! enddo
4297
4298	!jyg<
4299	!! DO k = i, nl + 1
4300	DO k = i, nl
4301	!>jyg
4302
4303	IF (iflag_mix/=0) THEN
4304	DO il = 1, ncum
4305	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
4306	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4307	cpinv = 1.0/cpn(il, i)
4308	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
4309	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
4310
4311
4312	END IF ! i
4313	END DO
4314	END IF
4315
4316	DO il = 1, ncum
4317	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
4318	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4319	cpinv = 1.0/cpn(il, i)
4320
4321	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
4322	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
4323	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
4324	END IF ! i and k
4325	END DO
4326	END DO
4327
4328	!AC! do j=1,ntra
4329	!AC! do k=i,nl+1
4330	!AC! do il=1,ncum
4331	!AC! if (i.le.inb(il) .and. k.le.inb(il)
4332	!AC! $ .and. iflag(il) .le. 1) then
4333	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4334	!AC! cpinv=1.0/cpn(il,i)
4335	!AC! if (cvflag_grav) then
4336	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
4337	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4338	!AC! else
4339	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
4340	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4341	!AC! endif
4342	!AC! endif ! i and k
4343	!AC! enddo
4344	!AC! enddo
4345	!AC! enddo
4346
4347	! sb: interface with the cloud parameterization: ! cld
4348
4349	DO k = i + 1, nl
4350	DO il = 1, ncum
4351	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
4352	! (saturated downdrafts resulting from mixing) ! cld
4353	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
4354	qdet(il,k,i) = qent(il,k,i) ! cld Louis : specific humidity in detraining water
4355	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4356	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4357	END IF ! cld
4358	END DO ! cld
4359	END DO ! cld
4360
4361	!ym BIG Warning : it seems that the k loop is missing !!!
4362	!ym Strong advice to check this
4363	!ym add a k loop temporary
4364
4365	! (particular case: no detraining level is found) ! cld
4366	! Verif merge Dynamico<<<<<<< .working
4367	DO il = 1, ncum ! cld
4368	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
4369	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
4370	!jyg< Bug correction 20180620
4371	! PROBLEM: Should not qent(il,i,i) be taken into account even if nent(il,i)/=0?
4372	!! qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4373	qdet(il,i,i) = qent(il,i,i) ! cld Louis : specific humidity in detraining water
4374	qtment(il, i) = qent(il,i,i) + qtment(il,i) ! cld
4375	!>jyg
4376	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4377	END IF ! cld
4378	END DO ! cld
4379	! Verif merge Dynamico =======
4380	! Verif merge Dynamico DO k = i + 1, nl
4381	! Verif merge Dynamico DO il = 1, ncum !ym k loop added ! cld
4382	! Verif merge Dynamico IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
4383	! Verif merge Dynamico qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
4384	! Verif merge Dynamico qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4385	! Verif merge Dynamico nqcond(il, i) = nqcond(il, i) + 1. ! cld
4386	! Verif merge Dynamico END IF ! cld
4387	! Verif merge Dynamico END DO
4388	! Verif merge Dynamico ENDDO ! cld
4389	! Verif merge Dynamico >>>>>>> .merge-right.r3413
4390
4391	DO il = 1, ncum ! cld
4392	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
4393	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
4394	qtment(il, i) = qtment(il,i)/nqcond(il, i) ! cld
4395	END IF ! cld
4396	END DO
4397
4398	!AC! do j=1,ntra
4399	!AC! do il=1,ncum
4400	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4401	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4402	!AC! cpinv=1.0/cpn(il,i)
4403	!AC!
4404	!AC! if (cvflag_grav) then
4405	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
4406	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
4407	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
4408	!AC! else
4409	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
4410	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
4411	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
4412	!AC! endif
4413	!AC! endif ! i
4414	!AC! enddo
4415	!AC! enddo
4416
4417
4418	500 END DO
4419
4420	!!!JYG<
4421	!!!Conservation de l'eau
4422	!! sumdq = 0.
4423	!! DO k = 1, nl
4424	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4425	!! END DO
4426	!! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4427	!!!JYG>
4428	! * move the detrainment at level inb down to level inb-1 *
4429	! * in such a way as to preserve the vertically *
4430	! * integrated enthalpy and water tendencies *
4431
4432	! Correction bug le 18-03-09
4433	DO il = 1, ncum
4434	IF (iflag(il)<=1) THEN
4435	ax = 0.01gravment(il, inb(il), inb(il))* &
4436	(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
4437	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
4438	ft(il, inb(il)) = ft(il, inb(il)) - ax
4439	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4440	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
4441
4442	bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
4443	(ph(il,inb(il))-ph(il,inb(il)+1))
4444	fr(il, inb(il)) = fr(il, inb(il)) - bx
4445	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4446	(ph(il,inb(il)-1)-ph(il,inb(il)))
4447
4448	cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
4449	(ph(il,inb(il))-ph(il,inb(il)+1))
4450	fu(il, inb(il)) = fu(il, inb(il)) - cx
4451	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4452	(ph(il,inb(il)-1)-ph(il,inb(il)))
4453
4454	dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
4455	(ph(il,inb(il))-ph(il,inb(il)+1))
4456	fv(il, inb(il)) = fv(il, inb(il)) - dx
4457	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4458	(ph(il,inb(il)-1)-ph(il,inb(il)))
4459	END IF !iflag
4460	END DO
4461
4462	!!!JYG<
4463	!!!Conservation de l'eau
4464	!! sumdq = 0.
4465	!! DO k = 1, nl
4466	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4467	!! END DO
4468	!! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4469	!!!JYG>
4470
4471	!AC! do j=1,ntra
4472	!AC! do il=1,ncum
4473	!AC! IF (iflag(il) .le. 1) THEN
4474	!AC! IF (cvflag_grav) then
4475	!AC! ex=0.01gravment(il,inb(il),inb(il))
4476	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
4477	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
4478	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
4479	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
4480	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
4481	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
4482	!AC! else
4483	!AC! ex=0.1*ment(il,inb(il),inb(il))
4484	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
4485	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
4486	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
4487	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
4488	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
4489	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
4490	!AC! ENDIF !cvflag grav
4491	!AC! ENDIF !iflag
4492	!AC! enddo
4493	!AC! enddo
4494
4495
4496	! * homogenize tendencies below cloud base *
4497
4498
4499	DO il = 1, ncum
4500	asum(il) = 0.0
4501	bsum(il) = 0.0
4502	csum(il) = 0.0
4503	dsum(il) = 0.0
4504	esum(il) = 0.0
4505	fsum(il) = 0.0
4506	gsum(il) = 0.0
4507	hsum(il) = 0.0
4508	END DO
4509
4510	!do i=1,nl
4511	!do il=1,ncum
4512	!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
4513	!enddo
4514	!enddo
4515
4516	DO i = 1, nl
4517	DO il = 1, ncum
4518	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
4519	!jyg Saturated part : use T profile
4520	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
4521	!jyg<20140311
4522	!Correction pour conserver l eau
4523	IF (ok_conserv_q) THEN
4524	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
4525	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
4526
4527	ELSE
4528	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
4529	(ph(il,i)-ph(il,i+1))
4530	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
4531	(ph(il,i)-ph(il,i+1))
4532	ENDIF ! (ok_conserv_q)
4533	!jyg>
4534	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
4535	!jyg Unsaturated part : use T_wake profile
4536	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
4537	!jyg<20140311
4538	!Correction pour conserver l eau
4539	IF (ok_conserv_q) THEN
4540	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
4541	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
4542	ELSE
4543	fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
4544	(ph(il,i)-ph(il,i+1))
4545	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
4546	(ph(il,i)-ph(il,i+1))
4547	ENDIF ! (ok_conserv_q)
4548	!jyg>
4549	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
4550	END IF
4551	END DO
4552	END DO
4553
4554	!!!! do 700 i=1,icb(il)-1
4555	IF (ok_homo_tend) THEN
4556	DO i = 1, nl
4557	DO il = 1, ncum
4558	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
4559	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
4560	fqd(il, i) = fsum(il)/gsum(il)
4561	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
4562	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
4563	END IF
4564	END DO
4565	END DO
4566	ENDIF
4567
4568	!jyg<
4569	!Conservation de l'eau
4570	!! sumdq = 0.
4571	!! DO k = 1, nl
4572	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4573	!! END DO
4574	!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4575	!jyg>
4576
4577
4578	! *** Check that moisture stays positive. If not, scale tendencies
4579	! in order to ensure moisture positivity
4580	DO il = 1, ncum
4581	alpha_qpos(il) = 1.
4582	IF (iflag(il)<=1) THEN
4583	IF (fr(il,1)<=0.) THEN
4584	alpha_qpos(il) = max(alpha_qpos(il), (-deltfr(il,1))/(s_wake(il)rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
4585	END IF
4586	END IF
4587	END DO
4588	DO i = 2, nl
4589	DO il = 1, ncum
4590	IF (iflag(il)<=1) THEN
4591	IF (fr(il,i)<=0.) THEN
4592	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
4593	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
4594	END IF
4595	END IF
4596	END DO
4597	END DO
4598	DO il = 1, ncum
4599	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
4600	alpha_qpos(il) = alpha_qpos(il)*1.1
4601	END IF
4602	END DO
4603	!
4604	IF (prt_level .GE. 5) THEN
4605	print *,' CV3_YIELD : alpha_qpos ',alpha_qpos(1)
4606	ENDIF
4607	!
4608	DO il = 1, ncum
4609	IF (iflag(il)<=1) THEN
4610	sigd(il) = sigd(il)/alpha_qpos(il)
4611	precip(il) = precip(il)/alpha_qpos(il)
4612	cbmf(il) = cbmf(il)/alpha_qpos(il)
4613	END IF
4614	END DO
4615	DO i = 1, nl
4616	DO il = 1, ncum
4617	IF (iflag(il)<=1) THEN
4618	fr(il, i) = fr(il, i)/alpha_qpos(il)
4619	ft(il, i) = ft(il, i)/alpha_qpos(il)
4620	fqd(il, i) = fqd(il, i)/alpha_qpos(il)
4621	ftd(il, i) = ftd(il, i)/alpha_qpos(il)
4622	fu(il, i) = fu(il, i)/alpha_qpos(il)
4623	fv(il, i) = fv(il, i)/alpha_qpos(il)
4624	m(il, i) = m(il, i)/alpha_qpos(il)
4625	mp(il, i) = mp(il, i)/alpha_qpos(il)
4626	Vprecip(il, i) = Vprecip(il, i)/alpha_qpos(il)
4627	Vprecipi(il, i) = Vprecipi(il, i)/alpha_qpos(il) ! jyg
4628	END IF
4629	END DO
4630	END DO
4631	!jyg<
4632	!-----------------------------------------------------------
4633	IF (ok_optim_yield) THEN !\|
4634	!-----------------------------------------------------------
4635	DO i = 1, nl
4636	DO il = 1, ncum
4637	IF (iflag(il)<=1) THEN
4638	upwd(il, i) = upwd(il, i)/alpha_qpos(il)
4639	dnwd(il, i) = dnwd(il, i)/alpha_qpos(il)
4640	END IF
4641	END DO
4642	END DO
4643	!-----------------------------------------------------------
4644	ENDIF !(ok_optim_yield) !\|
4645	!-----------------------------------------------------------
4646	!>jyg
4647	DO j = 1, nl !yor! inverted i and j loops
4648	DO i = 1, nl
4649	DO il = 1, ncum
4650	IF (iflag(il)<=1) THEN
4651	ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
4652	END IF
4653	END DO
4654	END DO
4655	END DO
4656
4657	!AC! DO j = 1,ntra
4658	!AC! DO i = 1,nl
4659	!AC! DO il = 1,ncum
4660	!AC! IF (iflag(il) .le. 1) THEN
4661	!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
4662	!AC! ENDIF
4663	!AC! ENDDO
4664	!AC! ENDDO
4665	!AC! ENDDO
4666
4667
4668	! * reset counter and return *
4669
4670	! Reset counter only for points actually convective (jyg)
4671	! In order take into account the possibility of changing the compression,
4672	! reset m, sig and w0 to zero for non-convecting points.
4673	DO il = 1, ncum
4674	IF (iflag(il) < 3) THEN
4675	sig(il, nd) = 2.0
4676	ENDIF
4677	END DO
4678
4679
4680	DO i = 1, nl
4681	DO il = 1, ncum
4682	dnwd0(il, i) = -mp(il, i)
4683	END DO
4684	END DO
4685	!jyg< (loops stop at nl)
4686	!! DO i = nl + 1, nd
4687	!! DO il = 1, ncum
4688	!! dnwd0(il, i) = 0.
4689	!! END DO
4690	!! END DO
4691	!>jyg
4692
4693
4694	!jyg<
4695	!-----------------------------------------------------------
4696	IF (.NOT.ok_optim_yield) THEN !\|
4697	!-----------------------------------------------------------
4698	DO i = 1, nl
4699	DO il = 1, ncum
4700	upwd(il, i) = 0.0
4701	dnwd(il, i) = 0.0
4702	END DO
4703	END DO
4704
4705	!! DO i = 1, nl ! useless; jyg
4706	!! DO il = 1, ncum ! useless; jyg
4707	!! IF (i>=icb(il) .AND. i<=inb(il)) THEN ! useless; jyg
4708	!! upwd(il, i) = 0.0 ! useless; jyg
4709	!! dnwd(il, i) = 0.0 ! useless; jyg
4710	!! END IF ! useless; jyg
4711	!! END DO ! useless; jyg
4712	!! END DO ! useless; jyg
4713
4714	DO i = 1, nl
4715	DO k = 1, nl
4716	DO il = 1, ncum
4717	up1(il, k, i) = 0.0
4718	dn1(il, k, i) = 0.0
4719	END DO
4720	END DO
4721	END DO
4722
4723	!yor! commented original
4724	! DO i = 1, nl
4725	! DO k = i, nl
4726	! DO n = 1, i - 1
4727	! DO il = 1, ncum
4728	! IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
4729	! up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
4730	! dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
4731	! END IF
4732	! END DO
4733	! END DO
4734	! END DO
4735	! END DO
4736	!yor! replaced with
4737	DO i = 1, nl
4738	DO k = i, nl
4739	DO n = 1, i - 1
4740	DO il = 1, ncum
4741	IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! as i always <= k
4742	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
4743	END IF
4744	END DO
4745	END DO
4746	END DO
4747	END DO
4748	DO i = 1, nl
4749	DO n = 1, i - 1
4750	DO k = i, nl
4751	DO il = 1, ncum
4752	IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! i always <= k
4753	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
4754	END IF
4755	END DO
4756	END DO
4757	END DO
4758	END DO
4759	!yor! end replace
4760
4761	DO i = 1, nl
4762	DO k = 1, nl
4763	DO il = 1, ncum
4764	IF (i>=icb(il)) THEN
4765	IF (k>=i .AND. k<=(inb(il))) THEN
4766	upwd(il, i) = upwd(il, i) + m(il, k)
4767	END IF
4768	ELSE
4769	IF (k<i) THEN
4770	upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
4771	END IF
4772	END IF
4773	! c print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
4774	END DO
4775	END DO
4776	END DO
4777
4778	DO i = 2, nl
4779	DO k = i, nl
4780	DO il = 1, ncum
4781	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
4782	IF (i<=inb(il) .AND. k<=inb(il)) THEN
4783	upwd(il, i) = upwd(il, i) + up1(il, k, i)
4784	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
4785	END IF
4786	! c print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
4787	END DO
4788	END DO
4789	END DO
4790
4791
4792	!!!! DO il=1,ncum
4793	!!!! do i=icb(il),inb(il)
4794	!!!!
4795	!!!! upwd(il,i)=0.0
4796	!!!! dnwd(il,i)=0.0
4797	!!!! do k=i,inb(il)
4798	!!!! up1=0.0
4799	!!!! dn1=0.0
4800	!!!! do n=1,i-1
4801	!!!! up1=up1+ment(il,n,k)
4802	!!!! dn1=dn1-ment(il,k,n)
4803	!!!! enddo
4804	!!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
4805	!!!! dnwd(il,i)=dnwd(il,i)+dn1
4806	!!!! enddo
4807	!!!! enddo
4808	!!!!
4809	!!!! ENDDO
4810
4811	!! DO i = 1, nlp
4812	!! DO il = 1, ncum
4813	!! ma(il, i) = 0
4814	!! END DO
4815	!! END DO
4816	!!
4817	!! DO i = 1, nl
4818	!! DO j = i, nl
4819	!! DO il = 1, ncum
4820	!! ma(il, i) = ma(il, i) + m(il, j)
4821	!! END DO
4822	!! END DO
4823	!! END DO
4824
4825	!jyg< (loops stop at nl)
4826	!! DO i = nl + 1, nd
4827	!! DO il = 1, ncum
4828	!! ma(il, i) = 0.
4829	!! END DO
4830	!! END DO
4831	!>jyg
4832
4833	!! DO i = 1, nl
4834	!! DO il = 1, ncum
4835	!! IF (i<=(icb(il)-1)) THEN
4836	!! ma(il, i) = 0
4837	!! END IF
4838	!! END DO
4839	!! END DO
4840
4841	!-----------------------------------------------------------
4842	ENDIF !(.NOT.ok_optim_yield) !\|
4843	!-----------------------------------------------------------
4844	!>jyg
4845
4846	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4847	! determination de la variation de flux ascendant entre
4848	! deux niveau non dilue mip
4849	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4850
4851	DO i = 1, nl
4852	DO il = 1, ncum
4853	mip(il, i) = m(il, i)
4854	END DO
4855	END DO
4856
4857	!jyg< (loops stop at nl)
4858	!! DO i = nl + 1, nd
4859	!! DO il = 1, ncum
4860	!! mip(il, i) = 0.
4861	!! END DO
4862	!! END DO
4863	!>jyg
4864
4865
4866	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4867	! icb represente de niveau ou se trouve la
4868	! base du nuage , et inb le top du nuage
4869	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4870
4871	!! DO i = 1, nd ! unused . jyg
4872	!! DO il = 1, ncum ! unused . jyg
4873	!! mke(il, i) = upwd(il, i) + dnwd(il, i) ! unused . jyg
4874	!! END DO ! unused . jyg
4875	!! END DO ! unused . jyg
4876
4877	!! DO i = 1, nd ! unused . jyg
4878	!! DO il = 1, ncum ! unused . jyg
4879	!! rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd(1.-rr(il,i))+rr(il,i)cpv) ! unused . jyg
4880	!! tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp ! unused . jyg
4881	!! tps(il, i) = tp(il, i) ! unused . jyg
4882	!! END DO ! unused . jyg
4883	!! END DO ! unused . jyg
4884
4885
4886	! * diagnose the in-cloud mixing ratio * ! cld
4887	! * of condensed water * ! cld
4888	!! cld
4889
4890	DO i = 1, nl+1 ! cld
4891	DO il = 1, ncum ! cld
4892	mac(il, i) = 0.0 ! cld
4893	wa(il, i) = 0.0 ! cld
4894	siga(il, i) = 0.0 ! cld
4895	sax(il, i) = 0.0 ! cld
4896	END DO ! cld
4897	END DO ! cld
4898
4899	DO i = minorig, nl ! cld
4900	DO k = i + 1, nl + 1 ! cld
4901	DO il = 1, ncum ! cld
4902	IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
4903	mac(il, i) = mac(il, i) + m(il, k) ! cld
4904	END IF ! cld
4905	END DO ! cld
4906	END DO ! cld
4907	END DO ! cld
4908
4909	DO i = 1, nl ! cld
4910	DO j = 1, i ! cld
4911	DO il = 1, ncum ! cld
4912	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
4913	.AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
4914	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
4915	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
4916	END IF ! cld
4917	END DO ! cld
4918	END DO ! cld
4919	END DO ! cld
4920
4921	DO i = 1, nl ! cld
4922	DO il = 1, ncum ! cld
4923	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
4924	.AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
4925	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
4926	END IF ! cld
4927	END DO ! cld
4928	END DO
4929	! cld
4930	DO i = 1, nl
4931
4932	! 14/01/15 AJ je remets les parties manquantes cf JYG
4933	! Initialize sument to 0
4934
4935	DO il = 1,ncum
4936	sument(il) = 0.
4937	ENDDO
4938
4939	! Sum mixed mass fluxes in sument
4940
4941	DO k = 1,nl
4942	DO il = 1,ncum
4943	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
4944	sument(il) =sument(il) + abs(ment(il,k,i))
4945	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
4946	ENDIF
4947	ENDDO ! il
4948	ENDDO ! k
4949
4950	! 14/01/15 AJ delta n'a rien � faire l�...
4951	DO il = 1, ncum ! cld
4952	!! IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
4953	!! siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
4954	!! rrdtvp(il, i)/p(il, i)/100. ! cld
4955	!!
4956	!! siga(il, i) = min(siga(il,i), 1.0) ! cld
4957	sigaq = 0.
4958	IF (wa(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
4959	siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
4960	rrdtvp(il, i)/p(il, i)/100. ! cld
4961	siga(il, i) = min(siga(il,i), 1.0) ! cld
4962	sigaq = siga(il,i)*qta(il,i-1) ! cld
4963	ENDIF
4964
4965	! IM cf. FH
4966	! 14/01/15 AJ ne correspond pas � ce qui a �t� cod� par JYG et SB
4967
4968	IF (iflag_clw==0) THEN ! cld
4969	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
4970	+(1.-siga(il,i))*qcond(il, i) ! cld
4971
4972
4973	sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1)) ! cld
4974	sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i)) ! cld
4975	!! qtc(il, i) = (siga(il,i)qta(il,i-1)+sigment(il,i)qtment(il,i)) & ! cld
4976	qtc(il, i) = (sigaq+sigment(il,i)*qtment(il,i)) & ! cld
4977	/(siga(il,i)+sigment(il,i)) ! cld
4978	sigt(il,i) = sigment(il, i) + siga(il, i)
4979
4980	! qtc(il, i) = siga(il,i)qta(il,i-1)+(1.-siga(il,i))qtment(il,i) ! cld
4981	! print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)
4982
4983	ELSE IF (iflag_clw==1) THEN ! cld
4984	qcondc(il, i) = qcond(il, i) ! cld
4985	qtc(il,i) = qtment(il,i) ! cld
4986	END IF ! cld
4987
4988	END DO ! cld
4989	END DO
4990	! print*,'cv3_yield fin'
4991
4992	RETURN
4993	END SUBROUTINE cv3_yield
4994
4995	!AC! et !RomP >>>
4996	SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
4997	ment, sigij, da, phi, phi2, d1a, dam, &
4998	ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
4999	icb, inb)
5000	USE lmdz_cv_ini, ONLY : nl,keep_bug_indices_cv3_tracer
5001	USE cvflag_mod_h
5002	USE ioipsl_getin_p_mod, ONLY : getin_p
5003	IMPLICIT NONE
5004
5005
5006	!inputs:
5007	!------
5008	INTEGER, INTENT (IN) :: ncum, nd, na, nloc, len
5009	INTEGER, DIMENSION (len), INTENT (IN) :: icb, inb
5010	REAL, DIMENSION (len, na, na), INTENT (IN) :: ment, sigij, elij
5011	REAL, DIMENSION (len, nd), INTENT (IN) :: clw
5012	REAL, DIMENSION (len, na), INTENT (IN) :: ep
5013	REAL, DIMENSION (len, nd+1), INTENT (IN) :: Vprecip
5014	!ouputs:
5015	!------
5016	REAL, DIMENSION (len, na, na), INTENT (OUT) :: phi, phi2, epmlmMm
5017	REAL, DIMENSION (len, na), INTENT (OUT) :: da, d1a, dam, eplaMm
5018	!
5019	!local variables:
5020	!---------------
5021	! variables pour tracer dans precip de l'AA et des mel
5022	INTEGER i, j, k
5023	REAL epm(nloc, na, na)
5024	!
5025	! variables d'Emanuel : du second indice au troisieme
5026	! ---> tab(i,k,j) -> de l origine k a l arrivee j
5027	! ment, sigij, elij
5028	! variables personnelles : du troisieme au second indice
5029	! ---> tab(i,j,k) -> de k a j
5030	! phi, phi2, epm, epmlmMm
5031
5032
5033	da(:, :) = 0.
5034	d1a(:, :) = 0.
5035	dam(:, :) = 0.
5036	epm(:, :, :) = 0.
5037	eplaMm(:, :) = 0.
5038	epmlmMm(:, :, :) = 0.
5039	phi(:, :, :) = 0.
5040	phi2(:, :, :) = 0.
5041
5042	! fraction deau condensee dans les melanges convertie en precip : epm
5043	! et eau condens�e pr�cipit�e dans masse d'air satur� : l_m*dM_m/dzdz.dzdz
5044	DO j = 1, nl
5045	DO k = 1, nl
5046	DO i = 1, ncum
5047	IF (k>=icb(i) .AND. k<=inb(i) .AND. &
5048	!!jyg j.ge.k.and.j.le.inb(i)) then
5049	!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
5050	j>k .AND. j<=inb(i)) THEN
5051	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
5052	!!
5053	epm(i, j, k) = max(epm(i,j,k), 0.0)
5054	END IF
5055	END DO
5056	END DO
5057	END DO
5058
5059
5060	DO j = 1, nl
5061	DO k = 1, nl
5062	DO i = 1, ncum
5063	IF (k>=icb(i) .AND. k<=inb(i)) THEN
5064	eplaMm(i, j) = eplamm(i, j) + &
5065	ep(i, j)clw(i, j)ment(i, j, k)*(1.-sigij(i,j,k))
5066	END IF
5067	END DO
5068	END DO
5069	END DO
5070
5071	DO j = 1, nl
5072	DO k = 1, j - 1
5073	DO i = 1, ncum
5074	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
5075	epmlmMm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
5076	END IF
5077	END DO
5078	END DO
5079	END DO
5080
5081	! matrices pour calculer la tendance des concentrations dans cvltr.F90
5082	DO j = 1, nl
5083	DO k = 1, nl
5084	DO i = 1, ncum
5085	da(i, j) = da(i, j) + (1.-sigij(i,k,j))*ment(i, k, j)
5086	phi(i, j, k) = sigij(i, k, j)*ment(i, k, j)
5087	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sigij(i,k,j))
5088	IF (k<=j) THEN
5089	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
5090	END IF
5091	END DO
5092	END DO
5093	END DO
5094
5095	IF (keep_bug_indices_cv3_tracer) THEN
5096	DO j = 1, nl
5097	DO k = 1, nl
5098	DO i = 1, ncum
5099	IF (k<=j) THEN
5100	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, k, j)(1.-ep(i,k))*(1.-sigij(i,k,j))
5101	END IF ! (k<=j)
5102	END DO
5103	END DO
5104	END DO
5105	ELSE ! (keep_bug_indices_cv3_tracer)
5106	DO j = 1, nl
5107	DO k = 1, nl
5108	DO i = 1, ncum
5109	IF (k<=j) THEN
5110	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, j, k)(1.-ep(i,k))*(1.-sigij(i,k,j))
5111	END IF ! (k<=j)
5112	END DO
5113	END DO
5114	END DO
5115	ENDIF ! (keep_bug_indices_cv3_tracer)
5116
5117	RETURN
5118	END SUBROUTINE cv3_tracer
5119	!AC! et !RomP <<<
5120
5121	SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
5122	iflag, &
5123	precip, sig, w0, &
5124	ft, fq, fu, fv, ftra, &
5125	Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
5126	epmax_diag, & ! epmax_cape
5127	iflag1, &
5128	precip1, sig1, w01, &
5129	ft1, fq1, fu1, fv1, ftra1, &
5130	Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1, &
5131	epmax_diag1) ! epmax_cape
5132	USE lmdz_cv_ini, ONLY : nl
5133	IMPLICIT NONE
5134
5135
5136	!inputs:
5137	INTEGER len, ncum, nd, ntra, nloc
5138	INTEGER idcum(nloc)
5139	INTEGER iflag(nloc)
5140	REAL precip(nloc)
5141	REAL sig(nloc, nd), w0(nloc, nd)
5142	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
5143	REAL ftra(nloc, nd, ntra)
5144	REAL ma(nloc, nd)
5145	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
5146	REAL qcondc(nloc, nd)
5147	REAL wd(nloc), cape(nloc)
5148	REAL epmax_diag(nloc)
5149
5150	!outputs:
5151	INTEGER iflag1(len)
5152	REAL precip1(len)
5153	REAL sig1(len, nd), w01(len, nd)
5154	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
5155	REAL ftra1(len, nd, ntra)
5156	REAL ma1(len, nd)
5157	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
5158	REAL qcondc1(nloc, nd)
5159	REAL wd1(nloc), cape1(nloc)
5160	REAL epmax_diag1(len) ! epmax_cape
5161
5162	!local variables:
5163	INTEGER i, k, j
5164
5165	DO i = 1, ncum
5166	precip1(idcum(i)) = precip(i)
5167	iflag1(idcum(i)) = iflag(i)
5168	wd1(idcum(i)) = wd(i)
5169	cape1(idcum(i)) = cape(i)
5170	epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
5171	END DO
5172
5173	DO k = 1, nl
5174	DO i = 1, ncum
5175	sig1(idcum(i), k) = sig(i, k)
5176	w01(idcum(i), k) = w0(i, k)
5177	ft1(idcum(i), k) = ft(i, k)
5178	fq1(idcum(i), k) = fq(i, k)
5179	fu1(idcum(i), k) = fu(i, k)
5180	fv1(idcum(i), k) = fv(i, k)
5181	ma1(idcum(i), k) = ma(i, k)
5182	upwd1(idcum(i), k) = upwd(i, k)
5183	dnwd1(idcum(i), k) = dnwd(i, k)
5184	dnwd01(idcum(i), k) = dnwd0(i, k)
5185	qcondc1(idcum(i), k) = qcondc(i, k)
5186	END DO
5187	END DO
5188
5189	DO i = 1, ncum
5190	sig1(idcum(i), nd) = sig(i, nd)
5191	END DO
5192
5193
5194	!AC! do 2100 j=1,ntra
5195	!AC!c oct3 do 2110 k=1,nl
5196	!AC! do 2110 k=1,nd ! oct3
5197	!AC! do 2120 i=1,ncum
5198	!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
5199	!AC! 2120 continue
5200	!AC! 2110 continue
5201	!AC! 2100 continue
5202	!
5203	RETURN
5204	END SUBROUTINE cv3_uncompress
5205
5206
5207	subroutine cv3_epmax_fn_cape(nloc,ncum,nd &
5208	, ep,hp,icb,inb,clw,nk,t,h,hnk,lv,lf,frac &
5209	, pbase, p, ph, tv, buoy, sig, w0,iflag &
5210	, epmax_diag)
5211	USE conema3_mod_h
5212	USE cvflag_mod_h
5213	USE lmdz_cv_ini, ONLY : nl,minorig,cpd,cpv
5214	implicit none
5215
5216	! On fait varier epmax en fn de la cape
5217	! Il faut donc recalculer ep, et hp qui a d�j� �t� calcul� et
5218	! qui en d�pend
5219	! Toutes les autres variables fn de ep sont calcul�es plus bas.
5220
5221
5222	! inputs:
5223	INTEGER, INTENT (IN) :: ncum, nd, nloc
5224	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
5225	REAL, DIMENSION (nloc), INTENT (IN) :: hnk,pbase
5226	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, lv, lf, tv, h
5227	REAL, DIMENSION (nloc, nd), INTENT (IN) :: clw, buoy,frac
5228	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig,w0
5229	INTEGER, DIMENSION (nloc), INTENT (IN) :: iflag(nloc)
5230	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
5231	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
5232	! inouts:
5233	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: ep,hp
5234	! outputs
5235	REAL, DIMENSION (nloc), INTENT (OUT) :: epmax_diag
5236
5237	! local
5238	integer i,k
5239	! real hp_bak(nloc,nd)
5240	! real ep_bak(nloc,nd)
5241	real m_loc(nloc,nd)
5242	real sig_loc(nloc,nd)
5243	real w0_loc(nloc,nd)
5244	integer iflag_loc(nloc)
5245	real cape(nloc)
5246
5247	if (coef_epmax_cape.gt.1e-12) then
5248
5249	! il faut calculer la cape: on fait un calcule simple car tant qu'on ne
5250	! connait pas ep, on ne connait pas les m�langes, ddfts etc... qui sont
5251	! necessaires au calcul de la cape dans la nouvelle physique
5252
5253	! write(,) 'cv3_routines check 4303'
5254	do i=1,ncum
5255	do k=1,nd
5256	sig_loc(i,k)=sig(i,k)
5257	w0_loc(i,k)=w0(i,k)
5258	iflag_loc(i)=iflag(i)
5259	! ep_bak(i,k)=ep(i,k)
5260	enddo ! do k=1,nd
5261	enddo !do i=1,ncum
5262
5263	! write(,) 'cv3_routines check 4311'
5264	! write(,) 'nl=',nl
5265	CALL cv3_closure(nloc, ncum, nd, icb, inb, & ! na->nd
5266	pbase, p, ph, tv, buoy, &
5267	sig_loc, w0_loc, cape, m_loc,iflag_loc)
5268
5269	! write(,) 'cv3_routines check 4316'
5270	! write(,) 'ep(1,:)=',ep(1,:)
5271	do i=1,ncum
5272	epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
5273	epmax_diag(i)=amax1(epmax_diag(i),0.0)
5274	! write(,) 'i,icb,inb,cape,epmax_diag=', &
5275	! i,icb(i),inb(i),cape(i),epmax_diag(i)
5276	do k=1,nl
5277	ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
5278	ep(i,k)=amax1(ep(i,k),0.0)
5279	ep(i,k)=amin1(ep(i,k),epmax_diag(i))
5280	enddo
5281	enddo
5282	! write(,) 'ep(1,:)=',ep(1,:)
5283
5284	!write(,) 'cv3_routines check 4326'
5285	! On recalcule hp:
5286	! do k=1,nl
5287	! do i=1,ncum
5288	! hp_bak(i,k)=hp(i,k)
5289	! enddo
5290	! enddo
5291	do k=1,nl
5292	do i=1,ncum
5293	hp(i,k)=h(i,k)
5294	enddo
5295	enddo
5296
5297	IF (cvflag_ice) THEN
5298
5299	do k=minorig+1,nl
5300	do i=1,ncum
5301	if((k.ge.icb(i)).and.(k.le.inb(i)))then
5302	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac(i,k)lf(i,k))* &
5303	ep(i, k)*clw(i, k)
5304	endif
5305	enddo
5306	enddo !do k=minorig+1,n
5307	ELSE !IF (cvflag_ice) THEN
5308
5309	DO k = minorig + 1, nl
5310	DO i = 1, ncum
5311	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
5312	hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
5313	endif
5314	enddo
5315	enddo !do k=minorig+1,n
5316
5317	ENDIF !IF (cvflag_ice) THEN
5318	!write(,) 'cv3_routines check 4345'
5319	! do i=1,ncum
5320	! do k=1,nl
5321	! if ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-1).or. &
5322	! ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-4).and. &
5323	! (ep(i,k)-ep_bak(i,k).lt.1e-4))) then
5324	! write(,) 'i,k=',i,k
5325	! write(,) 'coef_epmax_cape=',coef_epmax_cape
5326	! write(,) 'epmax_diag(i)=',epmax_diag(i)
5327	! write(,) 'ep(i,k)=',ep(i,k)
5328	! write(,) 'ep_bak(i,k)=',ep_bak(i,k)
5329	! write(,) 'hp(i,k)=',hp(i,k)
5330	! write(,) 'hp_bak(i,k)=',hp_bak(i,k)
5331	! write(,) 'h(i,k)=',h(i,k)
5332	! write(,) 'nk(i)=',nk(i)
5333	! write(,) 'h(i,nk(i))=',h(i,nk(i))
5334	! write(,) 'lv(i,k)=',lv(i,k)
5335	! write(,) 't(i,k)=',t(i,k)
5336	! write(,) 'clw(i,k)=',clw(i,k)
5337	! write(,) 'cpd,cpv=',cpd,cpv
5338	! stop
5339	! endif
5340	! enddo !do k=1,nl
5341	! enddo !do i=1,ncum
5342	endif !if (coef_epmax_cape.gt.1e-12) then
5343	!write(,) 'cv3_routines check 4367'
5344
5345	return
5346	end subroutine cv3_epmax_fn_cape
5347
5348	END MODULE cv3_routines_mod
5349

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

Download in other formats:

Original Format