Context Navigation

cv3_routines.F90 @ 2787

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

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

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

Download in other formats:

Original Format