Context Navigation

cv3_routines.F90 @ 5140

Last change on this file since 5140 was 5140, checked in by abarral, 3 months ago
Put comsoil.h, conema3.h, cvflag.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: 179.8 KB

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

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

Download in other formats:

Original Format