Context Navigation

cv3_routines.F90 @ 2594

Last change on this file since 2594 was 2594, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2545:2589 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 147.4 KB

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

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

Download in other formats:

Original Format