Context Navigation

cv3_routines.F90 @ 2459

Last change on this file since 2459 was 2459, checked in by jyg, 8 years ago
Fix to the previous commit (Yacine's optimisation was included by mistake).
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.4 KB

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

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

Download in other formats:

Original Format