Context Navigation

cv3_routines.F90 @ 2306

Last change on this file since 2306 was 2306, checked in by jyg, 9 years ago
Improved diagnostics: pmflxr and pmflxs are now the true vertical profiles of liquid and solid convective precipitation (previously they where merely diagnosed through a test on temperature). The convective scheme internal variable for ice precipitation is Vprecipi.
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: 130.6 KB

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

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

Download in other formats:

Original Format