Context Navigation

cv3_routines.F90 @ 2253

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

1/ Introduction of two variables in the ".def" files: (i) cvl_sig2feed is
the top of the convective feeding layer in sigma coordinates (D=0.97);
(ii) cvl_comp_threshold is the threshold fraction of convective points
below which compression occurs (D=1.).
2/ Corrections of various bugs revealed by the changes in compression:

correct bugs in cv3a_uncompress.F90 for 3 fields used for convective

scavenging.

add a reset to zero of "sig" and "w0" for non-convective points

(cva_driver.F90).

in cv3_routines.F90, correct bounds of a few loops in cv3_undilute2,

correct the reset of the no-convection counter in cv3_yield.

in phys_output_write_mod.F90, correct output of wdtrainA and wdtrainM.

3/ Improve declarations in various subroutines.

Modified files:

conema3.h
cv3param.h
cv3p1_closure.F90
conf_phys_m.F90
cv3a_compress.F90
phys_output_write_mod.F90
cv3_routines.F90
concvl.F90
cva_driver.F90
cv3a_uncompress.F90

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: 128.5 KB

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

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

Download in other formats:

Original Format