Context Navigation

cv3_routines.F90 @ 2686

Last change on this file since 2686 was 2671, checked in by jyg, 8 years ago
Bug fix in cv3_routine.F90 (cv3_unsat could create a precipitating downdraught when convection was off) and in wake.F90 and calwake.F90 (some array dimensions were incompatible and some initializations were missing)
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: 148.4 KB

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

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

Download in other formats:

Original Format