Context Navigation

cv3_routines.f90 @ 5300

Last change on this file since 5300 was 5299, checked in by abarral, 4 weeks ago
Turn cv3param.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 5299 2024-10-30 13:31:56Z 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	IMPLICIT NONE
1136
1137	! ---------------------------------------------------------------------
1138	! Purpose:
1139	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1140	! &
1141	! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
1142	! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
1143	! &
1144	! FIND THE LEVEL OF NEUTRAL BUOYANCY
1145
1146	! Main differences convect3/convect4:
1147	! - icbs (input) is the first level above LCL (may differ from icb)
1148	! - many minor differences in the iterations
1149	! - condensed water not removed from tvp in convect3
1150	! - vertical profile of buoyancy computed here (use of buoybase)
1151	! - the determination of inb is different
1152	! - no inb1, only inb in output
1153	! ---------------------------------------------------------------------
1154
1155	include "YOMCST2.h"
1156
1157	!inputs:
1158	INTEGER, INTENT (IN) :: ncum, nd, nloc
1159	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, icbs, nk
1160	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, q, qs, gz
1161	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
1162	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
1163	REAL, DIMENSION (nloc), INTENT (IN) :: tnk, qnk, gznk
1164	REAL, DIMENSION (nloc), INTENT (IN) :: hnk
1165	REAL, DIMENSION (nloc, nd), INTENT (IN) :: lv, lf, tv, h
1166	REAL, DIMENSION (nloc), INTENT (IN) :: pbase, buoybase, plcl
1167
1168	!input/outputs:
1169	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: tp, tvp, clw ! Input for k = 1, icb+1 (computed in cv3_undilute1)
1170	! Output above
1171	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
1172
1173	!outputs:
1174	INTEGER, DIMENSION (nloc), INTENT (OUT) :: inb
1175	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ep, sigp, hp
1176	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: buoy
1177	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: frac_a, frac_s
1178	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qpreca
1179	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qta
1180
1181	!local variables:
1182	INTEGER i, j, k
1183	REAL smallestreal
1184	REAL tg, qg, dqgdT, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
1185	REAL :: phinu2p
1186	REAL :: qhthreshold
1187	REAL :: als
1188	REAL :: qsat_new, snew
1189	REAL, DIMENSION (nloc,nd) :: qi
1190	REAL, DIMENSION (nloc,nd) :: ha ! moist static energy of adiabatic ascents
1191	! taking into account precip ejection
1192	REAL, DIMENSION (nloc,nd) :: hla ! liquid water static energy of adiabatic ascents
1193	! taking into account precip ejection
1194	REAL, DIMENSION (nloc,nd) :: qcld ! specific cloud water
1195	REAL, DIMENSION (nloc,nd) :: qhsat ! specific humidity at saturation
1196	REAL, DIMENSION (nloc,nd) :: dqhsatdT ! dqhsat/dT
1197	REAL, DIMENSION (nloc,nd) :: frac ! ice fraction function of envt temperature
1198	REAL, DIMENSION (nloc,nd) :: qps ! specific solid precipitation
1199	REAL, DIMENSION (nloc,nd) :: qpl ! specific liquid precipitation
1200	REAL, DIMENSION (nloc) :: ah0, cape, capem, byp
1201	LOGICAL, DIMENSION (nloc) :: lcape
1202	INTEGER, DIMENSION (nloc) :: iposit
1203	REAL :: denomm1
1204	REAL :: by, defrac, pden, tbis
1205	REAL :: fracg
1206	REAL :: deltap
1207	REAL, SAVE :: Tx, Tm
1208	DATA Tx/263.15/, Tm/243.15/
1209	!$OMP THREADPRIVATE(Tx, Tm)
1210	REAL :: aa, bb, dd, ddelta, discr
1211	REAL :: ff, fp
1212	REAL :: coefx, coefm, Zx, Zm, Ux, U, Um
1213
1214	IF (prt_level >= 10) THEN
1215	print *,'cv3_undilute2.0. icvflag_Tpa, t(1,k), q(1,k), qs(1,k) ', &
1216	icvflag_Tpa, (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
1217	ENDIF
1218	smallestreal=tiny(smallestreal)
1219
1220	! =====================================================================
1221	! --- SOME INITIALIZATIONS
1222	! =====================================================================
1223
1224	DO k = 1, nl
1225	DO i = 1, ncum
1226	qi(i, k) = 0.
1227	END DO
1228	END DO
1229
1230
1231	! =====================================================================
1232	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1233	! =====================================================================
1234
1235	! --- The procedure is to solve the equation.
1236	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
1237
1238	! * Calculate certain parcel quantities, including static energy *
1239
1240
1241	DO i = 1, ncum
1242	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))*tnk(i)+ &
1243	! debug qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
1244	qnk(i)(lv0-clmcpv(tnk(i)-273.15)) + gznk(i)
1245	END DO
1246	!
1247	! Ice fraction
1248	!
1249	IF (cvflag_ice) THEN
1250	DO k = minorig, nl
1251	DO i = 1, ncum
1252	frac(i, k) = (Tx - t(i,k))/(Tx - Tm)
1253	frac(i, k) = min(max(frac(i,k),0.0), 1.0)
1254	END DO
1255	END DO
1256	! Below cloud base, set ice fraction to cloud base value
1257	DO k = 1, nl
1258	DO i = 1, ncum
1259	IF (k<icb(i)) THEN
1260	frac(i,k) = frac(i,icb(i))
1261	END IF
1262	END DO
1263	END DO
1264	ELSE
1265	DO k = 1, nl
1266	DO i = 1, ncum
1267	frac(i,k) = 0.
1268	END DO
1269	END DO
1270	ENDIF ! (cvflag_ice)
1271
1272
1273	DO k = minorig, nl
1274	DO i = 1,ncum
1275	ha(i,k) = ah0(i)
1276	hla(i,k) = hnk(i)
1277	qta(i,k) = qnk(i)
1278	qpreca(i,k) = 0.
1279	frac_a(i,k) = 0.
1280	frac_s(i,k) = frac(i,k)
1281	qpl(i,k) = 0.
1282	qps(i,k) = 0.
1283	qhsat(i,k) = qs(i,k)
1284	qcld(i,k) = max(qta(i,k)-qhsat(i,k),0.)
1285	IF (k <= icb(i)+1) THEN
1286	qhsat(i,k) = qnk(i)-clw(i,k)
1287	qcld(i,k) = clw(i,k)
1288	ENDIF
1289	ENDDO
1290	ENDDO
1291
1292	!jyg<
1293	! =====================================================================
1294	! --- SET THE THE FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
1295	! =====================================================================
1296	DO k = 1, nl
1297	DO i = 1, ncum
1298	ep(i, k) = 0.0
1299	sigp(i, k) = spfac
1300	END DO
1301	END DO
1302	!>jyg
1303	!
1304
1305	! * Find lifted parcel quantities above cloud base *
1306
1307	!----------------------------------------------------------------------------
1308	!
1309	IF (icvflag_Tpa == 2) THEN
1310	!
1311	!----------------------------------------------------------------------------
1312	!
1313	DO k = minorig + 1, nl
1314	DO i = 1,ncum
1315	tp(i,k) = t(i,k)
1316	ENDDO
1317	!! alv = lv0 - clmcpv*(t(i,k)-273.15)
1318	!! alf = lf0 + clmci*(t(i,k)-273.15)
1319	!! als = alf + alv
1320	DO j = 1,4
1321	DO i = 1, ncum
1322	! ori if(k.ge.(icb(i)+1))then
1323	IF (k>=(icbs(i)+1)) THEN ! convect3
1324	tg = tp(i, k)
1325	IF (tg .gt. Tx) THEN
1326	es = 6.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
1327	qg = epses/(p(i,k)-es(1.-eps))
1328	ELSE
1329	esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1330	qg = epsesi/(p(i,k)-esi(1.-eps))
1331	ENDIF
1332	! Ice fraction
1333	ff = 0.
1334	fp = 1./(Tx - Tm)
1335	IF (tg < Tx) THEN
1336	IF (tg > Tm) THEN
1337	ff = (Tx - tg)*fp
1338	ELSE
1339	ff = 1.
1340	ENDIF ! (tg > Tm)
1341	ENDIF ! (tg < Tx)
1342	! Intermediate variables
1343	aa = cpd + (cl-cpd)qnk(i) + lv(i,k)lv(i,k)qg/(rrvtg*tg)
1344	ahg = (cpd + (cl-cpd)qnk(i))tg + lv(i,k)*qg - &
1345	lf(i,k)ff(qnk(i) - qg) + gz(i,k)
1346	dd = lf(i,k)lv(i,k)qg/(rrvtgtg)
1347	ddelta = lf(i,k)*(qnk(i) - qg)
1348	bb = aa + ddeltafp + ddfp*(Tx-tg)
1349	! Compute Zx and Zm
1350	coefx = aa
1351	coefm = aa + dd
1352	IF (tg .gt. Tx) THEN
1353	Zx = ahg + coefx*(Tx - tg)
1354	Zm = ahg - ddelta + coefm*(Tm - tg)
1355	ELSE
1356	IF (tg .gt. Tm) THEN
1357	Zx = ahg + (coefx +fpddelta)(Tx - Tg)
1358	Zm = ahg + (coefm +fpddelta)(Tm - Tg)
1359	ELSE
1360	Zx = ahg + ddelta + coefx*(Tx - tg)
1361	Zm = ahg + coefm*(Tm - tg)
1362	ENDIF ! (tg .gt. Tm)
1363	ENDIF ! (tg .gt. Tx)
1364	! Compute the masks Um, U, Ux
1365	Um = (sign(1., Zm-ah0(i))+1.)/2.
1366	Ux = (sign(1., ah0(i)-Zx)+1.)/2.
1367	U = (1. - Um)*(1. - Ux)
1368	! Compute the updated parcell temperature Tp : 3 cases depending on tg value
1369	IF (tg .gt. Tx) THEN
1370	discr = bbbb - 4ddfp(ah0(i) - ahg + ddeltafp(Tx-tg))
1371	Tp(i,k) = tg + &
1372	Um* (ah0(i) - ahg + ddelta) /(aa + dd) + &
1373	U 2(ah0(i) - ahg + ddeltafp(Tx-tg))/(bb + sqrt(discr)) + &
1374	Ux* (ah0(i) - ahg) /aa
1375	ELSEIF (tg .gt. Tm) THEN
1376	discr = bbbb - 4ddfp(ah0(i) - ahg)
1377	Tp(i,k) = tg + &
1378	Um* (ah0(i) - ahg + ddeltafp(tg-Tm))/(aa + dd) + &
1379	U 2(ah0(i) - ahg) /(bb + sqrt(discr)) + &
1380	Ux* (ah0(i) - ahg + ddeltafp(tg-Tx))/aa
1381	ELSE
1382	discr = bbbb - 4ddfp(ah0(i) - ahg + ddeltafp(Tm-tg))
1383	Tp(i,k) = tg + &
1384	Um* (ah0(i) - ahg) /(aa + dd) + &
1385	U 2(ah0(i) - ahg + ddeltafp(Tm-tg))/(bb + sqrt(discr)) + &
1386	Ux* (ah0(i) - ahg - ddelta) /aa
1387	ENDIF ! (tg .gt. Tx)
1388	!
1389	!! print *,' j, k, Um, U, Ux, aa, bb, discr, dd, ddelta ', j, k, Um, U, Ux, aa, bb, discr, dd, ddelta
1390	!! print *,' j, k, ah0(i), ahg, tg, qg, tp(i,k), ff ', j, k, ah0(i), ahg, tg, qg, tp(i,k), ff
1391	END IF ! (k>=(icbs(i)+1))
1392	END DO ! i = 1, ncum
1393	END DO ! j = 1,4
1394	DO i = 1, ncum
1395	IF (k>=(icbs(i)+1)) THEN ! convect3
1396	tg = tp(i, k)
1397	IF (tg .gt. Tx) THEN
1398	es = 6.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
1399	qg = epses/(p(i,k)-es(1.-eps))
1400	ELSE
1401	esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1402	qg = epsesi/(p(i,k)-esi(1.-eps))
1403	ENDIF
1404	clw(i, k) = qnk(i) - qg
1405	clw(i, k) = max(0.0, clw(i,k))
1406	tvp(i, k) = max(0., tp(i,k)*(1.+qg/eps-qnk(i)))
1407	! print*,tvp(i,k),'tvp'
1408	IF (clw(i,k)<1.E-11) THEN
1409	tp(i, k) = tv(i, k)
1410	tvp(i, k) = tv(i, k)
1411	END IF ! (clw(i,k)<1.E-11)
1412	END IF ! (k>=(icbs(i)+1))
1413	END DO ! i = 1, ncum
1414	END DO ! k = minorig + 1, nl
1415	!----------------------------------------------------------------------------
1416	!
1417	ELSE IF (icvflag_Tpa == 1) THEN ! (icvflag_Tpa == 2)
1418	!
1419	!----------------------------------------------------------------------------
1420	!
1421	DO k = minorig + 1, nl
1422	DO i = 1,ncum
1423	tp(i,k) = t(i,k)
1424	ENDDO
1425	!! alv = lv0 - clmcpv*(t(i,k)-273.15)
1426	!! alf = lf0 + clmci*(t(i,k)-273.15)
1427	!! als = alf + alv
1428	DO j = 1,4
1429	DO i = 1, ncum
1430	! ori if(k.ge.(icb(i)+1))then
1431	IF (k>=(icbs(i)+1)) THEN ! convect3
1432	tg = tp(i, k)
1433	IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
1434	es = 6.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
1435	qg = epses/(p(i,k)-es(1.-eps))
1436	dqgdT = lv(i,k)qg/(rrvtg*tg)
1437	ELSE
1438	esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1439	qg = epsesi/(p(i,k)-esi(1.-eps))
1440	dqgdT = (lv(i,k)+lf(i,k))qg/(rrvtg*tg)
1441	ENDIF
1442	IF (qsat_depends_on_qt) THEN
1443	dqgdT = dqgdT(1.-qta(i,k-1))/(1.-qg)*2
1444	qg = qg*(1.-qta(i,k-1))/(1.-qg)
1445	ENDIF
1446	ahg = (cpd + (cl-cpd)qta(i,k-1))tg + lv(i,k)*qg - &
1447	lf(i,k)frac(i,k)(qta(i,k-1) - qg) + gz(i,k)
1448	Tp(i,k) = tg + (ah0(i) - ahg)/ &
1449	(cpd + (cl-cpd)qta(i,k-1) + (lv(i,k)+frac(i,k)lf(i,k))*dqgdT)
1450	!! print *,'undilute2 iterations k, Tp(i,k), ah0(i), ahg ', &
1451	!! k, Tp(i,k), ah0(i), ahg
1452	END IF ! (k>=(icbs(i)+1))
1453	END DO ! i = 1, ncum
1454	END DO ! j = 1,4
1455	DO i = 1, ncum
1456	IF (k>=(icbs(i)+1)) THEN ! convect3
1457	tg = tp(i, k)
1458	IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
1459	es = 6.112exp(17.67(tg - 273.15)/(tg + 243.5 - 273.15))
1460	qg = epses/(p(i,k)-es(1.-eps))
1461	ELSE
1462	esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
1463	qg = epsesi/(p(i,k)-esi(1.-eps))
1464	ENDIF
1465	IF (qsat_depends_on_qt) THEN
1466	qg = qg*(1.-qta(i,k-1))/(1.-qg)
1467	ENDIF
1468	qhsat(i,k) = qg
1469	END IF ! (k>=(icbs(i)+1))
1470	END DO ! i = 1, ncum
1471	DO i = 1, ncum
1472	IF (k>=(icbs(i)+1)) THEN ! convect3
1473	clw(i, k) = qta(i,k-1) - qhsat(i,k)
1474	clw(i, k) = max(0.0, clw(i,k))
1475	tvp(i, k) = max(0., tp(i,k)*(1.+qhsat(i,k)/eps-qta(i,k-1)))
1476	! print*,tvp(i,k),'tvp'
1477	IF (clw(i,k)<1.E-11) THEN
1478	tp(i, k) = tv(i, k)
1479	tvp(i, k) = tv(i, k)
1480	END IF ! (clw(i,k)<1.E-11)
1481	END IF ! (k>=(icbs(i)+1))
1482	END DO ! i = 1, ncum
1483	!
1484	IF (cvflag_prec_eject) THEN
1485	DO i = 1, ncum
1486	IF (k>=(icbs(i)+1)) THEN ! convect3
1487	! Specific precipitation (liquid and solid) and ice content
1488	! before ejection of precipitation !!jygprl
1489	elacrit = elcrit*min(max(1.-(tp(i,k)-T0)/Tlcrit, 0.), 1.) !!jygprl
1490	!!!! qcld(i,k) = min(clw(i,k), elacrit) !!jygprl
1491	qhthreshold = elacrit*(1.-qta(i,k-1))/(1.-elacrit)
1492	qcld(i,k) = min(clw(i,k), qhthreshold) !!jygprl
1493	!!!! phinu2p = max(qhsat(i,k-1) + qcld(i,k-1) - (qhsat(i,k) + qcld(i,k)),0.) !!jygprl
1494	phinu2p = max(clw(i,k) - max(qta(i,k-1) - qhsat(i,k-1), qhthreshold), 0.)
1495	qpl(i,k) = qpl(i,k-1) + (1.-frac(i,k))*phinu2p !!jygprl
1496	qps(i,k) = qps(i,k-1) + frac(i,k) *phinu2p !!jygprl
1497	qi(i,k) = (1.-ejectliq)clw(i,k)frac(i,k) + & !!jygprl
1498	ejectliq(qps(i,k-1) + frac(i,k)(phinu2p+qcld(i,k))) !!jygprl
1499	!!
1500	! =====================================================================================
1501	! Ejection of precipitation from adiabatic ascents if requested (cvflag_prec_eject=True):
1502	! Compute the steps of total water (qta), of moist static energy (ha), of specific
1503	! precipitation (qpl and qps) and of specific cloud water (qcld) associated with precipitation
1504	! ejection.
1505	! =====================================================================================
1506	!
1507	! Verif
1508	qpreca(i,k) = ejectliqqpl(i,k) + ejecticeqps(i,k) !!jygprl
1509	frac_a(i,k) = ejectice*qps(i,k)/max(qpreca(i,k),smallestreal) !!jygprl
1510	frac_s(i,k) = (1.-ejectliq)*frac(i,k) + & !!jygprl
1511	ejectliq(1. - (qpl(i,k)+(1.-frac(i,k))qcld(i,k))/max(clw(i,k),smallestreal)) !!jygprl
1512	!
1513	denomm1 = 1./(1. - qpreca(i,k))
1514	!
1515	qta(i,k) = qta(i,k-1) - &
1516	qpreca(i,k)(1.-qta(i,k-1))denomm1
1517	ha(i,k) = ha(i,k-1) + &
1518	( qpreca(i,k)(-(1.-qta(i,k-1))(cl-cpd)*tp(i,k) + &
1519	lv(i,k)qhsat(i,k) - lf(i,k)(frac_s(i,k)*qcld(i,k)+qps(i,k))) + &
1520	lf(i,k)ejecticeqps(i,k))*denomm1
1521	hla(i,k) = hla(i,k-1) + &
1522	( qpreca(i,k)(-(1.-qta(i,k-1))(cpv-cpd)*tp(i,k) - &
1523	lv(i,k)((1.-frac_s(i,k))qcld(i,k)+qpl(i,k)) - &
1524	(lv(i,k)+lf(i,k))(frac_s(i,k)qcld(i,k)+qps(i,k))) + &
1525	lv(i,k)ejectliqqpl(i,k) + (lv(i,k)+lf(i,k))ejecticeqps(i,k))*denomm1
1526	qpl(i,k) = qpl(i,k)(1.-ejectliq)denomm1
1527	qps(i,k) = qps(i,k)(1.-ejectice)denomm1
1528	qcld(i,k) = qcld(i,k)*denomm1
1529	qhsat(i,k) = qhsat(i,k)*(1.-qta(i,k))/(1.-qta(i,k-1))
1530	END IF ! (k>=(icbs(i)+1))
1531	END DO ! i = 1, ncum
1532	ENDIF ! (cvflag_prec_eject)
1533	!
1534	END DO ! k = minorig + 1, nl
1535	!
1536	!----------------------------------------------------------------------------
1537	!
1538	ELSE IF (icvflag_Tpa == 0) THEN! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
1539	!
1540	!----------------------------------------------------------------------------
1541	!
1542	DO k = minorig + 1, nl
1543	DO i = 1, ncum
1544	! ori if(k.ge.(icb(i)+1))then
1545	IF (k>=(icbs(i)+1)) THEN ! convect3
1546	tg = t(i, k)
1547	qg = qs(i, k)
1548	! debug alv=lv0-clmcpv*(t(i,k)-t0)
1549	alv = lv0 - clmcpv*(t(i,k)-273.15)
1550
1551	! First iteration.
1552
1553	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1554	s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
1555	alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
1556	s = 1./s
1557	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1558	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1559	tg = tg + s*(ah0(i)-ahg)
1560	! ori tg=max(tg,35.0)
1561	! debug tc=tg-t0
1562	tc = tg - 273.15
1563	denom = 243.5 + tc
1564	denom = max(denom, 1.0) ! convect3
1565	! ori if(tc.ge.0.0)then
1566	es = 6.112exp(17.67tc/denom)
1567	! ori else
1568	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1569	! ori endif
1570	qg = epses/(p(i,k)-es(1.-eps))
1571
1572	! Second iteration.
1573
1574	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1575	! ori s=1./s
1576	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1577	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1578	tg = tg + s*(ah0(i)-ahg)
1579	! ori tg=max(tg,35.0)
1580	! debug tc=tg-t0
1581	tc = tg - 273.15
1582	denom = 243.5 + tc
1583	denom = max(denom, 1.0) ! convect3
1584	! ori if(tc.ge.0.0)then
1585	es = 6.112exp(17.67tc/denom)
1586	! ori else
1587	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1588	! ori endif
1589	qg = epses/(p(i,k)-es(1.-eps))
1590
1591	! debug alv=lv0-clmcpv*(t(i,k)-t0)
1592	alv = lv0 - clmcpv*(t(i,k)-273.15)
1593	! print*,'cpd dans convect2 ',cpd
1594	! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
1595	! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
1596
1597	! ori c approximation here:
1598	! ori tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
1599
1600	! convect3: no approximation:
1601	IF (cvflag_ice) THEN
1602	tp(i, k) = max(0., (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i)))
1603	ELSE
1604	tp(i, k) = (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
1605	END IF
1606
1607	clw(i, k) = qnk(i) - qg
1608	clw(i, k) = max(0.0, clw(i,k))
1609	rg = qg/(1.-qnk(i))
1610	! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
1611	! convect3: (qg utilise au lieu du vrai mixing ratio rg):
1612	tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
1613	IF (cvflag_ice) THEN
1614	IF (clw(i,k)<1.E-11) THEN
1615	tp(i, k) = tv(i, k)
1616	tvp(i, k) = tv(i, k)
1617	END IF
1618	END IF
1619	!jyg<
1620	!! END IF ! Endif moved to the end of the loop
1621	!>jyg
1622
1623	IF (cvflag_ice) THEN
1624	!CR:attention boucle en klon dans Icefrac
1625	! Call Icefrac(t,clw,qi,nl,nloc)
1626	IF (t(i,k)>263.15) THEN
1627	qi(i, k) = 0.
1628	ELSE
1629	IF (t(i,k)<243.15) THEN
1630	qi(i, k) = clw(i, k)
1631	ELSE
1632	fracg = (263.15-t(i,k))/20
1633	qi(i, k) = clw(i, k)*fracg
1634	END IF
1635	END IF
1636	!CR: fin test
1637	IF (t(i,k)<263.15) THEN
1638	!CR: on commente les calculs d'Arnaud car division par zero
1639	! nouveau calcul propose par JYG
1640	! alv=lv0-clmcpv*(t(i,k)-273.15)
1641	! alf=lf0-clmci*(t(i,k)-273.15)
1642	! tg=tp(i,k)
1643	! tc=tp(i,k)-273.15
1644	! denom=243.5+tc
1645	! do j=1,3
1646	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1647	! il faudra que esi vienne en argument de la convection
1648	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1649	! tbis=t(i,k)+(tp(i,k)-tg)
1650	! esi=exp(23.33086-(6111.72784/tbis) + &
1651	! 0.15215*log(tbis))
1652	! qsat_new=epsesi/(p(i,k)-esi(1.-eps))
1653	! snew=cpd(1.-qnk(i))+clqnk(i)+alvalvqsat_new/ &
1654	! (rrvtbistbis)
1655	! snew=1./snew
1656	! print*,esi,qsat_new,snew,'esi,qsat,snew'
1657	! tp(i,k)=tg+(alfqi(i,k)+alvqg(1.-(esi/es)))snew
1658	! print*,k,tp(i,k),qnk(i),'avec glace'
1659	! print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
1660	! enddo
1661
1662	alv = lv0 - clmcpv*(t(i,k)-273.15)
1663	alf = lf0 + clmci*(t(i,k)-273.15)
1664	als = alf + alv
1665	tg = tp(i, k)
1666	tp(i, k) = t(i, k)
1667	DO j = 1, 3
1668	esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
1669	qsat_new = epsesi/(p(i,k)-esi(1.-eps))
1670	snew = cpd(1.-qnk(i)) + clqnk(i) + alvalsqsat_new/ &
1671	(rrvtp(i,k)tp(i,k))
1672	snew = 1./snew
1673	! c print*,esi,qsat_new,snew,'esi,qsat,snew'
1674	tp(i, k) = tp(i, k) + &
1675	((cpd(1.-qnk(i))+clqnk(i))*(tg-tp(i,k)) + &
1676	alv(qg-qsat_new)+alfqi(i,k))*snew
1677	! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
1678	! 'k,tp,q,qt,qi avec glace'
1679	END DO
1680
1681	!CR:reprise du code AJ
1682	clw(i, k) = qnk(i) - qsat_new
1683	clw(i, k) = max(0.0, clw(i,k))
1684	tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
1685	! print*,tvp(i,k),'tvp'
1686	END IF
1687	IF (clw(i,k)<1.E-11) THEN
1688	tp(i, k) = tv(i, k)
1689	tvp(i, k) = tv(i, k)
1690	END IF
1691	END IF ! (cvflag_ice)
1692	!jyg<
1693	END IF ! (k>=(icbs(i)+1))
1694	!>jyg
1695	END DO
1696	END DO
1697
1698	!----------------------------------------------------------------------------
1699	!
1700	ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE (icvflag_Tpa == 0)
1701	!
1702	!----------------------------------------------------------------------------
1703	!
1704	! =====================================================================
1705	! --- SET THE PRECIPITATION EFFICIENCIES
1706	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
1707	! =====================================================================
1708	!
1709	IF (flag_epkeorig/=1) THEN
1710	DO k = 1, nl ! convect3
1711	DO i = 1, ncum
1712	!jyg<
1713	IF(k>=icb(i)) THEN
1714	!>jyg
1715	pden = ptcrit - pbcrit
1716	ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1717	ep(i, k) = max(ep(i,k), 0.0)
1718	ep(i, k) = min(ep(i,k), epmax)
1719	!! sigp(i, k) = spfac ! jyg
1720	ENDIF ! (k>=icb(i))
1721	END DO
1722	END DO
1723	ELSE
1724	DO k = 1, nl
1725	DO i = 1, ncum
1726	IF(k>=icb(i)) THEN
1727	!! IF (k>=(nk(i)+1)) THEN
1728	!>jyg
1729	tca = tp(i, k) - t0
1730	IF (tca>=0.0) THEN
1731	elacrit = elcrit
1732	ELSE
1733	elacrit = elcrit*(1.0-tca/tlcrit)
1734	END IF
1735	elacrit = max(elacrit, 0.0)
1736	ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
1737	ep(i, k) = max(ep(i,k), 0.0)
1738	ep(i, k) = min(ep(i,k), epmax)
1739	!! sigp(i, k) = spfac ! jyg
1740	END IF ! (k>=icb(i))
1741	END DO
1742	END DO
1743	END IF
1744	!
1745	! =========================================================================
1746	IF (prt_level >= 10) THEN
1747	print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
1748	(k, tp(1,k), tvp(1,k), k = 1,nl)
1749	ENDIF
1750	!
1751	! =====================================================================
1752	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1753	! --- VIRTUAL TEMPERATURE
1754	! =====================================================================
1755
1756	! dans convect3, tvp est calcule en une seule fois, et sans retirer
1757	! l'eau condensee (~> reversible CAPE)
1758
1759	! ori do 340 k=minorig+1,nl
1760	! ori do 330 i=1,ncum
1761	! ori if(k.ge.(icb(i)+1))then
1762	! ori tvp(i,k)=tvp(i,k)(1.0-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	! oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1765	! ori endif
1766	! ori 330 continue
1767	! ori 340 continue
1768
1769	! ori do 350 i=1,ncum
1770	! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1771	! ori 350 continue
1772
1773	DO i = 1, ncum ! convect3
1774	tp(i, nlp) = tp(i, nl) ! convect3
1775	END DO ! convect3
1776
1777	! =====================================================================
1778	! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1779	! =====================================================================
1780
1781	! -- this is for convect3 only:
1782
1783	! first estimate of buoyancy:
1784
1785	!jyg : k-loop outside i-loop (07042015)
1786	DO k = 1, nl
1787	DO i = 1, ncum
1788	buoy(i, k) = tvp(i, k) - tv(i, k)
1789	END DO
1790	END DO
1791
1792	! set buoyancy=buoybase for all levels below base
1793	! for safety, set buoy(icb)=buoybase
1794
1795	!jyg : k-loop outside i-loop (07042015)
1796	DO k = 1, nl
1797	DO i = 1, ncum
1798	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1799	buoy(i, k) = buoybase(i)
1800	END IF
1801	END DO
1802	END DO
1803	DO i = 1, ncum
1804	! buoy(icb(i),k)=buoybase(i)
1805	buoy(i, icb(i)) = buoybase(i)
1806	END DO
1807
1808	! -- end convect3
1809
1810	! =====================================================================
1811	! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1812	! --- LEVEL OF NEUTRAL BUOYANCY
1813	! =====================================================================
1814
1815	! -- this is for convect3 only:
1816
1817	DO i = 1, ncum
1818	inb(i) = nl - 1
1819	iposit(i) = nl
1820	END DO
1821
1822
1823	! -- iposit(i) = first level, above icb, with positive buoyancy
1824	DO k = 1, nl - 1
1825	DO i = 1, ncum
1826	IF (k>=icb(i) .AND. buoy(i,k)>0.) THEN
1827	iposit(i) = min(iposit(i), k)
1828	END IF
1829	END DO
1830	END DO
1831
1832	DO i = 1, ncum
1833	IF (iposit(i)==nl) THEN
1834	iposit(i) = icb(i)
1835	END IF
1836	END DO
1837
1838	DO k = 1, nl - 1
1839	DO i = 1, ncum
1840	IF ((k>=iposit(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1841	inb(i) = min(inb(i), k)
1842	END IF
1843	END DO
1844	END DO
1845
1846	!CR fix computation of inb
1847	!keep flag or modify in all cases?
1848	IF (iflag_mix_adiab.eq.1) THEN
1849	DO i = 1, ncum
1850	cape(i)=0.
1851	inb(i)=icb(i)+1
1852	ENDDO
1853
1854	DO k = 2, nl
1855	DO i = 1, ncum
1856	IF ((k>=iposit(i))) THEN
1857	deltap = min(plcl(i), ph(i,k-1)) - min(plcl(i), ph(i,k))
1858	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
1859	IF (cape(i).gt.0.) THEN
1860	inb(i) = max(inb(i), k)
1861	END IF
1862	ENDIF
1863	ENDDO
1864	ENDDO
1865
1866	! DO i = 1, ncum
1867	! print*,"inb",inb(i)
1868	! ENDDO
1869
1870	endif
1871
1872	! -- end convect3
1873
1874	! ori do 510 i=1,ncum
1875	! ori cape(i)=0.0
1876	! ori capem(i)=0.0
1877	! ori inb(i)=icb(i)+1
1878	! ori inb1(i)=inb(i)
1879	! ori 510 continue
1880
1881	! Originial Code
1882
1883	! do 530 k=minorig+1,nl-1
1884	! do 520 i=1,ncum
1885	! if(k.ge.(icb(i)+1))then
1886	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1887	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1888	! cape(i)=cape(i)+by
1889	! if(by.ge.0.0)inb1(i)=k+1
1890	! if(cape(i).gt.0.0)then
1891	! inb(i)=k+1
1892	! capem(i)=cape(i)
1893	! endif
1894	! endif
1895	!520 continue
1896	!530 continue
1897	! do 540 i=1,ncum
1898	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1899	! cape(i)=capem(i)+byp
1900	! defrac=capem(i)-cape(i)
1901	! defrac=max(defrac,0.001)
1902	! frac(i)=-cape(i)/defrac
1903	! frac(i)=min(frac(i),1.0)
1904	! frac(i)=max(frac(i),0.0)
1905	!540 continue
1906
1907	! K Emanuel fix
1908
1909	! call zilch(byp,ncum)
1910	! do 530 k=minorig+1,nl-1
1911	! do 520 i=1,ncum
1912	! if(k.ge.(icb(i)+1))then
1913	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1914	! cape(i)=cape(i)+by
1915	! if(by.ge.0.0)inb1(i)=k+1
1916	! if(cape(i).gt.0.0)then
1917	! inb(i)=k+1
1918	! capem(i)=cape(i)
1919	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1920	! endif
1921	! endif
1922	!520 continue
1923	!530 continue
1924	! do 540 i=1,ncum
1925	! inb(i)=max(inb(i),inb1(i))
1926	! cape(i)=capem(i)+byp(i)
1927	! defrac=capem(i)-cape(i)
1928	! defrac=max(defrac,0.001)
1929	! frac(i)=-cape(i)/defrac
1930	! frac(i)=min(frac(i),1.0)
1931	! frac(i)=max(frac(i),0.0)
1932	!540 continue
1933
1934	! J Teixeira fix
1935
1936	! ori call zilch(byp,ncum)
1937	! ori do 515 i=1,ncum
1938	! ori lcape(i)=.true.
1939	! ori 515 continue
1940	! ori do 530 k=minorig+1,nl-1
1941	! ori do 520 i=1,ncum
1942	! ori if(cape(i).lt.0.0)lcape(i)=.false.
1943	! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1944	! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1945	! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1946	! ori cape(i)=cape(i)+by
1947	! ori if(by.ge.0.0)inb1(i)=k+1
1948	! ori if(cape(i).gt.0.0)then
1949	! ori inb(i)=k+1
1950	! ori capem(i)=cape(i)
1951	! ori endif
1952	! ori endif
1953	! ori 520 continue
1954	! ori 530 continue
1955	! ori do 540 i=1,ncum
1956	! ori cape(i)=capem(i)+byp(i)
1957	! ori defrac=capem(i)-cape(i)
1958	! ori defrac=max(defrac,0.001)
1959	! ori frac(i)=-cape(i)/defrac
1960	! ori frac(i)=min(frac(i),1.0)
1961	! ori frac(i)=max(frac(i),0.0)
1962	! ori 540 continue
1963
1964	! --------------------------------------------------------------------
1965	! Prevent convection when top is too hot
1966	! --------------------------------------------------------------------
1967	DO i = 1,ncum
1968	IF (t(i,inb(i)) > T_top_max) iflag(i) = 10
1969	ENDDO
1970
1971	! =====================================================================
1972	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1973	! =====================================================================
1974
1975	DO k = 1, nl
1976	DO i = 1, ncum
1977	hp(i, k) = h(i, k)
1978	END DO
1979	END DO
1980
1981	!jyg : cvflag_ice test outside the loops (07042015)
1982	!
1983	IF (cvflag_ice) THEN
1984	!
1985	IF (cvflag_prec_eject) THEN
1986	!! DO k = minorig + 1, nl
1987	!! DO i = 1, ncum
1988	!! IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1989	!! frac_s(i,k) = qi(i,k)/max(clw(i,k),smallestreal)
1990	!! frac_s(i,k) = 1. - (qpl(i,k)+(1.-frac_s(i,k))*qcld(i,k))/max(clw(i,k),smallestreal)
1991	!! END IF
1992	!! END DO
1993	!! END DO
1994	ELSE ! (cvflag_prec_eject)
1995	DO k = minorig + 1, nl
1996	DO i = 1, ncum
1997	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1998	!jyg< frac computation moved to beginning of cv3_undilute2.
1999	! kept here for compatibility test with CMip6 version
2000	frac_s(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
2001	frac_s(i, k) = min(max(frac_s(i,k),0.0), 1.0)
2002	END IF
2003	END DO
2004	END DO
2005	ENDIF ! (cvflag_prec_eject) ELSE
2006	DO k = minorig + 1, nl
2007	DO i = 1, ncum
2008	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
2009	!! hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac_s(i,k)lf(i,k))* & !!jygprl
2010	!! ep(i, k)*clw(i, k) !!jygprl
2011	hp(i, k) = hla(i,k-1) + (lv(i,k)+(cpd-cpv)t(i,k)+frac_s(i,k)lf(i,k))* & !!jygprl
2012	ep(i, k)*clw(i, k) !!jygprl
2013	END IF
2014	END DO
2015	END DO
2016	!
2017	ELSE ! (cvflag_ice)
2018	!
2019	DO k = minorig + 1, nl
2020	DO i = 1, ncum
2021	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
2022	!jyg< (energy conservation tests)
2023	!! hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)tp(i,k))ep(i, k)*clw(i, k)
2024	!! hp(i, k) = ( hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k) ) / &
2025	!! (1. - ep(i,k)*clw(i,k))
2026	!! hp(i, k) = ( hnk(i) + (lv(i,k)+(cpd-cl)t(i,k))ep(i, k)*clw(i, k) ) / &
2027	!! (1. - ep(i,k)*clw(i,k))
2028	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k)
2029	END IF
2030	END DO
2031	END DO
2032	!
2033	END IF ! (cvflag_ice)
2034
2035	RETURN
2036	END SUBROUTINE cv3_undilute2
2037
2038	SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
2039	pbase, p, ph, tv, buoy, &
2040	sig, w0, cape, m, iflag)
2041	USE cv3param_mod_h
2042	USE cvthermo_mod_h
2043	USE cvflag_mod_h
2044	IMPLICIT NONE
2045
2046	! ===================================================================
2047	! --- CLOSURE OF CONVECT3
2048	!
2049	! vectorization: S. Bony
2050	! ===================================================================
2051
2052
2053	!input:
2054	INTEGER ncum, nd, nloc
2055	INTEGER icb(nloc), inb(nloc)
2056	REAL pbase(nloc)
2057	REAL p(nloc, nd), ph(nloc, nd+1)
2058	REAL tv(nloc, nd), buoy(nloc, nd)
2059
2060	!input/output:
2061	REAL sig(nloc, nd), w0(nloc, nd)
2062	INTEGER iflag(nloc)
2063
2064	!output:
2065	REAL cape(nloc)
2066	REAL m(nloc, nd)
2067
2068	!local variables:
2069	INTEGER i, j, k, icbmax
2070	REAL deltap, fac, w, amu
2071	REAL dtmin(nloc, nd), sigold(nloc, nd)
2072	REAL cbmflast(nloc)
2073
2074
2075	! -------------------------------------------------------
2076	! -- Initialization
2077	! -------------------------------------------------------
2078
2079	DO k = 1, nl
2080	DO i = 1, ncum
2081	m(i, k) = 0.0
2082	END DO
2083	END DO
2084
2085	! -------------------------------------------------------
2086	! -- Reset sig(i) and w0(i) for i>inb and i<icb
2087	! -------------------------------------------------------
2088
2089	! update sig and w0 above LNB:
2090
2091	DO k = 1, nl - 1
2092	DO i = 1, ncum
2093	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
2094	sig(i, k) = beta*sig(i, k) + &
2095	2.alphabuoy(i, inb(i))*abs(buoy(i,inb(i)))
2096	sig(i, k) = amax1(sig(i,k), 0.0)
2097	w0(i, k) = beta*w0(i, k)
2098	END IF
2099	END DO
2100	END DO
2101
2102	! compute icbmax:
2103
2104	icbmax = 2
2105	DO i = 1, ncum
2106	icbmax = max(icbmax, icb(i))
2107	END DO
2108
2109	! update sig and w0 below cloud base:
2110
2111	DO k = 1, icbmax
2112	DO i = 1, ncum
2113	IF (k<=icb(i)) THEN
2114	sig(i, k) = beta*sig(i, k) - &
2115	2.alphabuoy(i, icb(i))*buoy(i, icb(i))
2116	sig(i, k) = max(sig(i,k), 0.0)
2117	w0(i, k) = beta*w0(i, k)
2118	END IF
2119	END DO
2120	END DO
2121
2122	!! if(inb.lt.(nl-1))then
2123	!! do 85 i=inb+1,nl-1
2124	!! sig(i)=betasig(i)+2.alphabuoy(inb)
2125	!! 1 abs(buoy(inb))
2126	!! sig(i)=max(sig(i),0.0)
2127	!! w0(i)=beta*w0(i)
2128	!! 85 continue
2129	!! end if
2130
2131	!! do 87 i=1,icb
2132	!! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
2133	!! sig(i)=max(sig(i),0.0)
2134	!! w0(i)=beta*w0(i)
2135	!! 87 continue
2136
2137	! -------------------------------------------------------------
2138	! -- Reset fractional areas of updrafts and w0 at initial time
2139	! -- and after 10 time steps of no convection
2140	! -------------------------------------------------------------
2141
2142	DO k = 1, nl - 1
2143	DO i = 1, ncum
2144	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
2145	sig(i, k) = 0.0
2146	w0(i, k) = 0.0
2147	END IF
2148	END DO
2149	END DO
2150
2151	! -------------------------------------------------------------
2152	! -- Calculate convective available potential energy (cape),
2153	! -- vertical velocity (w), fractional area covered by
2154	! -- undilute updraft (sig), and updraft mass flux (m)
2155	! -------------------------------------------------------------
2156
2157	DO i = 1, ncum
2158	cape(i) = 0.0
2159	END DO
2160
2161	! compute dtmin (minimum buoyancy between ICB and given level k):
2162
2163	DO i = 1, ncum
2164	DO k = 1, nl
2165	dtmin(i, k) = 100.0
2166	END DO
2167	END DO
2168
2169	DO i = 1, ncum
2170	DO k = 1, nl
2171	DO j = minorig, nl
2172	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
2173	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
2174	END IF
2175	END DO
2176	END DO
2177	END DO
2178
2179	! the interval on which cape is computed starts at pbase :
2180
2181	DO k = 1, nl
2182	DO i = 1, ncum
2183
2184	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
2185
2186	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
2187	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
2188	cape(i) = amax1(0.0, cape(i))
2189	sigold(i, k) = sig(i, k)
2190
2191	! dtmin(i,k)=100.0
2192	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
2193	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
2194	! 97 continue
2195
2196	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
2197	sig(i, k) = max(sig(i,k), 0.0)
2198	sig(i, k) = amin1(sig(i,k), 0.01)
2199	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
2200	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
2201	amu = 0.5(sig(i,k)+sigold(i,k))w
2202	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
2203	w0(i, k) = w
2204	END IF
2205
2206	END DO
2207	END DO
2208
2209	DO i = 1, ncum
2210	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
2211	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))
2212	sig(i, icb(i)) = sig(i, icb(i)+1)
2213	sig(i, icb(i)-1) = sig(i, icb(i))
2214	END DO
2215
2216	! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
2217	! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
2218	! ccc the final mass flux (cbmflast) is greater than the target mass flux
2219	! ccc (cbmf) ??).
2220	! cc
2221	! c do i = 1,ncum
2222	! c cbmflast(i) = 0.
2223	! c enddo
2224	! cc
2225	! c do k= 1,nl
2226	! c do i = 1,ncum
2227	! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
2228	! c cbmflast(i) = cbmflast(i)+M(i,k)
2229	! c ENDIF
2230	! c enddo
2231	! c enddo
2232	! cc
2233	! c do i = 1,ncum
2234	! c IF (cbmflast(i) .lt. 1.e-6) THEN
2235	! c iflag(i) = 3
2236	! c ENDIF
2237	! c enddo
2238	! cc
2239	! c do k= 1,nl
2240	! c do i = 1,ncum
2241	! c IF (iflag(i) .ge. 3) THEN
2242	! c M(i,k) = 0.
2243	! c sig(i,k) = 0.
2244	! c w0(i,k) = 0.
2245	! c ENDIF
2246	! c enddo
2247	! c enddo
2248	! cc
2249	!! cape=0.0
2250	!! do 98 i=icb+1,inb
2251	!! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
2252	!! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
2253	!! dcape=rrdbuoy(i-1)deltap/p(i-1)
2254	!! dlnp=deltap/p(i-1)
2255	!! cape=max(0.0,cape)
2256	!! sigold=sig(i)
2257
2258	!! dtmin=100.0
2259	!! do 97 j=icb,i-1
2260	!! dtmin=amin1(dtmin,buoy(j))
2261	!! 97 continue
2262
2263	!! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
2264	!! sig(i)=max(sig(i),0.0)
2265	!! sig(i)=amin1(sig(i),0.01)
2266	!! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
2267	!! w=(1.-beta)facsqrt(cape)+beta*w0(i)
2268	!! amu=0.5(sig(i)+sigold)w
2269	!! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
2270	!! w0(i)=w
2271	!! 98 continue
2272	!! w0(icb)=0.5*w0(icb+1)
2273	!! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
2274	!! sig(icb)=sig(icb+1)
2275	!! sig(icb-1)=sig(icb)
2276
2277	RETURN
2278	END SUBROUTINE cv3_closure
2279
2280	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
2281	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
2282	unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
2283	ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
2284	USE cvflag_mod_h
2285	USE cvthermo_mod_h
2286	USE cv3param_mod_h
2287	IMPLICIT NONE
2288
2289	! ---------------------------------------------------------------------
2290	! a faire:
2291	! - vectorisation de la partie normalisation des flux (do 789...)
2292	! ---------------------------------------------------------------------
2293
2294	!inputs:
2295	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2296	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
2297	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
2298	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
2299	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2300	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
2301	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
2302	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
2303	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, h, hp
2304	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf, frac
2305	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp, ep, clw
2306	REAL, DIMENSION (nloc, na), INTENT (IN) :: m ! input of convect3
2307
2308	!outputs:
2309	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: ment, qent
2310	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
2311	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: sij, elij
2312	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
2313	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: ments, qents
2314	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
2315
2316	!local variables:
2317	INTEGER i, j, k, il, im, jm
2318	INTEGER num1, num2
2319	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
2320	REAL alt, smid, sjmin, sjmax, delp, delm
2321	REAL asij(nloc), smax(nloc), scrit(nloc)
2322	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
2323	REAL sigij(nloc, nd, nd)
2324	REAL wgh
2325	REAL zm(nloc, na)
2326	LOGICAL lwork(nloc)
2327
2328	! =====================================================================
2329	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
2330	! =====================================================================
2331
2332	! ori do 360 i=1,ncum*nlp
2333	DO j = 1, nl
2334	DO i = 1, ncum
2335	nent(i, j) = 0
2336	! in convect3, m is computed in cv3_closure
2337	! ori m(i,1)=0.0
2338	END DO
2339	END DO
2340
2341	! ori do 400 k=1,nlp
2342	! ori do 390 j=1,nlp
2343	DO j = 1, nl
2344	DO k = 1, nl
2345	DO i = 1, ncum
2346	qent(i, k, j) = rr(i, j)
2347	uent(i, k, j) = u(i, j)
2348	vent(i, k, j) = v(i, j)
2349	elij(i, k, j) = 0.0
2350	!ym ment(i,k,j)=0.0
2351	!ym sij(i,k,j)=0.0
2352	END DO
2353	END DO
2354	END DO
2355
2356	!ym
2357	ment(1:ncum, 1:nd, 1:nd) = 0.0
2358	sij(1:ncum, 1:nd, 1:nd) = 0.0
2359
2360	!AC! do k=1,ntra
2361	!AC! do j=1,nd ! instead nlp
2362	!AC! do i=1,nd ! instead nlp
2363	!AC! do il=1,ncum
2364	!AC! traent(il,i,j,k)=tra(il,j,k)
2365	!AC! enddo
2366	!AC! enddo
2367	!AC! enddo
2368	!AC! enddo
2369	zm(:, :) = 0.
2370
2371	! =====================================================================
2372	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
2373	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
2374	! --- FRACTION (sij)
2375	! =====================================================================
2376
2377	DO i = minorig + 1, nl
2378
2379	DO j = minorig, nl
2380	DO il = 1, ncum
2381	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
2382
2383	rti = qnk(il) - ep(il, i)*clw(il, i)
2384	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
2385
2386
2387	IF (cvflag_ice) THEN
2388	! print*,cvflag_ice,'cvflag_ice dans do 700'
2389	IF (t(il,j)<=263.15) THEN
2390	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
2391	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
2392	END IF
2393	END IF
2394
2395	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
2396	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
2397	dei = denom
2398	IF (abs(dei)<0.01) dei = 0.01
2399	sij(il, i, j) = anum/dei
2400	sij(il, i, i) = 1.0
2401	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
2402	altem = altem/bf2
2403	cwat = clw(il, j)*(1.-ep(il,j))
2404	stemp = sij(il, i, j)
2405	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
2406
2407	IF (cvflag_ice) THEN
2408	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
2409	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
2410	ELSE
2411	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
2412	denom = denom + lv(il, j)*(rr(il,i)-rti)
2413	END IF
2414
2415	IF (abs(denom)<0.01) denom = 0.01
2416	sij(il, i, j) = anum/denom
2417	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
2418	altem = altem - (bf2-1.)*cwat
2419	END IF
2420	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
2421	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
2422	uent(il, i, j) = sij(il, i, j)u(il, i) + (1.-sij(il,i,j))unk(il)
2423	vent(il, i, j) = sij(il, i, j)v(il, i) + (1.-sij(il,i,j))vnk(il)
2424	!!!! do k=1,ntra
2425	!!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2426	!!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2427	!!!! end do
2428	elij(il, i, j) = altem
2429	elij(il, i, j) = max(0.0, elij(il,i,j))
2430	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
2431	nent(il, i) = nent(il, i) + 1
2432	END IF
2433	sij(il, i, j) = max(0.0, sij(il,i,j))
2434	sij(il, i, j) = amin1(1.0, sij(il,i,j))
2435	END IF ! new
2436	END DO
2437	END DO
2438
2439	!AC! do k=1,ntra
2440	!AC! do j=minorig,nl
2441	!AC! do il=1,ncum
2442	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
2443	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
2444	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2445	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2446	!AC! endif
2447	!AC! enddo
2448	!AC! enddo
2449	!AC! enddo
2450
2451
2452	! * if no air can entrain at level i assume that updraft detrains *
2453	! * at that level and calculate detrained air flux and properties *
2454
2455
2456	! @ do 170 i=icb(il),inb(il)
2457
2458	DO il = 1, ncum
2459	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
2460	! @ if(nent(il,i).eq.0)then
2461	ment(il, i, i) = m(il, i)
2462	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2463	uent(il, i, i) = unk(il)
2464	vent(il, i, i) = vnk(il)
2465	elij(il, i, i) = clw(il, i)
2466	! MAF sij(il,i,i)=1.0
2467	sij(il, i, i) = 0.0
2468	END IF
2469	END DO
2470	END DO
2471
2472	!AC! do j=1,ntra
2473	!AC! do i=minorig+1,nl
2474	!AC! do il=1,ncum
2475	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
2476	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2477	!AC! endif
2478	!AC! enddo
2479	!AC! enddo
2480	!AC! enddo
2481
2482	DO j = minorig, nl
2483	DO i = minorig, nl
2484	DO il = 1, ncum
2485	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
2486	sigij(il, i, j) = sij(il, i, j)
2487	END IF
2488	END DO
2489	END DO
2490	END DO
2491	! @ enddo
2492
2493	! @170 continue
2494
2495	! =====================================================================
2496	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
2497	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
2498	! =====================================================================
2499
2500	CALL zilch(asum, nloc*nd)
2501	CALL zilch(csum, nloc*nd)
2502	CALL zilch(csum, nloc*nd)
2503
2504	DO il = 1, ncum
2505	lwork(il) = .FALSE.
2506	END DO
2507
2508	DO i = minorig + 1, nl
2509
2510	num1 = 0
2511	DO il = 1, ncum
2512	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
2513	END DO
2514	IF (num1<=0) GO TO 789
2515
2516
2517	DO il = 1, ncum
2518	IF (i>=icb(il) .AND. i<=inb(il)) THEN
2519	lwork(il) = (nent(il,i)/=0)
2520	qp = qnk(il) - ep(il, i)*clw(il, i)
2521
2522	IF (cvflag_ice) THEN
2523
2524	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
2525	(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
2526	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
2527	(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
2528	ELSE
2529
2530	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
2531	(cpv-cpd)t(il, i)(qp-rr(il,i))
2532	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
2533	(cpd-cpv)t(il, i)(rr(il,i)-qp)
2534	END IF
2535
2536	IF (abs(denom)<0.01) denom = 0.01
2537	scrit(il) = anum/denom
2538	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
2539	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
2540	smax(il) = 0.0
2541	asij(il) = 0.0
2542	END IF
2543	END DO
2544
2545	DO j = nl, minorig, -1
2546
2547	num2 = 0
2548	DO il = 1, ncum
2549	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2550	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2551	lwork(il)) num2 = num2 + 1
2552	END DO
2553	IF (num2<=0) GO TO 175
2554
2555	DO il = 1, ncum
2556	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2557	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2558	lwork(il)) THEN
2559
2560	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
2561	wgh = 1.0
2562	IF (j>i) THEN
2563	sjmax = max(sij(il,i,j+1), smax(il))
2564	sjmax = amin1(sjmax, scrit(il))
2565	smax(il) = max(sij(il,i,j), smax(il))
2566	sjmin = max(sij(il,i,j-1), smax(il))
2567	sjmin = amin1(sjmin, scrit(il))
2568	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
2569	smid = amin1(sij(il,i,j), scrit(il))
2570	ELSE
2571	sjmax = max(sij(il,i,j+1), scrit(il))
2572	smid = max(sij(il,i,j), scrit(il))
2573	sjmin = 0.0
2574	IF (j>1) sjmin = sij(il, i, j-1)
2575	sjmin = max(sjmin, scrit(il))
2576	END IF
2577	delp = abs(sjmax-smid)
2578	delm = abs(sjmin-smid)
2579	asij(il) = asij(il) + wgh*(delp+delm)
2580	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
2581	END IF
2582	END IF
2583	END DO
2584
2585	175 END DO
2586
2587	DO il = 1, ncum
2588	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2589	asij(il) = max(1.0E-16, asij(il))
2590	asij(il) = 1.0/asij(il)
2591	asum(il, i) = 0.0
2592	bsum(il, i) = 0.0
2593	csum(il, i) = 0.0
2594	END IF
2595	END DO
2596
2597	DO j = minorig, nl
2598	DO il = 1, ncum
2599	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2600	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2601	ment(il, i, j) = ment(il, i, j)*asij(il)
2602	END IF
2603	END DO
2604	END DO
2605
2606	DO j = minorig, nl
2607	DO il = 1, ncum
2608	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2609	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2610	asum(il, i) = asum(il, i) + ment(il, i, j)
2611	ment(il, i, j) = ment(il, i, j)*sig(il, j)
2612	bsum(il, i) = bsum(il, i) + ment(il, i, j)
2613	END IF
2614	END DO
2615	END DO
2616
2617	DO il = 1, ncum
2618	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2619	bsum(il, i) = max(bsum(il,i), 1.0E-16)
2620	bsum(il, i) = 1.0/bsum(il, i)
2621	END IF
2622	END DO
2623
2624	DO j = minorig, nl
2625	DO il = 1, ncum
2626	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2627	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2628	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
2629	END IF
2630	END DO
2631	END DO
2632
2633	DO j = minorig, nl
2634	DO il = 1, ncum
2635	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2636	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2637	csum(il, i) = csum(il, i) + ment(il, i, j)
2638	END IF
2639	END DO
2640	END DO
2641
2642	DO il = 1, ncum
2643	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2644	csum(il,i)<m(il,i)) THEN
2645	nent(il, i) = 0
2646	ment(il, i, i) = m(il, i)
2647	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2648	uent(il, i, i) = unk(il)
2649	vent(il, i, i) = vnk(il)
2650	elij(il, i, i) = clw(il, i)
2651	! MAF sij(il,i,i)=1.0
2652	sij(il, i, i) = 0.0
2653	END IF
2654	END DO ! il
2655
2656	!AC! do j=1,ntra
2657	!AC! do il=1,ncum
2658	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
2659	!AC! : .and. csum(il,i).lt.m(il,i) ) then
2660	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2661	!AC! endif
2662	!AC! enddo
2663	!AC! enddo
2664	789 END DO
2665
2666	! MAF: renormalisation de MENT
2667	CALL zilch(zm, nloc*na)
2668	DO jm = 1, nl
2669	DO im = 1, nl
2670	DO il = 1, ncum
2671	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
2672	END DO
2673	END DO
2674	END DO
2675
2676	DO jm = 1, nl
2677	DO im = 1, nl
2678	DO il = 1, ncum
2679	IF (zm(il,im)/=0.) THEN
2680	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
2681	END IF
2682	END DO
2683	END DO
2684	END DO
2685
2686	DO jm = 1, nl
2687	DO im = 1, nl
2688	DO il = 1, ncum
2689	qents(il, im, jm) = qent(il, im, jm)
2690	ments(il, im, jm) = ment(il, im, jm)
2691	END DO
2692	END DO
2693	END DO
2694
2695	RETURN
2696	END SUBROUTINE cv3_mixing
2697
2698	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
2699	t, rr, rs, gz, u, v, tra, p, ph, &
2700	th, tv, lv, lf, cpn, ep, sigp, clw, frac_s, qpreca, frac_a, qta , & !!jygprl
2701	m, ment, elij, delt, plcl, coef_clos, &
2702	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
2703	faci, b, sigd, &
2704	wdtrainA, wdtrainS, wdtrainM) ! RomP
2705	USE cv3param_mod_h
2706	USE cvthermo_mod_h
2707	USE cvflag_mod_h
2708	USE print_control_mod, ONLY: prt_level, lunout
2709	IMPLICIT NONE
2710
2711	include "nuage.h"
2712
2713	!inputs:
2714	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2715	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
2716	REAL, INTENT(IN) :: delt
2717	REAL, DIMENSION (nloc), INTENT (IN) :: plcl
2718	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
2719	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz
2720	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
2721	REAL, DIMENSION (nloc, nd, ntra), INTENT(IN) :: tra
2722	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
2723	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2724	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, sigp, clw !adiab ascent shedding
2725	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_s !ice fraction in adiab ascent shedding !!jygprl
2726	REAL, DIMENSION (nloc, na), INTENT (IN) :: qpreca !adiab ascent precip !!jygprl
2727	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_a !ice fraction in adiab ascent precip !!jygprl
2728	REAL, DIMENSION (nloc, na), INTENT (IN) :: qta !adiab ascent specific total water !!jygprl
2729	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tv, lv, cpn
2730	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
2731	REAL, DIMENSION (nloc, na), INTENT (IN) :: m
2732	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: ment, elij
2733	REAL, DIMENSION (nloc), INTENT (IN) :: coef_clos
2734
2735	!input/output
2736	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag(nloc)
2737
2738	!outputs:
2739	REAL, DIMENSION (nloc, na), INTENT (OUT) :: mp, rp, up, vp
2740	REAL, DIMENSION (nloc, na), INTENT (OUT) :: water, evap, wt
2741	REAL, DIMENSION (nloc, na), INTENT (OUT) :: ice, fondue
2742	REAL, DIMENSION (nloc, na), INTENT (OUT) :: faci ! ice fraction in precipitation
2743	REAL, DIMENSION (nloc, na, ntra), INTENT (OUT) :: trap
2744	REAL, DIMENSION (nloc, na), INTENT (OUT) :: b
2745	REAL, DIMENSION (nloc), INTENT (OUT) :: sigd
2746	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
2747	! de l ascendance adiabatique et des flux melanges Pa et Pm.
2748	! Distinction des wdtrain
2749	! Pa = wdtrainA Pm = wdtrainM
2750	REAL, DIMENSION (nloc, na), INTENT (OUT) :: wdtrainA, wdtrainS, wdtrainM
2751
2752	!local variables
2753	INTEGER i, j, k, il, num1, ndp1
2754	REAL smallestreal
2755	REAL tinv, delti, coef
2756	REAL awat, afac, afac1, afac2, bfac
2757	REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
2758	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
2759	REAL ampmax, thaw
2760	REAL tevap(nloc)
2761	REAL, DIMENSION (nloc, na) :: lvcp, lfcp
2762	REAL, DIMENSION (nloc, na) :: h, hm
2763	REAL, DIMENSION (nloc, na) :: ma
2764	REAL, DIMENSION (nloc, na) :: frac ! ice fraction in precipitation source
2765	REAL, DIMENSION (nloc, na) :: fraci ! provisionnal ice fraction in precipitation
2766	REAL, DIMENSION (nloc, na) :: prec
2767	REAL wdtrain(nloc)
2768	LOGICAL lwork(nloc), mplus(nloc)
2769
2770
2771	! ------------------------------------------------------
2772	IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1)
2773
2774	smallestreal=tiny(smallestreal)
2775
2776	! =============================
2777	! --- INITIALIZE OUTPUT ARRAYS
2778	! =============================
2779	! (loops up to nl+1)
2780	mp(:,:) = 0.
2781	rp(:,:) = 0.
2782	up(:,:) = 0.
2783	vp(:,:) = 0.
2784	water(:,:) = 0.
2785	evap(:,:) = 0.
2786	wt(:,:) = 0.
2787	ice(:,:) = 0.
2788	fondue(:,:) = 0.
2789	faci(:,:) = 0.
2790	b(:,:) = 0.
2791	sigd(:) = 0.
2792	!! RomP >>>
2793	wdtrainA(:,:) = 0.
2794	wdtrainS(:,:) = 0.
2795	wdtrainM(:,:) = 0.
2796	!! RomP <<<
2797
2798	DO i = 1, nlp
2799	DO il = 1, ncum
2800	rp(il, i) = rr(il, i)
2801	up(il, i) = u(il, i)
2802	vp(il, i) = v(il, i)
2803	wt(il, i) = 0.001
2804	END DO
2805	END DO
2806
2807	! *** Set the fractionnal area sigd of precipitating downdraughts
2808	DO il = 1, ncum
2809	sigd(il) = sigdz*coef_clos(il)
2810	END DO
2811
2812	! =====================================================================
2813	! --- INITIALIZE VARIOUS ARRAYS AND PARAMETERS USED IN THE COMPUTATIONS
2814	! =====================================================================
2815	! (loops up to nl+1)
2816
2817	delti = 1./delt
2818	tinv = 1./3.
2819
2820	DO i = 1, nlp
2821	DO il = 1, ncum
2822	frac(il, i) = 0.0
2823	fraci(il, i) = 0.0
2824	prec(il, i) = 0.0
2825	lvcp(il, i) = lv(il, i)/cpn(il, i)
2826	lfcp(il, i) = lf(il, i)/cpn(il, i)
2827	END DO
2828	END DO
2829
2830	!AC! do k=1,ntra
2831	!AC! do i=1,nd
2832	!AC! do il=1,ncum
2833	!AC! trap(il,i,k)=tra(il,i,k)
2834	!AC! enddo
2835	!AC! enddo
2836	!AC! enddo
2837
2838	! * check whether ep(inb)=0, if so, skip precipitating *
2839	! * downdraft calculation *
2840
2841
2842	DO il = 1, ncum
2843	!! lwork(il)=.TRUE.
2844	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
2845	!jyg<
2846	!! lwork(il) = ep(il, inb(il)) >= 0.0001
2847	lwork(il) = ep(il, inb(il)) >= 0.0001 .AND. iflag(il) <= 2
2848	END DO
2849
2850	!
2851	! Get adiabatic ascent mass flux
2852	!
2853	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2854	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
2855	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2856	!!! Warning : this option leads to water conservation violation
2857	!!! Expert only
2858	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2859	DO il = 1, ncum
2860	ma(il, nlp) = 0.
2861	ma(il, 1) = 0.
2862	END DO
2863
2864	DO i = nl, 2, -1
2865	DO il = 1, ncum
2866	ma(il, i) = ma(il, i+1)*(1.-qta(il,i))/(1.-qta(il,i-1)) + m(il, i)
2867	END DO
2868	END DO
2869	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2870	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
2871	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2872	DO il = 1, ncum
2873	ma(il, nlp) = 0.
2874	ma(il, 1) = 0.
2875	END DO
2876
2877	DO i = nl, 2, -1
2878	DO il = 1, ncum
2879	ma(il, i) = ma(il, i+1) + m(il, i)
2880	END DO
2881	END DO
2882
2883	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
2884	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2885
2886	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2887	!
2888	! * begin downdraft loop *
2889	!
2890	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2891
2892	DO i = nl + 1, 1, -1
2893
2894	num1 = 0
2895	DO il = 1, ncum
2896	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
2897	END DO
2898	IF (num1<=0) GO TO 400
2899
2900	CALL zilch(wdtrain, ncum)
2901
2902
2903	! * integrate liquid water equation to find condensed water *
2904	! * and condensed water flux *
2905	!
2906	!
2907	! * calculate detrained precipitation *
2908
2909
2910	DO il = 1, ncum
2911	IF (i<=inb(il) .AND. lwork(il)) THEN
2912	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
2913	wdtrainS(il, i) = wdtrain(il)/grav ! Ps jyg
2914	!! wdtrainA(il, i) = wdtrain(il)/grav ! Ps RomP
2915	END IF
2916	END DO
2917
2918	IF (i>1) THEN
2919	DO j = 1, i - 1
2920	DO il = 1, ncum
2921	IF (i<=inb(il) .AND. lwork(il)) THEN
2922	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
2923	awat = max(awat, 0.0)
2924	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
2925	wdtrainM(il, i) = wdtrain(il)/grav - wdtrainS(il, i) ! Pm jyg
2926	!! wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP
2927	END IF
2928	END DO
2929	END DO
2930	END IF
2931
2932	IF (cvflag_prec_eject) THEN
2933	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2934	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
2935	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2936	!!! Warning : this option leads to water conservation violation
2937	!!! Expert only
2938	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2939	IF ( i > 1) THEN
2940	DO il = 1, ncum
2941	IF (i<=inb(il) .AND. lwork(il)) THEN
2942	wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))/(1. - qta(il, i-1)) ! Pa jygprl
2943	wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
2944	END IF
2945	END DO
2946	ENDIF ! ( i > 1)
2947	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2948	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
2949	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2950	IF ( i > 1) THEN
2951	DO il = 1, ncum
2952	IF (i<=inb(il) .AND. lwork(il)) THEN
2953	wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i)) ! Pa jygprl
2954	wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
2955	END IF
2956	END DO
2957	ENDIF ! ( i > 1)
2958
2959	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
2960	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2961	ENDIF ! (cvflag_prec_eject)
2962
2963
2964	! * find rain water and evaporation using provisional *
2965	! * estimates of rp(i)and rp(i-1) *
2966
2967
2968	IF (cvflag_ice) THEN !!jygprl
2969	IF (cvflag_prec_eject) THEN
2970	DO il = 1, ncum !!jygprl
2971	IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl
2972	frac(il, i) = (frac_a(il,i)wdtrainA(il,i)+frac_s(il,i)(wdtrainS(il,i)+wdtrainM(il,i))) / & !!jygprl
2973	max(wdtrainA(il,i)+wdtrainS(il,i)+wdtrainM(il,i),smallestreal) !!jygprl
2974	fraci(il, i) = frac(il, i) !!jygprl
2975	END IF !!jygprl
2976	END DO !!jygprl
2977	ELSE ! (cvflag_prec_eject)
2978	DO il = 1, ncum !!jygprl
2979	IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl
2980	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2981	IF (keepbug_ice_frac) THEN
2982	frac(il, i) = frac_s(il, i)
2983	! Ice fraction computed again here as a function of the temperature seen by unsaturated downdraughts
2984	! (i.e. the cold pool temperature) for compatibility with earlier versions.
2985	fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
2986	fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
2987	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2988	ELSE ! (keepbug_ice_frac)
2989	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2990	frac(il, i) = frac_s(il, i)
2991	fraci(il, i) = frac(il, i) !!jygprl
2992	ENDIF ! (keepbug_ice_frac)
2993	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2994	END IF !!jygprl
2995	END DO !!jygprl
2996	ENDIF ! (cvflag_prec_eject)
2997	END IF !!jygprl
2998
2999
3000	DO il = 1, ncum
3001	IF (i<=inb(il) .AND. lwork(il)) THEN
3002
3003	wt(il, i) = 45.0
3004
3005	IF (i<inb(il)) THEN
3006	rp(il, i) = rp(il, i+1) + &
3007	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
3008	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
3009	END IF
3010	rp(il, i) = max(rp(il,i), 0.0)
3011	rp(il, i) = amin1(rp(il,i), rs(il,i))
3012	rp(il, inb(il)) = rr(il, inb(il))
3013
3014	IF (i==1) THEN
3015	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
3016	IF (cvflag_ice) THEN
3017	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
3018	END IF
3019	ELSE
3020	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)
3021	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
3022	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
3023	rp(il, i-1) = max(rp(il,i-1), 0.0)
3024	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
3025	afac2 = p(il, i-1)(rs(il,i-1)-rp(il,i-1))/(1.0E4+2000.0p(il,i-1)*rs(il,i-1))
3026	afac = 0.5*(afac1+afac2)
3027	END IF
3028	IF (i==inb(il)) afac = 0.0
3029	afac = max(afac, 0.0)
3030	bfac = 1./(sigd(il)*wt(il,i))
3031
3032	!
3033	IF (prt_level >= 20) THEN
3034	Print*, 'cv3_unsat after provisional rp estimate: rp, afac, bfac ', &
3035	i, rp(1, i), afac,bfac
3036	ENDIF
3037	!
3038	!JYG1
3039	! cc sigt=1.0
3040	! cc if(i.ge.icb)sigt=sigp(i)
3041	! prise en compte de la variation progressive de sigt dans
3042	! les couches icb et icb-1:
3043	! pour plcl<ph(i+1), pr1=0 & pr2=1
3044	! pour plcl>ph(i), pr1=1 & pr2=0
3045	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
3046	! sur le nuage, et pr2 est la proportion sous la base du
3047	! nuage.
3048	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
3049	pr1 = max(0., min(1.,pr1))
3050	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
3051	pr2 = max(0., min(1.,pr2))
3052	sigt = sigp(il, i)*pr1 + pr2
3053	!JYG2
3054
3055	!JYG----
3056	! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3057	! c6 = water(il,i+1) + wdtrain(il)*bfac
3058	! c6 = prec(il,i+1) + wdtrain(il)*bfac
3059	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
3060	! evap(il,i)=sigtafacrevap
3061	! water(il,i)=revap*revap
3062	! prec(il,i)=revap*revap
3063	!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
3064	!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
3065	!!---end jyg---
3066
3067	! --------retour � la formulation originale d'Emanuel.
3068	IF (cvflag_ice) THEN
3069
3070	! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3071	! c6=prec(il,i+1)+bfac*wdtrain(il) &
3072	! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
3073	! if(c6.gt.0.0)then
3074	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
3075
3076	!JAM Attention: evap=sigt*E
3077	! Modification: evap devient l'�vaporation en milieu de couche
3078	! car n�cessaire dans cv3_yield
3079	! Du coup, il faut modifier pas mal d'�quations...
3080	! et l'expression de afac qui devient afac1
3081	! revap=sqrt((prec(i+1)+prec(i))/2)
3082
3083	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
3084	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
3085	! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
3086	! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
3087	! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
3088	IF (c6>b6*b6+1.E-20) THEN
3089	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
3090	ELSE
3091	revap = (-b6+sqrt(b6b6+4.c6))/2.
3092	END IF
3093	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
3094	! print*,prec(il,i),'neige'
3095
3096	!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
3097	! c evap(il,i)=sigtafacrevap
3098	! ce qui n'est pas correct. Dans cv_routines, la formulation a �t� modifiee.
3099	! Ici,l'evaporation evap est simplement calculee par l'equation de
3100	! conservation.
3101	! prec(il,i)=revap*revap
3102	! else
3103	!JYG---- Correction : si c6 <= 0, water(il,i)=0.
3104	! prec(il,i)=0.
3105	! endif
3106
3107	!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
3108	! moins [tt ce qui sort de la couche i]
3109	! print *, 'evap avec ice'
3110	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
3111	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
3112	!
3113	IF (prt_level >= 20) THEN
3114	Print*, 'cv3_unsat after evap computation: wdtrain, sigd, wt, prec(i+1),prec(i) ', &
3115	i, wdtrain(1), sigd(1), wt(1,i), prec(1,i+1),prec(1,i)
3116	ENDIF
3117	!
3118
3119	!jyg<
3120	d6 = prec(il,i)-prec(il,i+1)
3121
3122	!! d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
3123	!! e6 = bfac*wdtrain(il)
3124	!! f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
3125	!>jyg
3126	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
3127	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
3128	thaw = min(max(thaw,0.0), 1.0)
3129	!jyg<
3130	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
3131	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
3132	water(il, i) = min(prec(il,i), max(water(il,i), 0.))
3133	ice(il, i) = min(prec(il,i), max(ice(il,i), 0.))
3134
3135	!! water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
3136	!! water(il, i) = max(water(il,i), 0.)
3137	!! ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
3138	!! ice(il, i) = max(ice(il,i), 0.)
3139	!>jyg
3140	fondue(il, i) = ice(il, i)*thaw
3141	water(il, i) = water(il, i) + fondue(il, i)
3142	ice(il, i) = ice(il, i) - fondue(il, i)
3143
3144	IF (water(il,i)+ice(il,i)<1.E-30) THEN
3145	faci(il, i) = 0.
3146	ELSE
3147	faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
3148	END IF
3149
3150	! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
3151	! water(il,i)=max(water(il,i),0.)
3152	! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
3153	! ice(il,i)=max(ice(il,i),0.)
3154	! fondue(il,i)=ice(il,i)*thaw
3155	! water(il,i)=water(il,i)+fondue(il,i)
3156	! ice(il,i)=ice(il,i)-fondue(il,i)
3157
3158	! if((water(il,i)+ice(il,i)).lt.1.e-30)then
3159	! faci(il,i)=0.
3160	! else
3161	! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
3162	! endif
3163
3164	ELSE
3165	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
3166	c6 = water(il, i+1) + bfac*wdtrain(il) - &
3167	50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
3168	IF (c6>0.0) THEN
3169	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
3170	water(il, i) = revap*revap
3171	ELSE
3172	water(il, i) = 0.
3173	END IF
3174	! print *, 'evap sans ice'
3175	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
3176	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
3177
3178	END IF
3179	END IF !(i.le.inb(il) .and. lwork(il))
3180	END DO
3181	! ----------------------------------------------------------------
3182
3183	! cc
3184	! * calculate precipitating downdraft mass flux under *
3185	! * hydrostatic approximation *
3186
3187	DO il = 1, ncum
3188	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
3189
3190	tevap(il) = max(0.0, evap(il,i))
3191	delth = max(0.001, (th(il,i)-th(il,i-1)))
3192	IF (cvflag_ice) THEN
3193	IF (cvflag_grav) THEN
3194	mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
3195	(p(il,i-1)-p(il,i))/delth + &
3196	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
3197	(p(il,i-1)-p(il,i))/delth + &
3198	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
3199	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
3200	ELSE
3201	mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
3202	(p(il,i-1)-p(il,i))/delth + &
3203	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
3204	(p(il,i-1)-p(il,i))/delth + &
3205	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
3206	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
3207
3208	END IF
3209	ELSE
3210	IF (cvflag_grav) THEN
3211	mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
3212	(p(il,i-1)-p(il,i))/delth
3213	ELSE
3214	mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
3215	(p(il,i-1)-p(il,i))/delth
3216	END IF
3217
3218	END IF
3219
3220	END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
3221	IF (prt_level .GE. 20) THEN
3222	PRINT *,'cv3_unsat, mp hydrostatic ', i, mp(il,i)
3223	ENDIF
3224	END DO
3225	! ----------------------------------------------------------------
3226
3227	! * if hydrostatic assumption fails, *
3228	! * solve cubic difference equation for downdraft theta *
3229	! * and mass flux from two simultaneous differential eqns *
3230
3231	DO il = 1, ncum
3232	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
3233
3234	amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
3235	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
3236	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
3237
3238	IF (amp2>(0.1*amfac)) THEN
3239	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
3240	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
3241	(lvcp(il,i)sigd(il)th(il,i))
3242	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
3243
3244	IF (cvflag_ice) THEN
3245	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
3246	50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
3247	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
3248	ELSE
3249
3250	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
3251	50.(p(il,i-1)-p(il,i))xf*tevap(il)
3252	END IF
3253
3254	fac2 = 1.0
3255	IF (bf<0.0) fac2 = -1.0
3256	bf = abs(bf)
3257	ur = 0.25bfbf - afafaftinvtinv*tinv
3258	IF (ur>=0.0) THEN
3259	sru = sqrt(ur)
3260	fac = 1.0
3261	IF ((0.5*bf-sru)<0.0) fac = -1.0
3262	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
3263	fac(abs(0.5bf-sru))**tinv
3264	ELSE
3265	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
3266	IF (fac2<0.0) d = 3.14159 - d
3267	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
3268	END IF
3269	mp(il, i) = max(0.0, mp(il,i))
3270	IF (prt_level .GE. 20) THEN
3271	PRINT *,'cv3_unsat, mp cubic ', i, mp(il,i)
3272	ENDIF
3273
3274	IF (cvflag_ice) THEN
3275	IF (cvflag_grav) THEN
3276	!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
3277	! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
3278	! Et il faut bien revoir les facteurs 100.
3279	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i)) &
3280	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
3281	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
3282	(ph(il,i)-ph(il,i+1))) / &
3283	(mp(il,i)+sigd(il)*0.1) - &
3284	10.0(th(il,i)-th(il,i-1))t(il, i) / &
3285	(lvcp(il,i)sigd(il)th(il,i))
3286	ELSE
3287	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
3288	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
3289	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
3290	(ph(il,i)-ph(il,i+1))) / &
3291	(mp(il,i)+sigd(il)*0.1) - &
3292	10.0(th(il,i)-th(il,i-1))t(il, i) / &
3293	(lvcp(il,i)sigd(il)th(il,i))
3294	END IF
3295	ELSE
3296	IF (cvflag_grav) THEN
3297	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
3298	(mp(il,i)+sigd(il)*0.1) - &
3299	10.0(th(il,i)-th(il,i-1))t(il, i) / &
3300	(lvcp(il,i)sigd(il)th(il,i))
3301	ELSE
3302	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
3303	(mp(il,i)+sigd(il)*0.1) - &
3304	10.0(th(il,i)-th(il,i-1))t(il, i) / &
3305	(lvcp(il,i)sigd(il)th(il,i))
3306	END IF
3307	END IF
3308	b(il, i-1) = max(b(il,i-1), 0.0)
3309
3310	END IF !(amp2.gt.(0.1*amfac))
3311
3312	!jyg< This part shifted 10 lines farther
3313	!!! * limit magnitude of mp(i) to meet cfl condition *
3314	!!
3315	!! ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
3316	!! amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
3317	!! ampmax = min(ampmax, amp2)
3318	!! mp(il, i) = min(mp(il,i), ampmax)
3319	!>jyg
3320
3321	! * force mp to decrease linearly to zero *
3322	! * between cloud base and the surface *
3323
3324
3325	! c if(p(il,i).gt.p(il,icb(il)))then
3326	! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
3327	! c endif
3328	IF (ph(il,i)>0.9*plcl(il)) THEN
3329	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
3330	END IF
3331
3332	!jyg< Shifted part
3333	! * limit magnitude of mp(i) to meet cfl condition *
3334
3335	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
3336	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
3337	ampmax = min(ampmax, amp2)
3338	mp(il, i) = min(mp(il,i), ampmax)
3339	!>jyg
3340
3341	END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
3342	END DO
3343	! ----------------------------------------------------------------
3344	!
3345	IF (prt_level >= 20) THEN
3346	Print*, 'cv3_unsat after mp computation: mp, b(i), b(i-1) ', &
3347	i, mp(1, i), b(1,i), b(1,max(i-1,1))
3348	ENDIF
3349	!
3350
3351	! * find mixing ratio of precipitating downdraft *
3352
3353	DO il = 1, ncum
3354	IF (i<inb(il) .AND. lwork(il)) THEN
3355	mplus(il) = mp(il, i) > mp(il, i+1)
3356	END IF ! (i.lt.inb(il) .and. lwork(il))
3357	END DO
3358
3359	DO il = 1, ncum
3360	IF (i<inb(il) .AND. lwork(il)) THEN
3361
3362	rp(il, i) = rr(il, i)
3363
3364	IF (mplus(il)) THEN
3365
3366	IF (cvflag_grav) THEN
3367	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
3368	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
3369	ELSE
3370	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
3371	5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
3372	END IF
3373	rp(il, i) = rp(il, i)/mp(il, i)
3374	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
3375	up(il, i) = up(il, i)/mp(il, i)
3376	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
3377	vp(il, i) = vp(il, i)/mp(il, i)
3378
3379	ELSE ! if (mplus(il))
3380
3381	IF (mp(il,i+1)>1.0E-16) THEN
3382	IF (cvflag_grav) THEN
3383	rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
3384	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
3385	ELSE
3386	rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
3387	(evap(il,i+1)+evap(il,i))/mp(il, i+1)
3388	END IF
3389	up(il, i) = up(il, i+1)
3390	vp(il, i) = vp(il, i+1)
3391	END IF ! (mp(il,i+1).gt.1.0e-16)
3392	END IF ! (mplus(il)) else if (.not.mplus(il))
3393
3394	rp(il, i) = amin1(rp(il,i), rs(il,i))
3395	rp(il, i) = max(rp(il,i), 0.0)
3396
3397	END IF ! (i.lt.inb(il) .and. lwork(il))
3398	END DO
3399	! ----------------------------------------------------------------
3400
3401	! * find tracer concentrations in precipitating downdraft *
3402
3403	!AC! do j=1,ntra
3404	!AC! do il = 1,ncum
3405	!AC! if (i.lt.inb(il) .and. lwork(il)) then
3406	!AC!c
3407	!AC! if(mplus(il))then
3408	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
3409	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
3410	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
3411	!AC! else ! if (mplus(il))
3412	!AC! if(mp(il,i+1).gt.1.0e-16)then
3413	!AC! trap(il,i,j)=trap(il,i+1,j)
3414	!AC! endif
3415	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
3416	!AC!c
3417	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
3418	!AC! enddo
3419	!AC! end do
3420
3421	400 END DO
3422	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3423
3424	! * end of downdraft loop *
3425
3426	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3427
3428
3429	RETURN
3430
3431	END SUBROUTINE cv3_unsat
3432
3433	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
3434	icb, inb, delt, &
3435	t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
3436	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
3437	ep, clw, qpreca, m, tp, mp, rp, up, vp, trap, &
3438	wt, water, ice, evap, fondue, faci, b, sigd, &
3439	ment, qent, hent, iflag_mix, uent, vent, &
3440	nent, elij, traent, sig, &
3441	tv, tvp, wghti, &
3442	iflag, precip, Vprecip, Vprecipi, & ! jyg: Vprecipi
3443	ft, fr, fr_comp, fu, fv, ftra, & ! jyg
3444	cbmf, upwd, dnwd, dnwd0, ma, mip, &
3445	!! tls, tps, ! useless . jyg
3446	qcondc, wd, &
3447	ftd, fqd, qta, qtc, sigt, detrain, tau_cld_cv, coefw_cld_cv)
3448
3449	USE cv3param_mod_h
3450	USE conema3_mod_h
3451	USE print_control_mod, ONLY: lunout, prt_level
3452	USE add_phys_tend_mod, only : fl_cor_ebil
3453	USE cvflag_mod_h
3454	USE cvthermo_mod_h
3455	IMPLICIT NONE
3456
3457
3458	!inputs:
3459	INTEGER, INTENT (IN) :: iflag_mix
3460	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
3461	LOGICAL, INTENT (IN) :: ok_conserv_q
3462	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
3463	REAL, INTENT (IN) :: delt
3464	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, u, v
3465	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t_wake, rr_wake
3466	REAL, DIMENSION (nloc), INTENT (IN) :: s_wake
3467	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra
3468	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
3469	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
3470	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz, h, hp
3471	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tp
3472	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, cpn, ep, clw
3473	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
3474	REAL, DIMENSION (nloc, na), INTENT (IN) :: rp, up
3475	REAL, DIMENSION (nloc, na), INTENT (IN) :: vp
3476	REAL, DIMENSION (nloc, nd), INTENT (IN) :: wt
3477	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: trap
3478	REAL, DIMENSION (nloc, na), INTENT (IN) :: water, evap, b
3479	REAL, DIMENSION (nloc, na), INTENT (IN) :: fondue, faci, ice
3480	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: qent, uent
3481	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: hent
3482	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: vent, elij
3483	INTEGER, DIMENSION (nloc, nd), INTENT (IN) :: nent
3484	REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent
3485	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti
3486	REAL, DIMENSION (nloc, nd), INTENT (IN) :: qta
3487	REAL, DIMENSION (nloc, na),INTENT(IN) :: qpreca
3488	REAL, INTENT(IN) :: tau_cld_cv, coefw_cld_cv
3489	!
3490	!input/output:
3491	REAL, DIMENSION (nloc, na), INTENT (INOUT) :: m, mp
3492	REAL, DIMENSION (nloc, na, na), INTENT (INOUT) :: ment
3493	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
3494	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig
3495	REAL, DIMENSION (nloc), INTENT (INOUT) :: sigd
3496	!
3497	!outputs:
3498	REAL, DIMENSION (nloc), INTENT (OUT) :: precip
3499	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ft, fr, fu, fv , fr_comp
3500	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ftd, fqd
3501	REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT) :: ftra
3502	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: upwd, dnwd, ma
3503	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: dnwd0, mip
3504	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecip
3505	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecipi
3506	!! REAL tls(nloc, nd), tps(nloc, nd) ! useless . jyg
3507	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qcondc ! cld
3508	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qtc, sigt ! cld
3509	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: detrain ! Louis : pour le calcul de Klein du terme de variance qui détraine dans lenvironnement
3510	REAL, DIMENSION (nloc), INTENT (OUT) :: wd ! gust
3511	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf
3512	!
3513	!local variables:
3514	INTEGER :: i, k, il, n, j, num1
3515	REAL :: rat, delti
3516	REAL :: ax, bx, cx, dx, ex
3517	REAL :: cpinv, rdcp, dpinv
3518	REAL :: sigaq
3519	REAL, DIMENSION (nloc) :: awat
3520	REAL, DIMENSION (nloc, nd) :: lvcp, lfcp ! , mke ! unused . jyg
3521	REAL, DIMENSION (nloc) :: am, work, ad, amp1
3522	!! real up1(nloc), dn1(nloc)
3523	REAL, DIMENSION (nloc, nd, nd) :: up1, dn1
3524	!jyg<
3525	REAL, DIMENSION (nloc, nd) :: up_to, up_from
3526	REAL, DIMENSION (nloc, nd) :: dn_to, dn_from
3527	!>jyg
3528	REAL, DIMENSION (nloc) :: asum, bsum, csum, dsum
3529	REAL, DIMENSION (nloc) :: esum, fsum, gsum, hsum
3530	REAL, DIMENSION (nloc, nd) :: th_wake
3531	REAL, DIMENSION (nloc) :: alpha_qpos, alpha_qpos1
3532	REAL, DIMENSION (nloc, nd) :: qcond, nqcond, wa ! cld
3533	REAL, DIMENSION (nloc, nd) :: siga, sax, mac ! cld
3534	REAL, DIMENSION (nloc) :: sument
3535	REAL, DIMENSION (nloc, nd) :: sigment, qtment ! cld
3536	REAL, DIMENSION (nloc, nd, nd) :: qdet
3537	REAL sumdq !jyg
3538	!
3539	! -------------------------------------------------------------
3540
3541	! initialization:
3542
3543	delti = 1.0/delt
3544	! print*,'cv3_yield initialisation delt', delt
3545
3546	DO il = 1, ncum
3547	precip(il) = 0.0
3548	wd(il) = 0.0 ! gust
3549	END DO
3550
3551	! Fluxes are on a staggered grid : loops extend up to nl+1
3552	DO i = 1, nlp
3553	DO il = 1, ncum
3554	Vprecip(il, i) = 0.0
3555	Vprecipi(il, i) = 0.0 ! jyg
3556	upwd(il, i) = 0.0
3557	dnwd(il, i) = 0.0
3558	dnwd0(il, i) = 0.0
3559	mip(il, i) = 0.0
3560	END DO
3561	END DO
3562	DO i = 1, nl
3563	DO il = 1, ncum
3564	ft(il, i) = 0.0
3565	fr(il, i) = 0.0
3566	fr_comp(il,i) = 0.0
3567	fu(il, i) = 0.0
3568	fv(il, i) = 0.0
3569	ftd(il, i) = 0.0
3570	fqd(il, i) = 0.0
3571	qcondc(il, i) = 0.0 ! cld
3572	qcond(il, i) = 0.0 ! cld
3573	qtc(il, i) = 0.0 ! cld
3574	qtment(il, i) = 0.0 ! cld
3575	sigment(il, i) = 0.0 ! cld
3576	sigt(il, i) = 0.0 ! cld
3577	qdet(il,i,:) = 0.0 ! cld
3578	detrain(il, i) = 0.0 ! cld
3579	nqcond(il, i) = 0.0 ! cld
3580	END DO
3581	END DO
3582	! print*,'cv3_yield initialisation 2'
3583	!AC! do j=1,ntra
3584	!AC! do i=1,nd
3585	!AC! do il=1,ncum
3586	!AC! ftra(il,i,j)=0.0
3587	!AC! enddo
3588	!AC! enddo
3589	!AC! enddo
3590	! print*,'cv3_yield initialisation 3'
3591	DO i = 1, nl
3592	DO il = 1, ncum
3593	lvcp(il, i) = lv(il, i)/cpn(il, i)
3594	lfcp(il, i) = lf(il, i)/cpn(il, i)
3595	END DO
3596	END DO
3597
3598
3599
3600	! * calculate surface precipitation in mm/day *
3601
3602	DO il = 1, ncum
3603	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
3604	IF (cvflag_ice) THEN
3605	precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
3606	86400.1000./(rowl*grav)
3607	ELSE
3608	precip(il) = wt(il, 1)sigd(il)water(il, 1) &
3609	86400.1000./(rowl*grav)
3610	END IF
3611	END IF
3612	END DO
3613	! print*,'cv3_yield apres calcul precip'
3614
3615
3616	! === calculate vertical profile of precipitation in kg/m2/s ===
3617
3618	DO i = 1, nl
3619	DO il = 1, ncum
3620	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
3621	IF (cvflag_ice) THEN
3622	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
3623	Vprecipi(il, i) = wt(il, i)sigd(il)ice(il,i)/grav ! jyg
3624	ELSE
3625	Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
3626	Vprecipi(il, i) = 0. ! jyg
3627	END IF
3628	END IF
3629	END DO
3630	END DO
3631
3632
3633	! * Calculate downdraft velocity scale *
3634	! * NE PAS UTILISER POUR L'INSTANT *
3635
3636	!! do il=1,ncum
3637	!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
3638	!! /(sigd(il)*p(il,icb(il)))
3639	!! enddo
3640
3641
3642	! * calculate tendencies of lowest level potential temperature *
3643	! * and mixing ratio *
3644
3645	DO il = 1, ncum
3646	work(il) = 1.0/(ph(il,1)-ph(il,2))
3647	cbmf(il) = 0.0
3648	END DO
3649
3650	! - Adiabatic ascent mass flux "ma" and cloud base mass flux "cbmf"
3651	!-----------------------------------------------------------------
3652	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3653	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
3654	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3655	!!! Warning : this option leads to water conservation violation
3656	!!! Expert only
3657	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3658	DO il = 1, ncum
3659	ma(il, nlp) = 0.
3660	ma(il, 1) = 0.
3661	END DO
3662	DO k = nl, 2, -1
3663	DO il = 1, ncum
3664	ma(il, k) = ma(il, k+1)*(1.-qta(il, k))/(1.-qta(il, k-1)) + m(il, k)
3665	cbmf(il) = max(cbmf(il), ma(il,k))
3666	END DO
3667	END DO
3668	DO k = 2,nl
3669	DO il = 1, ncum
3670	IF (k <icb(il)) THEN
3671	ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
3672	ENDIF
3673	END DO
3674	END DO
3675	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3676	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
3677	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3678	!! Line kept for compatibility with earlier versions
3679	DO k = 2, nl
3680	DO il = 1, ncum
3681	IF (k>=icb(il)) THEN
3682	cbmf(il) = cbmf(il) + m(il, k)
3683	END IF
3684	END DO
3685	END DO
3686
3687	DO il = 1, ncum
3688	ma(il, nlp) = 0.
3689	ma(il, 1) = 0.
3690	END DO
3691	DO k = nl, 2, -1
3692	DO il = 1, ncum
3693	ma(il, k) = ma(il, k+1) + m(il, k)
3694	END DO
3695	END DO
3696	DO k = 2,nl
3697	DO il = 1, ncum
3698	IF (k <icb(il)) THEN
3699	ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
3700	ENDIF
3701	END DO
3702	END DO
3703
3704	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
3705	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3706	!
3707	! print*,'cv3_yield avant ft'
3708	! am is the part of cbmf taken from the first level
3709	DO il = 1, ncum
3710	am(il) = cbmf(il)*wghti(il, 1)
3711	END DO
3712
3713	DO il = 1, ncum
3714	IF (iflag(il)<=1) THEN
3715	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
3716	!JYG Correction pour conserver l'eau
3717	! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
3718	IF (cvflag_ice) THEN
3719	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
3720	lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
3721	lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
3722	(100.*(ph(il,1)-ph(il,2))) !precip
3723	ELSE
3724	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
3725	END IF
3726
3727	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
3728
3729	IF (cvflag_ice) THEN
3730	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3731	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
3732	0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
3733	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3734	ELSE
3735	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3736	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3737	END IF
3738
3739	ftd(il, 1) = ft(il, 1) ! fin precip
3740
3741	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
3742	!jyg<
3743	IF (fl_cor_ebil >= 2) THEN
3744	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3745	((t(il,2)-t(il,1))*cpn(il,2)+gz(il,2)-gz(il,1))/cpn(il,1)
3746	ELSE
3747	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3748	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
3749	ENDIF
3750	!>jyg
3751	END IF ! iflag
3752	END DO
3753
3754
3755	DO j = 2, nl
3756	IF (iflag_mix>0) THEN
3757	DO il = 1, ncum
3758	! FH WARNING a modifier :
3759	cpinv = 0.
3760	! cpinv=1.0/cpn(il,1)
3761	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3762	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
3763	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
3764	END IF ! j
3765	END DO
3766	END IF
3767	END DO
3768	! fin sature
3769
3770
3771	DO il = 1, ncum
3772	IF (iflag(il)<=1) THEN
3773	!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
3774	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
3775	sigd(il)*evap(il, 1)
3776	!!! sigd(il)0.5(evap(il,1)+evap(il,2))
3777
3778	fqd(il, 1) = fr(il, 1) !precip
3779
3780	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
3781
3782	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
3783	am(il)*(u(il,2)-u(il,1)))
3784	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
3785	am(il)*(v(il,2)-v(il,1)))
3786	END IF ! iflag
3787	END DO ! il
3788
3789
3790	!AC! do j=1,ntra
3791	!AC! do il=1,ncum
3792	!AC! if (iflag(il) .le. 1) then
3793	!AC! if (cvflag_grav) then
3794	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
3795	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3796	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3797	!AC! else
3798	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
3799	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3800	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3801	!AC! endif
3802	!AC! endif ! iflag
3803	!AC! enddo
3804	!AC! enddo
3805
3806	DO j = 2, nl
3807	DO il = 1, ncum
3808	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3809	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3810	fr_comp(il,1) = fr_comp(il,1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3811	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
3812	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
3813	END IF ! j
3814	END DO
3815	END DO
3816
3817	!AC! do k=1,ntra
3818	!AC! do j=2,nl
3819	!AC! do il=1,ncum
3820	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
3821	!AC!
3822	!AC! if (cvflag_grav) then
3823	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
3824	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3825	!AC! else
3826	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
3827	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3828	!AC! endif
3829	!AC!
3830	!AC! endif
3831	!AC! enddo
3832	!AC! enddo
3833	!AC! enddo
3834	! print*,'cv3_yield apres ft'
3835
3836	!jyg<
3837	!-----------------------------------------------------------
3838	IF (ok_optim_yield) THEN !\|
3839	!-----------------------------------------------------------
3840	!
3841	!* *
3842	!* Compute convective mass fluxes upwd and dnwd *
3843
3844	!
3845	! =================================================
3846	! upward fluxes \|
3847	! ------------------------------------------------
3848	!
3849	upwd(:,:) = 0.
3850	up_to(:,:) = 0.
3851	up_from(:,:) = 0.
3852	!
3853	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3854	IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
3855	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3856	!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
3857	!! is taken into account.
3858	!! WARNING : in the present version, taking into account the mass-flux decrease due to
3859	!! precipitation ejection leads to water conservation violation.
3860	!
3861	! - Upward mass flux of mixed draughts
3862	!---------------------------------------
3863	DO i = 2, nl
3864	DO j = 1, i-1
3865	DO il = 1, ncum
3866	IF (i<=inb(il)) THEN
3867	up_to(il,i) = up_to(il,i) + ment(il,j,i)
3868	ENDIF
3869	ENDDO
3870	ENDDO
3871	ENDDO
3872	!
3873	DO j = 3, nl
3874	DO i = 2, j-1
3875	DO il = 1, ncum
3876	IF (j<=inb(il)) THEN
3877	up_from(il,i) = up_from(il,i) + ment(il,i,j)
3878	ENDIF
3879	ENDDO
3880	ENDDO
3881	ENDDO
3882	!
3883	! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
3884	!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
3885	!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
3886	!
3887	DO i = 2, nlp
3888	DO il = 1, ncum
3889	IF (i<=inb(il)+1) THEN
3890	upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
3891	ENDIF
3892	ENDDO
3893	ENDDO
3894	!
3895	! - Total upward mass flux
3896	!---------------------------
3897	DO i = 2, nlp
3898	DO il = 1, ncum
3899	IF (i<=inb(il)+1) THEN
3900	upwd(il,i) = upwd(il,i) + ma(il,i)
3901	ENDIF
3902	ENDDO
3903	ENDDO
3904	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3905	ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
3906	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3907	!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
3908	!! is not taken into account.
3909	!
3910	! - Upward mass flux
3911	!-------------------
3912	DO i = 2, nl
3913	DO il = 1, ncum
3914	IF (i<=inb(il)) THEN
3915	up_to(il,i) = m(il,i)
3916	ENDIF
3917	ENDDO
3918	DO j = 1, i-1
3919	DO il = 1, ncum
3920	IF (i<=inb(il)) THEN
3921	up_to(il,i) = up_to(il,i) + ment(il,j,i)
3922	ENDIF
3923	ENDDO
3924	ENDDO
3925	ENDDO
3926	!
3927	DO i = 1, nl
3928	DO il = 1, ncum
3929	IF (i<=inb(il)) THEN
3930	up_from(il,i) = cbmf(il)*wghti(il,i)
3931	ENDIF
3932	ENDDO
3933	ENDDO
3934	!
3935	DO j = 3, nl
3936	DO i = 2, j-1
3937	DO il = 1, ncum
3938	IF (j<=inb(il)) THEN
3939	up_from(il,i) = up_from(il,i) + ment(il,i,j)
3940	ENDIF
3941	ENDDO
3942	ENDDO
3943	ENDDO
3944	!
3945	! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
3946	!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
3947	!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
3948	!
3949	DO i = 2, nlp
3950	DO il = 1, ncum
3951	IF (i<=inb(il)+1) THEN
3952	upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
3953	ENDIF
3954	ENDDO
3955	ENDDO
3956
3957
3958	ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
3959	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3960
3961	!
3962	! =================================================
3963	! downward fluxes \|
3964	! ------------------------------------------------
3965	dnwd(:,:) = 0.
3966	dn_to(:,:) = 0.
3967	dn_from(:,:) = 0.
3968	DO i = 1, nl
3969	DO j = i+1, nl
3970	DO il = 1, ncum
3971	IF (j<=inb(il)) THEN
3972	!! dn_to(il,i) = dn_to(il,i) + ment(il,j,i) !jyg,20220202
3973	dn_to(il,i) = dn_to(il,i) - ment(il,j,i)
3974	ENDIF
3975	ENDDO
3976	ENDDO
3977	ENDDO
3978	!
3979	DO j = 1, nl
3980	DO i = j+1, nl
3981	DO il = 1, ncum
3982	IF (i<=inb(il)) THEN
3983	!! dn_from(il,i) = dn_from(il,i) + ment(il,i,j) !jyg,20220202
3984	dn_from(il,i) = dn_from(il,i) - ment(il,i,j)
3985	ENDIF
3986	ENDDO
3987	ENDDO
3988	ENDDO
3989	!
3990	! The difference between dnwd(il,i) and dnwd(il,i+1) is due to downdrafts ending in layer
3991	!(i) (theses drafts cross interface (i+1) but not interface(i)) and to downdrafts
3992	!starting from layer (i) (theses drafts cross interface (i) but not interface(i+1)):
3993	!
3994	DO i = nl-1, 1, -1
3995	DO il = 1, ncum
3996	!! dnwd(il,i) = max(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i)) !jyg,20220202
3997	dnwd(il,i) = min(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i))
3998	ENDDO
3999	ENDDO
4000	! =================================================
4001	!
4002	!-----------------------------------------------------------
4003	ENDIF !(ok_optim_yield) !\|
4004	!-----------------------------------------------------------
4005	!>jyg
4006
4007	! * calculate tendencies of potential temperature and mixing ratio *
4008	! * at levels above the lowest level *
4009
4010	! * first find the net saturated updraft and downdraft mass fluxes *
4011	! * through each level *
4012
4013
4014	!jyg<
4015	!! DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
4016	DO i = 2, nl
4017	!>jyg
4018
4019	num1 = 0
4020	DO il = 1, ncum
4021	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
4022	END DO
4023	IF (num1<=0) GO TO 500
4024
4025	!
4026	!jyg<
4027	!-----------------------------------------------------------
4028	IF (ok_optim_yield) THEN !\|
4029	!-----------------------------------------------------------
4030	DO il = 1, ncum
4031	amp1(il) = upwd(il,i+1)
4032	ad(il) = dnwd(il,i)
4033	ENDDO
4034	!-----------------------------------------------------------
4035	ELSE !(ok_optim_yield) !\|
4036	!-----------------------------------------------------------
4037	!>jyg
4038	DO il = 1,ncum
4039	amp1(il) = 0.
4040	ad(il) = 0.
4041	ENDDO
4042
4043	DO k = 1, nl + 1
4044	DO il = 1, ncum
4045	IF (i>=icb(il)) THEN
4046	IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
4047	amp1(il) = amp1(il) + m(il, k)
4048	END IF
4049	ELSE
4050	! AMP1 is the part of cbmf taken from layers I and lower
4051	IF (k<=i) THEN
4052	amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
4053	END IF
4054	END IF
4055	END DO
4056	END DO
4057
4058	DO j = i + 1, nl + 1
4059	DO k = 1, i
4060	!yor! reverted j and k loops
4061	DO il = 1, ncum
4062	!yor! IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN ! the second condition implies the first !
4063	IF (j<=(inb(il)+1)) THEN
4064	amp1(il) = amp1(il) + ment(il, k, j)
4065	END IF
4066	END DO
4067	END DO
4068	END DO
4069
4070	DO k = 1, i - 1
4071	!jyg<
4072	!! DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
4073	DO j = i, nl
4074	!>jyg
4075	DO il = 1, ncum
4076	!yor! IF (i<=inb(il) .AND. j<=inb(il)) THEN ! the second condition implies the 1st !
4077	IF (j<=inb(il)) THEN
4078	ad(il) = ad(il) + ment(il, j, k)
4079	END IF
4080	END DO
4081	END DO
4082	END DO
4083	!
4084	!-----------------------------------------------------------
4085	ENDIF !(ok_optim_yield) !\|
4086	!-----------------------------------------------------------
4087	!
4088	!! print *,'yield, i, amp1, ad', i, amp1(1), ad(1)
4089
4090	DO il = 1, ncum
4091	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4092	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4093	cpinv = 1.0/cpn(il, i)
4094
4095	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
4096	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
4097
4098	! precip
4099	! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
4100	IF (cvflag_ice) THEN
4101	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
4102	sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
4103	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
4104	ELSE
4105	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
4106	END IF
4107
4108	rat = cpn(il, i-1)*cpinv
4109
4110	ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
4111	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
4112	IF (cvflag_ice) THEN
4113	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
4114	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
4115	0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
4116	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
4117	ELSE
4118	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
4119	(t_wake(il,i+1)-t_wake(il,i))dpinv &
4120	cpinv
4121	END IF
4122
4123	ftd(il, i) = ft(il, i)
4124	! fin precip
4125
4126	! sature
4127	!jyg<
4128	IF (fl_cor_ebil >= 2) THEN
4129	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
4130	( amp1(il)( (t(il,i+1)-t(il,i))cpn(il,i+1) + gz(il,i+1)-gz(il,i))*cpinv - &
4131	ad(il)( (t(il,i)-t(il,i-1))cpn(il,i-1) + gz(il,i)-gz(il,i-1))*cpinv)
4132	ELSE
4133	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
4134	(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
4135	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
4136	ENDIF
4137	!>jyg
4138
4139
4140	IF (iflag_mix==0) THEN
4141	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
4142	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
4143	END IF
4144	!
4145	! sb: on ne fait pas encore la correction permettant de mieux
4146	! conserver l'eau:
4147	!JYG: correction permettant de mieux conserver l'eau:
4148	! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
4149	fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
4150	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
4151	fqd(il, i) = fr(il, i) ! precip
4152
4153	fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
4154	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
4155	fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
4156	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
4157
4158
4159	fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
4160	ad(il)*(rr(il,i)-rr(il,i-1)))
4161	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
4162	ad(il)*(u(il,i)-u(il,i-1)))
4163	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
4164	ad(il)*(v(il,i)-v(il,i-1)))
4165
4166	END IF ! i
4167	END DO
4168
4169	!AC! do k=1,ntra
4170	!AC! do il=1,ncum
4171	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4172	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4173	!AC! cpinv=1.0/cpn(il,i)
4174	!AC! if (cvflag_grav) then
4175	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
4176	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
4177	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
4178	!AC! else
4179	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
4180	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
4181	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
4182	!AC! endif
4183	!AC! endif
4184	!AC! enddo
4185	!AC! enddo
4186
4187	DO k = 1, i - 1
4188
4189	DO il = 1, ncum
4190	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
4191	awat(il) = max(awat(il), 0.0)
4192	END DO
4193
4194	IF (iflag_mix/=0) THEN
4195	DO il = 1, ncum
4196	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4197	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4198	cpinv = 1.0/cpn(il, i)
4199	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
4200	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
4201	!
4202	!
4203	END IF ! i
4204	END DO
4205	END IF
4206
4207	DO il = 1, ncum
4208	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
4209	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4210	cpinv = 1.0/cpn(il, i)
4211	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
4212	(qent(il,k,i)-awat(il)-rr(il,i))
4213	fr_comp(il,i) = fr_comp(il,i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-awat(il)-rr(il,i))
4214	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
4215	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
4216
4217	! (saturated updrafts resulting from mixing) ! cld
4218	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
4219	qdet(il,k,i) = (qent(il,k,i)-awat(il)) ! cld Louis : specific humidity in detraining water
4220	qtment(il, i) = qtment(il, i) + qent(il,k,i) ! cld
4221	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4222	END IF ! i
4223	END DO
4224	END DO
4225
4226	!AC! do j=1,ntra
4227	!AC! do k=1,i-1
4228	!AC! do il=1,ncum
4229	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4230	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4231	!AC! cpinv=1.0/cpn(il,i)
4232	!AC! if (cvflag_grav) then
4233	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
4234	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4235	!AC! else
4236	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
4237	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4238	!AC! endif
4239	!AC! endif
4240	!AC! enddo
4241	!AC! enddo
4242	!AC! enddo
4243
4244	!jyg<
4245	!! DO k = i, nl + 1
4246	DO k = i, nl
4247	!>jyg
4248
4249	IF (iflag_mix/=0) THEN
4250	DO il = 1, ncum
4251	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
4252	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4253	cpinv = 1.0/cpn(il, i)
4254	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
4255	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
4256
4257
4258	END IF ! i
4259	END DO
4260	END IF
4261
4262	DO il = 1, ncum
4263	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
4264	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
4265	cpinv = 1.0/cpn(il, i)
4266
4267	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
4268	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
4269	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
4270	END IF ! i and k
4271	END DO
4272	END DO
4273
4274	!AC! do j=1,ntra
4275	!AC! do k=i,nl+1
4276	!AC! do il=1,ncum
4277	!AC! if (i.le.inb(il) .and. k.le.inb(il)
4278	!AC! $ .and. iflag(il) .le. 1) then
4279	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4280	!AC! cpinv=1.0/cpn(il,i)
4281	!AC! if (cvflag_grav) then
4282	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
4283	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4284	!AC! else
4285	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
4286	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
4287	!AC! endif
4288	!AC! endif ! i and k
4289	!AC! enddo
4290	!AC! enddo
4291	!AC! enddo
4292
4293	! sb: interface with the cloud parameterization: ! cld
4294
4295	DO k = i + 1, nl
4296	DO il = 1, ncum
4297	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
4298	! (saturated downdrafts resulting from mixing) ! cld
4299	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
4300	qdet(il,k,i) = qent(il,k,i) ! cld Louis : specific humidity in detraining water
4301	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4302	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4303	END IF ! cld
4304	END DO ! cld
4305	END DO ! cld
4306
4307	!ym BIG Warning : it seems that the k loop is missing !!!
4308	!ym Strong advice to check this
4309	!ym add a k loop temporary
4310
4311	! (particular case: no detraining level is found) ! cld
4312	! Verif merge Dynamico<<<<<<< .working
4313	DO il = 1, ncum ! cld
4314	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
4315	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
4316	!jyg< Bug correction 20180620
4317	! PROBLEM: Should not qent(il,i,i) be taken into account even if nent(il,i)/=0?
4318	!! qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4319	qdet(il,i,i) = qent(il,i,i) ! cld Louis : specific humidity in detraining water
4320	qtment(il, i) = qent(il,i,i) + qtment(il,i) ! cld
4321	!>jyg
4322	nqcond(il, i) = nqcond(il, i) + 1. ! cld
4323	END IF ! cld
4324	END DO ! cld
4325	! Verif merge Dynamico =======
4326	! Verif merge Dynamico DO k = i + 1, nl
4327	! Verif merge Dynamico DO il = 1, ncum !ym k loop added ! cld
4328	! Verif merge Dynamico IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
4329	! Verif merge Dynamico qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
4330	! Verif merge Dynamico qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
4331	! Verif merge Dynamico nqcond(il, i) = nqcond(il, i) + 1. ! cld
4332	! Verif merge Dynamico END IF ! cld
4333	! Verif merge Dynamico END DO
4334	! Verif merge Dynamico ENDDO ! cld
4335	! Verif merge Dynamico >>>>>>> .merge-right.r3413
4336
4337	DO il = 1, ncum ! cld
4338	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
4339	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
4340	qtment(il, i) = qtment(il,i)/nqcond(il, i) ! cld
4341	END IF ! cld
4342	END DO
4343
4344	!AC! do j=1,ntra
4345	!AC! do il=1,ncum
4346	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
4347	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
4348	!AC! cpinv=1.0/cpn(il,i)
4349	!AC!
4350	!AC! if (cvflag_grav) then
4351	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
4352	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
4353	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
4354	!AC! else
4355	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
4356	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
4357	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
4358	!AC! endif
4359	!AC! endif ! i
4360	!AC! enddo
4361	!AC! enddo
4362
4363
4364	500 END DO
4365
4366	!JYG<
4367	!Conservation de l'eau
4368	! sumdq = 0.
4369	! DO k = 1, nl
4370	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4371	! END DO
4372	! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4373	!JYG>
4374	! * move the detrainment at level inb down to level inb-1 *
4375	! * in such a way as to preserve the vertically *
4376	! * integrated enthalpy and water tendencies *
4377
4378	! Correction bug le 18-03-09
4379	DO il = 1, ncum
4380	IF (iflag(il)<=1) THEN
4381	ax = 0.01gravment(il, inb(il), inb(il))* &
4382	(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
4383	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
4384	ft(il, inb(il)) = ft(il, inb(il)) - ax
4385	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4386	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
4387
4388	bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
4389	(ph(il,inb(il))-ph(il,inb(il)+1))
4390	fr(il, inb(il)) = fr(il, inb(il)) - bx
4391	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4392	(ph(il,inb(il)-1)-ph(il,inb(il)))
4393
4394	cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
4395	(ph(il,inb(il))-ph(il,inb(il)+1))
4396	fu(il, inb(il)) = fu(il, inb(il)) - cx
4397	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4398	(ph(il,inb(il)-1)-ph(il,inb(il)))
4399
4400	dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
4401	(ph(il,inb(il))-ph(il,inb(il)+1))
4402	fv(il, inb(il)) = fv(il, inb(il)) - dx
4403	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
4404	(ph(il,inb(il)-1)-ph(il,inb(il)))
4405	END IF !iflag
4406	END DO
4407
4408	!JYG<
4409	!Conservation de l'eau
4410	! sumdq = 0.
4411	! DO k = 1, nl
4412	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4413	! END DO
4414	! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4415	!JYG>
4416
4417	!AC! do j=1,ntra
4418	!AC! do il=1,ncum
4419	!AC! IF (iflag(il) .le. 1) THEN
4420	!AC! IF (cvflag_grav) then
4421	!AC! ex=0.01gravment(il,inb(il),inb(il))
4422	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
4423	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
4424	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
4425	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
4426	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
4427	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
4428	!AC! else
4429	!AC! ex=0.1*ment(il,inb(il),inb(il))
4430	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
4431	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
4432	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
4433	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
4434	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
4435	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
4436	!AC! ENDIF !cvflag grav
4437	!AC! ENDIF !iflag
4438	!AC! enddo
4439	!AC! enddo
4440
4441
4442	! * homogenize tendencies below cloud base *
4443
4444
4445	DO il = 1, ncum
4446	asum(il) = 0.0
4447	bsum(il) = 0.0
4448	csum(il) = 0.0
4449	dsum(il) = 0.0
4450	esum(il) = 0.0
4451	fsum(il) = 0.0
4452	gsum(il) = 0.0
4453	hsum(il) = 0.0
4454	END DO
4455
4456	!do i=1,nl
4457	!do il=1,ncum
4458	!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
4459	!enddo
4460	!enddo
4461
4462	DO i = 1, nl
4463	DO il = 1, ncum
4464	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
4465	!jyg Saturated part : use T profile
4466	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
4467	!jyg<20140311
4468	!Correction pour conserver l eau
4469	IF (ok_conserv_q) THEN
4470	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
4471	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
4472
4473	ELSE
4474	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
4475	(ph(il,i)-ph(il,i+1))
4476	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
4477	(ph(il,i)-ph(il,i+1))
4478	ENDIF ! (ok_conserv_q)
4479	!jyg>
4480	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
4481	!jyg Unsaturated part : use T_wake profile
4482	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
4483	!jyg<20140311
4484	!Correction pour conserver l eau
4485	IF (ok_conserv_q) THEN
4486	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
4487	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
4488	ELSE
4489	fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
4490	(ph(il,i)-ph(il,i+1))
4491	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
4492	(ph(il,i)-ph(il,i+1))
4493	ENDIF ! (ok_conserv_q)
4494	!jyg>
4495	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
4496	END IF
4497	END DO
4498	END DO
4499
4500	!!!! do 700 i=1,icb(il)-1
4501	IF (ok_homo_tend) THEN
4502	DO i = 1, nl
4503	DO il = 1, ncum
4504	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
4505	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
4506	fqd(il, i) = fsum(il)/gsum(il)
4507	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
4508	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
4509	END IF
4510	END DO
4511	END DO
4512	ENDIF
4513
4514	!jyg<
4515	!Conservation de l'eau
4516	!! sumdq = 0.
4517	!! DO k = 1, nl
4518	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
4519	!! END DO
4520	!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
4521	!jyg>
4522
4523
4524	! *** Check that moisture stays positive. If not, scale tendencies
4525	! in order to ensure moisture positivity
4526	DO il = 1, ncum
4527	alpha_qpos(il) = 1.
4528	IF (iflag(il)<=1) THEN
4529	IF (fr(il,1)<=0.) THEN
4530	alpha_qpos(il) = max(alpha_qpos(il), (-deltfr(il,1))/(s_wake(il)rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
4531	END IF
4532	END IF
4533	END DO
4534	DO i = 2, nl
4535	DO il = 1, ncum
4536	IF (iflag(il)<=1) THEN
4537	IF (fr(il,i)<=0.) THEN
4538	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
4539	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
4540	END IF
4541	END IF
4542	END DO
4543	END DO
4544	DO il = 1, ncum
4545	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
4546	alpha_qpos(il) = alpha_qpos(il)*1.1
4547	END IF
4548	END DO
4549	!
4550	IF (prt_level .GE. 5) THEN
4551	print *,' CV3_YIELD : alpha_qpos ',alpha_qpos(1)
4552	ENDIF
4553	!
4554	DO il = 1, ncum
4555	IF (iflag(il)<=1) THEN
4556	sigd(il) = sigd(il)/alpha_qpos(il)
4557	precip(il) = precip(il)/alpha_qpos(il)
4558	cbmf(il) = cbmf(il)/alpha_qpos(il)
4559	END IF
4560	END DO
4561	DO i = 1, nl
4562	DO il = 1, ncum
4563	IF (iflag(il)<=1) THEN
4564	fr(il, i) = fr(il, i)/alpha_qpos(il)
4565	ft(il, i) = ft(il, i)/alpha_qpos(il)
4566	fqd(il, i) = fqd(il, i)/alpha_qpos(il)
4567	ftd(il, i) = ftd(il, i)/alpha_qpos(il)
4568	fu(il, i) = fu(il, i)/alpha_qpos(il)
4569	fv(il, i) = fv(il, i)/alpha_qpos(il)
4570	m(il, i) = m(il, i)/alpha_qpos(il)
4571	mp(il, i) = mp(il, i)/alpha_qpos(il)
4572	Vprecip(il, i) = Vprecip(il, i)/alpha_qpos(il)
4573	Vprecipi(il, i) = Vprecipi(il, i)/alpha_qpos(il) ! jyg
4574	END IF
4575	END DO
4576	END DO
4577	!jyg<
4578	!-----------------------------------------------------------
4579	IF (ok_optim_yield) THEN !\|
4580	!-----------------------------------------------------------
4581	DO i = 1, nl
4582	DO il = 1, ncum
4583	IF (iflag(il)<=1) THEN
4584	upwd(il, i) = upwd(il, i)/alpha_qpos(il)
4585	dnwd(il, i) = dnwd(il, i)/alpha_qpos(il)
4586	END IF
4587	END DO
4588	END DO
4589	!-----------------------------------------------------------
4590	ENDIF !(ok_optim_yield) !\|
4591	!-----------------------------------------------------------
4592	!>jyg
4593	DO j = 1, nl !yor! inverted i and j loops
4594	DO i = 1, nl
4595	DO il = 1, ncum
4596	IF (iflag(il)<=1) THEN
4597	ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
4598	END IF
4599	END DO
4600	END DO
4601	END DO
4602
4603	!AC! DO j = 1,ntra
4604	!AC! DO i = 1,nl
4605	!AC! DO il = 1,ncum
4606	!AC! IF (iflag(il) .le. 1) THEN
4607	!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
4608	!AC! ENDIF
4609	!AC! ENDDO
4610	!AC! ENDDO
4611	!AC! ENDDO
4612
4613
4614	! * reset counter and return *
4615
4616	! Reset counter only for points actually convective (jyg)
4617	! In order take into account the possibility of changing the compression,
4618	! reset m, sig and w0 to zero for non-convecting points.
4619	DO il = 1, ncum
4620	IF (iflag(il) < 3) THEN
4621	sig(il, nd) = 2.0
4622	ENDIF
4623	END DO
4624
4625
4626	DO i = 1, nl
4627	DO il = 1, ncum
4628	dnwd0(il, i) = -mp(il, i)
4629	END DO
4630	END DO
4631	!jyg< (loops stop at nl)
4632	!! DO i = nl + 1, nd
4633	!! DO il = 1, ncum
4634	!! dnwd0(il, i) = 0.
4635	!! END DO
4636	!! END DO
4637	!>jyg
4638
4639
4640	!jyg<
4641	!-----------------------------------------------------------
4642	IF (.NOT.ok_optim_yield) THEN !\|
4643	!-----------------------------------------------------------
4644	DO i = 1, nl
4645	DO il = 1, ncum
4646	upwd(il, i) = 0.0
4647	dnwd(il, i) = 0.0
4648	END DO
4649	END DO
4650
4651	!! DO i = 1, nl ! useless; jyg
4652	!! DO il = 1, ncum ! useless; jyg
4653	!! IF (i>=icb(il) .AND. i<=inb(il)) THEN ! useless; jyg
4654	!! upwd(il, i) = 0.0 ! useless; jyg
4655	!! dnwd(il, i) = 0.0 ! useless; jyg
4656	!! END IF ! useless; jyg
4657	!! END DO ! useless; jyg
4658	!! END DO ! useless; jyg
4659
4660	DO i = 1, nl
4661	DO k = 1, nl
4662	DO il = 1, ncum
4663	up1(il, k, i) = 0.0
4664	dn1(il, k, i) = 0.0
4665	END DO
4666	END DO
4667	END DO
4668
4669	!yor! commented original
4670	! DO i = 1, nl
4671	! DO k = i, nl
4672	! DO n = 1, i - 1
4673	! DO il = 1, ncum
4674	! IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
4675	! up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
4676	! dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
4677	! END IF
4678	! END DO
4679	! END DO
4680	! END DO
4681	! END DO
4682	!yor! replaced with
4683	DO i = 1, nl
4684	DO k = i, nl
4685	DO n = 1, i - 1
4686	DO il = 1, ncum
4687	IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! as i always <= k
4688	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
4689	END IF
4690	END DO
4691	END DO
4692	END DO
4693	END DO
4694	DO i = 1, nl
4695	DO n = 1, i - 1
4696	DO k = i, nl
4697	DO il = 1, ncum
4698	IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! i always <= k
4699	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
4700	END IF
4701	END DO
4702	END DO
4703	END DO
4704	END DO
4705	!yor! end replace
4706
4707	DO i = 1, nl
4708	DO k = 1, nl
4709	DO il = 1, ncum
4710	IF (i>=icb(il)) THEN
4711	IF (k>=i .AND. k<=(inb(il))) THEN
4712	upwd(il, i) = upwd(il, i) + m(il, k)
4713	END IF
4714	ELSE
4715	IF (k<i) THEN
4716	upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
4717	END IF
4718	END IF
4719	! c print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
4720	END DO
4721	END DO
4722	END DO
4723
4724	DO i = 2, nl
4725	DO k = i, nl
4726	DO il = 1, ncum
4727	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
4728	IF (i<=inb(il) .AND. k<=inb(il)) THEN
4729	upwd(il, i) = upwd(il, i) + up1(il, k, i)
4730	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
4731	END IF
4732	! c print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
4733	END DO
4734	END DO
4735	END DO
4736
4737
4738	!!!! DO il=1,ncum
4739	!!!! do i=icb(il),inb(il)
4740	!!!!
4741	!!!! upwd(il,i)=0.0
4742	!!!! dnwd(il,i)=0.0
4743	!!!! do k=i,inb(il)
4744	!!!! up1=0.0
4745	!!!! dn1=0.0
4746	!!!! do n=1,i-1
4747	!!!! up1=up1+ment(il,n,k)
4748	!!!! dn1=dn1-ment(il,k,n)
4749	!!!! enddo
4750	!!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
4751	!!!! dnwd(il,i)=dnwd(il,i)+dn1
4752	!!!! enddo
4753	!!!! enddo
4754	!!!!
4755	!!!! ENDDO
4756
4757	!! DO i = 1, nlp
4758	!! DO il = 1, ncum
4759	!! ma(il, i) = 0
4760	!! END DO
4761	!! END DO
4762	!!
4763	!! DO i = 1, nl
4764	!! DO j = i, nl
4765	!! DO il = 1, ncum
4766	!! ma(il, i) = ma(il, i) + m(il, j)
4767	!! END DO
4768	!! END DO
4769	!! END DO
4770
4771	!jyg< (loops stop at nl)
4772	!! DO i = nl + 1, nd
4773	!! DO il = 1, ncum
4774	!! ma(il, i) = 0.
4775	!! END DO
4776	!! END DO
4777	!>jyg
4778
4779	!! DO i = 1, nl
4780	!! DO il = 1, ncum
4781	!! IF (i<=(icb(il)-1)) THEN
4782	!! ma(il, i) = 0
4783	!! END IF
4784	!! END DO
4785	!! END DO
4786
4787	!-----------------------------------------------------------
4788	ENDIF !(.NOT.ok_optim_yield) !\|
4789	!-----------------------------------------------------------
4790	!>jyg
4791
4792	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4793	! determination de la variation de flux ascendant entre
4794	! deux niveau non dilue mip
4795	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4796
4797	DO i = 1, nl
4798	DO il = 1, ncum
4799	mip(il, i) = m(il, i)
4800	END DO
4801	END DO
4802
4803	!jyg< (loops stop at nl)
4804	!! DO i = nl + 1, nd
4805	!! DO il = 1, ncum
4806	!! mip(il, i) = 0.
4807	!! END DO
4808	!! END DO
4809	!>jyg
4810
4811
4812	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4813	! icb represente de niveau ou se trouve la
4814	! base du nuage , et inb le top du nuage
4815	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4816
4817	!! DO i = 1, nd ! unused . jyg
4818	!! DO il = 1, ncum ! unused . jyg
4819	!! mke(il, i) = upwd(il, i) + dnwd(il, i) ! unused . jyg
4820	!! END DO ! unused . jyg
4821	!! END DO ! unused . jyg
4822
4823	!! DO i = 1, nd ! unused . jyg
4824	!! DO il = 1, ncum ! unused . jyg
4825	!! rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd(1.-rr(il,i))+rr(il,i)cpv) ! unused . jyg
4826	!! tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp ! unused . jyg
4827	!! tps(il, i) = tp(il, i) ! unused . jyg
4828	!! END DO ! unused . jyg
4829	!! END DO ! unused . jyg
4830
4831
4832	! * diagnose the in-cloud mixing ratio * ! cld
4833	! * of condensed water * ! cld
4834	!! cld
4835
4836	DO i = 1, nl+1 ! cld
4837	DO il = 1, ncum ! cld
4838	mac(il, i) = 0.0 ! cld
4839	wa(il, i) = 0.0 ! cld
4840	siga(il, i) = 0.0 ! cld
4841	sax(il, i) = 0.0 ! cld
4842	END DO ! cld
4843	END DO ! cld
4844
4845	DO i = minorig, nl ! cld
4846	DO k = i + 1, nl + 1 ! cld
4847	DO il = 1, ncum ! cld
4848	IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
4849	mac(il, i) = mac(il, i) + m(il, k) ! cld
4850	END IF ! cld
4851	END DO ! cld
4852	END DO ! cld
4853	END DO ! cld
4854
4855	DO i = 1, nl ! cld
4856	DO j = 1, i ! cld
4857	DO il = 1, ncum ! cld
4858	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
4859	.AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
4860	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
4861	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
4862	END IF ! cld
4863	END DO ! cld
4864	END DO ! cld
4865	END DO ! cld
4866
4867	DO i = 1, nl ! cld
4868	DO il = 1, ncum ! cld
4869	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
4870	.AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
4871	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
4872	END IF ! cld
4873	END DO ! cld
4874	END DO
4875	! cld
4876	DO i = 1, nl
4877
4878	! 14/01/15 AJ je remets les parties manquantes cf JYG
4879	! Initialize sument to 0
4880
4881	DO il = 1,ncum
4882	sument(il) = 0.
4883	ENDDO
4884
4885	! Sum mixed mass fluxes in sument
4886
4887	DO k = 1,nl
4888	DO il = 1,ncum
4889	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
4890	sument(il) =sument(il) + abs(ment(il,k,i))
4891	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
4892	ENDIF
4893	ENDDO ! il
4894	ENDDO ! k
4895
4896	! 14/01/15 AJ delta n'a rien � faire l�...
4897	DO il = 1, ncum ! cld
4898	!! IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
4899	!! siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
4900	!! rrdtvp(il, i)/p(il, i)/100. ! cld
4901	!!
4902	!! siga(il, i) = min(siga(il,i), 1.0) ! cld
4903	sigaq = 0.
4904	IF (wa(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
4905	siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
4906	rrdtvp(il, i)/p(il, i)/100. ! cld
4907	siga(il, i) = min(siga(il,i), 1.0) ! cld
4908	sigaq = siga(il,i)*qta(il,i-1) ! cld
4909	ENDIF
4910
4911	! IM cf. FH
4912	! 14/01/15 AJ ne correspond pas � ce qui a �t� cod� par JYG et SB
4913
4914	IF (iflag_clw==0) THEN ! cld
4915	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
4916	+(1.-siga(il,i))*qcond(il, i) ! cld
4917
4918
4919	sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1)) ! cld
4920	sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i)) ! cld
4921	!! qtc(il, i) = (siga(il,i)qta(il,i-1)+sigment(il,i)qtment(il,i)) & ! cld
4922	qtc(il, i) = (sigaq+sigment(il,i)*qtment(il,i)) & ! cld
4923	/(siga(il,i)+sigment(il,i)) ! cld
4924	sigt(il,i) = sigment(il, i) + siga(il, i)
4925
4926	! qtc(il, i) = siga(il,i)qta(il,i-1)+(1.-siga(il,i))qtment(il,i) ! cld
4927	! print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)
4928
4929	ELSE IF (iflag_clw==1) THEN ! cld
4930	qcondc(il, i) = qcond(il, i) ! cld
4931	qtc(il,i) = qtment(il,i) ! cld
4932	END IF ! cld
4933
4934	END DO ! cld
4935	END DO
4936	! print*,'cv3_yield fin'
4937
4938	RETURN
4939	END SUBROUTINE cv3_yield
4940
4941	!AC! et !RomP >>>
4942	SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
4943	ment, sigij, da, phi, phi2, d1a, dam, &
4944	ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
4945	icb, inb)
4946	USE cv3param_mod_h
4947	USE cvthermo_mod_h
4948	USE cvflag_mod_h
4949	IMPLICIT NONE
4950
4951
4952	!inputs:
4953	INTEGER, INTENT (IN) :: ncum, nd, na, nloc, len
4954	INTEGER, DIMENSION (len), INTENT (IN) :: icb, inb
4955	REAL, DIMENSION (len, na, na), INTENT (IN) :: ment, sigij, elij
4956	REAL, DIMENSION (len, nd), INTENT (IN) :: clw
4957	REAL, DIMENSION (len, na), INTENT (IN) :: ep
4958	REAL, DIMENSION (len, nd+1), INTENT (IN) :: Vprecip
4959	!ouputs:
4960	REAL, DIMENSION (len, na, na), INTENT (OUT) :: phi, phi2, epmlmMm
4961	REAL, DIMENSION (len, na), INTENT (OUT) :: da, d1a, dam, eplaMm
4962	!
4963	! variables pour tracer dans precip de l'AA et des mel
4964	!local variables:
4965	INTEGER i, j, k
4966	REAL epm(nloc, na, na)
4967
4968	! variables d'Emanuel : du second indice au troisieme
4969	! ---> tab(i,k,j) -> de l origine k a l arrivee j
4970	! ment, sigij, elij
4971	! variables personnelles : du troisieme au second indice
4972	! ---> tab(i,j,k) -> de k a j
4973	! phi, phi2
4974
4975	! initialisations
4976
4977	da(:, :) = 0.
4978	d1a(:, :) = 0.
4979	dam(:, :) = 0.
4980	epm(:, :, :) = 0.
4981	eplaMm(:, :) = 0.
4982	epmlmMm(:, :, :) = 0.
4983	phi(:, :, :) = 0.
4984	phi2(:, :, :) = 0.
4985
4986	! fraction deau condensee dans les melanges convertie en precip : epm
4987	! et eau condens�e pr�cipit�e dans masse d'air satur� : l_m*dM_m/dzdz.dzdz
4988	DO j = 1, nl
4989	DO k = 1, nl
4990	DO i = 1, ncum
4991	IF (k>=icb(i) .AND. k<=inb(i) .AND. &
4992	!!jyg j.ge.k.and.j.le.inb(i)) then
4993	!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
4994	j>k .AND. j<=inb(i)) THEN
4995	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
4996	!!
4997	epm(i, j, k) = max(epm(i,j,k), 0.0)
4998	END IF
4999	END DO
5000	END DO
5001	END DO
5002
5003
5004	DO j = 1, nl
5005	DO k = 1, nl
5006	DO i = 1, ncum
5007	IF (k>=icb(i) .AND. k<=inb(i)) THEN
5008	eplaMm(i, j) = eplamm(i, j) + &
5009	ep(i, j)clw(i, j)ment(i, j, k)*(1.-sigij(i,j,k))
5010	END IF
5011	END DO
5012	END DO
5013	END DO
5014
5015	DO j = 1, nl
5016	DO k = 1, j - 1
5017	DO i = 1, ncum
5018	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
5019	epmlmMm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
5020	END IF
5021	END DO
5022	END DO
5023	END DO
5024
5025	! matrices pour calculer la tendance des concentrations dans cvltr.F90
5026	DO j = 1, nl
5027	DO k = 1, nl
5028	DO i = 1, ncum
5029	da(i, j) = da(i, j) + (1.-sigij(i,k,j))*ment(i, k, j)
5030	phi(i, j, k) = sigij(i, k, j)*ment(i, k, j)
5031	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sigij(i,k,j))
5032	IF (k<=j) THEN
5033	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, k, j)(1.-ep(i,k))*(1.-sigij(i,k,j))
5034	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
5035	END IF
5036	END DO
5037	END DO
5038	END DO
5039
5040	RETURN
5041	END SUBROUTINE cv3_tracer
5042	!AC! et !RomP <<<
5043
5044	SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
5045	iflag, &
5046	precip, sig, w0, &
5047	ft, fq, fu, fv, ftra, &
5048	Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
5049	epmax_diag, & ! epmax_cape
5050	iflag1, &
5051	precip1, sig1, w01, &
5052	ft1, fq1, fu1, fv1, ftra1, &
5053	Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1, &
5054	epmax_diag1) ! epmax_cape
5055	USE cv3param_mod_h
5056	IMPLICIT NONE
5057
5058
5059	!inputs:
5060	INTEGER len, ncum, nd, ntra, nloc
5061	INTEGER idcum(nloc)
5062	INTEGER iflag(nloc)
5063	REAL precip(nloc)
5064	REAL sig(nloc, nd), w0(nloc, nd)
5065	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
5066	REAL ftra(nloc, nd, ntra)
5067	REAL ma(nloc, nd)
5068	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
5069	REAL qcondc(nloc, nd)
5070	REAL wd(nloc), cape(nloc)
5071	REAL epmax_diag(nloc)
5072
5073	!outputs:
5074	INTEGER iflag1(len)
5075	REAL precip1(len)
5076	REAL sig1(len, nd), w01(len, nd)
5077	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
5078	REAL ftra1(len, nd, ntra)
5079	REAL ma1(len, nd)
5080	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
5081	REAL qcondc1(nloc, nd)
5082	REAL wd1(nloc), cape1(nloc)
5083	REAL epmax_diag1(len) ! epmax_cape
5084
5085	!local variables:
5086	INTEGER i, k, j
5087
5088	DO i = 1, ncum
5089	precip1(idcum(i)) = precip(i)
5090	iflag1(idcum(i)) = iflag(i)
5091	wd1(idcum(i)) = wd(i)
5092	cape1(idcum(i)) = cape(i)
5093	epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
5094	END DO
5095
5096	DO k = 1, nl
5097	DO i = 1, ncum
5098	sig1(idcum(i), k) = sig(i, k)
5099	w01(idcum(i), k) = w0(i, k)
5100	ft1(idcum(i), k) = ft(i, k)
5101	fq1(idcum(i), k) = fq(i, k)
5102	fu1(idcum(i), k) = fu(i, k)
5103	fv1(idcum(i), k) = fv(i, k)
5104	ma1(idcum(i), k) = ma(i, k)
5105	upwd1(idcum(i), k) = upwd(i, k)
5106	dnwd1(idcum(i), k) = dnwd(i, k)
5107	dnwd01(idcum(i), k) = dnwd0(i, k)
5108	qcondc1(idcum(i), k) = qcondc(i, k)
5109	END DO
5110	END DO
5111
5112	DO i = 1, ncum
5113	sig1(idcum(i), nd) = sig(i, nd)
5114	END DO
5115
5116
5117	!AC! do 2100 j=1,ntra
5118	!AC!c oct3 do 2110 k=1,nl
5119	!AC! do 2110 k=1,nd ! oct3
5120	!AC! do 2120 i=1,ncum
5121	!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
5122	!AC! 2120 continue
5123	!AC! 2110 continue
5124	!AC! 2100 continue
5125	!
5126	RETURN
5127	END SUBROUTINE cv3_uncompress
5128
5129
5130	subroutine cv3_epmax_fn_cape(nloc,ncum,nd &
5131	, ep,hp,icb,inb,clw,nk,t,h,hnk,lv,lf,frac &
5132	, pbase, p, ph, tv, buoy, sig, w0,iflag &
5133	, epmax_diag)
5134	USE cv3param_mod_h
5135	USE conema3_mod_h
5136	USE cvflag_mod_h
5137	USE cvthermo_mod_h
5138	implicit none
5139
5140	! On fait varier epmax en fn de la cape
5141	! Il faut donc recalculer ep, et hp qui a d�j� �t� calcul� et
5142	! qui en d�pend
5143	! Toutes les autres variables fn de ep sont calcul�es plus bas.
5144
5145
5146	! inputs:
5147	INTEGER, INTENT (IN) :: ncum, nd, nloc
5148	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
5149	REAL, DIMENSION (nloc), INTENT (IN) :: hnk,pbase
5150	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, lv, lf, tv, h
5151	REAL, DIMENSION (nloc, nd), INTENT (IN) :: clw, buoy,frac
5152	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig,w0
5153	INTEGER, DIMENSION (nloc), INTENT (IN) :: iflag(nloc)
5154	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
5155	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
5156	! inouts:
5157	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: ep,hp
5158	! outputs
5159	REAL, DIMENSION (nloc), INTENT (OUT) :: epmax_diag
5160
5161	! local
5162	integer i,k
5163	! real hp_bak(nloc,nd)
5164	! real ep_bak(nloc,nd)
5165	real m_loc(nloc,nd)
5166	real sig_loc(nloc,nd)
5167	real w0_loc(nloc,nd)
5168	integer iflag_loc(nloc)
5169	real cape(nloc)
5170
5171	if (coef_epmax_cape.gt.1e-12) then
5172
5173	! il faut calculer la cape: on fait un calcule simple car tant qu'on ne
5174	! connait pas ep, on ne connait pas les m�langes, ddfts etc... qui sont
5175	! necessaires au calcul de la cape dans la nouvelle physique
5176
5177	! write(,) 'cv3_routines check 4303'
5178	do i=1,ncum
5179	do k=1,nd
5180	sig_loc(i,k)=sig(i,k)
5181	w0_loc(i,k)=w0(i,k)
5182	iflag_loc(i)=iflag(i)
5183	! ep_bak(i,k)=ep(i,k)
5184	enddo ! do k=1,nd
5185	enddo !do i=1,ncum
5186
5187	! write(,) 'cv3_routines check 4311'
5188	! write(,) 'nl=',nl
5189	CALL cv3_closure(nloc, ncum, nd, icb, inb, & ! na->nd
5190	pbase, p, ph, tv, buoy, &
5191	sig_loc, w0_loc, cape, m_loc,iflag_loc)
5192
5193	! write(,) 'cv3_routines check 4316'
5194	! write(,) 'ep(1,:)=',ep(1,:)
5195	do i=1,ncum
5196	epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
5197	epmax_diag(i)=amax1(epmax_diag(i),0.0)
5198	! write(,) 'i,icb,inb,cape,epmax_diag=', &
5199	! i,icb(i),inb(i),cape(i),epmax_diag(i)
5200	do k=1,nl
5201	ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
5202	ep(i,k)=amax1(ep(i,k),0.0)
5203	ep(i,k)=amin1(ep(i,k),epmax_diag(i))
5204	enddo
5205	enddo
5206	! write(,) 'ep(1,:)=',ep(1,:)
5207
5208	!write(,) 'cv3_routines check 4326'
5209	! On recalcule hp:
5210	! do k=1,nl
5211	! do i=1,ncum
5212	! hp_bak(i,k)=hp(i,k)
5213	! enddo
5214	! enddo
5215	do k=1,nl
5216	do i=1,ncum
5217	hp(i,k)=h(i,k)
5218	enddo
5219	enddo
5220
5221	IF (cvflag_ice) THEN
5222
5223	do k=minorig+1,nl
5224	do i=1,ncum
5225	if((k.ge.icb(i)).and.(k.le.inb(i)))then
5226	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac(i,k)lf(i,k))* &
5227	ep(i, k)*clw(i, k)
5228	endif
5229	enddo
5230	enddo !do k=minorig+1,n
5231	ELSE !IF (cvflag_ice) THEN
5232
5233	DO k = minorig + 1, nl
5234	DO i = 1, ncum
5235	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
5236	hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
5237	endif
5238	enddo
5239	enddo !do k=minorig+1,n
5240
5241	ENDIF !IF (cvflag_ice) THEN
5242	!write(,) 'cv3_routines check 4345'
5243	! do i=1,ncum
5244	! do k=1,nl
5245	! if ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-1).or. &
5246	! ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-4).and. &
5247	! (ep(i,k)-ep_bak(i,k).lt.1e-4))) then
5248	! write(,) 'i,k=',i,k
5249	! write(,) 'coef_epmax_cape=',coef_epmax_cape
5250	! write(,) 'epmax_diag(i)=',epmax_diag(i)
5251	! write(,) 'ep(i,k)=',ep(i,k)
5252	! write(,) 'ep_bak(i,k)=',ep_bak(i,k)
5253	! write(,) 'hp(i,k)=',hp(i,k)
5254	! write(,) 'hp_bak(i,k)=',hp_bak(i,k)
5255	! write(,) 'h(i,k)=',h(i,k)
5256	! write(,) 'nk(i)=',nk(i)
5257	! write(,) 'h(i,nk(i))=',h(i,nk(i))
5258	! write(,) 'lv(i,k)=',lv(i,k)
5259	! write(,) 't(i,k)=',t(i,k)
5260	! write(,) 'clw(i,k)=',clw(i,k)
5261	! write(,) 'cpd,cpv=',cpd,cpv
5262	! stop
5263	! endif
5264	! enddo !do k=1,nl
5265	! enddo !do i=1,ncum
5266	endif !if (coef_epmax_cape.gt.1e-12) then
5267	!write(,) 'cv3_routines check 4367'
5268
5269	return
5270	end subroutine cv3_epmax_fn_cape
5271
5272
5273

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

Download in other formats:

Original Format