Context Navigation

cv3_routines.F90 @ 3492

Last change on this file since 3492 was 3492, checked in by jyg, 5 years ago

New computations of the lifted parcel in Emanuel's
convective scheme, controlled by two new flags:

Integer flag icvflag_Tpa:

0 -> Ice fraction in the adiabatic ascents =

function of environment tempertaure.
Computation in two steps: the first is ice
free; the second adds ice (Arnaud Jam
scheme)

1 -> Ice fraction in the adiabatic ascents =

function of environment tempertaure.
Computation in one single step.

2 -> Ice fraction in the adiabatic ascents =

function of ascent temperature.
Computation in one single step.

Logical flag qsat_depends_on_qt:

When True, specific hunidity at saturation

is a function of the amount of condsensate.

Else it is equal to qsat.

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: 162.5 KB

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

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

Download in other formats:

Original Format