Context Navigation

cv3_routines.f90 @ 5840

Last change on this file since 5840 was 5840, checked in by jyg, 2 months ago
Getting rid of tracer arrays within cva_driver. Lot of 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: 178.5 KB

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

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

Download in other formats:

Original Format