Context Navigation

cv3_routines.F90 @ 5218

Last change on this file since 5218 was 5160, checked in by abarral, 5 months ago
Put .h into modules
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: 185.4 KB

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

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

Download in other formats:

Original Format