Context Navigation

cv3_routines.F90 @ 3479

Last change on this file since 3479 was 3435, checked in by Laurent Fairhead, 5 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

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

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

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

Download in other formats:

Original Format