Context Navigation

cv3_routines.F90 @ 2518

Last change on this file since 2518 was 2508, checked in by jyg, 9 years ago

Optimization of cv3_yield (in cv3_routines.F90):
suppression of the triple vertical loops. The
optimized code is activated by:

ok_optim_yield=y

which changes numerical results.

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 2508 2016-05-07 17:12:33Z jbmadeleine $
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) = max(icb(i), 2) !convect3
586	icb1(i) = min(icb(i), nl) !convect3
587	! if icb is below LCL, start loop at ICB+1:
588	! (icbs est le premier niveau au-dessus du LCL)
589	icbs(i) = icb1(i) !convect3
590	IF (plcl(i)<p(i,icb1(i))) THEN
591	icbs(i) = min(icbs(i)+1, nl) !convect3
592	END IF
593	END DO !convect3
594
595	DO i = 1, len !convect3
596	ticb(i) = t(i, icbs(i)) !convect3
597	gzicb(i) = gz(i, icbs(i)) !convect3
598	qsicb(i) = qs(i, icbs(i)) !convect3
599	END DO !convect3
600
601
602	! Re-compute icbsmax (icbsmax2): !convect3
603	! !convect3
604	icbsmax2 = 2 !convect3
605	DO i = 1, len !convect3
606	icbsmax2 = max(icbsmax2, icbs(i)) !convect3
607	END DO !convect3
608
609	! initialization outputs:
610
611	DO k = 1, icbsmax2 ! convect3
612	DO i = 1, len ! convect3
613	tp(i, k) = 0.0 ! convect3
614	tvp(i, k) = 0.0 ! convect3
615	clw(i, k) = 0.0 ! convect3
616	END DO ! convect3
617	END DO ! convect3
618
619	! tp and tvp below cloud base:
620
621	DO k = minorig, icbsmax2 - 1
622	DO i = 1, len
623	tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
624	tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
625	END DO
626	END DO
627
628	! * Find lifted parcel quantities above cloud base *
629
630	DO i = 1, len
631	tg = ticb(i)
632	! ori qg=qs(i,icb(i))
633	qg = qsicb(i) ! convect3
634	! debug alv=lv0-clmcpv*(ticb(i)-t0)
635	alv = lv0 - clmcpv*(ticb(i)-273.15)
636
637	! First iteration.
638
639	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
640	s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
641	alvalvqg/(rrvticb(i)ticb(i)) ! convect3
642	s = 1./s
643	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
644	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
645	tg = tg + s*(ah0(i)-ahg)
646	! ori tg=max(tg,35.0)
647	! debug tc=tg-t0
648	tc = tg - 273.15
649	denom = 243.5 + tc
650	denom = max(denom, 1.0) ! convect3
651	! ori if(tc.ge.0.0)then
652	es = 6.112exp(17.67tc/denom)
653	! ori else
654	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
655	! ori endif
656	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
657	qg = epses/(p(i,icbs(i))-es(1.-eps))
658
659	! Second iteration.
660
661
662	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
663	! ori s=1./s
664	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
665	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
666	tg = tg + s*(ah0(i)-ahg)
667	! ori tg=max(tg,35.0)
668	! debug tc=tg-t0
669	tc = tg - 273.15
670	denom = 243.5 + tc
671	denom = max(denom, 1.0) ! convect3
672	! ori if(tc.ge.0.0)then
673	es = 6.112exp(17.67tc/denom)
674	! ori else
675	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
676	! ori end if
677	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
678	qg = epses/(p(i,icbs(i))-es(1.-eps))
679
680	alv = lv0 - clmcpv*(ticb(i)-273.15)
681
682	! ori c approximation here:
683	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
684	! ori & -gz(i,icb(i))-alv*qg)/cpd
685
686	! convect3: no approximation:
687	tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alvqg)/(cpd+(cl-cpd)qnk(i))
688
689	! ori clw(i,icb(i))=qnk(i)-qg
690	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
691	clw(i, icbs(i)) = qnk(i) - qg
692	clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
693
694	rg = qg/(1.-qnk(i))
695	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
696	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
697	tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
698
699	END DO
700
701	! ori do 380 k=minorig,icbsmax2
702	! ori do 370 i=1,len
703	! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
704	! ori 370 continue
705	! ori 380 continue
706
707
708	! -- The following is only for convect3:
709
710	! * icbs is the first level above the LCL:
711	! if plcl<p(icb), then icbs=icb+1
712	! if plcl>p(icb), then icbs=icb
713
714	! * the routine above computes tvp from minorig to icbs (included).
715
716	! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
717	! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
718
719	! * therefore, in the case icbs=icb, we compute tvp at level icb+1
720	! (tvp at other levels will be computed in cv3_undilute2.F)
721
722
723	DO i = 1, len
724	ticb(i) = t(i, icb(i)+1)
725	gzicb(i) = gz(i, icb(i)+1)
726	qsicb(i) = qs(i, icb(i)+1)
727	END DO
728
729	DO i = 1, len
730	tg = ticb(i)
731	qg = qsicb(i) ! convect3
732	! debug alv=lv0-clmcpv*(ticb(i)-t0)
733	alv = lv0 - clmcpv*(ticb(i)-273.15)
734
735	! First iteration.
736
737	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
738	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
739	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
740	s = 1./s
741	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
742	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
743	tg = tg + s*(ah0(i)-ahg)
744	! ori tg=max(tg,35.0)
745	! debug tc=tg-t0
746	tc = tg - 273.15
747	denom = 243.5 + tc
748	denom = max(denom, 1.0) ! convect3
749	! ori if(tc.ge.0.0)then
750	es = 6.112exp(17.67tc/denom)
751	! ori else
752	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
753	! ori endif
754	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
755	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
756
757	! Second iteration.
758
759
760	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
761	! ori s=1./s
762	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
763	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
764	tg = tg + s*(ah0(i)-ahg)
765	! ori tg=max(tg,35.0)
766	! debug tc=tg-t0
767	tc = tg - 273.15
768	denom = 243.5 + tc
769	denom = max(denom, 1.0) ! convect3
770	! ori if(tc.ge.0.0)then
771	es = 6.112exp(17.67tc/denom)
772	! ori else
773	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
774	! ori end if
775	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
776	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
777
778	alv = lv0 - clmcpv*(ticb(i)-273.15)
779
780	! ori c approximation here:
781	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
782	! ori & -gz(i,icb(i))-alv*qg)/cpd
783
784	! convect3: no approximation:
785	tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alvqg)/(cpd+(cl-cpd)qnk(i))
786
787	! ori clw(i,icb(i))=qnk(i)-qg
788	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
789	clw(i, icb(i)+1) = qnk(i) - qg
790	clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
791
792	rg = qg/(1.-qnk(i))
793	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
794	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
795	tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
796
797	END DO
798
799	RETURN
800	END SUBROUTINE cv3_undilute1
801
802	SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, &
803	pbase, buoybase, iflag, sig, w0)
804	IMPLICIT NONE
805
806	! -------------------------------------------------------------------
807	! --- TRIGGERING
808
809	! - computes the cloud base
810	! - triggering (crude in this version)
811	! - relaxation of sig and w0 when no convection
812
813	! Caution1: if no convection, we set iflag=4
814	! (it used to be 0 in convect3)
815
816	! Caution2: at this stage, tvp (and thus buoy) are know up
817	! through icb only!
818	! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
819	! -------------------------------------------------------------------
820
821	include "cv3param.h"
822
823	! input:
824	INTEGER len, nd
825	INTEGER icb(len)
826	REAL plcl(len), p(len, nd)
827	REAL th(len, nd), tv(len, nd), tvp(len, nd)
828	REAL thnk(len)
829
830	! output:
831	REAL pbase(len), buoybase(len)
832
833	! input AND output:
834	INTEGER iflag(len)
835	REAL sig(len, nd), w0(len, nd)
836
837	! local variables:
838	INTEGER i, k
839	REAL tvpbase, tvbase, tdif, ath, ath1
840
841
842	! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
843
844	DO i = 1, len
845	pbase(i) = plcl(i) + dpbase
846	tvpbase = tvp(i, icb(i)) *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
847	tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i)) /(p(i,icb(i))-p(i,icb(i)+1))
848	tvbase = tv(i, icb(i)) *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
849	tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i)) /(p(i,icb(i))-p(i,icb(i)+1))
850	buoybase(i) = tvpbase - tvbase
851	END DO
852
853
854	! * make sure that column is dry adiabatic between the surface *
855	! * and cloud base, and that lifted air is positively buoyant *
856	! * at cloud base *
857	! * if not, return to calling program after resetting *
858	! * sig(i) and w0(i) *
859
860
861	! oct3 do 200 i=1,len
862	! oct3
863	! oct3 tdif = buoybase(i)
864	! oct3 ath1 = th(i,1)
865	! oct3 ath = th(i,icb(i)-1) - dttrig
866	! oct3
867	! oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
868	! oct3 do 60 k=1,nl
869	! oct3 sig(i,k) = betasig(i,k) - 2.alphatdiftdif
870	! oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
871	! oct3 w0(i,k) = beta*w0(i,k)
872	! oct3 60 continue
873	! oct3 iflag(i)=4 ! pour version vectorisee
874	! oct3c convect3 iflag(i)=0
875	! oct3cccc return
876	! oct3 endif
877	! oct3
878	! oct3200 continue
879
880	! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
881
882	DO k = 1, nl
883	DO i = 1, len
884
885	tdif = buoybase(i)
886	ath1 = thnk(i)
887	ath = th(i, icb(i)-1) - dttrig
888
889	IF (tdif<dtcrit .OR. ath>ath1) THEN
890	sig(i, k) = betasig(i, k) - 2.alphatdiftdif
891	sig(i, k) = amax1(sig(i,k), 0.0)
892	w0(i, k) = beta*w0(i, k)
893	iflag(i) = 4 ! pour version vectorisee
894	! convect3 iflag(i)=0
895	END IF
896
897	END DO
898	END DO
899
900	! fin oct3 --
901
902	RETURN
903	END SUBROUTINE cv3_trigger
904
905	SUBROUTINE cv3_compress(len, nloc, ncum, nd, ntra, &
906	iflag1, nk1, icb1, icbs1, &
907	plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, &
908	t1, q1, qs1, u1, v1, gz1, th1, &
909	tra1, &
910	h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, &
911	sig1, w01, &
912	iflag, nk, icb, icbs, &
913	plcl, tnk, qnk, gznk, pbase, buoybase, &
914	t, q, qs, u, v, gz, th, &
915	tra, &
916	h, lv, cpn, p, ph, tv, tp, tvp, clw, &
917	sig, w0)
918	USE print_control_mod, ONLY: lunout
919	IMPLICIT NONE
920
921	include "cv3param.h"
922
923	!inputs:
924	INTEGER len, ncum, nd, ntra, nloc
925	INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
926	REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
927	REAL pbase1(len), buoybase1(len)
928	REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
929	REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
930	REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
931	REAL tvp1(len, nd), clw1(len, nd)
932	REAL th1(len, nd)
933	REAL sig1(len, nd), w01(len, nd)
934	REAL tra1(len, nd, ntra)
935
936	!outputs:
937	! en fait, on a nloc=len pour l'instant (cf cv_driver)
938	INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
939	REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
940	REAL pbase(nloc), buoybase(nloc)
941	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
942	REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
943	REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
944	REAL tvp(nloc, nd), clw(nloc, nd)
945	REAL th(nloc, nd)
946	REAL sig(nloc, nd), w0(nloc, nd)
947	REAL tra(nloc, nd, ntra)
948
949	!local variables:
950	INTEGER i, k, nn, j
951
952	CHARACTER (LEN=20) :: modname = 'cv3_compress'
953	CHARACTER (LEN=80) :: abort_message
954
955	DO k = 1, nl + 1
956	nn = 0
957	DO i = 1, len
958	IF (iflag1(i)==0) THEN
959	nn = nn + 1
960	sig(nn, k) = sig1(i, k)
961	w0(nn, k) = w01(i, k)
962	t(nn, k) = t1(i, k)
963	q(nn, k) = q1(i, k)
964	qs(nn, k) = qs1(i, k)
965	u(nn, k) = u1(i, k)
966	v(nn, k) = v1(i, k)
967	gz(nn, k) = gz1(i, k)
968	h(nn, k) = h1(i, k)
969	lv(nn, k) = lv1(i, k)
970	cpn(nn, k) = cpn1(i, k)
971	p(nn, k) = p1(i, k)
972	ph(nn, k) = ph1(i, k)
973	tv(nn, k) = tv1(i, k)
974	tp(nn, k) = tp1(i, k)
975	tvp(nn, k) = tvp1(i, k)
976	clw(nn, k) = clw1(i, k)
977	th(nn, k) = th1(i, k)
978	END IF
979	END DO
980	END DO
981
982	!AC! do 121 j=1,ntra
983	!AC!ccccc do 111 k=1,nl+1
984	!AC! do 111 k=1,nd
985	!AC! nn=0
986	!AC! do 101 i=1,len
987	!AC! if(iflag1(i).eq.0)then
988	!AC! nn=nn+1
989	!AC! tra(nn,k,j)=tra1(i,k,j)
990	!AC! endif
991	!AC! 101 continue
992	!AC! 111 continue
993	!AC! 121 continue
994
995	IF (nn/=ncum) THEN
996	WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
997	abort_message = ''
998	CALL abort_physic(modname, abort_message, 1)
999	END IF
1000
1001	nn = 0
1002	DO i = 1, len
1003	IF (iflag1(i)==0) THEN
1004	nn = nn + 1
1005	pbase(nn) = pbase1(i)
1006	buoybase(nn) = buoybase1(i)
1007	plcl(nn) = plcl1(i)
1008	tnk(nn) = tnk1(i)
1009	qnk(nn) = qnk1(i)
1010	gznk(nn) = gznk1(i)
1011	nk(nn) = nk1(i)
1012	icb(nn) = icb1(i)
1013	icbs(nn) = icbs1(i)
1014	iflag(nn) = iflag1(i)
1015	END IF
1016	END DO
1017
1018	RETURN
1019	END SUBROUTINE cv3_compress
1020
1021	SUBROUTINE icefrac(t, clw, qi, nl, len)
1022	IMPLICIT NONE
1023
1024
1025	!JAM--------------------------------------------------------------------
1026	! Calcul de la quantité d'eau sous forme de glace
1027	! --------------------------------------------------------------------
1028	INTEGER nl, len
1029	REAL qi(len, nl)
1030	REAL t(len, nl), clw(len, nl)
1031	REAL fracg
1032	INTEGER k, i
1033
1034	DO k = 3, nl
1035	DO i = 1, len
1036	IF (t(i,k)>263.15) THEN
1037	qi(i, k) = 0.
1038	ELSE
1039	IF (t(i,k)<243.15) THEN
1040	qi(i, k) = clw(i, k)
1041	ELSE
1042	fracg = (263.15-t(i,k))/20
1043	qi(i, k) = clw(i, k)*fracg
1044	END IF
1045	END IF
1046	! print*,t(i,k),qi(i,k),'temp,testglace'
1047	END DO
1048	END DO
1049
1050	RETURN
1051
1052	END SUBROUTINE icefrac
1053
1054	SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, &
1055	tnk, qnk, gznk, hnk, t, q, qs, gz, &
1056	p, ph, h, tv, lv, lf, pbase, buoybase, plcl, &
1057	inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
1058	IMPLICIT NONE
1059
1060	! ---------------------------------------------------------------------
1061	! Purpose:
1062	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1063	! &
1064	! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
1065	! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
1066	! &
1067	! FIND THE LEVEL OF NEUTRAL BUOYANCY
1068
1069	! Main differences convect3/convect4:
1070	! - icbs (input) is the first level above LCL (may differ from icb)
1071	! - many minor differences in the iterations
1072	! - condensed water not removed from tvp in convect3
1073	! - vertical profile of buoyancy computed here (use of buoybase)
1074	! - the determination of inb is different
1075	! - no inb1, only inb in output
1076	! ---------------------------------------------------------------------
1077
1078	include "cvthermo.h"
1079	include "cv3param.h"
1080	include "conema3.h"
1081	include "cvflag.h"
1082	include "YOMCST2.h"
1083
1084	!inputs:
1085	INTEGER, INTENT (IN) :: ncum, nd, nloc
1086	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, icbs, nk
1087	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, q, qs, gz
1088	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
1089	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
1090	REAL, DIMENSION (nloc), INTENT (IN) :: tnk, qnk, gznk
1091	REAL, DIMENSION (nloc), INTENT (IN) :: hnk
1092	REAL, DIMENSION (nloc, nd), INTENT (IN) :: lv, lf, tv, h
1093	REAL, DIMENSION (nloc), INTENT (IN) :: pbase, buoybase, plcl
1094
1095	!input/outputs:
1096	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: tp, tvp, clw ! Input for k = 1, icb+1 (computed in cv3_undilute1)
1097	! Output above
1098
1099	!outputs:
1100	INTEGER, DIMENSION (nloc), INTENT (OUT) :: inb
1101	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ep, sigp, hp
1102	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: buoy
1103	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: frac
1104
1105	!local variables:
1106	INTEGER i, j, k
1107	REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
1108	REAL als
1109	REAL qsat_new, snew, qi(nloc, nd)
1110	REAL by, defrac, pden, tbis
1111	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
1112	LOGICAL lcape(nloc)
1113	INTEGER iposit(nloc)
1114	REAL fracg
1115	REAL deltap
1116
1117	! =====================================================================
1118	! --- SOME INITIALIZATIONS
1119	! =====================================================================
1120
1121	DO k = 1, nl
1122	DO i = 1, ncum
1123	qi(i, k) = 0.
1124	END DO
1125	END DO
1126
1127	! =====================================================================
1128	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1129	! =====================================================================
1130
1131	! --- The procedure is to solve the equation.
1132	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
1133
1134	! * Calculate certain parcel quantities, including static energy *
1135
1136
1137	DO i = 1, ncum
1138	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))*tnk(i)+ &
1139	! debug qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
1140	qnk(i)(lv0-clmcpv(tnk(i)-273.15)) + gznk(i)
1141	END DO
1142
1143
1144	! * Find lifted parcel quantities above cloud base *
1145
1146
1147	DO k = minorig + 1, nl
1148	DO i = 1, ncum
1149	! ori if(k.ge.(icb(i)+1))then
1150	IF (k>=(icbs(i)+1)) THEN ! convect3
1151	tg = t(i, k)
1152	qg = qs(i, k)
1153	! debug alv=lv0-clmcpv*(t(i,k)-t0)
1154	alv = lv0 - clmcpv*(t(i,k)-273.15)
1155
1156	! First iteration.
1157
1158	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1159	s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
1160	alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
1161	s = 1./s
1162	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1163	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1164	tg = tg + s*(ah0(i)-ahg)
1165	! ori tg=max(tg,35.0)
1166	! debug tc=tg-t0
1167	tc = tg - 273.15
1168	denom = 243.5 + tc
1169	denom = max(denom, 1.0) ! convect3
1170	! ori if(tc.ge.0.0)then
1171	es = 6.112exp(17.67tc/denom)
1172	! ori else
1173	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1174	! ori endif
1175	qg = epses/(p(i,k)-es(1.-eps))
1176
1177	! Second iteration.
1178
1179	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1180	! ori s=1./s
1181	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1182	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1183	tg = tg + s*(ah0(i)-ahg)
1184	! ori tg=max(tg,35.0)
1185	! debug tc=tg-t0
1186	tc = tg - 273.15
1187	denom = 243.5 + tc
1188	denom = max(denom, 1.0) ! convect3
1189	! ori if(tc.ge.0.0)then
1190	es = 6.112exp(17.67tc/denom)
1191	! ori else
1192	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1193	! ori endif
1194	qg = epses/(p(i,k)-es(1.-eps))
1195
1196	! debug alv=lv0-clmcpv*(t(i,k)-t0)
1197	alv = lv0 - clmcpv*(t(i,k)-273.15)
1198	! print*,'cpd dans convect2 ',cpd
1199	! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
1200	! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
1201
1202	! ori c approximation here:
1203	! ori tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
1204
1205	! convect3: no approximation:
1206	IF (cvflag_ice) THEN
1207	tp(i, k) = max(0., (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i)))
1208	ELSE
1209	tp(i, k) = (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
1210	END IF
1211
1212	clw(i, k) = qnk(i) - qg
1213	clw(i, k) = max(0.0, clw(i,k))
1214	rg = qg/(1.-qnk(i))
1215	! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
1216	! convect3: (qg utilise au lieu du vrai mixing ratio rg):
1217	tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
1218	IF (cvflag_ice) THEN
1219	IF (clw(i,k)<1.E-11) THEN
1220	tp(i, k) = tv(i, k)
1221	tvp(i, k) = tv(i, k)
1222	END IF
1223	END IF
1224	!jyg<
1225	!! END IF ! Endif moved to the end of the loop
1226	!>jyg
1227
1228	IF (cvflag_ice) THEN
1229	!CR:attention boucle en klon dans Icefrac
1230	! Call Icefrac(t,clw,qi,nl,nloc)
1231	IF (t(i,k)>263.15) THEN
1232	qi(i, k) = 0.
1233	ELSE
1234	IF (t(i,k)<243.15) THEN
1235	qi(i, k) = clw(i, k)
1236	ELSE
1237	fracg = (263.15-t(i,k))/20
1238	qi(i, k) = clw(i, k)*fracg
1239	END IF
1240	END IF
1241	!CR: fin test
1242	IF (t(i,k)<263.15) THEN
1243	!CR: on commente les calculs d'Arnaud car division par zero
1244	! nouveau calcul propose par JYG
1245	! alv=lv0-clmcpv*(t(i,k)-273.15)
1246	! alf=lf0-clmci*(t(i,k)-273.15)
1247	! tg=tp(i,k)
1248	! tc=tp(i,k)-273.15
1249	! denom=243.5+tc
1250	! do j=1,3
1251	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1252	! il faudra que esi vienne en argument de la convection
1253	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1254	! tbis=t(i,k)+(tp(i,k)-tg)
1255	! esi=exp(23.33086-(6111.72784/tbis) + &
1256	! 0.15215*log(tbis))
1257	! qsat_new=epsesi/(p(i,k)-esi(1.-eps))
1258	! snew=cpd(1.-qnk(i))+clqnk(i)+alvalvqsat_new/ &
1259	! (rrvtbistbis)
1260	! snew=1./snew
1261	! print*,esi,qsat_new,snew,'esi,qsat,snew'
1262	! tp(i,k)=tg+(alfqi(i,k)+alvqg(1.-(esi/es)))snew
1263	! print*,k,tp(i,k),qnk(i),'avec glace'
1264	! print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
1265	! enddo
1266
1267	alv = lv0 - clmcpv*(t(i,k)-273.15)
1268	alf = lf0 + clmci*(t(i,k)-273.15)
1269	als = alf + alv
1270	tg = tp(i, k)
1271	tp(i, k) = t(i, k)
1272	DO j = 1, 3
1273	esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
1274	qsat_new = epsesi/(p(i,k)-esi(1.-eps))
1275	snew = cpd(1.-qnk(i)) + clqnk(i) + alvalsqsat_new/ &
1276	(rrvtp(i,k)tp(i,k))
1277	snew = 1./snew
1278	! c print*,esi,qsat_new,snew,'esi,qsat,snew'
1279	tp(i, k) = tp(i, k) + &
1280	((cpd(1.-qnk(i))+clqnk(i))*(tg-tp(i,k)) + &
1281	alv(qg-qsat_new)+alfqi(i,k))*snew
1282	! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
1283	! 'k,tp,q,qt,qi avec glace'
1284	END DO
1285
1286	!CR:reprise du code AJ
1287	clw(i, k) = qnk(i) - qsat_new
1288	clw(i, k) = max(0.0, clw(i,k))
1289	tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
1290	! print*,tvp(i,k),'tvp'
1291	END IF
1292	IF (clw(i,k)<1.E-11) THEN
1293	tp(i, k) = tv(i, k)
1294	tvp(i, k) = tv(i, k)
1295	END IF
1296	END IF ! (cvflag_ice)
1297	!jyg<
1298	END IF ! (k>=(icbs(i)+1))
1299	!>jyg
1300	END DO
1301	END DO
1302
1303	! =====================================================================
1304	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
1305	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
1306	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
1307	! =====================================================================
1308	!
1309	!jyg<
1310	DO k = 1, nl
1311	DO i = 1, ncum
1312	ep(i, k) = 0.0
1313	sigp(i, k) = spfac
1314	END DO
1315	END DO
1316	!>jyg
1317	!
1318	IF (flag_epkeorig/=1) THEN
1319	DO k = 1, nl ! convect3
1320	DO i = 1, ncum
1321	!jyg<
1322	IF(k>=icb(i)) THEN
1323	!>jyg
1324	pden = ptcrit - pbcrit
1325	ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1326	ep(i, k) = max(ep(i,k), 0.0)
1327	ep(i, k) = min(ep(i,k), epmax)
1328	!! sigp(i, k) = spfac ! jyg
1329	ENDIF ! (k>=icb(i))
1330	END DO
1331	END DO
1332	ELSE
1333	DO k = 1, nl
1334	DO i = 1, ncum
1335	IF(k>=icb(i)) THEN
1336	!! IF (k>=(nk(i)+1)) THEN
1337	!>jyg
1338	tca = tp(i, k) - t0
1339	IF (tca>=0.0) THEN
1340	elacrit = elcrit
1341	ELSE
1342	elacrit = elcrit*(1.0-tca/tlcrit)
1343	END IF
1344	elacrit = max(elacrit, 0.0)
1345	ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
1346	ep(i, k) = max(ep(i,k), 0.0)
1347	ep(i, k) = min(ep(i,k), epmax)
1348	!! sigp(i, k) = spfac ! jyg
1349	END IF ! (k>=icb(i))
1350	END DO
1351	END DO
1352	END IF
1353	!
1354	! =====================================================================
1355	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1356	! --- VIRTUAL TEMPERATURE
1357	! =====================================================================
1358
1359	! dans convect3, tvp est calcule en une seule fois, et sans retirer
1360	! l'eau condensee (~> reversible CAPE)
1361
1362	! ori do 340 k=minorig+1,nl
1363	! ori do 330 i=1,ncum
1364	! ori if(k.ge.(icb(i)+1))then
1365	! ori tvp(i,k)=tvp(i,k)(1.0-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	! oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1368	! ori endif
1369	! ori 330 continue
1370	! ori 340 continue
1371
1372	! ori do 350 i=1,ncum
1373	! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1374	! ori 350 continue
1375
1376	DO i = 1, ncum ! convect3
1377	tp(i, nlp) = tp(i, nl) ! convect3
1378	END DO ! convect3
1379
1380	! =====================================================================
1381	! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1382	! =====================================================================
1383
1384	! -- this is for convect3 only:
1385
1386	! first estimate of buoyancy:
1387
1388	!jyg : k-loop outside i-loop (07042015)
1389	DO k = 1, nl
1390	DO i = 1, ncum
1391	buoy(i, k) = tvp(i, k) - tv(i, k)
1392	END DO
1393	END DO
1394
1395	! set buoyancy=buoybase for all levels below base
1396	! for safety, set buoy(icb)=buoybase
1397
1398	!jyg : k-loop outside i-loop (07042015)
1399	DO k = 1, nl
1400	DO i = 1, ncum
1401	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1402	buoy(i, k) = buoybase(i)
1403	END IF
1404	END DO
1405	END DO
1406	DO i = 1, ncum
1407	! buoy(icb(i),k)=buoybase(i)
1408	buoy(i, icb(i)) = buoybase(i)
1409	END DO
1410
1411	! -- end convect3
1412
1413	! =====================================================================
1414	! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1415	! --- LEVEL OF NEUTRAL BUOYANCY
1416	! =====================================================================
1417
1418	! -- this is for convect3 only:
1419
1420	DO i = 1, ncum
1421	inb(i) = nl - 1
1422	iposit(i) = nl
1423	END DO
1424
1425
1426	! -- iposit(i) = first level, above icb, with positive buoyancy
1427	DO k = 1, nl - 1
1428	DO i = 1, ncum
1429	IF (k>=icb(i) .AND. buoy(i,k)>0.) THEN
1430	iposit(i) = min(iposit(i), k)
1431	END IF
1432	END DO
1433	END DO
1434
1435	DO i = 1, ncum
1436	IF (iposit(i)==nl) THEN
1437	iposit(i) = icb(i)
1438	END IF
1439	END DO
1440
1441	DO k = 1, nl - 1
1442	DO i = 1, ncum
1443	IF ((k>=iposit(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1444	inb(i) = min(inb(i), k)
1445	END IF
1446	END DO
1447	END DO
1448
1449	!CR fix computation of inb
1450	!keep flag or modify in all cases?
1451	IF (iflag_mix_adiab.eq.1) THEN
1452	DO i = 1, ncum
1453	cape(i)=0.
1454	inb(i)=icb(i)+1
1455	ENDDO
1456
1457	DO k = 2, nl
1458	DO i = 1, ncum
1459	IF ((k>=iposit(i))) THEN
1460	deltap = min(plcl(i), ph(i,k-1)) - min(plcl(i), ph(i,k))
1461	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
1462	IF (cape(i).gt.0.) THEN
1463	inb(i) = max(inb(i), k)
1464	END IF
1465	ENDIF
1466	ENDDO
1467	ENDDO
1468
1469	! DO i = 1, ncum
1470	! print*,"inb",inb(i)
1471	! ENDDO
1472
1473	endif
1474
1475	! -- end convect3
1476
1477	! ori do 510 i=1,ncum
1478	! ori cape(i)=0.0
1479	! ori capem(i)=0.0
1480	! ori inb(i)=icb(i)+1
1481	! ori inb1(i)=inb(i)
1482	! ori 510 continue
1483
1484	! Originial Code
1485
1486	! do 530 k=minorig+1,nl-1
1487	! do 520 i=1,ncum
1488	! if(k.ge.(icb(i)+1))then
1489	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1490	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1491	! cape(i)=cape(i)+by
1492	! if(by.ge.0.0)inb1(i)=k+1
1493	! if(cape(i).gt.0.0)then
1494	! inb(i)=k+1
1495	! capem(i)=cape(i)
1496	! endif
1497	! endif
1498	!520 continue
1499	!530 continue
1500	! do 540 i=1,ncum
1501	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1502	! cape(i)=capem(i)+byp
1503	! defrac=capem(i)-cape(i)
1504	! defrac=max(defrac,0.001)
1505	! frac(i)=-cape(i)/defrac
1506	! frac(i)=min(frac(i),1.0)
1507	! frac(i)=max(frac(i),0.0)
1508	!540 continue
1509
1510	! K Emanuel fix
1511
1512	! call zilch(byp,ncum)
1513	! do 530 k=minorig+1,nl-1
1514	! do 520 i=1,ncum
1515	! if(k.ge.(icb(i)+1))then
1516	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1517	! cape(i)=cape(i)+by
1518	! if(by.ge.0.0)inb1(i)=k+1
1519	! if(cape(i).gt.0.0)then
1520	! inb(i)=k+1
1521	! capem(i)=cape(i)
1522	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1523	! endif
1524	! endif
1525	!520 continue
1526	!530 continue
1527	! do 540 i=1,ncum
1528	! inb(i)=max(inb(i),inb1(i))
1529	! cape(i)=capem(i)+byp(i)
1530	! defrac=capem(i)-cape(i)
1531	! defrac=max(defrac,0.001)
1532	! frac(i)=-cape(i)/defrac
1533	! frac(i)=min(frac(i),1.0)
1534	! frac(i)=max(frac(i),0.0)
1535	!540 continue
1536
1537	! J Teixeira fix
1538
1539	! ori call zilch(byp,ncum)
1540	! ori do 515 i=1,ncum
1541	! ori lcape(i)=.true.
1542	! ori 515 continue
1543	! ori do 530 k=minorig+1,nl-1
1544	! ori do 520 i=1,ncum
1545	! ori if(cape(i).lt.0.0)lcape(i)=.false.
1546	! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1547	! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1548	! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1549	! ori cape(i)=cape(i)+by
1550	! ori if(by.ge.0.0)inb1(i)=k+1
1551	! ori if(cape(i).gt.0.0)then
1552	! ori inb(i)=k+1
1553	! ori capem(i)=cape(i)
1554	! ori endif
1555	! ori endif
1556	! ori 520 continue
1557	! ori 530 continue
1558	! ori do 540 i=1,ncum
1559	! ori cape(i)=capem(i)+byp(i)
1560	! ori defrac=capem(i)-cape(i)
1561	! ori defrac=max(defrac,0.001)
1562	! ori frac(i)=-cape(i)/defrac
1563	! ori frac(i)=min(frac(i),1.0)
1564	! ori frac(i)=max(frac(i),0.0)
1565	! ori 540 continue
1566
1567	! =====================================================================
1568	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1569	! =====================================================================
1570
1571	DO k = 1, nl
1572	DO i = 1, ncum
1573	hp(i, k) = h(i, k)
1574	END DO
1575	END DO
1576
1577	!jyg : cvflag_ice test outside the loops (07042015)
1578	!
1579	IF (cvflag_ice) THEN
1580	!
1581	DO k = minorig + 1, nl
1582	DO i = 1, ncum
1583	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1584	frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
1585	frac(i, k) = min(max(frac(i,k),0.0), 1.0)
1586	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac(i,k)lf(i,k))* &
1587	ep(i, k)*clw(i, k)
1588	END IF
1589	END DO
1590	END DO
1591	! Below cloud base, set ice fraction to cloud base value
1592	DO k = 1, nl
1593	DO i = 1, ncum
1594	IF (k<icb(i)) THEN
1595	frac(i,k) = frac(i,icb(i))
1596	END IF
1597	END DO
1598	END DO
1599	!
1600	ELSE
1601	!
1602	DO k = minorig + 1, nl
1603	DO i = 1, ncum
1604	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1605	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k)
1606	END IF
1607	END DO
1608	END DO
1609	!
1610	END IF ! (cvflag_ice)
1611
1612	RETURN
1613	END SUBROUTINE cv3_undilute2
1614
1615	SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
1616	pbase, p, ph, tv, buoy, &
1617	sig, w0, cape, m, iflag)
1618	IMPLICIT NONE
1619
1620	! ===================================================================
1621	! --- CLOSURE OF CONVECT3
1622	!
1623	! vectorization: S. Bony
1624	! ===================================================================
1625
1626	include "cvthermo.h"
1627	include "cv3param.h"
1628
1629	!input:
1630	INTEGER ncum, nd, nloc
1631	INTEGER icb(nloc), inb(nloc)
1632	REAL pbase(nloc)
1633	REAL p(nloc, nd), ph(nloc, nd+1)
1634	REAL tv(nloc, nd), buoy(nloc, nd)
1635
1636	!input/output:
1637	REAL sig(nloc, nd), w0(nloc, nd)
1638	INTEGER iflag(nloc)
1639
1640	!output:
1641	REAL cape(nloc)
1642	REAL m(nloc, nd)
1643
1644	!local variables:
1645	INTEGER i, j, k, icbmax
1646	REAL deltap, fac, w, amu
1647	REAL dtmin(nloc, nd), sigold(nloc, nd)
1648	REAL cbmflast(nloc)
1649
1650
1651	! -------------------------------------------------------
1652	! -- Initialization
1653	! -------------------------------------------------------
1654
1655	DO k = 1, nl
1656	DO i = 1, ncum
1657	m(i, k) = 0.0
1658	END DO
1659	END DO
1660
1661	! -------------------------------------------------------
1662	! -- Reset sig(i) and w0(i) for i>inb and i<icb
1663	! -------------------------------------------------------
1664
1665	! update sig and w0 above LNB:
1666
1667	DO k = 1, nl - 1
1668	DO i = 1, ncum
1669	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
1670	sig(i, k) = beta*sig(i, k) + &
1671	2.alphabuoy(i, inb(i))*abs(buoy(i,inb(i)))
1672	sig(i, k) = amax1(sig(i,k), 0.0)
1673	w0(i, k) = beta*w0(i, k)
1674	END IF
1675	END DO
1676	END DO
1677
1678	! compute icbmax:
1679
1680	icbmax = 2
1681	DO i = 1, ncum
1682	icbmax = max(icbmax, icb(i))
1683	END DO
1684
1685	! update sig and w0 below cloud base:
1686
1687	DO k = 1, icbmax
1688	DO i = 1, ncum
1689	IF (k<=icb(i)) THEN
1690	sig(i, k) = beta*sig(i, k) - &
1691	2.alphabuoy(i, icb(i))*buoy(i, icb(i))
1692	sig(i, k) = max(sig(i,k), 0.0)
1693	w0(i, k) = beta*w0(i, k)
1694	END IF
1695	END DO
1696	END DO
1697
1698	!! if(inb.lt.(nl-1))then
1699	!! do 85 i=inb+1,nl-1
1700	!! sig(i)=betasig(i)+2.alphabuoy(inb)
1701	!! 1 abs(buoy(inb))
1702	!! sig(i)=max(sig(i),0.0)
1703	!! w0(i)=beta*w0(i)
1704	!! 85 continue
1705	!! end if
1706
1707	!! do 87 i=1,icb
1708	!! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
1709	!! sig(i)=max(sig(i),0.0)
1710	!! w0(i)=beta*w0(i)
1711	!! 87 continue
1712
1713	! -------------------------------------------------------------
1714	! -- Reset fractional areas of updrafts and w0 at initial time
1715	! -- and after 10 time steps of no convection
1716	! -------------------------------------------------------------
1717
1718	DO k = 1, nl - 1
1719	DO i = 1, ncum
1720	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
1721	sig(i, k) = 0.0
1722	w0(i, k) = 0.0
1723	END IF
1724	END DO
1725	END DO
1726
1727	! -------------------------------------------------------------
1728	! -- Calculate convective available potential energy (cape),
1729	! -- vertical velocity (w), fractional area covered by
1730	! -- undilute updraft (sig), and updraft mass flux (m)
1731	! -------------------------------------------------------------
1732
1733	DO i = 1, ncum
1734	cape(i) = 0.0
1735	END DO
1736
1737	! compute dtmin (minimum buoyancy between ICB and given level k):
1738
1739	DO i = 1, ncum
1740	DO k = 1, nl
1741	dtmin(i, k) = 100.0
1742	END DO
1743	END DO
1744
1745	DO i = 1, ncum
1746	DO k = 1, nl
1747	DO j = minorig, nl
1748	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
1749	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
1750	END IF
1751	END DO
1752	END DO
1753	END DO
1754
1755	! the interval on which cape is computed starts at pbase :
1756
1757	DO k = 1, nl
1758	DO i = 1, ncum
1759
1760	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
1761
1762	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
1763	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
1764	cape(i) = amax1(0.0, cape(i))
1765	sigold(i, k) = sig(i, k)
1766
1767	! dtmin(i,k)=100.0
1768	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
1769	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1770	! 97 continue
1771
1772	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
1773	sig(i, k) = max(sig(i,k), 0.0)
1774	sig(i, k) = amin1(sig(i,k), 0.01)
1775	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
1776	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
1777	amu = 0.5(sig(i,k)+sigold(i,k))w
1778	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
1779	w0(i, k) = w
1780	END IF
1781
1782	END DO
1783	END DO
1784
1785	DO i = 1, ncum
1786	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
1787	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))
1788	sig(i, icb(i)) = sig(i, icb(i)+1)
1789	sig(i, icb(i)-1) = sig(i, icb(i))
1790	END DO
1791
1792	! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
1793	! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
1794	! ccc the final mass flux (cbmflast) is greater than the target mass flux
1795	! ccc (cbmf) ??).
1796	! cc
1797	! c do i = 1,ncum
1798	! c cbmflast(i) = 0.
1799	! c enddo
1800	! cc
1801	! c do k= 1,nl
1802	! c do i = 1,ncum
1803	! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
1804	! c cbmflast(i) = cbmflast(i)+M(i,k)
1805	! c ENDIF
1806	! c enddo
1807	! c enddo
1808	! cc
1809	! c do i = 1,ncum
1810	! c IF (cbmflast(i) .lt. 1.e-6) THEN
1811	! c iflag(i) = 3
1812	! c ENDIF
1813	! c enddo
1814	! cc
1815	! c do k= 1,nl
1816	! c do i = 1,ncum
1817	! c IF (iflag(i) .ge. 3) THEN
1818	! c M(i,k) = 0.
1819	! c sig(i,k) = 0.
1820	! c w0(i,k) = 0.
1821	! c ENDIF
1822	! c enddo
1823	! c enddo
1824	! cc
1825	!! cape=0.0
1826	!! do 98 i=icb+1,inb
1827	!! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1828	!! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
1829	!! dcape=rrdbuoy(i-1)deltap/p(i-1)
1830	!! dlnp=deltap/p(i-1)
1831	!! cape=max(0.0,cape)
1832	!! sigold=sig(i)
1833
1834	!! dtmin=100.0
1835	!! do 97 j=icb,i-1
1836	!! dtmin=amin1(dtmin,buoy(j))
1837	!! 97 continue
1838
1839	!! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
1840	!! sig(i)=max(sig(i),0.0)
1841	!! sig(i)=amin1(sig(i),0.01)
1842	!! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1843	!! w=(1.-beta)facsqrt(cape)+beta*w0(i)
1844	!! amu=0.5(sig(i)+sigold)w
1845	!! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
1846	!! w0(i)=w
1847	!! 98 continue
1848	!! w0(icb)=0.5*w0(icb+1)
1849	!! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1850	!! sig(icb)=sig(icb+1)
1851	!! sig(icb-1)=sig(icb)
1852
1853	RETURN
1854	END SUBROUTINE cv3_closure
1855
1856	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
1857	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
1858	unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
1859	ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
1860	IMPLICIT NONE
1861
1862	! ---------------------------------------------------------------------
1863	! a faire:
1864	! - vectorisation de la partie normalisation des flux (do 789...)
1865	! ---------------------------------------------------------------------
1866
1867	include "cvthermo.h"
1868	include "cv3param.h"
1869	include "cvflag.h"
1870
1871	!inputs:
1872	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
1873	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
1874	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
1875	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
1876	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
1877	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
1878	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
1879	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
1880	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, h, hp
1881	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf, frac
1882	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp, ep, clw
1883	REAL, DIMENSION (nloc, na), INTENT (IN) :: m ! input of convect3
1884
1885	!outputs:
1886	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: ment, qent
1887	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
1888	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: sij, elij
1889	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
1890	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: ments, qents
1891	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
1892
1893	!local variables:
1894	INTEGER i, j, k, il, im, jm
1895	INTEGER num1, num2
1896	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1897	REAL alt, smid, sjmin, sjmax, delp, delm
1898	REAL asij(nloc), smax(nloc), scrit(nloc)
1899	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
1900	REAL sigij(nloc, nd, nd)
1901	REAL wgh
1902	REAL zm(nloc, na)
1903	LOGICAL lwork(nloc)
1904
1905	! =====================================================================
1906	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1907	! =====================================================================
1908
1909	! ori do 360 i=1,ncum*nlp
1910	DO j = 1, nl
1911	DO i = 1, ncum
1912	nent(i, j) = 0
1913	! in convect3, m is computed in cv3_closure
1914	! ori m(i,1)=0.0
1915	END DO
1916	END DO
1917
1918	! ori do 400 k=1,nlp
1919	! ori do 390 j=1,nlp
1920	DO j = 1, nl
1921	DO k = 1, nl
1922	DO i = 1, ncum
1923	qent(i, k, j) = rr(i, j)
1924	uent(i, k, j) = u(i, j)
1925	vent(i, k, j) = v(i, j)
1926	elij(i, k, j) = 0.0
1927	!ym ment(i,k,j)=0.0
1928	!ym sij(i,k,j)=0.0
1929	END DO
1930	END DO
1931	END DO
1932
1933	!ym
1934	ment(1:ncum, 1:nd, 1:nd) = 0.0
1935	sij(1:ncum, 1:nd, 1:nd) = 0.0
1936
1937	!AC! do k=1,ntra
1938	!AC! do j=1,nd ! instead nlp
1939	!AC! do i=1,nd ! instead nlp
1940	!AC! do il=1,ncum
1941	!AC! traent(il,i,j,k)=tra(il,j,k)
1942	!AC! enddo
1943	!AC! enddo
1944	!AC! enddo
1945	!AC! enddo
1946	zm(:, :) = 0.
1947
1948	! =====================================================================
1949	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1950	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1951	! --- FRACTION (sij)
1952	! =====================================================================
1953
1954	DO i = minorig + 1, nl
1955
1956	DO j = minorig, nl
1957	DO il = 1, ncum
1958	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
1959
1960	rti = qnk(il) - ep(il, i)*clw(il, i)
1961	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
1962
1963
1964	IF (cvflag_ice) THEN
1965	! print*,cvflag_ice,'cvflag_ice dans do 700'
1966	IF (t(il,j)<=263.15) THEN
1967	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
1968	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
1969	END IF
1970	END IF
1971
1972	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
1973	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
1974	dei = denom
1975	IF (abs(dei)<0.01) dei = 0.01
1976	sij(il, i, j) = anum/dei
1977	sij(il, i, i) = 1.0
1978	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
1979	altem = altem/bf2
1980	cwat = clw(il, j)*(1.-ep(il,j))
1981	stemp = sij(il, i, j)
1982	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
1983
1984	IF (cvflag_ice) THEN
1985	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
1986	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
1987	ELSE
1988	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
1989	denom = denom + lv(il, j)*(rr(il,i)-rti)
1990	END IF
1991
1992	IF (abs(denom)<0.01) denom = 0.01
1993	sij(il, i, j) = anum/denom
1994	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
1995	altem = altem - (bf2-1.)*cwat
1996	END IF
1997	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
1998	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
1999	uent(il, i, j) = sij(il, i, j)u(il, i) + (1.-sij(il,i,j))unk(il)
2000	vent(il, i, j) = sij(il, i, j)v(il, i) + (1.-sij(il,i,j))vnk(il)
2001	!!!! do k=1,ntra
2002	!!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2003	!!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2004	!!!! end do
2005	elij(il, i, j) = altem
2006	elij(il, i, j) = max(0.0, elij(il,i,j))
2007	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
2008	nent(il, i) = nent(il, i) + 1
2009	END IF
2010	sij(il, i, j) = max(0.0, sij(il,i,j))
2011	sij(il, i, j) = amin1(1.0, sij(il,i,j))
2012	END IF ! new
2013	END DO
2014	END DO
2015
2016	!AC! do k=1,ntra
2017	!AC! do j=minorig,nl
2018	!AC! do il=1,ncum
2019	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
2020	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
2021	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2022	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2023	!AC! endif
2024	!AC! enddo
2025	!AC! enddo
2026	!AC! enddo
2027
2028
2029	! * if no air can entrain at level i assume that updraft detrains *
2030	! * at that level and calculate detrained air flux and properties *
2031
2032
2033	! @ do 170 i=icb(il),inb(il)
2034
2035	DO il = 1, ncum
2036	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
2037	! @ if(nent(il,i).eq.0)then
2038	ment(il, i, i) = m(il, i)
2039	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2040	uent(il, i, i) = unk(il)
2041	vent(il, i, i) = vnk(il)
2042	elij(il, i, i) = clw(il, i)
2043	! MAF sij(il,i,i)=1.0
2044	sij(il, i, i) = 0.0
2045	END IF
2046	END DO
2047	END DO
2048
2049	!AC! do j=1,ntra
2050	!AC! do i=minorig+1,nl
2051	!AC! do il=1,ncum
2052	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
2053	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2054	!AC! endif
2055	!AC! enddo
2056	!AC! enddo
2057	!AC! enddo
2058
2059	DO j = minorig, nl
2060	DO i = minorig, nl
2061	DO il = 1, ncum
2062	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
2063	sigij(il, i, j) = sij(il, i, j)
2064	END IF
2065	END DO
2066	END DO
2067	END DO
2068	! @ enddo
2069
2070	! @170 continue
2071
2072	! =====================================================================
2073	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
2074	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
2075	! =====================================================================
2076
2077	CALL zilch(asum, nloc*nd)
2078	CALL zilch(csum, nloc*nd)
2079	CALL zilch(csum, nloc*nd)
2080
2081	DO il = 1, ncum
2082	lwork(il) = .FALSE.
2083	END DO
2084
2085	DO i = minorig + 1, nl
2086
2087	num1 = 0
2088	DO il = 1, ncum
2089	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
2090	END DO
2091	IF (num1<=0) GO TO 789
2092
2093
2094	DO il = 1, ncum
2095	IF (i>=icb(il) .AND. i<=inb(il)) THEN
2096	lwork(il) = (nent(il,i)/=0)
2097	qp = qnk(il) - ep(il, i)*clw(il, i)
2098
2099	IF (cvflag_ice) THEN
2100
2101	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
2102	(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
2103	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
2104	(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
2105	ELSE
2106
2107	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
2108	(cpv-cpd)t(il, i)(qp-rr(il,i))
2109	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
2110	(cpd-cpv)t(il, i)(rr(il,i)-qp)
2111	END IF
2112
2113	IF (abs(denom)<0.01) denom = 0.01
2114	scrit(il) = anum/denom
2115	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
2116	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
2117	smax(il) = 0.0
2118	asij(il) = 0.0
2119	END IF
2120	END DO
2121
2122	DO j = nl, minorig, -1
2123
2124	num2 = 0
2125	DO il = 1, ncum
2126	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2127	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2128	lwork(il)) num2 = num2 + 1
2129	END DO
2130	IF (num2<=0) GO TO 175
2131
2132	DO il = 1, ncum
2133	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2134	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2135	lwork(il)) THEN
2136
2137	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
2138	wgh = 1.0
2139	IF (j>i) THEN
2140	sjmax = max(sij(il,i,j+1), smax(il))
2141	sjmax = amin1(sjmax, scrit(il))
2142	smax(il) = max(sij(il,i,j), smax(il))
2143	sjmin = max(sij(il,i,j-1), smax(il))
2144	sjmin = amin1(sjmin, scrit(il))
2145	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
2146	smid = amin1(sij(il,i,j), scrit(il))
2147	ELSE
2148	sjmax = max(sij(il,i,j+1), scrit(il))
2149	smid = max(sij(il,i,j), scrit(il))
2150	sjmin = 0.0
2151	IF (j>1) sjmin = sij(il, i, j-1)
2152	sjmin = max(sjmin, scrit(il))
2153	END IF
2154	delp = abs(sjmax-smid)
2155	delm = abs(sjmin-smid)
2156	asij(il) = asij(il) + wgh*(delp+delm)
2157	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
2158	END IF
2159	END IF
2160	END DO
2161
2162	175 END DO
2163
2164	DO il = 1, ncum
2165	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2166	asij(il) = max(1.0E-16, asij(il))
2167	asij(il) = 1.0/asij(il)
2168	asum(il, i) = 0.0
2169	bsum(il, i) = 0.0
2170	csum(il, i) = 0.0
2171	END IF
2172	END DO
2173
2174	DO j = minorig, nl
2175	DO il = 1, ncum
2176	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2177	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2178	ment(il, i, j) = ment(il, i, j)*asij(il)
2179	END IF
2180	END DO
2181	END DO
2182
2183	DO j = minorig, nl
2184	DO il = 1, ncum
2185	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2186	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2187	asum(il, i) = asum(il, i) + ment(il, i, j)
2188	ment(il, i, j) = ment(il, i, j)*sig(il, j)
2189	bsum(il, i) = bsum(il, i) + ment(il, i, j)
2190	END IF
2191	END DO
2192	END DO
2193
2194	DO il = 1, ncum
2195	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2196	bsum(il, i) = max(bsum(il,i), 1.0E-16)
2197	bsum(il, i) = 1.0/bsum(il, i)
2198	END IF
2199	END DO
2200
2201	DO j = minorig, nl
2202	DO il = 1, ncum
2203	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2204	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2205	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
2206	END IF
2207	END DO
2208	END DO
2209
2210	DO j = minorig, nl
2211	DO il = 1, ncum
2212	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2213	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2214	csum(il, i) = csum(il, i) + ment(il, i, j)
2215	END IF
2216	END DO
2217	END DO
2218
2219	DO il = 1, ncum
2220	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2221	csum(il,i)<m(il,i)) THEN
2222	nent(il, i) = 0
2223	ment(il, i, i) = m(il, i)
2224	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2225	uent(il, i, i) = unk(il)
2226	vent(il, i, i) = vnk(il)
2227	elij(il, i, i) = clw(il, i)
2228	! MAF sij(il,i,i)=1.0
2229	sij(il, i, i) = 0.0
2230	END IF
2231	END DO ! il
2232
2233	!AC! do j=1,ntra
2234	!AC! do il=1,ncum
2235	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
2236	!AC! : .and. csum(il,i).lt.m(il,i) ) then
2237	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2238	!AC! endif
2239	!AC! enddo
2240	!AC! enddo
2241	789 END DO
2242
2243	! MAF: renormalisation de MENT
2244	CALL zilch(zm, nloc*na)
2245	DO jm = 1, nl
2246	DO im = 1, nl
2247	DO il = 1, ncum
2248	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
2249	END DO
2250	END DO
2251	END DO
2252
2253	DO jm = 1, nl
2254	DO im = 1, nl
2255	DO il = 1, ncum
2256	IF (zm(il,im)/=0.) THEN
2257	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
2258	END IF
2259	END DO
2260	END DO
2261	END DO
2262
2263	DO jm = 1, nl
2264	DO im = 1, nl
2265	DO il = 1, ncum
2266	qents(il, im, jm) = qent(il, im, jm)
2267	ments(il, im, jm) = ment(il, im, jm)
2268	END DO
2269	END DO
2270	END DO
2271
2272	RETURN
2273	END SUBROUTINE cv3_mixing
2274
2275	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
2276	t, rr, rs, gz, u, v, tra, p, ph, &
2277	th, tv, lv, lf, cpn, ep, sigp, clw, &
2278	m, ment, elij, delt, plcl, coef_clos, &
2279	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
2280	faci, b, sigd, &
2281	wdtrainA, wdtrainM) ! RomP
2282	USE print_control_mod, ONLY: prt_level, lunout
2283	IMPLICIT NONE
2284
2285
2286	include "cvthermo.h"
2287	include "cv3param.h"
2288	include "cvflag.h"
2289	include "nuage.h"
2290
2291	!inputs:
2292	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2293	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
2294	REAL, INTENT(IN) :: delt
2295	REAL, DIMENSION (nloc), INTENT (IN) :: plcl
2296	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
2297	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz
2298	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
2299	REAL tra(nloc, nd, ntra)
2300	REAL p(nloc, nd), ph(nloc, nd+1)
2301	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, sigp, clw
2302	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tv, lv, cpn
2303	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
2304	REAL, DIMENSION (nloc, na), INTENT (IN) :: m
2305	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: ment, elij
2306	REAL, DIMENSION (nloc), INTENT (IN) :: coef_clos
2307
2308	!input/output
2309	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag(nloc)
2310
2311	!outputs:
2312	REAL, DIMENSION (nloc, na), INTENT (OUT) :: mp, rp, up, vp
2313	REAL, DIMENSION (nloc, na), INTENT (OUT) :: water, evap, wt
2314	REAL, DIMENSION (nloc, na), INTENT (OUT) :: ice, fondue, faci
2315	REAL, DIMENSION (nloc, na, ntra), INTENT (OUT) :: trap
2316	REAL, DIMENSION (nloc, na), INTENT (OUT) :: b
2317	REAL, DIMENSION (nloc), INTENT (OUT) :: sigd
2318	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
2319	! de l ascendance adiabatique et des flux melanges Pa et Pm.
2320	! Distinction des wdtrain
2321	! Pa = wdtrainA Pm = wdtrainM
2322	REAL, DIMENSION (nloc, na), INTENT (OUT) :: wdtrainA, wdtrainM
2323
2324	!local variables
2325	INTEGER i, j, k, il, num1, ndp1
2326	REAL tinv, delti, coef
2327	REAL awat, afac, afac1, afac2, bfac
2328	REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
2329	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
2330	REAL ampmax, thaw
2331	REAL tevap(nloc)
2332	REAL lvcp(nloc, na), lfcp(nloc, na)
2333	REAL h(nloc, na), hm(nloc, na)
2334	REAL frac(nloc, na)
2335	REAL fraci(nloc, na), prec(nloc, na)
2336	REAL wdtrain(nloc)
2337	LOGICAL lwork(nloc), mplus(nloc)
2338
2339
2340	! ------------------------------------------------------
2341
2342	! =============================
2343	! --- INITIALIZE OUTPUT ARRAYS
2344	! =============================
2345	! (loops up to nl+1)
2346
2347	DO i = 1, nlp
2348	DO il = 1, ncum
2349	mp(il, i) = 0.0
2350	rp(il, i) = rr(il, i)
2351	up(il, i) = u(il, i)
2352	vp(il, i) = v(il, i)
2353	wt(il, i) = 0.001
2354	water(il, i) = 0.0
2355	faci(il, i) = 0.0
2356	ice(il, i) = 0.0
2357	fondue(il, i) = 0.0
2358	evap(il, i) = 0.0
2359	b(il, i) = 0.0
2360	END DO
2361	END DO
2362	!! RomP >>>
2363	DO i = 1, nlp
2364	DO il = 1, ncum
2365	wdtrainA(il, i) = 0.0
2366	wdtrainM(il, i) = 0.0
2367	END DO
2368	END DO
2369	!! RomP <<<
2370
2371	! *** Set the fractionnal area sigd of precipitating downdraughts
2372	DO il = 1, ncum
2373	sigd(il) = sigdz*coef_clos(il)
2374	END DO
2375
2376	! =====================================================================
2377	! --- INITIALIZE VARIOUS ARRAYS AND PARAMETERS USED IN THE COMPUTATIONS
2378	! =====================================================================
2379	! (loops up to nl+1)
2380
2381	delti = 1./delt
2382	tinv = 1./3.
2383
2384	DO i = 1, nlp
2385	DO il = 1, ncum
2386	frac(il, i) = 0.0
2387	fraci(il, i) = 0.0
2388	prec(il, i) = 0.0
2389	lvcp(il, i) = lv(il, i)/cpn(il, i)
2390	lfcp(il, i) = lf(il, i)/cpn(il, i)
2391	END DO
2392	END DO
2393
2394	!AC! do k=1,ntra
2395	!AC! do i=1,nd
2396	!AC! do il=1,ncum
2397	!AC! trap(il,i,k)=tra(il,i,k)
2398	!AC! enddo
2399	!AC! enddo
2400	!AC! enddo
2401
2402	! * check whether ep(inb)=0, if so, skip precipitating *
2403	! * downdraft calculation *
2404
2405
2406	DO il = 1, ncum
2407	!! lwork(il)=.TRUE.
2408	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
2409	lwork(il) = ep(il, inb(il)) >= 0.0001
2410	END DO
2411
2412
2413	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2414	!
2415	! * begin downdraft loop *
2416	!
2417	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2418
2419	DO i = nl + 1, 1, -1
2420
2421	num1 = 0
2422	DO il = 1, ncum
2423	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
2424	END DO
2425	IF (num1<=0) GO TO 400
2426
2427	CALL zilch(wdtrain, ncum)
2428
2429
2430	! * integrate liquid water equation to find condensed water *
2431	! * and condensed water flux *
2432	!
2433	!
2434	! * calculate detrained precipitation *
2435
2436	DO il = 1, ncum
2437	IF (i<=inb(il) .AND. lwork(il)) THEN
2438	IF (cvflag_grav) THEN
2439	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
2440	wdtrainA(il, i) = wdtrain(il)/grav ! Pa RomP
2441	ELSE
2442	wdtrain(il) = 10.0ep(il, i)m(il, i)*clw(il, i)
2443	wdtrainA(il, i) = wdtrain(il)/10. ! Pa RomP
2444	END IF
2445	END IF
2446	END DO
2447
2448	IF (i>1) THEN
2449	DO j = 1, i - 1
2450	DO il = 1, ncum
2451	IF (i<=inb(il) .AND. lwork(il)) THEN
2452	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
2453	awat = max(awat, 0.0)
2454	IF (cvflag_grav) THEN
2455	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
2456	wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP
2457	ELSE
2458	wdtrain(il) = wdtrain(il) + 10.0awatment(il, j, i)
2459	wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i) ! Pm RomP
2460	END IF
2461	END IF
2462	END DO
2463	END DO
2464	END IF
2465
2466
2467	! * find rain water and evaporation using provisional *
2468	! * estimates of rp(i)and rp(i-1) *
2469
2470
2471	DO il = 1, ncum
2472	IF (i<=inb(il) .AND. lwork(il)) THEN
2473
2474	wt(il, i) = 45.0
2475
2476	IF (cvflag_ice) THEN
2477	frac(il, inb(il)) = 1. - (t(il,inb(il))-243.15)/(263.15-243.15)
2478	frac(il, inb(il)) = min(max(frac(il,inb(il)),0.), 1.)
2479	fraci(il, inb(il)) = frac(il, inb(il))
2480	ELSE
2481	CONTINUE
2482	END IF
2483
2484	IF (i<inb(il)) THEN
2485
2486	IF (cvflag_ice) THEN
2487	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
2488	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
2489	thaw = min(max(thaw,0.0), 1.0)
2490	frac(il, i) = frac(il, i)*(1.-thaw)
2491	ELSE
2492	CONTINUE
2493	END IF
2494
2495	rp(il, i) = rp(il, i+1) + &
2496	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
2497	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
2498	END IF
2499	fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
2500	fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
2501	rp(il, i) = max(rp(il,i), 0.0)
2502	rp(il, i) = amin1(rp(il,i), rs(il,i))
2503	rp(il, inb(il)) = rr(il, inb(il))
2504
2505	IF (i==1) THEN
2506	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
2507	IF (cvflag_ice) THEN
2508	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
2509	END IF
2510	ELSE
2511	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)
2512	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
2513	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
2514	rp(il, i-1) = max(rp(il,i-1), 0.0)
2515	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
2516	afac2 = p(il, i-1)(rs(il,i-1)-rp(il,i-1))/(1.0E4+2000.0p(il,i-1)*rs(il,i-1))
2517	afac = 0.5*(afac1+afac2)
2518	END IF
2519	IF (i==inb(il)) afac = 0.0
2520	afac = max(afac, 0.0)
2521	bfac = 1./(sigd(il)*wt(il,i))
2522
2523	!
2524	IF (prt_level >= 20) THEN
2525	Print*, 'cv3_unsat after provisional rp estimate: rp, afac, bfac ', &
2526	i, rp(1, i), afac,bfac
2527	ENDIF
2528	!
2529	!JYG1
2530	! cc sigt=1.0
2531	! cc if(i.ge.icb)sigt=sigp(i)
2532	! prise en compte de la variation progressive de sigt dans
2533	! les couches icb et icb-1:
2534	! pour plcl<ph(i+1), pr1=0 & pr2=1
2535	! pour plcl>ph(i), pr1=1 & pr2=0
2536	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2537	! sur le nuage, et pr2 est la proportion sous la base du
2538	! nuage.
2539	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2540	pr1 = max(0., min(1.,pr1))
2541	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2542	pr2 = max(0., min(1.,pr2))
2543	sigt = sigp(il, i)*pr1 + pr2
2544	!JYG2
2545
2546	!JYG----
2547	! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2548	! c6 = water(il,i+1) + wdtrain(il)*bfac
2549	! c6 = prec(il,i+1) + wdtrain(il)*bfac
2550	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
2551	! evap(il,i)=sigtafacrevap
2552	! water(il,i)=revap*revap
2553	! prec(il,i)=revap*revap
2554	!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
2555	!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
2556	!!---end jyg---
2557
2558	! --------retour à la formulation originale d'Emanuel.
2559	IF (cvflag_ice) THEN
2560
2561	! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2562	! c6=prec(il,i+1)+bfac*wdtrain(il) &
2563	! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2564	! if(c6.gt.0.0)then
2565	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
2566
2567	!JAM Attention: evap=sigt*E
2568	! Modification: evap devient l'évaporation en milieu de couche
2569	! car nécessaire dans cv3_yield
2570	! Du coup, il faut modifier pas mal d'équations...
2571	! et l'expression de afac qui devient afac1
2572	! revap=sqrt((prec(i+1)+prec(i))/2)
2573
2574	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
2575	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
2576	! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
2577	! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
2578	! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
2579	IF (c6>b6*b6+1.E-20) THEN
2580	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
2581	ELSE
2582	revap = (-b6+sqrt(b6b6+4.c6))/2.
2583	END IF
2584	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
2585	! print*,prec(il,i),'neige'
2586
2587	!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
2588	! c evap(il,i)=sigtafacrevap
2589	! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
2590	! Ici,l'evaporation evap est simplement calculee par l'equation de
2591	! conservation.
2592	! prec(il,i)=revap*revap
2593	! else
2594	!JYG---- Correction : si c6 <= 0, water(il,i)=0.
2595	! prec(il,i)=0.
2596	! endif
2597
2598	!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
2599	! moins [tt ce qui sort de la couche i]
2600	! print *, 'evap avec ice'
2601	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
2602	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2603	!
2604	IF (prt_level >= 20) THEN
2605	Print*, 'cv3_unsat after evap computation: wdtrain, sigd, wt, prec(i+1),prec(i) ', &
2606	i, wdtrain(1), sigd(1), wt(1,i), prec(1,i+1),prec(1,i)
2607	ENDIF
2608	!
2609
2610	!jyg<
2611	d6 = prec(il,i)-prec(il,i+1)
2612
2613	!! d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
2614	!! e6 = bfac*wdtrain(il)
2615	!! f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
2616	!>jyg
2617	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
2618	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
2619	thaw = min(max(thaw,0.0), 1.0)
2620	!jyg<
2621	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
2622	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
2623	water(il, i) = min(prec(il,i), max(water(il,i), 0.))
2624	ice(il, i) = min(prec(il,i), max(ice(il,i), 0.))
2625
2626	!! water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
2627	!! water(il, i) = max(water(il,i), 0.)
2628	!! ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
2629	!! ice(il, i) = max(ice(il,i), 0.)
2630	!>jyg
2631	fondue(il, i) = ice(il, i)*thaw
2632	water(il, i) = water(il, i) + fondue(il, i)
2633	ice(il, i) = ice(il, i) - fondue(il, i)
2634
2635	IF (water(il,i)+ice(il,i)<1.E-30) THEN
2636	faci(il, i) = 0.
2637	ELSE
2638	faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
2639	END IF
2640
2641	! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
2642	! water(il,i)=max(water(il,i),0.)
2643	! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
2644	! ice(il,i)=max(ice(il,i),0.)
2645	! fondue(il,i)=ice(il,i)*thaw
2646	! water(il,i)=water(il,i)+fondue(il,i)
2647	! ice(il,i)=ice(il,i)-fondue(il,i)
2648
2649	! if((water(il,i)+ice(il,i)).lt.1.e-30)then
2650	! faci(il,i)=0.
2651	! else
2652	! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
2653	! endif
2654
2655	ELSE
2656	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2657	c6 = water(il, i+1) + bfac*wdtrain(il) - &
2658	50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
2659	IF (c6>0.0) THEN
2660	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
2661	water(il, i) = revap*revap
2662	ELSE
2663	water(il, i) = 0.
2664	END IF
2665	! print *, 'evap sans ice'
2666	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
2667	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2668
2669	END IF
2670	END IF !(i.le.inb(il) .and. lwork(il))
2671	END DO
2672	! ----------------------------------------------------------------
2673
2674	! cc
2675	! * calculate precipitating downdraft mass flux under *
2676	! * hydrostatic approximation *
2677
2678	DO il = 1, ncum
2679	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
2680
2681	tevap(il) = max(0.0, evap(il,i))
2682	delth = max(0.001, (th(il,i)-th(il,i-1)))
2683	IF (cvflag_ice) THEN
2684	IF (cvflag_grav) THEN
2685	mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
2686	(p(il,i-1)-p(il,i))/delth + &
2687	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
2688	(p(il,i-1)-p(il,i))/delth + &
2689	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
2690	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
2691	ELSE
2692	mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
2693	(p(il,i-1)-p(il,i))/delth + &
2694	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
2695	(p(il,i-1)-p(il,i))/delth + &
2696	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
2697	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
2698
2699	END IF
2700	ELSE
2701	IF (cvflag_grav) THEN
2702	mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
2703	(p(il,i-1)-p(il,i))/delth
2704	ELSE
2705	mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
2706	(p(il,i-1)-p(il,i))/delth
2707	END IF
2708
2709	END IF
2710
2711	END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2712	END DO
2713	! ----------------------------------------------------------------
2714
2715	! * if hydrostatic assumption fails, *
2716	! * solve cubic difference equation for downdraft theta *
2717	! * and mass flux from two simultaneous differential eqns *
2718
2719	DO il = 1, ncum
2720	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
2721
2722	amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
2723	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2724	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2725
2726	IF (amp2>(0.1*amfac)) THEN
2727	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
2728	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
2729	(lvcp(il,i)sigd(il)th(il,i))
2730	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
2731
2732	IF (cvflag_ice) THEN
2733	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2734	50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
2735	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
2736	ELSE
2737
2738	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2739	50.(p(il,i-1)-p(il,i))xf*tevap(il)
2740	END IF
2741
2742	fac2 = 1.0
2743	IF (bf<0.0) fac2 = -1.0
2744	bf = abs(bf)
2745	ur = 0.25bfbf - afafaftinvtinv*tinv
2746	IF (ur>=0.0) THEN
2747	sru = sqrt(ur)
2748	fac = 1.0
2749	IF ((0.5*bf-sru)<0.0) fac = -1.0
2750	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
2751	fac(abs(0.5bf-sru))**tinv
2752	ELSE
2753	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
2754	IF (fac2<0.0) d = 3.14159 - d
2755	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
2756	END IF
2757	mp(il, i) = max(0.0, mp(il,i))
2758
2759	IF (cvflag_ice) THEN
2760	IF (cvflag_grav) THEN
2761	!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2762	! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
2763	! Et il faut bien revoir les facteurs 100.
2764	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i)) &
2765	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
2766	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
2767	(ph(il,i)-ph(il,i+1))) / &
2768	(mp(il,i)+sigd(il)*0.1) - &
2769	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2770	(lvcp(il,i)sigd(il)th(il,i))
2771	ELSE
2772	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
2773	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
2774	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
2775	(ph(il,i)-ph(il,i+1))) / &
2776	(mp(il,i)+sigd(il)*0.1) - &
2777	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2778	(lvcp(il,i)sigd(il)th(il,i))
2779	END IF
2780	ELSE
2781	IF (cvflag_grav) THEN
2782	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
2783	(mp(il,i)+sigd(il)*0.1) - &
2784	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2785	(lvcp(il,i)sigd(il)th(il,i))
2786	ELSE
2787	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
2788	(mp(il,i)+sigd(il)*0.1) - &
2789	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2790	(lvcp(il,i)sigd(il)th(il,i))
2791	END IF
2792	END IF
2793	b(il, i-1) = max(b(il,i-1), 0.0)
2794
2795	END IF !(amp2.gt.(0.1*amfac))
2796
2797	! * limit magnitude of mp(i) to meet cfl condition *
2798
2799	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
2800	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
2801	ampmax = min(ampmax, amp2)
2802	mp(il, i) = min(mp(il,i), ampmax)
2803
2804	! * force mp to decrease linearly to zero *
2805	! * between cloud base and the surface *
2806
2807
2808	! c if(p(il,i).gt.p(il,icb(il)))then
2809	! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2810	! c endif
2811	IF (ph(il,i)>0.9*plcl(il)) THEN
2812	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
2813	END IF
2814
2815	END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2816	END DO
2817	! ----------------------------------------------------------------
2818	!
2819	IF (prt_level >= 20) THEN
2820	Print*, 'cv3_unsat after mp computation: mp, b(i), b(i-1) ', &
2821	i, mp(1, i), b(1,i), b(1,max(i-1,1))
2822	ENDIF
2823	!
2824
2825	! * find mixing ratio of precipitating downdraft *
2826
2827	DO il = 1, ncum
2828	IF (i<inb(il) .AND. lwork(il)) THEN
2829	mplus(il) = mp(il, i) > mp(il, i+1)
2830	END IF ! (i.lt.inb(il) .and. lwork(il))
2831	END DO
2832
2833	DO il = 1, ncum
2834	IF (i<inb(il) .AND. lwork(il)) THEN
2835
2836	rp(il, i) = rr(il, i)
2837
2838	IF (mplus(il)) THEN
2839
2840	IF (cvflag_grav) THEN
2841	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
2842	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2843	ELSE
2844	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
2845	5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2846	END IF
2847	rp(il, i) = rp(il, i)/mp(il, i)
2848	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
2849	up(il, i) = up(il, i)/mp(il, i)
2850	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
2851	vp(il, i) = vp(il, i)/mp(il, i)
2852
2853	ELSE ! if (mplus(il))
2854
2855	IF (mp(il,i+1)>1.0E-16) THEN
2856	IF (cvflag_grav) THEN
2857	rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
2858	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2859	ELSE
2860	rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
2861	(evap(il,i+1)+evap(il,i))/mp(il, i+1)
2862	END IF
2863	up(il, i) = up(il, i+1)
2864	vp(il, i) = vp(il, i+1)
2865	END IF ! (mp(il,i+1).gt.1.0e-16)
2866	END IF ! (mplus(il)) else if (.not.mplus(il))
2867
2868	rp(il, i) = amin1(rp(il,i), rs(il,i))
2869	rp(il, i) = max(rp(il,i), 0.0)
2870
2871	END IF ! (i.lt.inb(il) .and. lwork(il))
2872	END DO
2873	! ----------------------------------------------------------------
2874
2875	! * find tracer concentrations in precipitating downdraft *
2876
2877	!AC! do j=1,ntra
2878	!AC! do il = 1,ncum
2879	!AC! if (i.lt.inb(il) .and. lwork(il)) then
2880	!AC!c
2881	!AC! if(mplus(il))then
2882	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2883	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2884	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2885	!AC! else ! if (mplus(il))
2886	!AC! if(mp(il,i+1).gt.1.0e-16)then
2887	!AC! trap(il,i,j)=trap(il,i+1,j)
2888	!AC! endif
2889	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
2890	!AC!c
2891	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
2892	!AC! enddo
2893	!AC! end do
2894
2895	400 END DO
2896	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2897
2898	! * end of downdraft loop *
2899
2900	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2901
2902
2903	RETURN
2904	END SUBROUTINE cv3_unsat
2905
2906	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
2907	icb, inb, delt, &
2908	t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
2909	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
2910	ep, clw, m, tp, mp, rp, up, vp, trap, &
2911	wt, water, ice, evap, fondue, faci, b, sigd, &
2912	ment, qent, hent, iflag_mix, uent, vent, &
2913	nent, elij, traent, sig, &
2914	tv, tvp, wghti, &
2915	iflag, precip, Vprecip, Vprecipi, & ! jyg: Vprecipi
2916	ft, fr, fu, fv, ftra, & ! jyg
2917	cbmf, upwd, dnwd, dnwd0, ma, mip, &
2918	!! tls, tps, ! useless . jyg
2919	qcondc, wd, &
2920	ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
2921
2922	IMPLICIT NONE
2923
2924	include "cvthermo.h"
2925	include "cv3param.h"
2926	include "cvflag.h"
2927	include "conema3.h"
2928
2929	!inputs:
2930	INTEGER, INTENT (IN) :: iflag_mix
2931	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2932	LOGICAL, INTENT (IN) :: ok_conserv_q
2933	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
2934	REAL, INTENT (IN) :: delt
2935	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, u, v
2936	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t_wake, rr_wake
2937	REAL, DIMENSION (nloc), INTENT (IN) :: s_wake
2938	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra
2939	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
2940	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2941	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz, h, hp
2942	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tp
2943	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, cpn, ep, clw
2944	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
2945	REAL, DIMENSION (nloc, na), INTENT (IN) :: rp, up
2946	REAL, DIMENSION (nloc, na), INTENT (IN) :: vp
2947	REAL, DIMENSION (nloc, nd), INTENT (IN) :: wt
2948	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: trap
2949	REAL, DIMENSION (nloc, na), INTENT (IN) :: water, evap, b
2950	REAL, DIMENSION (nloc, na), INTENT (IN) :: fondue, faci, ice
2951	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: qent, uent
2952	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: hent
2953	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: vent, elij
2954	INTEGER, DIMENSION (nloc, nd), INTENT (IN) :: nent
2955	REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent
2956	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti
2957	REAL,INTENT(IN) :: tau_cld_cv, coefw_cld_cv
2958	!
2959	!input/output:
2960	REAL, DIMENSION (nloc, na), INTENT (INOUT) :: m, mp
2961	REAL, DIMENSION (nloc, na, na), INTENT (INOUT) :: ment
2962	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
2963	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig
2964	REAL, DIMENSION (nloc), INTENT (INOUT) :: sigd
2965	!
2966	!outputs:
2967	REAL, DIMENSION (nloc), INTENT (OUT) :: precip
2968	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ft, fr, fu, fv
2969	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ftd, fqd
2970	REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT) :: ftra
2971	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: upwd, dnwd, ma
2972	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: dnwd0, mip
2973	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecip
2974	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecipi
2975	!! REAL tls(nloc, nd), tps(nloc, nd) ! useless . jyg
2976	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qcondc ! cld
2977	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qtc, sigt ! cld
2978	REAL, DIMENSION (nloc), INTENT (OUT) :: wd ! gust
2979	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf
2980	!
2981	!local variables:
2982	INTEGER :: i, k, il, n, j, num1
2983	REAL :: rat, delti
2984	REAL :: ax, bx, cx, dx, ex
2985	REAL :: cpinv, rdcp, dpinv
2986	REAL, DIMENSION (nloc) :: awat
2987	REAL, DIMENSION (nloc, nd) :: lvcp, lfcp ! , mke ! unused . jyg
2988	REAL, DIMENSION (nloc) :: am, work, ad, amp1
2989	!! real up1(nloc), dn1(nloc)
2990	REAL, DIMENSION (nloc, nd, nd) :: up1, dn1
2991	!jyg<
2992	REAL, DIMENSION (nloc, nd) :: up_to, up_from
2993	REAL, DIMENSION (nloc, nd) :: dn_to, dn_from
2994	!>jyg
2995	REAL, DIMENSION (nloc) :: asum, bsum, csum, dsum
2996	REAL, DIMENSION (nloc) :: esum, fsum, gsum, hsum
2997	REAL, DIMENSION (nloc, nd) :: th_wake
2998	REAL, DIMENSION (nloc) :: alpha_qpos, alpha_qpos1
2999	REAL, DIMENSION (nloc, nd) :: qcond, nqcond, wa ! cld
3000	REAL, DIMENSION (nloc, nd) :: siga, sax, mac ! cld
3001	REAL, DIMENSION (nloc) :: sument
3002	REAL, DIMENSION (nloc, nd) :: sigment, qtment ! cld
3003	REAL, DIMENSION (nloc) :: qnk
3004	REAL sumdq !jyg
3005	!
3006	! -------------------------------------------------------------
3007
3008	! initialization:
3009
3010	delti = 1.0/delt
3011	! print*,'cv3_yield initialisation delt', delt
3012
3013	DO il = 1, ncum
3014	precip(il) = 0.0
3015	wd(il) = 0.0 ! gust
3016	END DO
3017
3018	! Fluxes are on a staggered grid : loops extend up to nl+1
3019	DO i = 1, nlp
3020	DO il = 1, ncum
3021	Vprecip(il, i) = 0.0
3022	Vprecipi(il, i) = 0.0 ! jyg
3023	upwd(il, i) = 0.0
3024	dnwd(il, i) = 0.0
3025	dnwd0(il, i) = 0.0
3026	mip(il, i) = 0.0
3027	END DO
3028	END DO
3029	DO i = 1, nl
3030	DO il = 1, ncum
3031	ft(il, i) = 0.0
3032	fr(il, i) = 0.0
3033	fu(il, i) = 0.0
3034	fv(il, i) = 0.0
3035	ftd(il, i) = 0.0
3036	fqd(il, i) = 0.0
3037	qcondc(il, i) = 0.0 ! cld
3038	qcond(il, i) = 0.0 ! cld
3039	qtc(il, i) = 0.0 ! cld
3040	qtment(il, i) = 0.0 ! cld
3041	sigment(il, i) = 0.0 ! cld
3042	sigt(il, i) = 0.0 ! cld
3043	nqcond(il, i) = 0.0 ! cld
3044	END DO
3045	END DO
3046	! print*,'cv3_yield initialisation 2'
3047	!AC! do j=1,ntra
3048	!AC! do i=1,nd
3049	!AC! do il=1,ncum
3050	!AC! ftra(il,i,j)=0.0
3051	!AC! enddo
3052	!AC! enddo
3053	!AC! enddo
3054	! print*,'cv3_yield initialisation 3'
3055	DO i = 1, nl
3056	DO il = 1, ncum
3057	lvcp(il, i) = lv(il, i)/cpn(il, i)
3058	lfcp(il, i) = lf(il, i)/cpn(il, i)
3059	END DO
3060	END DO
3061
3062
3063
3064	! * calculate surface precipitation in mm/day *
3065
3066	DO il = 1, ncum
3067	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
3068	IF (cvflag_ice) THEN
3069	precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
3070	86400.1000./(rowl*grav)
3071	ELSE
3072	precip(il) = wt(il, 1)sigd(il)water(il, 1) &
3073	86400.1000./(rowl*grav)
3074	END IF
3075	END IF
3076	END DO
3077	! print*,'cv3_yield apres calcul precip'
3078
3079
3080	! === calculate vertical profile of precipitation in kg/m2/s ===
3081
3082	DO i = 1, nl
3083	DO il = 1, ncum
3084	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
3085	IF (cvflag_ice) THEN
3086	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
3087	Vprecipi(il, i) = wt(il, i)sigd(il)ice(il,i)/grav ! jyg
3088	ELSE
3089	Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
3090	Vprecipi(il, i) = 0. ! jyg
3091	END IF
3092	END IF
3093	END DO
3094	END DO
3095
3096
3097	! * Calculate downdraft velocity scale *
3098	! * NE PAS UTILISER POUR L'INSTANT *
3099
3100	!! do il=1,ncum
3101	!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
3102	!! /(sigd(il)*p(il,icb(il)))
3103	!! enddo
3104
3105
3106	! * calculate tendencies of lowest level potential temperature *
3107	! * and mixing ratio *
3108
3109	DO il = 1, ncum
3110	work(il) = 1.0/(ph(il,1)-ph(il,2))
3111	cbmf(il) = 0.0
3112	END DO
3113
3114	DO k = 2, nl
3115	DO il = 1, ncum
3116	IF (k>=icb(il)) THEN
3117	cbmf(il) = cbmf(il) + m(il, k)
3118	END IF
3119	END DO
3120	END DO
3121
3122	! print*,'cv3_yield avant ft'
3123	! am is the part of cbmf taken from the first level
3124	DO il = 1, ncum
3125	am(il) = cbmf(il)*wghti(il, 1)
3126	END DO
3127
3128	DO il = 1, ncum
3129	IF (iflag(il)<=1) THEN
3130	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
3131	!JYG Correction pour conserver l'eau
3132	! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
3133	IF (cvflag_ice) THEN
3134	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
3135	lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
3136	lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
3137	(100.*(ph(il,1)-ph(il,2))) !precip
3138	ELSE
3139	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
3140	END IF
3141
3142	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
3143
3144	IF (cvflag_ice) THEN
3145	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3146	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
3147	0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
3148	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3149	ELSE
3150	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3151	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3152	END IF
3153
3154	ftd(il, 1) = ft(il, 1) ! fin precip
3155
3156	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
3157	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3158	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
3159	END IF ! iflag
3160	END DO
3161
3162
3163	DO j = 2, nl
3164	IF (iflag_mix>0) THEN
3165	DO il = 1, ncum
3166	! FH WARNING a modifier :
3167	cpinv = 0.
3168	! cpinv=1.0/cpn(il,1)
3169	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3170	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
3171	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
3172	END IF ! j
3173	END DO
3174	END IF
3175	END DO
3176	! fin sature
3177
3178
3179	DO il = 1, ncum
3180	IF (iflag(il)<=1) THEN
3181	!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
3182	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
3183	sigd(il)*evap(il, 1)
3184	!!! sigd(il)0.5(evap(il,1)+evap(il,2))
3185
3186	fqd(il, 1) = fr(il, 1) !precip
3187
3188	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
3189
3190	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
3191	am(il)*(u(il,2)-u(il,1)))
3192	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
3193	am(il)*(v(il,2)-v(il,1)))
3194	END IF ! iflag
3195	END DO ! il
3196
3197
3198	!AC! do j=1,ntra
3199	!AC! do il=1,ncum
3200	!AC! if (iflag(il) .le. 1) then
3201	!AC! if (cvflag_grav) then
3202	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
3203	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3204	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3205	!AC! else
3206	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
3207	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3208	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3209	!AC! endif
3210	!AC! endif ! iflag
3211	!AC! enddo
3212	!AC! enddo
3213
3214	DO j = 2, nl
3215	DO il = 1, ncum
3216	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3217	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3218	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
3219	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
3220	END IF ! j
3221	END DO
3222	END DO
3223
3224	!AC! do k=1,ntra
3225	!AC! do j=2,nl
3226	!AC! do il=1,ncum
3227	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
3228	!AC!
3229	!AC! if (cvflag_grav) then
3230	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
3231	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3232	!AC! else
3233	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
3234	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3235	!AC! endif
3236	!AC!
3237	!AC! endif
3238	!AC! enddo
3239	!AC! enddo
3240	!AC! enddo
3241	! print*,'cv3_yield apres ft'
3242
3243	!jyg<
3244	!-----------------------------------------------------------
3245	IF (ok_optim_yield) THEN !\|
3246	!-----------------------------------------------------------
3247	!
3248	!* *
3249	!* Compute convective mass fluxes upwd and dnwd *
3250
3251	upwd(:,:) = 0.
3252	up_to(:,:) = 0.
3253	up_from(:,:) = 0.
3254	dnwd(:,:) = 0.
3255	dn_to(:,:) = 0.
3256	dn_from(:,:) = 0.
3257	!
3258	! =================================================
3259	! upward fluxes \|
3260	! ------------------------------------------------
3261	DO i = 2, nl
3262	DO il = 1, ncum
3263	IF (i<=inb(il)) THEN
3264	up_to(il,i) = m(il,i)
3265	ENDIF
3266	ENDDO
3267	DO j = 1, i-1
3268	DO il = 1, ncum
3269	IF (i<=inb(il)) THEN
3270	up_to(il,i) = up_to(il,i) + ment(il,j,i)
3271	ENDIF
3272	ENDDO
3273	ENDDO
3274	ENDDO
3275	!
3276	DO i = 1, nl
3277	DO il = 1, ncum
3278	IF (i<=inb(il)) THEN
3279	up_from(il,i) = cbmf(il)*wghti(il,i)
3280	ENDIF
3281	ENDDO
3282	ENDDO
3283	!!DO i = 2, nl
3284	!! DO j = i+1, nl !! Permuter les boucles i et j
3285	DO j = 3, nl
3286	DO i = 2, j-1
3287	DO il = 1, ncum
3288	IF (j<=inb(il)) THEN
3289	up_from(il,i) = up_from(il,i) + ment(il,i,j)
3290	ENDIF
3291	ENDDO
3292	ENDDO
3293	ENDDO
3294	!
3295	! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
3296	!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
3297	!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
3298	!
3299	DO i = 2, nlp
3300	DO il = 1, ncum
3301	upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
3302	ENDDO
3303	ENDDO
3304	!
3305	! =================================================
3306	! downward fluxes \|
3307	! ------------------------------------------------
3308	DO i = 1, nl
3309	DO j = i+1, nl
3310	DO il = 1, ncum
3311	IF (j<=inb(il)) THEN
3312	dn_to(il,i) = dn_to(il,i) + ment(il,j,i)
3313	ENDIF
3314	ENDDO
3315	ENDDO
3316	ENDDO
3317	!
3318	!!DO i = 2, nl
3319	!! DO j = 1, i-1 !! Permuter les boucles i et j
3320	DO j = 1, nl
3321	DO i = j+1, nl
3322	DO il = 1, ncum
3323	IF (i<=inb(il)) THEN
3324	dn_from(il,i) = dn_from(il,i) + ment(il,i,j)
3325	ENDIF
3326	ENDDO
3327	ENDDO
3328	ENDDO
3329	!
3330	! The difference between dnwd(il,i) and dnwd(il,i+1) is due to downdrafts ending in layer
3331	!(i) (theses drafts cross interface (i+1) but not interface(i)) and to downdrafts
3332	!starting from layer (i) (theses drafts cross interface (i) but not interface(i+1)):
3333	!
3334	DO i = nl-1, 1, -1
3335	DO il = 1, ncum
3336	dnwd(il,i) = max(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i))
3337	ENDDO
3338	ENDDO
3339	! =================================================
3340	!
3341	!-----------------------------------------------------------
3342	ENDIF !(ok_optim_yield) !\|
3343	!-----------------------------------------------------------
3344	!>jyg
3345
3346	! * calculate tendencies of potential temperature and mixing ratio *
3347	! * at levels above the lowest level *
3348
3349	! * first find the net saturated updraft and downdraft mass fluxes *
3350	! * through each level *
3351
3352
3353	!jyg<
3354	!! DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
3355	DO i = 2, nl
3356	!>jyg
3357
3358	num1 = 0
3359	DO il = 1, ncum
3360	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
3361	END DO
3362	IF (num1<=0) GO TO 500
3363
3364	!
3365	!jyg<
3366	!-----------------------------------------------------------
3367	IF (ok_optim_yield) THEN !\|
3368	!-----------------------------------------------------------
3369	DO il = 1, ncum
3370	amp1(il) = upwd(il,i+1)
3371	ad(il) = dnwd(il,i)
3372	ENDDO
3373	!-----------------------------------------------------------
3374	ELSE !(ok_optim_yield) !\|
3375	!-----------------------------------------------------------
3376	!>jyg
3377	DO il = 1,ncum
3378	amp1(il) = 0.
3379	ad(il) = 0.
3380	ENDDO
3381
3382	DO k = 1, nl + 1
3383	DO il = 1, ncum
3384	IF (i>=icb(il)) THEN
3385	IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
3386	amp1(il) = amp1(il) + m(il, k)
3387	END IF
3388	ELSE
3389	! AMP1 is the part of cbmf taken from layers I and lower
3390	IF (k<=i) THEN
3391	amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
3392	END IF
3393	END IF
3394	END DO
3395	END DO
3396
3397	DO j = i + 1, nl + 1
3398	DO k = 1, i
3399	!yor! reverted j and k loops
3400	DO il = 1, ncum
3401	!yor! IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN ! the second condition implies the first !
3402	IF (j<=(inb(il)+1)) THEN
3403	amp1(il) = amp1(il) + ment(il, k, j)
3404	END IF
3405	END DO
3406	END DO
3407	END DO
3408
3409	DO k = 1, i - 1
3410	!jyg<
3411	!! DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
3412	DO j = i, nl
3413	!>jyg
3414	DO il = 1, ncum
3415	!yor! IF (i<=inb(il) .AND. j<=inb(il)) THEN ! the second condition implies the 1st !
3416	IF (j<=inb(il)) THEN
3417	ad(il) = ad(il) + ment(il, j, k)
3418	END IF
3419	END DO
3420	END DO
3421	END DO
3422	!
3423	!-----------------------------------------------------------
3424	ENDIF !(ok_optim_yield) !\|
3425	!-----------------------------------------------------------
3426	!
3427	!! print *,'yield, i, amp1, ad', i, amp1(1), ad(1)
3428
3429	DO il = 1, ncum
3430	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3431	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3432	cpinv = 1.0/cpn(il, i)
3433
3434	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
3435	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
3436
3437	! precip
3438	! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
3439	IF (cvflag_ice) THEN
3440	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
3441	sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
3442	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
3443	ELSE
3444	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
3445	END IF
3446
3447	rat = cpn(il, i-1)*cpinv
3448
3449	ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
3450	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
3451	IF (cvflag_ice) THEN
3452	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
3453	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
3454	0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
3455	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
3456	ELSE
3457	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
3458	(t_wake(il,i+1)-t_wake(il,i))dpinv &
3459	cpinv
3460	END IF
3461
3462	ftd(il, i) = ft(il, i)
3463	! fin precip
3464
3465	! sature
3466	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
3467	(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
3468	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
3469
3470
3471	IF (iflag_mix==0) THEN
3472	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
3473	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
3474	END IF
3475
3476
3477
3478	! sb: on ne fait pas encore la correction permettant de mieux
3479	! conserver l'eau:
3480	!JYG: correction permettant de mieux conserver l'eau:
3481	! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
3482	fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
3483	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
3484	fqd(il, i) = fr(il, i) ! precip
3485
3486	fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
3487	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
3488	fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
3489	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
3490
3491
3492	fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
3493	ad(il)*(rr(il,i)-rr(il,i-1)))
3494	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
3495	ad(il)*(u(il,i)-u(il,i-1)))
3496	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
3497	ad(il)*(v(il,i)-v(il,i-1)))
3498
3499	END IF ! i
3500	END DO
3501
3502	!AC! do k=1,ntra
3503	!AC! do il=1,ncum
3504	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3505	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3506	!AC! cpinv=1.0/cpn(il,i)
3507	!AC! if (cvflag_grav) then
3508	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
3509	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
3510	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
3511	!AC! else
3512	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
3513	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
3514	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
3515	!AC! endif
3516	!AC! endif
3517	!AC! enddo
3518	!AC! enddo
3519
3520	DO k = 1, i - 1
3521
3522	DO il = 1, ncum
3523	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
3524	awat(il) = max(awat(il), 0.0)
3525	END DO
3526
3527	IF (iflag_mix/=0) THEN
3528	DO il = 1, ncum
3529	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3530	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3531	cpinv = 1.0/cpn(il, i)
3532	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
3533	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
3534	!
3535	!
3536	END IF ! i
3537	END DO
3538	END IF
3539
3540	DO il = 1, ncum
3541	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3542	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3543	cpinv = 1.0/cpn(il, i)
3544	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
3545	(qent(il,k,i)-awat(il)-rr(il,i))
3546	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
3547	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
3548
3549	! (saturated updrafts resulting from mixing) ! cld
3550	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
3551	qtment(il, i) = qtment(il, i) + qent(il,k,i) ! cld
3552	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3553	END IF ! i
3554	END DO
3555	END DO
3556
3557	!AC! do j=1,ntra
3558	!AC! do k=1,i-1
3559	!AC! do il=1,ncum
3560	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3561	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3562	!AC! cpinv=1.0/cpn(il,i)
3563	!AC! if (cvflag_grav) then
3564	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3565	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3566	!AC! else
3567	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3568	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3569	!AC! endif
3570	!AC! endif
3571	!AC! enddo
3572	!AC! enddo
3573	!AC! enddo
3574
3575	!jyg<
3576	!! DO k = i, nl + 1
3577	DO k = i, nl
3578	!>jyg
3579
3580	IF (iflag_mix/=0) THEN
3581	DO il = 1, ncum
3582	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
3583	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3584	cpinv = 1.0/cpn(il, i)
3585	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
3586	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
3587
3588
3589	END IF ! i
3590	END DO
3591	END IF
3592
3593	DO il = 1, ncum
3594	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
3595	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3596	cpinv = 1.0/cpn(il, i)
3597
3598	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
3599	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
3600	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
3601	END IF ! i and k
3602	END DO
3603	END DO
3604
3605	!AC! do j=1,ntra
3606	!AC! do k=i,nl+1
3607	!AC! do il=1,ncum
3608	!AC! if (i.le.inb(il) .and. k.le.inb(il)
3609	!AC! $ .and. iflag(il) .le. 1) then
3610	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3611	!AC! cpinv=1.0/cpn(il,i)
3612	!AC! if (cvflag_grav) then
3613	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3614	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3615	!AC! else
3616	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3617	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3618	!AC! endif
3619	!AC! endif ! i and k
3620	!AC! enddo
3621	!AC! enddo
3622	!AC! enddo
3623
3624	! sb: interface with the cloud parameterization: ! cld
3625
3626	DO k = i + 1, nl
3627	DO il = 1, ncum
3628	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
3629	! (saturated downdrafts resulting from mixing) ! cld
3630	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
3631	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
3632	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3633	END IF ! cld
3634	END DO ! cld
3635	END DO ! cld
3636
3637	! (particular case: no detraining level is found) ! cld
3638	DO il = 1, ncum ! cld
3639	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
3640	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
3641	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
3642	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3643	END IF ! cld
3644	END DO ! cld
3645
3646	DO il = 1, ncum ! cld
3647	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
3648	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
3649	qtment(il, i) = qtment(il,i)/nqcond(il, i) ! cld
3650	END IF ! cld
3651	END DO
3652
3653	!AC! do j=1,ntra
3654	!AC! do il=1,ncum
3655	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3656	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3657	!AC! cpinv=1.0/cpn(il,i)
3658	!AC!
3659	!AC! if (cvflag_grav) then
3660	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
3661	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3662	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3663	!AC! else
3664	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
3665	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3666	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3667	!AC! endif
3668	!AC! endif ! i
3669	!AC! enddo
3670	!AC! enddo
3671
3672
3673	500 END DO
3674
3675	!JYG<
3676	!Conservation de l'eau
3677	! sumdq = 0.
3678	! DO k = 1, nl
3679	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3680	! END DO
3681	! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3682	!JYG>
3683	! * move the detrainment at level inb down to level inb-1 *
3684	! * in such a way as to preserve the vertically *
3685	! * integrated enthalpy and water tendencies *
3686
3687	! Correction bug le 18-03-09
3688	DO il = 1, ncum
3689	IF (iflag(il)<=1) THEN
3690	ax = 0.01gravment(il, inb(il), inb(il))* &
3691	(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
3692	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
3693	ft(il, inb(il)) = ft(il, inb(il)) - ax
3694	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3695	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
3696
3697	bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
3698	(ph(il,inb(il))-ph(il,inb(il)+1))
3699	fr(il, inb(il)) = fr(il, inb(il)) - bx
3700	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3701	(ph(il,inb(il)-1)-ph(il,inb(il)))
3702
3703	cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
3704	(ph(il,inb(il))-ph(il,inb(il)+1))
3705	fu(il, inb(il)) = fu(il, inb(il)) - cx
3706	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3707	(ph(il,inb(il)-1)-ph(il,inb(il)))
3708
3709	dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
3710	(ph(il,inb(il))-ph(il,inb(il)+1))
3711	fv(il, inb(il)) = fv(il, inb(il)) - dx
3712	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3713	(ph(il,inb(il)-1)-ph(il,inb(il)))
3714	END IF !iflag
3715	END DO
3716
3717	!JYG<
3718	!Conservation de l'eau
3719	! sumdq = 0.
3720	! DO k = 1, nl
3721	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3722	! END DO
3723	! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3724	!JYG>
3725
3726	!AC! do j=1,ntra
3727	!AC! do il=1,ncum
3728	!AC! IF (iflag(il) .le. 1) THEN
3729	!AC! IF (cvflag_grav) then
3730	!AC! ex=0.01gravment(il,inb(il),inb(il))
3731	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3732	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3733	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3734	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3735	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3736	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3737	!AC! else
3738	!AC! ex=0.1*ment(il,inb(il),inb(il))
3739	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3740	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3741	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3742	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3743	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3744	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3745	!AC! ENDIF !cvflag grav
3746	!AC! ENDIF !iflag
3747	!AC! enddo
3748	!AC! enddo
3749
3750
3751	! * homogenize tendencies below cloud base *
3752
3753
3754	DO il = 1, ncum
3755	asum(il) = 0.0
3756	bsum(il) = 0.0
3757	csum(il) = 0.0
3758	dsum(il) = 0.0
3759	esum(il) = 0.0
3760	fsum(il) = 0.0
3761	gsum(il) = 0.0
3762	hsum(il) = 0.0
3763	END DO
3764
3765	!do i=1,nl
3766	!do il=1,ncum
3767	!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
3768	!enddo
3769	!enddo
3770
3771	DO i = 1, nl
3772	DO il = 1, ncum
3773	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
3774	!jyg Saturated part : use T profile
3775	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
3776	!jyg<20140311
3777	!Correction pour conserver l eau
3778	IF (ok_conserv_q) THEN
3779	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
3780	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
3781
3782	ELSE
3783	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
3784	(ph(il,i)-ph(il,i+1))
3785	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
3786	(ph(il,i)-ph(il,i+1))
3787	ENDIF ! (ok_conserv_q)
3788	!jyg>
3789	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
3790	!jyg Unsaturated part : use T_wake profile
3791	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
3792	!jyg<20140311
3793	!Correction pour conserver l eau
3794	IF (ok_conserv_q) THEN
3795	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
3796	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
3797	ELSE
3798	fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
3799	(ph(il,i)-ph(il,i+1))
3800	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
3801	(ph(il,i)-ph(il,i+1))
3802	ENDIF ! (ok_conserv_q)
3803	!jyg>
3804	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
3805	END IF
3806	END DO
3807	END DO
3808
3809	!!!! do 700 i=1,icb(il)-1
3810	DO i = 1, nl
3811	DO il = 1, ncum
3812	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
3813	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
3814	fqd(il, i) = fsum(il)/gsum(il)
3815	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
3816	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
3817	END IF
3818	END DO
3819	END DO
3820
3821	!jyg<
3822	!Conservation de l'eau
3823	!! sumdq = 0.
3824	!! DO k = 1, nl
3825	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3826	!! END DO
3827	!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3828	!jyg>
3829
3830
3831	! *** Check that moisture stays positive. If not, scale tendencies
3832	! in order to ensure moisture positivity
3833	DO il = 1, ncum
3834	alpha_qpos(il) = 1.
3835	IF (iflag(il)<=1) THEN
3836	IF (fr(il,1)<=0.) THEN
3837	alpha_qpos(il) = max(alpha_qpos(il), (-deltfr(il,1))/(s_wake(il)rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
3838	END IF
3839	END IF
3840	END DO
3841	DO i = 2, nl
3842	DO il = 1, ncum
3843	IF (iflag(il)<=1) THEN
3844	IF (fr(il,i)<=0.) THEN
3845	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
3846	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
3847	END IF
3848	END IF
3849	END DO
3850	END DO
3851	DO il = 1, ncum
3852	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
3853	alpha_qpos(il) = alpha_qpos(il)*1.1
3854	END IF
3855	END DO
3856	!
3857	! print *,' YIELD : alpha_qpos ',alpha_qpos(1)
3858	!
3859	DO il = 1, ncum
3860	IF (iflag(il)<=1) THEN
3861	sigd(il) = sigd(il)/alpha_qpos(il)
3862	precip(il) = precip(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