Context Navigation

cv3_routines.F90 @ 3496

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

Implementation of the ejection of liquid precipitation from the adiabatic ascents.
New flags:
+cvflag_prec_eject: logical

n -> old code, y -> new code

+ejectliq: real; possible values 0. & 1.

-> no liquid precipitation is ejected
-> all liquid precipitation is ejected

+ejectice: real; any value between 0. and 1.

fraction of solid precipitation ejected at each level

Note that the adiabatic ascent mass flux decrease due to precipitation ejection is not taken into account.

Attempts to do it led to water conservation violation.

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

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

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

Download in other formats:

Original Format