Context Navigation

cv3_routines.F90 @ 2435

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

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

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

Download in other formats:

Original Format