Context Navigation

cv3_routines.f90 @ 5898

Last change on this file since 5898 was 5843, checked in by jyg, 3 months ago
Getting rid of "ments" and "qents" arrays within cva_driver. Several comments to be cleared later.
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: 179.2 KB

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

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

Download in other formats:

Original Format