Context Navigation

cv3_routines.f90 @ 5692

Last change on this file since 5692 was 5692, checked in by yann meurdesoif, 5 weeks ago

Convection GPU porting : set convection subroutines into module

Files will be renamed later to *_mod.f90

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

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

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

Download in other formats:

Original Format