Context Navigation

cv3_routines.f90 @ 5698

Last change on this file since 5698 was 5698, checked in by yann meurdesoif, 4 weeks ago
Convection GPU porting : fix for previous commit YM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 183.1 KB

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

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

Download in other formats:

Original Format