Context Navigation

cv3_routines.F90 @ 5306

Last change on this file since 5306 was 4727, checked in by idelkadi, 13 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
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.1 KB

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

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

Download in other formats:

Original Format