Context Navigation

cv3_routines.F90 @ 3183

Last change on this file since 3183 was 2605, checked in by jyg, 8 years ago
Adding some initialization in cv3_routines.F90. Should help guarantee 1+1=2.
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: 147.8 KB

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

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

Download in other formats:

Original Format