Context Navigation

cv3_routines.F90 @ 2759

Last change on this file since 2759 was 2759, checked in by jyg, 8 years ago
Small improvements to convection scheme (cv3_unsat) and to wake scheme.
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.2 KB

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

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

Download in other formats:

Original Format