Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

cv3_routines.f90

Last change on this file was 5305, checked in by abarral, 42 hours ago
Turn YOMCST2.h.h into module
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: 180.0 KB

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

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

Download in other formats:

Original Format