Context Navigation

cv3_routines.F90 @ 3623

Last change on this file since 3623 was 3614, checked in by jyg, 5 years ago
Bug fix: Missing a max( . . . ,0.) in the computation of the conversion of cloud water in precipitation (subroutine cv3_undilute2 in cv3_routines.F90).
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: 178.6 KB

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

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

Download in other formats:

Original Format