Context Navigation

cv3_routines.F90 @ 2376

Last change on this file since 2376 was 2376, checked in by jyg, 9 years ago
Add initialization of ice fraction below cloud base in cv3_undilute2 in cv3_routines.F90.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 132.3 KB

Line
1
2	! $Id: cv3_routines.F90 2376 2015-10-20 09:57:39Z 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	USE print_control_mod, ONLY: lunout
894	IMPLICIT NONE
895
896	include "cv3param.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_physic(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	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: frac
1077
1078	!local variables:
1079	INTEGER i, j, k
1080	REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
1081	REAL als
1082	REAL qsat_new, snew, qi(nloc, nd)
1083	REAL by, defrac, pden, tbis
1084	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
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	! Below cloud base, set ice fraction to cloud base value
1538	DO k = 1, nl
1539	DO i = 1, ncum
1540	IF (k<icb(i)) THEN
1541	frac(i,k) = frac(i,icb(i))
1542	END IF
1543	END DO
1544	END DO
1545	!
1546	ELSE
1547	!
1548	DO k = minorig + 1, nl
1549	DO i = 1, ncum
1550	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1551	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k)
1552	END IF
1553	END DO
1554	END DO
1555	!
1556	END IF ! (cvflag_ice)
1557
1558	RETURN
1559	END SUBROUTINE cv3_undilute2
1560
1561	SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
1562	pbase, p, ph, tv, buoy, &
1563	sig, w0, cape, m, iflag)
1564	IMPLICIT NONE
1565
1566	! ===================================================================
1567	! --- CLOSURE OF CONVECT3
1568	!
1569	! vectorization: S. Bony
1570	! ===================================================================
1571
1572	include "cvthermo.h"
1573	include "cv3param.h"
1574
1575	!input:
1576	INTEGER ncum, nd, nloc
1577	INTEGER icb(nloc), inb(nloc)
1578	REAL pbase(nloc)
1579	REAL p(nloc, nd), ph(nloc, nd+1)
1580	REAL tv(nloc, nd), buoy(nloc, nd)
1581
1582	!input/output:
1583	REAL sig(nloc, nd), w0(nloc, nd)
1584	INTEGER iflag(nloc)
1585
1586	!output:
1587	REAL cape(nloc)
1588	REAL m(nloc, nd)
1589
1590	!local variables:
1591	INTEGER i, j, k, icbmax
1592	REAL deltap, fac, w, amu
1593	REAL dtmin(nloc, nd), sigold(nloc, nd)
1594	REAL cbmflast(nloc)
1595
1596
1597	! -------------------------------------------------------
1598	! -- Initialization
1599	! -------------------------------------------------------
1600
1601	DO k = 1, nl
1602	DO i = 1, ncum
1603	m(i, k) = 0.0
1604	END DO
1605	END DO
1606
1607	! -------------------------------------------------------
1608	! -- Reset sig(i) and w0(i) for i>inb and i<icb
1609	! -------------------------------------------------------
1610
1611	! update sig and w0 above LNB:
1612
1613	DO k = 1, nl - 1
1614	DO i = 1, ncum
1615	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
1616	sig(i, k) = beta*sig(i, k) + &
1617	2.alphabuoy(i, inb(i))*abs(buoy(i,inb(i)))
1618	sig(i, k) = amax1(sig(i,k), 0.0)
1619	w0(i, k) = beta*w0(i, k)
1620	END IF
1621	END DO
1622	END DO
1623
1624	! compute icbmax:
1625
1626	icbmax = 2
1627	DO i = 1, ncum
1628	icbmax = max(icbmax, icb(i))
1629	END DO
1630
1631	! update sig and w0 below cloud base:
1632
1633	DO k = 1, icbmax
1634	DO i = 1, ncum
1635	IF (k<=icb(i)) THEN
1636	sig(i, k) = beta*sig(i, k) - &
1637	2.alphabuoy(i, icb(i))*buoy(i, icb(i))
1638	sig(i, k) = max(sig(i,k), 0.0)
1639	w0(i, k) = beta*w0(i, k)
1640	END IF
1641	END DO
1642	END DO
1643
1644	!! if(inb.lt.(nl-1))then
1645	!! do 85 i=inb+1,nl-1
1646	!! sig(i)=betasig(i)+2.alphabuoy(inb)
1647	!! 1 abs(buoy(inb))
1648	!! sig(i)=max(sig(i),0.0)
1649	!! w0(i)=beta*w0(i)
1650	!! 85 continue
1651	!! end if
1652
1653	!! do 87 i=1,icb
1654	!! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
1655	!! sig(i)=max(sig(i),0.0)
1656	!! w0(i)=beta*w0(i)
1657	!! 87 continue
1658
1659	! -------------------------------------------------------------
1660	! -- Reset fractional areas of updrafts and w0 at initial time
1661	! -- and after 10 time steps of no convection
1662	! -------------------------------------------------------------
1663
1664	DO k = 1, nl - 1
1665	DO i = 1, ncum
1666	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
1667	sig(i, k) = 0.0
1668	w0(i, k) = 0.0
1669	END IF
1670	END DO
1671	END DO
1672
1673	! -------------------------------------------------------------
1674	! -- Calculate convective available potential energy (cape),
1675	! -- vertical velocity (w), fractional area covered by
1676	! -- undilute updraft (sig), and updraft mass flux (m)
1677	! -------------------------------------------------------------
1678
1679	DO i = 1, ncum
1680	cape(i) = 0.0
1681	END DO
1682
1683	! compute dtmin (minimum buoyancy between ICB and given level k):
1684
1685	DO i = 1, ncum
1686	DO k = 1, nl
1687	dtmin(i, k) = 100.0
1688	END DO
1689	END DO
1690
1691	DO i = 1, ncum
1692	DO k = 1, nl
1693	DO j = minorig, nl
1694	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
1695	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
1696	END IF
1697	END DO
1698	END DO
1699	END DO
1700
1701	! the interval on which cape is computed starts at pbase :
1702
1703	DO k = 1, nl
1704	DO i = 1, ncum
1705
1706	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
1707
1708	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
1709	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
1710	cape(i) = amax1(0.0, cape(i))
1711	sigold(i, k) = sig(i, k)
1712
1713	! dtmin(i,k)=100.0
1714	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
1715	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1716	! 97 continue
1717
1718	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
1719	sig(i, k) = max(sig(i,k), 0.0)
1720	sig(i, k) = amin1(sig(i,k), 0.01)
1721	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
1722	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
1723	amu = 0.5(sig(i,k)+sigold(i,k))w
1724	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
1725	w0(i, k) = w
1726	END IF
1727
1728	END DO
1729	END DO
1730
1731	DO i = 1, ncum
1732	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
1733	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))
1734	sig(i, icb(i)) = sig(i, icb(i)+1)
1735	sig(i, icb(i)-1) = sig(i, icb(i))
1736	END DO
1737
1738	! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
1739	! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
1740	! ccc the final mass flux (cbmflast) is greater than the target mass flux
1741	! ccc (cbmf) ??).
1742	! cc
1743	! c do i = 1,ncum
1744	! c cbmflast(i) = 0.
1745	! c enddo
1746	! cc
1747	! c do k= 1,nl
1748	! c do i = 1,ncum
1749	! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
1750	! c cbmflast(i) = cbmflast(i)+M(i,k)
1751	! c ENDIF
1752	! c enddo
1753	! c enddo
1754	! cc
1755	! c do i = 1,ncum
1756	! c IF (cbmflast(i) .lt. 1.e-6) THEN
1757	! c iflag(i) = 3
1758	! c ENDIF
1759	! c enddo
1760	! cc
1761	! c do k= 1,nl
1762	! c do i = 1,ncum
1763	! c IF (iflag(i) .ge. 3) THEN
1764	! c M(i,k) = 0.
1765	! c sig(i,k) = 0.
1766	! c w0(i,k) = 0.
1767	! c ENDIF
1768	! c enddo
1769	! c enddo
1770	! cc
1771	!! cape=0.0
1772	!! do 98 i=icb+1,inb
1773	!! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1774	!! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
1775	!! dcape=rrdbuoy(i-1)deltap/p(i-1)
1776	!! dlnp=deltap/p(i-1)
1777	!! cape=max(0.0,cape)
1778	!! sigold=sig(i)
1779
1780	!! dtmin=100.0
1781	!! do 97 j=icb,i-1
1782	!! dtmin=amin1(dtmin,buoy(j))
1783	!! 97 continue
1784
1785	!! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
1786	!! sig(i)=max(sig(i),0.0)
1787	!! sig(i)=amin1(sig(i),0.01)
1788	!! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1789	!! w=(1.-beta)facsqrt(cape)+beta*w0(i)
1790	!! amu=0.5(sig(i)+sigold)w
1791	!! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
1792	!! w0(i)=w
1793	!! 98 continue
1794	!! w0(icb)=0.5*w0(icb+1)
1795	!! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1796	!! sig(icb)=sig(icb+1)
1797	!! sig(icb-1)=sig(icb)
1798
1799	RETURN
1800	END SUBROUTINE cv3_closure
1801
1802	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
1803	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
1804	unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
1805	ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
1806	IMPLICIT NONE
1807
1808	! ---------------------------------------------------------------------
1809	! a faire:
1810	! - vectorisation de la partie normalisation des flux (do 789...)
1811	! ---------------------------------------------------------------------
1812
1813	include "cvthermo.h"
1814	include "cv3param.h"
1815	include "cvflag.h"
1816
1817	!inputs:
1818	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
1819	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
1820	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
1821	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
1822	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
1823	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
1824	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
1825	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
1826	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, h, hp
1827	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf, frac
1828	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp, ep, clw
1829	REAL, DIMENSION (nloc, na), INTENT (IN) :: m ! input of convect3
1830
1831	!outputs:
1832	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: ment, qent
1833	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
1834	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: sij, elij
1835	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
1836	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: ments, qents
1837	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
1838
1839	!local variables:
1840	INTEGER i, j, k, il, im, jm
1841	INTEGER num1, num2
1842	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1843	REAL alt, smid, sjmin, sjmax, delp, delm
1844	REAL asij(nloc), smax(nloc), scrit(nloc)
1845	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
1846	REAL sigij(nloc, nd, nd)
1847	REAL wgh
1848	REAL zm(nloc, na)
1849	LOGICAL lwork(nloc)
1850
1851	! =====================================================================
1852	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1853	! =====================================================================
1854
1855	! ori do 360 i=1,ncum*nlp
1856	DO j = 1, nl
1857	DO i = 1, ncum
1858	nent(i, j) = 0
1859	! in convect3, m is computed in cv3_closure
1860	! ori m(i,1)=0.0
1861	END DO
1862	END DO
1863
1864	! ori do 400 k=1,nlp
1865	! ori do 390 j=1,nlp
1866	DO j = 1, nl
1867	DO k = 1, nl
1868	DO i = 1, ncum
1869	qent(i, k, j) = rr(i, j)
1870	uent(i, k, j) = u(i, j)
1871	vent(i, k, j) = v(i, j)
1872	elij(i, k, j) = 0.0
1873	!ym ment(i,k,j)=0.0
1874	!ym sij(i,k,j)=0.0
1875	END DO
1876	END DO
1877	END DO
1878
1879	!ym
1880	ment(1:ncum, 1:nd, 1:nd) = 0.0
1881	sij(1:ncum, 1:nd, 1:nd) = 0.0
1882
1883	!AC! do k=1,ntra
1884	!AC! do j=1,nd ! instead nlp
1885	!AC! do i=1,nd ! instead nlp
1886	!AC! do il=1,ncum
1887	!AC! traent(il,i,j,k)=tra(il,j,k)
1888	!AC! enddo
1889	!AC! enddo
1890	!AC! enddo
1891	!AC! enddo
1892	zm(:, :) = 0.
1893
1894	! =====================================================================
1895	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1896	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1897	! --- FRACTION (sij)
1898	! =====================================================================
1899
1900	DO i = minorig + 1, nl
1901
1902	DO j = minorig, nl
1903	DO il = 1, ncum
1904	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
1905
1906	rti = qnk(il) - ep(il, i)*clw(il, i)
1907	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
1908
1909
1910	IF (cvflag_ice) THEN
1911	! print*,cvflag_ice,'cvflag_ice dans do 700'
1912	IF (t(il,j)<=263.15) THEN
1913	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
1914	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
1915	END IF
1916	END IF
1917
1918	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
1919	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
1920	dei = denom
1921	IF (abs(dei)<0.01) dei = 0.01
1922	sij(il, i, j) = anum/dei
1923	sij(il, i, i) = 1.0
1924	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
1925	altem = altem/bf2
1926	cwat = clw(il, j)*(1.-ep(il,j))
1927	stemp = sij(il, i, j)
1928	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
1929
1930	IF (cvflag_ice) THEN
1931	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
1932	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
1933	ELSE
1934	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
1935	denom = denom + lv(il, j)*(rr(il,i)-rti)
1936	END IF
1937
1938	IF (abs(denom)<0.01) denom = 0.01
1939	sij(il, i, j) = anum/denom
1940	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
1941	altem = altem - (bf2-1.)*cwat
1942	END IF
1943	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
1944	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
1945	uent(il, i, j) = sij(il, i, j)u(il, i) + (1.-sij(il,i,j))unk(il)
1946	vent(il, i, j) = sij(il, i, j)v(il, i) + (1.-sij(il,i,j))vnk(il)
1947	!!!! do k=1,ntra
1948	!!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1949	!!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1950	!!!! end do
1951	elij(il, i, j) = altem
1952	elij(il, i, j) = max(0.0, elij(il,i,j))
1953	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
1954	nent(il, i) = nent(il, i) + 1
1955	END IF
1956	sij(il, i, j) = max(0.0, sij(il,i,j))
1957	sij(il, i, j) = amin1(1.0, sij(il,i,j))
1958	END IF ! new
1959	END DO
1960	END DO
1961
1962	!AC! do k=1,ntra
1963	!AC! do j=minorig,nl
1964	!AC! do il=1,ncum
1965	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1966	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1967	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1968	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1969	!AC! endif
1970	!AC! enddo
1971	!AC! enddo
1972	!AC! enddo
1973
1974
1975	! * if no air can entrain at level i assume that updraft detrains *
1976	! * at that level and calculate detrained air flux and properties *
1977
1978
1979	! @ do 170 i=icb(il),inb(il)
1980
1981	DO il = 1, ncum
1982	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
1983	! @ if(nent(il,i).eq.0)then
1984	ment(il, i, i) = m(il, i)
1985	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
1986	uent(il, i, i) = unk(il)
1987	vent(il, i, i) = vnk(il)
1988	elij(il, i, i) = clw(il, i)
1989	! MAF sij(il,i,i)=1.0
1990	sij(il, i, i) = 0.0
1991	END IF
1992	END DO
1993	END DO
1994
1995	!AC! do j=1,ntra
1996	!AC! do i=minorig+1,nl
1997	!AC! do il=1,ncum
1998	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1999	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2000	!AC! endif
2001	!AC! enddo
2002	!AC! enddo
2003	!AC! enddo
2004
2005	DO j = minorig, nl
2006	DO i = minorig, nl
2007	DO il = 1, ncum
2008	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
2009	sigij(il, i, j) = sij(il, i, j)
2010	END IF
2011	END DO
2012	END DO
2013	END DO
2014	! @ enddo
2015
2016	! @170 continue
2017
2018	! =====================================================================
2019	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
2020	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
2021	! =====================================================================
2022
2023	CALL zilch(asum, nloc*nd)
2024	CALL zilch(csum, nloc*nd)
2025	CALL zilch(csum, nloc*nd)
2026
2027	DO il = 1, ncum
2028	lwork(il) = .FALSE.
2029	END DO
2030
2031	DO i = minorig + 1, nl
2032
2033	num1 = 0
2034	DO il = 1, ncum
2035	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
2036	END DO
2037	IF (num1<=0) GO TO 789
2038
2039
2040	DO il = 1, ncum
2041	IF (i>=icb(il) .AND. i<=inb(il)) THEN
2042	lwork(il) = (nent(il,i)/=0)
2043	qp = qnk(il) - ep(il, i)*clw(il, i)
2044
2045	IF (cvflag_ice) THEN
2046
2047	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
2048	(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
2049	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
2050	(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
2051	ELSE
2052
2053	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
2054	(cpv-cpd)t(il, i)(qp-rr(il,i))
2055	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
2056	(cpd-cpv)t(il, i)(rr(il,i)-qp)
2057	END IF
2058
2059	IF (abs(denom)<0.01) denom = 0.01
2060	scrit(il) = anum/denom
2061	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
2062	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
2063	smax(il) = 0.0
2064	asij(il) = 0.0
2065	END IF
2066	END DO
2067
2068	DO j = nl, minorig, -1
2069
2070	num2 = 0
2071	DO il = 1, ncum
2072	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2073	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2074	lwork(il)) num2 = num2 + 1
2075	END DO
2076	IF (num2<=0) GO TO 175
2077
2078	DO il = 1, ncum
2079	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2080	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2081	lwork(il)) THEN
2082
2083	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
2084	wgh = 1.0
2085	IF (j>i) THEN
2086	sjmax = max(sij(il,i,j+1), smax(il))
2087	sjmax = amin1(sjmax, scrit(il))
2088	smax(il) = max(sij(il,i,j), smax(il))
2089	sjmin = max(sij(il,i,j-1), smax(il))
2090	sjmin = amin1(sjmin, scrit(il))
2091	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
2092	smid = amin1(sij(il,i,j), scrit(il))
2093	ELSE
2094	sjmax = max(sij(il,i,j+1), scrit(il))
2095	smid = max(sij(il,i,j), scrit(il))
2096	sjmin = 0.0
2097	IF (j>1) sjmin = sij(il, i, j-1)
2098	sjmin = max(sjmin, scrit(il))
2099	END IF
2100	delp = abs(sjmax-smid)
2101	delm = abs(sjmin-smid)
2102	asij(il) = asij(il) + wgh*(delp+delm)
2103	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
2104	END IF
2105	END IF
2106	END DO
2107
2108	175 END DO
2109
2110	DO il = 1, ncum
2111	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2112	asij(il) = max(1.0E-16, asij(il))
2113	asij(il) = 1.0/asij(il)
2114	asum(il, i) = 0.0
2115	bsum(il, i) = 0.0
2116	csum(il, i) = 0.0
2117	END IF
2118	END DO
2119
2120	DO j = minorig, nl
2121	DO il = 1, ncum
2122	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2123	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2124	ment(il, i, j) = ment(il, i, j)*asij(il)
2125	END IF
2126	END DO
2127	END DO
2128
2129	DO j = minorig, nl
2130	DO il = 1, ncum
2131	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2132	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2133	asum(il, i) = asum(il, i) + ment(il, i, j)
2134	ment(il, i, j) = ment(il, i, j)*sig(il, j)
2135	bsum(il, i) = bsum(il, i) + ment(il, i, j)
2136	END IF
2137	END DO
2138	END DO
2139
2140	DO il = 1, ncum
2141	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2142	bsum(il, i) = max(bsum(il,i), 1.0E-16)
2143	bsum(il, i) = 1.0/bsum(il, i)
2144	END IF
2145	END DO
2146
2147	DO j = minorig, nl
2148	DO il = 1, ncum
2149	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2150	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2151	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
2152	END IF
2153	END DO
2154	END DO
2155
2156	DO j = minorig, nl
2157	DO il = 1, ncum
2158	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2159	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2160	csum(il, i) = csum(il, i) + ment(il, i, j)
2161	END IF
2162	END DO
2163	END DO
2164
2165	DO il = 1, ncum
2166	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2167	csum(il,i)<m(il,i)) THEN
2168	nent(il, i) = 0
2169	ment(il, i, i) = m(il, i)
2170	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2171	uent(il, i, i) = unk(il)
2172	vent(il, i, i) = vnk(il)
2173	elij(il, i, i) = clw(il, i)
2174	! MAF sij(il,i,i)=1.0
2175	sij(il, i, i) = 0.0
2176	END IF
2177	END DO ! il
2178
2179	!AC! do j=1,ntra
2180	!AC! do il=1,ncum
2181	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
2182	!AC! : .and. csum(il,i).lt.m(il,i) ) then
2183	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2184	!AC! endif
2185	!AC! enddo
2186	!AC! enddo
2187	789 END DO
2188
2189	! MAF: renormalisation de MENT
2190	CALL zilch(zm, nloc*na)
2191	DO jm = 1, nd
2192	DO im = 1, nd
2193	DO il = 1, ncum
2194	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
2195	END DO
2196	END DO
2197	END DO
2198
2199	DO jm = 1, nd
2200	DO im = 1, nd
2201	DO il = 1, ncum
2202	IF (zm(il,im)/=0.) THEN
2203	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
2204	END IF
2205	END DO
2206	END DO
2207	END DO
2208
2209	DO jm = 1, nd
2210	DO im = 1, nd
2211	DO il = 1, ncum
2212	qents(il, im, jm) = qent(il, im, jm)
2213	ments(il, im, jm) = ment(il, im, jm)
2214	END DO
2215	END DO
2216	END DO
2217
2218	RETURN
2219	END SUBROUTINE cv3_mixing
2220
2221	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
2222	t, rr, rs, gz, u, v, tra, p, ph, &
2223	th, tv, lv, lf, cpn, ep, sigp, clw, &
2224	m, ment, elij, delt, plcl, coef_clos, &
2225	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
2226	faci, b, sigd, &
2227	wdtrainA, wdtrainM) ! RomP
2228	IMPLICIT NONE
2229
2230
2231	include "cvthermo.h"
2232	include "cv3param.h"
2233	include "cvflag.h"
2234	include "nuage.h"
2235
2236	!inputs:
2237	INTEGER ncum, nd, na, ntra, nloc
2238	INTEGER icb(nloc), inb(nloc)
2239	REAL delt, plcl(nloc)
2240	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd), gz(nloc, na)
2241	REAL u(nloc, nd), v(nloc, nd)
2242	REAL tra(nloc, nd, ntra)
2243	REAL p(nloc, nd), ph(nloc, nd+1)
2244	REAL ep(nloc, na), sigp(nloc, na), clw(nloc, na)
2245	REAL th(nloc, na), tv(nloc, na), lv(nloc, na), cpn(nloc, na)
2246	REAL lf(nloc, na)
2247	REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
2248	REAL coef_clos(nloc)
2249
2250	!input/output
2251	INTEGER iflag(nloc)
2252
2253	!outputs:
2254	REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
2255	REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
2256	REAL ice(nloc, na), fondue(nloc, na), faci(nloc, na)
2257	REAL trap(nloc, na, ntra)
2258	REAL b(nloc, na), sigd(nloc)
2259	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
2260	! de l ascendance adiabatique et des flux melanges Pa et Pm.
2261	! Distinction des wdtrain
2262	! Pa = wdtrainA Pm = wdtrainM
2263	REAL wdtrainA(nloc, na), wdtrainM(nloc, na)
2264
2265	!local variables
2266	INTEGER i, j, k, il, num1, ndp1
2267	REAL tinv, delti, coef
2268	REAL awat, afac, afac1, afac2, bfac
2269	REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
2270	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
2271	REAL ampmax, thaw
2272	REAL tevap(nloc)
2273	REAL lvcp(nloc, na), lfcp(nloc, na)
2274	REAL h(nloc, na), hm(nloc, na)
2275	REAL frac(nloc, na)
2276	REAL fraci(nloc, na), prec(nloc, na)
2277	REAL wdtrain(nloc)
2278	LOGICAL lwork(nloc), mplus(nloc)
2279
2280
2281	! ------------------------------------------------------
2282
2283	delti = 1./delt
2284	tinv = 1./3.
2285
2286	mp(:, :) = 0.
2287
2288	DO i = 1, nl
2289	DO il = 1, ncum
2290	mp(il, i) = 0.0
2291	rp(il, i) = rr(il, i)
2292	up(il, i) = u(il, i)
2293	vp(il, i) = v(il, i)
2294	wt(il, i) = 0.001
2295	water(il, i) = 0.0
2296	frac(il, i) = 0.0
2297	faci(il, i) = 0.0
2298	fraci(il, i) = 0.0
2299	ice(il, i) = 0.0
2300	prec(il, i) = 0.0
2301	fondue(il, i) = 0.0
2302	evap(il, i) = 0.0
2303	b(il, i) = 0.0
2304	lvcp(il, i) = lv(il, i)/cpn(il, i)
2305	lfcp(il, i) = lf(il, i)/cpn(il, i)
2306	END DO
2307	END DO
2308	!AC! do k=1,ntra
2309	!AC! do i=1,nd
2310	!AC! do il=1,ncum
2311	!AC! trap(il,i,k)=tra(il,i,k)
2312	!AC! enddo
2313	!AC! enddo
2314	!AC! enddo
2315	!! RomP >>>
2316	DO i = 1, nd
2317	DO il = 1, ncum
2318	wdtrainA(il, i) = 0.0
2319	wdtrainM(il, i) = 0.0
2320	END DO
2321	END DO
2322	!! RomP <<<
2323
2324	! * check whether ep(inb)=0, if so, skip precipitating *
2325	! * downdraft calculation *
2326
2327
2328	DO il = 1, ncum
2329	!! lwork(il)=.TRUE.
2330	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
2331	lwork(il) = ep(il, inb(il)) >= 0.0001
2332	END DO
2333
2334	! *** Set the fractionnal area sigd of precipitating downdraughts
2335	DO il = 1, ncum
2336	sigd(il) = sigdz*coef_clos(il)
2337	END DO
2338
2339
2340	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2341	!
2342	! * begin downdraft loop *
2343	!
2344	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2345
2346	DO i = nl + 1, 1, -1
2347
2348	num1 = 0
2349	DO il = 1, ncum
2350	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
2351	END DO
2352	IF (num1<=0) GO TO 400
2353
2354	CALL zilch(wdtrain, ncum)
2355
2356
2357	! * integrate liquid water equation to find condensed water *
2358	! * and condensed water flux *
2359	!
2360	!
2361	! * calculate detrained precipitation *
2362
2363	DO il = 1, ncum
2364	IF (i<=inb(il) .AND. lwork(il)) THEN
2365	IF (cvflag_grav) THEN
2366	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
2367	wdtrainA(il, i) = wdtrain(il)/grav ! Pa RomP
2368	ELSE
2369	wdtrain(il) = 10.0ep(il, i)m(il, i)*clw(il, i)
2370	wdtrainA(il, i) = wdtrain(il)/10. ! Pa RomP
2371	END IF
2372	END IF
2373	END DO
2374
2375	IF (i>1) THEN
2376	DO j = 1, i - 1
2377	DO il = 1, ncum
2378	IF (i<=inb(il) .AND. lwork(il)) THEN
2379	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
2380	awat = max(awat, 0.0)
2381	IF (cvflag_grav) THEN
2382	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
2383	wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP
2384	ELSE
2385	wdtrain(il) = wdtrain(il) + 10.0awatment(il, j, i)
2386	wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i) ! Pm RomP
2387	END IF
2388	END IF
2389	END DO
2390	END DO
2391	END IF
2392
2393
2394	! * find rain water and evaporation using provisional *
2395	! * estimates of rp(i)and rp(i-1) *
2396
2397
2398	DO il = 1, ncum
2399	IF (i<=inb(il) .AND. lwork(il)) THEN
2400
2401	wt(il, i) = 45.0
2402
2403	IF (cvflag_ice) THEN
2404	frac(il, inb(il)) = 1. - (t(il,inb(il))-243.15)/(263.15-243.15)
2405	frac(il, inb(il)) = min(max(frac(il,inb(il)),0.), 1.)
2406	fraci(il, inb(il)) = frac(il, inb(il))
2407	ELSE
2408	CONTINUE
2409	END IF
2410
2411	IF (i<inb(il)) THEN
2412
2413	IF (cvflag_ice) THEN
2414	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
2415	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
2416	thaw = min(max(thaw,0.0), 1.0)
2417	frac(il, i) = frac(il, i)*(1.-thaw)
2418	ELSE
2419	CONTINUE
2420	END IF
2421
2422	rp(il, i) = rp(il, i+1) + &
2423	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
2424	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
2425	END IF
2426	fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
2427	fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
2428	rp(il, i) = max(rp(il,i), 0.0)
2429	rp(il, i) = amin1(rp(il,i), rs(il,i))
2430	rp(il, inb(il)) = rr(il, inb(il))
2431
2432	IF (i==1) THEN
2433	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
2434	IF (cvflag_ice) THEN
2435	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
2436	END IF
2437	ELSE
2438	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)
2439	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
2440	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
2441	rp(il, i-1) = max(rp(il,i-1), 0.0)
2442	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
2443	afac2 = p(il, i-1)(rs(il,i-1)-rp(il,i-1))/(1.0E4+2000.0p(il,i-1)*rs(il,i-1))
2444	afac = 0.5*(afac1+afac2)
2445	END IF
2446	IF (i==inb(il)) afac = 0.0
2447	afac = max(afac, 0.0)
2448	bfac = 1./(sigd(il)*wt(il,i))
2449
2450	!JYG1
2451	! cc sigt=1.0
2452	! cc if(i.ge.icb)sigt=sigp(i)
2453	! prise en compte de la variation progressive de sigt dans
2454	! les couches icb et icb-1:
2455	! pour plcl<ph(i+1), pr1=0 & pr2=1
2456	! pour plcl>ph(i), pr1=1 & pr2=0
2457	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2458	! sur le nuage, et pr2 est la proportion sous la base du
2459	! nuage.
2460	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2461	pr1 = max(0., min(1.,pr1))
2462	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2463	pr2 = max(0., min(1.,pr2))
2464	sigt = sigp(il, i)*pr1 + pr2
2465	!JYG2
2466
2467	!JYG----
2468	! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2469	! c6 = water(il,i+1) + wdtrain(il)*bfac
2470	! c6 = prec(il,i+1) + wdtrain(il)*bfac
2471	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
2472	! evap(il,i)=sigtafacrevap
2473	! water(il,i)=revap*revap
2474	! prec(il,i)=revap*revap
2475	!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
2476	!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
2477	!!---end jyg---
2478
2479	! --------retour à la formulation originale d'Emanuel.
2480	IF (cvflag_ice) THEN
2481
2482	! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2483	! c6=prec(il,i+1)+bfac*wdtrain(il) &
2484	! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2485	! if(c6.gt.0.0)then
2486	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
2487
2488	!JAM Attention: evap=sigt*E
2489	! Modification: evap devient l'évaporation en milieu de couche
2490	! car nécessaire dans cv3_yield
2491	! Du coup, il faut modifier pas mal d'équations...
2492	! et l'expression de afac qui devient afac1
2493	! revap=sqrt((prec(i+1)+prec(i))/2)
2494
2495	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
2496	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
2497	! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
2498	! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
2499	! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
2500	IF (c6>b6*b6+1.E-20) THEN
2501	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
2502	ELSE
2503	revap = (-b6+sqrt(b6b6+4.c6))/2.
2504	END IF
2505	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
2506	! print*,prec(il,i),'neige'
2507
2508	!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
2509	! c evap(il,i)=sigtafacrevap
2510	! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
2511	! Ici,l'evaporation evap est simplement calculee par l'equation de
2512	! conservation.
2513	! prec(il,i)=revap*revap
2514	! else
2515	!JYG---- Correction : si c6 <= 0, water(il,i)=0.
2516	! prec(il,i)=0.
2517	! endif
2518
2519	!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
2520	! moins [tt ce qui sort de la couche i]
2521	! print *, 'evap avec ice'
2522	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
2523	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2524
2525	d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
2526	e6 = bfac*wdtrain(il)
2527	f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
2528
2529	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
2530	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
2531	thaw = min(max(thaw,0.0), 1.0)
2532	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
2533	water(il, i) = max(water(il,i), 0.)
2534	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
2535	ice(il, i) = max(ice(il,i), 0.)
2536	fondue(il, i) = ice(il, i)*thaw
2537	water(il, i) = water(il, i) + fondue(il, i)
2538	ice(il, i) = ice(il, i) - fondue(il, i)
2539
2540	IF (water(il,i)+ice(il,i)<1.E-30) THEN
2541	faci(il, i) = 0.
2542	ELSE
2543	faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
2544	END IF
2545
2546	! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
2547	! water(il,i)=max(water(il,i),0.)
2548	! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
2549	! ice(il,i)=max(ice(il,i),0.)
2550	! fondue(il,i)=ice(il,i)*thaw
2551	! water(il,i)=water(il,i)+fondue(il,i)
2552	! ice(il,i)=ice(il,i)-fondue(il,i)
2553
2554	! if((water(il,i)+ice(il,i)).lt.1.e-30)then
2555	! faci(il,i)=0.
2556	! else
2557	! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
2558	! endif
2559
2560	ELSE
2561	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2562	c6 = water(il, i+1) + bfac*wdtrain(il) - &
2563	50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
2564	IF (c6>0.0) THEN
2565	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
2566	water(il, i) = revap*revap
2567	ELSE
2568	water(il, i) = 0.
2569	END IF
2570	! print *, 'evap sans ice'
2571	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
2572	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2573
2574	END IF
2575	END IF !(i.le.inb(il) .and. lwork(il))
2576	END DO
2577	! ----------------------------------------------------------------
2578
2579	! cc
2580	! * calculate precipitating downdraft mass flux under *
2581	! * hydrostatic approximation *
2582
2583	DO il = 1, ncum
2584	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
2585
2586	tevap(il) = max(0.0, evap(il,i))
2587	delth = max(0.001, (th(il,i)-th(il,i-1)))
2588	IF (cvflag_ice) THEN
2589	IF (cvflag_grav) THEN
2590	mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
2591	(p(il,i-1)-p(il,i))/delth + &
2592	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
2593	(p(il,i-1)-p(il,i))/delth + &
2594	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
2595	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
2596	ELSE
2597	mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
2598	(p(il,i-1)-p(il,i))/delth + &
2599	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
2600	(p(il,i-1)-p(il,i))/delth + &
2601	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
2602	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
2603
2604	END IF
2605	ELSE
2606	IF (cvflag_grav) THEN
2607	mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
2608	(p(il,i-1)-p(il,i))/delth
2609	ELSE
2610	mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
2611	(p(il,i-1)-p(il,i))/delth
2612	END IF
2613
2614	END IF
2615
2616	END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2617	END DO
2618	! ----------------------------------------------------------------
2619
2620	! * if hydrostatic assumption fails, *
2621	! * solve cubic difference equation for downdraft theta *
2622	! * and mass flux from two simultaneous differential eqns *
2623
2624	DO il = 1, ncum
2625	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
2626
2627	amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
2628	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2629	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2630
2631	IF (amp2>(0.1*amfac)) THEN
2632	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
2633	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
2634	(lvcp(il,i)sigd(il)th(il,i))
2635	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
2636
2637	IF (cvflag_ice) THEN
2638	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2639	50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
2640	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
2641	ELSE
2642
2643	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2644	50.(p(il,i-1)-p(il,i))xf*tevap(il)
2645	END IF
2646
2647	fac2 = 1.0
2648	IF (bf<0.0) fac2 = -1.0
2649	bf = abs(bf)
2650	ur = 0.25bfbf - afafaftinvtinv*tinv
2651	IF (ur>=0.0) THEN
2652	sru = sqrt(ur)
2653	fac = 1.0
2654	IF ((0.5*bf-sru)<0.0) fac = -1.0
2655	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
2656	fac(abs(0.5bf-sru))**tinv
2657	ELSE
2658	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
2659	IF (fac2<0.0) d = 3.14159 - d
2660	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
2661	END IF
2662	mp(il, i) = max(0.0, mp(il,i))
2663
2664	IF (cvflag_ice) THEN
2665	IF (cvflag_grav) THEN
2666	!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2667	! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
2668	! Et il faut bien revoir les facteurs 100.
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	ELSE
2677	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
2678	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
2679	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
2680	(ph(il,i)-ph(il,i+1))) / &
2681	(mp(il,i)+sigd(il)*0.1) - &
2682	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2683	(lvcp(il,i)sigd(il)th(il,i))
2684	END IF
2685	ELSE
2686	IF (cvflag_grav) THEN
2687	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
2688	(mp(il,i)+sigd(il)*0.1) - &
2689	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2690	(lvcp(il,i)sigd(il)th(il,i))
2691	ELSE
2692	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
2693	(mp(il,i)+sigd(il)*0.1) - &
2694	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2695	(lvcp(il,i)sigd(il)th(il,i))
2696	END IF
2697	END IF
2698	b(il, i-1) = max(b(il,i-1), 0.0)
2699
2700	END IF !(amp2.gt.(0.1*amfac))
2701
2702	! * limit magnitude of mp(i) to meet cfl condition *
2703
2704	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
2705	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
2706	ampmax = min(ampmax, amp2)
2707	mp(il, i) = min(mp(il,i), ampmax)
2708
2709	! * force mp to decrease linearly to zero *
2710	! * between cloud base and the surface *
2711
2712
2713	! c if(p(il,i).gt.p(il,icb(il)))then
2714	! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2715	! c endif
2716	IF (ph(il,i)>0.9*plcl(il)) THEN
2717	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
2718	END IF
2719
2720	END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2721	END DO
2722	! ----------------------------------------------------------------
2723
2724	! * find mixing ratio of precipitating downdraft *
2725
2726	DO il = 1, ncum
2727	IF (i<inb(il) .AND. lwork(il)) THEN
2728	mplus(il) = mp(il, i) > mp(il, i+1)
2729	END IF ! (i.lt.inb(il) .and. lwork(il))
2730	END DO
2731
2732	DO il = 1, ncum
2733	IF (i<inb(il) .AND. lwork(il)) THEN
2734
2735	rp(il, i) = rr(il, i)
2736
2737	IF (mplus(il)) THEN
2738
2739	IF (cvflag_grav) THEN
2740	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
2741	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2742	ELSE
2743	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
2744	5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2745	END IF
2746	rp(il, i) = rp(il, i)/mp(il, i)
2747	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
2748	up(il, i) = up(il, i)/mp(il, i)
2749	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
2750	vp(il, i) = vp(il, i)/mp(il, i)
2751
2752	ELSE ! if (mplus(il))
2753
2754	IF (mp(il,i+1)>1.0E-16) THEN
2755	IF (cvflag_grav) THEN
2756	rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
2757	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2758	ELSE
2759	rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
2760	(evap(il,i+1)+evap(il,i))/mp(il, i+1)
2761	END IF
2762	up(il, i) = up(il, i+1)
2763	vp(il, i) = vp(il, i+1)
2764	END IF ! (mp(il,i+1).gt.1.0e-16)
2765	END IF ! (mplus(il)) else if (.not.mplus(il))
2766
2767	rp(il, i) = amin1(rp(il,i), rs(il,i))
2768	rp(il, i) = max(rp(il,i), 0.0)
2769
2770	END IF ! (i.lt.inb(il) .and. lwork(il))
2771	END DO
2772	! ----------------------------------------------------------------
2773
2774	! * find tracer concentrations in precipitating downdraft *
2775
2776	!AC! do j=1,ntra
2777	!AC! do il = 1,ncum
2778	!AC! if (i.lt.inb(il) .and. lwork(il)) then
2779	!AC!c
2780	!AC! if(mplus(il))then
2781	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2782	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2783	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2784	!AC! else ! if (mplus(il))
2785	!AC! if(mp(il,i+1).gt.1.0e-16)then
2786	!AC! trap(il,i,j)=trap(il,i+1,j)
2787	!AC! endif
2788	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
2789	!AC!c
2790	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
2791	!AC! enddo
2792	!AC! end do
2793
2794	400 END DO
2795	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2796
2797	! * end of downdraft loop *
2798
2799	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2800
2801
2802	RETURN
2803	END SUBROUTINE cv3_unsat
2804
2805	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
2806	icb, inb, delt, &
2807	t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
2808	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
2809	ep, clw, m, tp, mp, rp, up, vp, trap, &
2810	wt, water, ice, evap, fondue, faci, b, sigd, &
2811	ment, qent, hent, iflag_mix, uent, vent, &
2812	nent, elij, traent, sig, &
2813	tv, tvp, wghti, &
2814	iflag, precip, Vprecip, Vprecipi, & ! jyg: Vprecipi
2815	ft, fr, fu, fv, ftra, & ! jyg
2816	cbmf, upwd, dnwd, dnwd0, ma, mip, &
2817	!! tls, tps, ! useless . jyg
2818	qcondc, wd, &
2819	ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
2820
2821	IMPLICIT NONE
2822
2823	include "cvthermo.h"
2824	include "cv3param.h"
2825	include "cvflag.h"
2826	include "conema3.h"
2827
2828	!inputs:
2829	INTEGER, INTENT (IN) :: iflag_mix
2830	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2831	LOGICAL, INTENT (IN) :: ok_conserv_q
2832	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
2833	REAL, INTENT (IN) :: delt
2834	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, u, v
2835	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t_wake, rr_wake
2836	REAL, DIMENSION (nloc), INTENT (IN) :: s_wake
2837	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra
2838	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
2839	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2840	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz, h, hp
2841	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tp
2842	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, cpn, ep, clw
2843	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
2844	REAL, DIMENSION (nloc, na), INTENT (IN) :: rp, up
2845	REAL, DIMENSION (nloc, na), INTENT (IN) :: vp
2846	REAL, DIMENSION (nloc, nd), INTENT (IN) :: wt
2847	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: trap
2848	REAL, DIMENSION (nloc, na), INTENT (IN) :: water, evap, b
2849	REAL, DIMENSION (nloc, na), INTENT (IN) :: fondue, faci, ice
2850	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: qent, uent
2851	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: hent
2852	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: vent, elij
2853	INTEGER, DIMENSION (nloc, nd), INTENT (IN) :: nent
2854	REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent
2855	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti
2856	REAL,INTENT(IN) :: tau_cld_cv, coefw_cld_cv
2857	!
2858	!input/output:
2859	REAL, DIMENSION (nloc, na), INTENT (INOUT) :: m, mp
2860	REAL, DIMENSION (nloc, na, na), INTENT (INOUT) :: ment
2861	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
2862	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig
2863	REAL, DIMENSION (nloc), INTENT (INOUT) :: sigd
2864	!
2865	!outputs:
2866	REAL, DIMENSION (nloc), INTENT (OUT) :: precip
2867	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ft, fr, fu, fv
2868	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ftd, fqd
2869	REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT) :: ftra
2870	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: upwd, dnwd, ma
2871	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: dnwd0, mip
2872	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecip
2873	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecipi
2874	!! REAL tls(nloc, nd), tps(nloc, nd) ! useless . jyg
2875	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qcondc ! cld
2876	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qtc, sigt ! cld
2877	REAL, DIMENSION (nloc), INTENT (OUT) :: wd ! gust
2878	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf
2879	!
2880	!local variables:
2881	INTEGER i, k, il, n, j, num1
2882	REAL rat, delti
2883	REAL ax, bx, cx, dx, ex
2884	REAL cpinv, rdcp, dpinv
2885	REAL awat(nloc)
2886	REAL lvcp(nloc, na), lfcp(nloc, na) ! , mke(nloc, na) ! unused . jyg
2887	REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
2888	!! real up1(nloc), dn1(nloc)
2889	REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
2890	REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2891	REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
2892	REAL th_wake(nloc, nd)
2893	REAL alpha_qpos(nloc), alpha_qpos1(nloc)
2894	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
2895	REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
2896	REAL sument(nloc), sigment(nloc,nd), qtment(nloc,nd) ! cld
2897	REAL qnk(nloc)
2898	REAL sumdq !jyg
2899	!
2900	! -------------------------------------------------------------
2901
2902	! initialization:
2903
2904	delti = 1.0/delt
2905	! print*,'cv3_yield initialisation delt', delt
2906	!
2907	DO il = 1, ncum
2908	precip(il) = 0.0
2909	Vprecip(il, nd+1) = 0.0
2910	Vprecipi(il, nd+1) = 0.0 ! jyg: Vprecipi
2911	wd(il) = 0.0 ! gust
2912	END DO
2913
2914	DO i = 1, nd
2915	DO il = 1, ncum
2916	Vprecip(il, i) = 0.0
2917	Vprecipi(il, i) = 0.0 ! jyg
2918	ft(il, i) = 0.0
2919	fr(il, i) = 0.0
2920	fu(il, i) = 0.0
2921	fv(il, i) = 0.0
2922	upwd(il, i) = 0.0
2923	dnwd(il, i) = 0.0
2924	dnwd0(il, i) = 0.0
2925	mip(il, i) = 0.0
2926	ftd(il, i) = 0.0
2927	fqd(il, i) = 0.0
2928	qcondc(il, i) = 0.0 ! cld
2929	qcond(il, i) = 0.0 ! cld
2930	qtc(il, i) = 0.0 ! cld
2931	qtment(il, i) = 0.0 ! cld
2932	sigment(il, i) = 0.0 ! cld
2933	sigt(il, i) = 0.0 ! cld
2934	nqcond(il, i) = 0.0 ! cld
2935	END DO
2936	END DO
2937	! print*,'cv3_yield initialisation 2'
2938	!AC! do j=1,ntra
2939	!AC! do i=1,nd
2940	!AC! do il=1,ncum
2941	!AC! ftra(il,i,j)=0.0
2942	!AC! enddo
2943	!AC! enddo
2944	!AC! enddo
2945	! print*,'cv3_yield initialisation 3'
2946	DO i = 1, nl
2947	DO il = 1, ncum
2948	lvcp(il, i) = lv(il, i)/cpn(il, i)
2949	lfcp(il, i) = lf(il, i)/cpn(il, i)
2950	END DO
2951	END DO
2952
2953
2954
2955	! * calculate surface precipitation in mm/day *
2956
2957	DO il = 1, ncum
2958	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
2959	IF (cvflag_ice) THEN
2960	precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
2961	86400.1000./(rowl*grav)
2962	ELSE
2963	precip(il) = wt(il, 1)sigd(il)water(il, 1) &
2964	86400.1000./(rowl*grav)
2965	END IF
2966	END IF
2967	END DO
2968	! print*,'cv3_yield apres calcul precip'
2969
2970
2971	! === calculate vertical profile of precipitation in kg/m2/s ===
2972
2973	DO i = 1, nl
2974	DO il = 1, ncum
2975	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
2976	IF (cvflag_ice) THEN
2977	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
2978	Vprecipi(il, i) = wt(il, i)sigd(il)ice(il,i)/grav ! jyg
2979	ELSE
2980	Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
2981	Vprecipi(il, i) = 0. ! jyg
2982	END IF
2983	END IF
2984	END DO
2985	END DO
2986
2987
2988	! * Calculate downdraft velocity scale *
2989	! * NE PAS UTILISER POUR L'INSTANT *
2990
2991	!! do il=1,ncum
2992	!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
2993	!! /(sigd(il)*p(il,icb(il)))
2994	!! enddo
2995
2996
2997	! * calculate tendencies of lowest level potential temperature *
2998	! * and mixing ratio *
2999
3000	DO il = 1, ncum
3001	work(il) = 1.0/(ph(il,1)-ph(il,2))
3002	cbmf(il) = 0.0
3003	END DO
3004
3005	DO k = 2, nl
3006	DO il = 1, ncum
3007	IF (k>=icb(il)) THEN
3008	cbmf(il) = cbmf(il) + m(il, k)
3009	END IF
3010	END DO
3011	END DO
3012
3013	! print*,'cv3_yield avant ft'
3014	! am is the part of cbmf taken from the first level
3015	DO il = 1, ncum
3016	am(il) = cbmf(il)*wghti(il, 1)
3017	END DO
3018
3019	DO il = 1, ncum
3020	IF (iflag(il)<=1) THEN
3021	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
3022	!JYG Correction pour conserver l'eau
3023	! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
3024	IF (cvflag_ice) THEN
3025	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
3026	lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
3027	lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
3028	(100.*(ph(il,1)-ph(il,2))) !precip
3029	ELSE
3030	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
3031	END IF
3032
3033	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
3034
3035	IF (cvflag_ice) THEN
3036	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3037	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
3038	0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
3039	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3040	ELSE
3041	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3042	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3043	END IF
3044
3045	ftd(il, 1) = ft(il, 1) ! fin precip
3046
3047	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
3048	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3049	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
3050	END IF ! iflag
3051	END DO
3052
3053
3054	DO j = 2, nl
3055	IF (iflag_mix>0) THEN
3056	DO il = 1, ncum
3057	! FH WARNING a modifier :
3058	cpinv = 0.
3059	! cpinv=1.0/cpn(il,1)
3060	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3061	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
3062	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
3063	END IF ! j
3064	END DO
3065	END IF
3066	END DO
3067	! fin sature
3068
3069
3070	DO il = 1, ncum
3071	IF (iflag(il)<=1) THEN
3072	!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
3073	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
3074	sigd(il)*evap(il, 1)
3075	!!! sigd(il)0.5(evap(il,1)+evap(il,2))
3076
3077	fqd(il, 1) = fr(il, 1) !precip
3078
3079	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
3080
3081	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
3082	am(il)*(u(il,2)-u(il,1)))
3083	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
3084	am(il)*(v(il,2)-v(il,1)))
3085	END IF ! iflag
3086	END DO ! il
3087
3088
3089	!AC! do j=1,ntra
3090	!AC! do il=1,ncum
3091	!AC! if (iflag(il) .le. 1) then
3092	!AC! if (cvflag_grav) then
3093	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
3094	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3095	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3096	!AC! else
3097	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
3098	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3099	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3100	!AC! endif
3101	!AC! endif ! iflag
3102	!AC! enddo
3103	!AC! enddo
3104
3105	DO j = 2, nl
3106	DO il = 1, ncum
3107	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3108	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3109	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
3110	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
3111	END IF ! j
3112	END DO
3113	END DO
3114
3115	!AC! do k=1,ntra
3116	!AC! do j=2,nl
3117	!AC! do il=1,ncum
3118	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
3119	!AC!
3120	!AC! if (cvflag_grav) then
3121	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
3122	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3123	!AC! else
3124	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
3125	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3126	!AC! endif
3127	!AC!
3128	!AC! endif
3129	!AC! enddo
3130	!AC! enddo
3131	!AC! enddo
3132	! print*,'cv3_yield apres ft'
3133
3134	! * calculate tendencies of potential temperature and mixing ratio *
3135	! * at levels above the lowest level *
3136
3137	! * first find the net saturated updraft and downdraft mass fluxes *
3138	! * through each level *
3139
3140
3141	DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
3142
3143	num1 = 0
3144	DO il = 1, ncum
3145	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
3146	END DO
3147	IF (num1<=0) GO TO 500
3148
3149	CALL zilch(amp1, ncum)
3150	CALL zilch(ad, ncum)
3151
3152	DO k = 1, nl + 1
3153	DO il = 1, ncum
3154	IF (i>=icb(il)) THEN
3155	IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
3156	amp1(il) = amp1(il) + m(il, k)
3157	END IF
3158	ELSE
3159	! AMP1 is the part of cbmf taken from layers I and lower
3160	IF (k<=i) THEN
3161	amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
3162	END IF
3163	END IF
3164	END DO
3165	END DO
3166
3167	DO k = 1, i
3168	DO j = i + 1, nl + 1
3169	DO il = 1, ncum
3170	IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
3171	amp1(il) = amp1(il) + ment(il, k, j)
3172	END IF
3173	END DO
3174	END DO
3175	END DO
3176
3177	DO k = 1, i - 1
3178	DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
3179	DO il = 1, ncum
3180	IF (i<=inb(il) .AND. j<=inb(il)) THEN
3181	ad(il) = ad(il) + ment(il, j, k)
3182	END IF
3183	END DO
3184	END DO
3185	END DO
3186
3187	DO il = 1, ncum
3188	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3189	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3190	cpinv = 1.0/cpn(il, i)
3191
3192	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
3193	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
3194
3195	! precip
3196	! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
3197	IF (cvflag_ice) THEN
3198	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
3199	sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
3200	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
3201	ELSE
3202	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
3203	END IF
3204
3205	rat = cpn(il, i-1)*cpinv
3206
3207	ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
3208	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
3209	IF (cvflag_ice) THEN
3210	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
3211	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
3212	0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
3213	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
3214	ELSE
3215	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
3216	(t_wake(il,i+1)-t_wake(il,i))dpinv &
3217	cpinv
3218	END IF
3219
3220	ftd(il, i) = ft(il, i)
3221	! fin precip
3222
3223	! sature
3224	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
3225	(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
3226	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
3227
3228
3229	IF (iflag_mix==0) THEN
3230	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
3231	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
3232	END IF
3233
3234
3235
3236	! sb: on ne fait pas encore la correction permettant de mieux
3237	! conserver l'eau:
3238	!JYG: correction permettant de mieux conserver l'eau:
3239	! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
3240	fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
3241	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
3242	fqd(il, i) = fr(il, i) ! precip
3243
3244	fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
3245	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
3246	fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
3247	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
3248
3249
3250	fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
3251	ad(il)*(rr(il,i)-rr(il,i-1)))
3252	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
3253	ad(il)*(u(il,i)-u(il,i-1)))
3254	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
3255	ad(il)*(v(il,i)-v(il,i-1)))
3256
3257	END IF ! i
3258	END DO
3259
3260	!AC! do k=1,ntra
3261	!AC! do il=1,ncum
3262	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3263	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3264	!AC! cpinv=1.0/cpn(il,i)
3265	!AC! if (cvflag_grav) then
3266	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
3267	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
3268	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
3269	!AC! else
3270	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
3271	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
3272	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
3273	!AC! endif
3274	!AC! endif
3275	!AC! enddo
3276	!AC! enddo
3277
3278	DO k = 1, i - 1
3279
3280	DO il = 1, ncum
3281	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
3282	awat(il) = max(awat(il), 0.0)
3283	END DO
3284
3285	IF (iflag_mix/=0) THEN
3286	DO il = 1, ncum
3287	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3288	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3289	cpinv = 1.0/cpn(il, i)
3290	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
3291	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
3292	!
3293	!
3294	END IF ! i
3295	END DO
3296	END IF
3297
3298	DO il = 1, ncum
3299	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3300	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3301	cpinv = 1.0/cpn(il, i)
3302	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
3303	(qent(il,k,i)-awat(il)-rr(il,i))
3304	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
3305	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
3306
3307	! (saturated updrafts resulting from mixing) ! cld
3308	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
3309	qtment(il, i) = qtment(il, i) + qent(il,k,i) ! cld
3310	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3311	END IF ! i
3312	END DO
3313	END DO
3314
3315	!AC! do j=1,ntra
3316	!AC! do k=1,i-1
3317	!AC! do il=1,ncum
3318	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3319	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3320	!AC! cpinv=1.0/cpn(il,i)
3321	!AC! if (cvflag_grav) then
3322	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3323	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3324	!AC! else
3325	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3326	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3327	!AC! endif
3328	!AC! endif
3329	!AC! enddo
3330	!AC! enddo
3331	!AC! enddo
3332
3333	DO k = i, nl + 1
3334
3335	IF (iflag_mix/=0) THEN
3336	DO il = 1, ncum
3337	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
3338	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3339	cpinv = 1.0/cpn(il, i)
3340	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
3341	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
3342
3343
3344	END IF ! i
3345	END DO
3346	END IF
3347
3348	DO il = 1, ncum
3349	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
3350	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3351	cpinv = 1.0/cpn(il, i)
3352
3353	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
3354	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
3355	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
3356	END IF ! i and k
3357	END DO
3358	END DO
3359
3360	!AC! do j=1,ntra
3361	!AC! do k=i,nl+1
3362	!AC! do il=1,ncum
3363	!AC! if (i.le.inb(il) .and. k.le.inb(il)
3364	!AC! $ .and. iflag(il) .le. 1) then
3365	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3366	!AC! cpinv=1.0/cpn(il,i)
3367	!AC! if (cvflag_grav) then
3368	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3369	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3370	!AC! else
3371	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3372	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3373	!AC! endif
3374	!AC! endif ! i and k
3375	!AC! enddo
3376	!AC! enddo
3377	!AC! enddo
3378
3379	! sb: interface with the cloud parameterization: ! cld
3380
3381	DO k = i + 1, nl
3382	DO il = 1, ncum
3383	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
3384	! (saturated downdrafts resulting from mixing) ! cld
3385	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
3386	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
3387	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3388	END IF ! cld
3389	END DO ! cld
3390	END DO ! cld
3391
3392	! (particular case: no detraining level is found) ! cld
3393	DO il = 1, ncum ! cld
3394	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
3395	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
3396	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
3397	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3398	END IF ! cld
3399	END DO ! cld
3400
3401	DO il = 1, ncum ! cld
3402	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
3403	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
3404	qtment(il, i) = qtment(il,i)/nqcond(il, i) ! cld
3405	END IF ! cld
3406	END DO
3407
3408	!AC! do j=1,ntra
3409	!AC! do il=1,ncum
3410	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3411	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3412	!AC! cpinv=1.0/cpn(il,i)
3413	!AC!
3414	!AC! if (cvflag_grav) then
3415	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
3416	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3417	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3418	!AC! else
3419	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
3420	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3421	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3422	!AC! endif
3423	!AC! endif ! i
3424	!AC! enddo
3425	!AC! enddo
3426
3427
3428	500 END DO
3429
3430	!JYG<
3431	!Conservation de l'eau
3432	! sumdq = 0.
3433	! DO k = 1, nl
3434	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3435	! END DO
3436	! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3437	!JYG>
3438	! * move the detrainment at level inb down to level inb-1 *
3439	! * in such a way as to preserve the vertically *
3440	! * integrated enthalpy and water tendencies *
3441
3442	! Correction bug le 18-03-09
3443	DO il = 1, ncum
3444	IF (iflag(il)<=1) THEN
3445	ax = 0.01gravment(il, inb(il), inb(il))* &
3446	(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
3447	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
3448	ft(il, inb(il)) = ft(il, inb(il)) - ax
3449	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3450	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
3451
3452	bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
3453	(ph(il,inb(il))-ph(il,inb(il)+1))
3454	fr(il, inb(il)) = fr(il, inb(il)) - bx
3455	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3456	(ph(il,inb(il)-1)-ph(il,inb(il)))
3457
3458	cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
3459	(ph(il,inb(il))-ph(il,inb(il)+1))
3460	fu(il, inb(il)) = fu(il, inb(il)) - cx
3461	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3462	(ph(il,inb(il)-1)-ph(il,inb(il)))
3463
3464	dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
3465	(ph(il,inb(il))-ph(il,inb(il)+1))
3466	fv(il, inb(il)) = fv(il, inb(il)) - dx
3467	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3468	(ph(il,inb(il)-1)-ph(il,inb(il)))
3469	END IF !iflag
3470	END DO
3471
3472	!JYG<
3473	!Conservation de l'eau
3474	! sumdq = 0.
3475	! DO k = 1, nl
3476	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3477	! END DO
3478	! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3479	!JYG>
3480
3481	!AC! do j=1,ntra
3482	!AC! do il=1,ncum
3483	!AC! IF (iflag(il) .le. 1) THEN
3484	!AC! IF (cvflag_grav) then
3485	!AC! ex=0.01gravment(il,inb(il),inb(il))
3486	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3487	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3488	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3489	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3490	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3491	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3492	!AC! else
3493	!AC! ex=0.1*ment(il,inb(il),inb(il))
3494	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3495	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3496	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3497	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3498	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3499	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3500	!AC! ENDIF !cvflag grav
3501	!AC! ENDIF !iflag
3502	!AC! enddo
3503	!AC! enddo
3504
3505
3506	! * homogenize tendencies below cloud base *
3507
3508
3509	DO il = 1, ncum
3510	asum(il) = 0.0
3511	bsum(il) = 0.0
3512	csum(il) = 0.0
3513	dsum(il) = 0.0
3514	esum(il) = 0.0
3515	fsum(il) = 0.0
3516	gsum(il) = 0.0
3517	hsum(il) = 0.0
3518	END DO
3519
3520	!do i=1,nl
3521	!do il=1,ncum
3522	!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
3523	!enddo
3524	!enddo
3525
3526	DO i = 1, nl
3527	DO il = 1, ncum
3528	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
3529	!jyg Saturated part : use T profile
3530	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
3531	!jyg<20140311
3532	!Correction pour conserver l eau
3533	IF (ok_conserv_q) THEN
3534	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
3535	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
3536
3537	ELSE
3538	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
3539	(ph(il,i)-ph(il,i+1))
3540	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
3541	(ph(il,i)-ph(il,i+1))
3542	ENDIF ! (ok_conserv_q)
3543	!jyg>
3544	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
3545	!jyg Unsaturated part : use T_wake profile
3546	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
3547	!jyg<20140311
3548	!Correction pour conserver l eau
3549	IF (ok_conserv_q) THEN
3550	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
3551	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
3552	ELSE
3553	fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
3554	(ph(il,i)-ph(il,i+1))
3555	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
3556	(ph(il,i)-ph(il,i+1))
3557	ENDIF ! (ok_conserv_q)
3558	!jyg>
3559	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
3560	END IF
3561	END DO
3562	END DO
3563
3564	!!!! do 700 i=1,icb(il)-1
3565	DO i = 1, nl
3566	DO il = 1, ncum
3567	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
3568	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
3569	fqd(il, i) = fsum(il)/gsum(il)
3570	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
3571	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
3572	END IF
3573	END DO
3574	END DO
3575
3576	!jyg<
3577	!Conservation de l'eau
3578	!! sumdq = 0.
3579	!! DO k = 1, nl
3580	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3581	!! END DO
3582	!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3583	!jyg>
3584
3585
3586	! *** Check that moisture stays positive. If not, scale tendencies
3587	! in order to ensure moisture positivity
3588	DO il = 1, ncum
3589	alpha_qpos(il) = 1.
3590	IF (iflag(il)<=1) THEN
3591	IF (fr(il,1)<=0.) THEN
3592	alpha_qpos(il) = max(alpha_qpos(il), (-deltfr(il,1))/(s_wake(il)rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
3593	END IF
3594	END IF
3595	END DO
3596	DO i = 2, nl
3597	DO il = 1, ncum
3598	IF (iflag(il)<=1) THEN
3599	IF (fr(il,i)<=0.) THEN
3600	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
3601	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
3602	END IF
3603	END IF
3604	END DO
3605	END DO
3606	DO il = 1, ncum
3607	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
3608	alpha_qpos(il) = alpha_qpos(il)*1.1
3609	END IF
3610	END DO
3611	!
3612	! print *,' YIELD : alpha_qpos ',alpha_qpos(1)
3613	!
3614	DO il = 1, ncum
3615	IF (iflag(il)<=1) THEN
3616	sigd(il) = sigd(il)/alpha_qpos(il)
3617	precip(il) = precip(il)/alpha_qpos(il)
3618	END IF
3619	END DO
3620	DO i = 1, nl
3621	DO il = 1, ncum
3622	IF (iflag(il)<=1) THEN
3623	fr(il, i) = fr(il, i)/alpha_qpos(il)
3624	ft(il, i) = ft(il, i)/alpha_qpos(il)
3625	fqd(il, i) = fqd(il, i)/alpha_qpos(il)
3626	ftd(il, i) = ftd(il, i)/alpha_qpos(il)
3627	fu(il, i) = fu(il, i)/alpha_qpos(il)
3628	fv(il, i) = fv(il, i)/alpha_qpos(il)
3629	m(il, i) = m(il, i)/alpha_qpos(il)
3630	mp(il, i) = mp(il, i)/alpha_qpos(il)
3631	Vprecip(il, i) = Vprecip(il, i)/alpha_qpos(il)
3632	Vprecipi(il, i) = Vprecipi(il, i)/alpha_qpos(il) ! jyg
3633	END IF
3634	END DO
3635	END DO
3636	DO i = 1, nl
3637	DO j = 1, nl
3638	DO il = 1, ncum
3639	IF (iflag(il)<=1) THEN
3640	ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
3641	END IF
3642	END DO
3643	END DO
3644	END DO
3645
3646	!AC! DO j = 1,ntra
3647	!AC! DO i = 1,nl
3648	!AC! DO il = 1,ncum
3649	!AC! IF (iflag(il) .le. 1) THEN
3650	!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
3651	!AC! ENDIF
3652	!AC! ENDDO
3653	!AC! ENDDO
3654	!AC! ENDDO
3655
3656
3657	! * reset counter and return *
3658
3659	! Reset counter only for points actually convective (jyg)
3660	! In order take into account the possibility of changing the compression,
3661	! reset m, sig and w0 to zero for non-convecting points.
3662	DO il = 1, ncum
3663	IF (iflag(il) < 3) THEN
3664	sig(il, nd) = 2.0
3665	ENDIF
3666	END DO
3667
3668
3669	DO i = 1, nd
3670	DO il = 1, ncum
3671	upwd(il, i) = 0.0
3672	dnwd(il, i) = 0.0
3673	END DO
3674	END DO
3675
3676	DO i = 1, nl
3677	DO il = 1, ncum
3678	dnwd0(il, i) = -mp(il, i)
3679	END DO
3680	END DO
3681	DO i = nl + 1, nd
3682	DO il = 1, ncum
3683	dnwd0(il, i) = 0.
3684	END DO
3685	END DO
3686
3687
3688	DO i = 1, nl
3689	DO il = 1, ncum
3690	IF (i>=icb(il) .AND. i<=inb(il)) THEN
3691	upwd(il, i) = 0.0
3692	dnwd(il, i) = 0.0
3693	END IF
3694	END DO
3695	END DO
3696
3697	DO i = 1, nl
3698	DO k = 1, nl
3699	DO il = 1, ncum
3700	up1(il, k, i) = 0.0
3701	dn1(il, k, i) = 0.0
3702	END DO
3703	END DO
3704	END DO
3705
3706	DO i = 1, nl
3707	DO k = i, nl
3708	DO n = 1, i - 1
3709	DO il = 1, ncum
3710	IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
3711	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
3712	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
3713	END IF
3714	END DO
3715	END DO
3716	END DO
3717	END DO
3718
3719	DO i = 1, nl
3720	DO k = 1, nl
3721	DO il = 1, ncum
3722	IF (i>=icb(il)) THEN
3723	IF (k>=i .AND. k<=(inb(il))) THEN
3724	upwd(il, i) = upwd(il, i) + m(il, k)
3725	END IF
3726	ELSE
3727	IF (k<i) THEN
3728	upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
3729	END IF
3730	END IF
3731	! c print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
3732	END DO
3733	END DO
3734	END DO
3735
3736	DO i = 2, nl
3737	DO k = i, nl
3738	DO il = 1, ncum
3739	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
3740	IF (i<=inb(il) .AND. k<=inb(il)) THEN
3741	upwd(il, i) = upwd(il, i) + up1(il, k, i)
3742	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
3743	END IF
3744	! c print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
3745	END DO
3746	END DO
3747	END DO
3748
3749
3750	!!!! DO il=1,ncum
3751	!!!! do i=icb(il),inb(il)
3752	!!!!
3753	!!!! upwd(il,i)=0.0
3754	!!!! dnwd(il,i)=0.0
3755	!!!! do k=i,inb(il)
3756	!!!! up1=0.0
3757	!!!! dn1=0.0
3758	!!!! do n=1,i-1
3759	!!!! up1=up1+ment(il,n,k)
3760	!!!! dn1=dn1-ment(il,k,n)
3761	!!!! enddo
3762	!!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
3763	!!!! dnwd(il,i)=dnwd(il,i)+dn1
3764	!!!! enddo
3765	!!!! enddo
3766	!!!!
3767	!!!! ENDDO
3768
3769	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3770	! determination de la variation de flux ascendant entre
3771	! deux niveau non dilue mip
3772	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3773
3774	DO i = 1, nl
3775	DO il = 1, ncum
3776	mip(il, i) = m(il, i)
3777	END DO
3778	END DO
3779
3780	DO i = nl + 1, nd
3781	DO il = 1, ncum
3782	mip(il, i) = 0.
3783	END DO
3784	END DO
3785
3786	DO i = 1, nd
3787	DO il = 1, ncum
3788	ma(il, i) = 0
3789	END DO
3790	END DO
3791
3792	DO i = 1, nl
3793	DO j = i, nl
3794	DO il = 1, ncum
3795	ma(il, i) = ma(il, i) + m(il, j)
3796	END DO
3797	END DO
3798	END DO
3799
3800	DO i = nl + 1, nd
3801	DO il = 1, ncum
3802	ma(il, i) = 0.
3803	END DO
3804	END DO
3805
3806	DO i = 1, nl
3807	DO il = 1, ncum
3808	IF (i<=(icb(il)-1)) THEN
3809	ma(il, i) = 0
3810	END IF
3811	END DO
3812	END DO
3813
3814	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3815	! icb represente de niveau ou se trouve la
3816	! base du nuage , et inb le top du nuage
3817	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3818
3819	!! DO i = 1, nd ! unused . jyg
3820	!! DO il = 1, ncum ! unused . jyg
3821	!! mke(il, i) = upwd(il, i) + dnwd(il, i) ! unused . jyg
3822	!! END DO ! unused . jyg
3823	!! END DO ! unused . jyg
3824
3825	!! DO i = 1, nd ! unused . jyg
3826	!! DO il = 1, ncum ! unused . jyg
3827	!! rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd(1.-rr(il,i))+rr(il,i)cpv) ! unused . jyg
3828	!! tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp ! unused . jyg
3829	!! tps(il, i) = tp(il, i) ! unused . jyg
3830	!! END DO ! unused . jyg
3831	!! END DO ! unused . jyg
3832
3833
3834	! * diagnose the in-cloud mixing ratio * ! cld
3835	! * of condensed water * ! cld
3836	!! cld
3837
3838	DO i = 1, nd ! cld
3839	DO il = 1, ncum ! cld
3840	mac(il, i) = 0.0 ! cld
3841	wa(il, i) = 0.0 ! cld
3842	siga(il, i) = 0.0 ! cld
3843	sax(il, i) = 0.0 ! cld
3844	END DO ! cld
3845	END DO ! cld
3846
3847	DO i = minorig, nl ! cld
3848	DO k = i + 1, nl + 1 ! cld
3849	DO il = 1, ncum ! cld
3850	IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
3851	mac(il, i) = mac(il, i) + m(il, k) ! cld
3852	END IF ! cld
3853	END DO ! cld
3854	END DO ! cld
3855	END DO ! cld
3856
3857	DO i = 1, nl ! cld
3858	DO j = 1, i ! cld
3859	DO il = 1, ncum ! cld
3860	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
3861	.AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
3862	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
3863	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
3864	END IF ! cld
3865	END DO ! cld
3866	END DO ! cld
3867	END DO ! cld
3868
3869	DO i = 1, nl ! cld
3870	DO il = 1, ncum ! cld
3871	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
3872	.AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
3873	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
3874	END IF ! cld
3875	END DO ! cld
3876	END DO
3877	! cld
3878	DO i = 1, nl
3879
3880	! 14/01/15 AJ je remets les parties manquantes cf JYG
3881	! Initialize sument to 0
3882
3883	DO il = 1,ncum
3884	sument(il) = 0.
3885	ENDDO
3886
3887	! Sum mixed mass fluxes in sument
3888
3889	DO k = 1,nl
3890	DO il = 1,ncum
3891	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
3892	sument(il) =sument(il) + abs(ment(il,k,i))
3893	ENDIF
3894	ENDDO ! il
3895	ENDDO ! k
3896
3897	! 14/01/15 AJ delta n'a rien à faire là...
3898	DO il = 1, ncum ! cld
3899	IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
3900	siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
3901	rrdtvp(il, i)/p(il, i)/100. ! cld
3902
3903	siga(il, i) = min(siga(il,i), 1.0) ! cld
3904
3905	! IM cf. FH
3906	! 14/01/15 AJ ne correspond pas à ce qui a été codé par JYG et SB
3907
3908	IF (iflag_clw==0) THEN ! cld
3909	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
3910	+(1.-siga(il,i))*qcond(il, i) ! cld
3911
3912
3913	sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1)) ! cld
3914	sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i)) ! cld
3915	qtc(il, i) = (siga(il,i)qnk(il)+sigment(il,i)qtment(il,i)) & ! cld
3916	/(siga(il,i)+sigment(il,i)) ! cld
3917	sigt(il,i) = sigment(il, i) + siga(il, i)
3918
3919	! qtc(il, i) = siga(il,i)qnk(il)+(1.-siga(il,i))qtment(il,i) ! cld
3920	! print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)
3921
3922	ELSE IF (iflag_clw==1) THEN ! cld
3923	qcondc(il, i) = qcond(il, i) ! cld
3924	qtc(il,i) = qtment(il,i) ! cld
3925	END IF ! cld
3926
3927	END DO ! cld
3928	END DO
3929	! print*,'cv3_yield fin'
3930
3931	RETURN
3932	END SUBROUTINE cv3_yield
3933
3934	!AC! et !RomP >>>
3935	SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
3936	ment, sigij, da, phi, phi2, d1a, dam, &
3937	ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
3938	icb, inb)
3939	IMPLICIT NONE
3940
3941	include "cv3param.h"
3942
3943	!inputs:
3944	INTEGER ncum, nd, na, nloc, len
3945	REAL ment(nloc, na, na), sigij(nloc, na, na)
3946	REAL clw(nloc, nd), elij(nloc, na, na)
3947	REAL ep(nloc, na)
3948	INTEGER icb(nloc), inb(nloc)
3949	REAL Vprecip(nloc, nd+1)
3950	!ouputs:
3951	REAL da(nloc, na), phi(nloc, na, na)
3952	REAL phi2(nloc, na, na)
3953	REAL d1a(nloc, na), dam(nloc, na)
3954	REAL epmlmMm(nloc, na, na), eplaMm(nloc, na)
3955	! variables pour tracer dans precip de l'AA et des mel
3956	!local variables:
3957	INTEGER i, j, k
3958	REAL epm(nloc, na, na)
3959
3960	! variables d'Emanuel : du second indice au troisieme
3961	! ---> tab(i,k,j) -> de l origine k a l arrivee j
3962	! ment, sigij, elij
3963	! variables personnelles : du troisieme au second indice
3964	! ---> tab(i,j,k) -> de k a j
3965	! phi, phi2
3966
3967	! initialisations
3968
3969	da(:, :) = 0.
3970	d1a(:, :) = 0.
3971	dam(:, :) = 0.
3972	epm(:, :, :) = 0.
3973	eplaMm(:, :) = 0.
3974	epmlmMm(:, :, :) = 0.
3975	phi(:, :, :) = 0.
3976	phi2(:, :, :) = 0.
3977
3978	! fraction deau condensee dans les melanges convertie en precip : epm
3979	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3980	DO j = 1, na
3981	DO k = 1, na
3982	DO i = 1, ncum
3983	IF (k>=icb(i) .AND. k<=inb(i) .AND. &
3984	!!jyg j.ge.k.and.j.le.inb(i)) then
3985	!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3986	j>k .AND. j<=inb(i)) THEN
3987	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
3988	!!
3989	epm(i, j, k) = max(epm(i,j,k), 0.0)
3990	END IF
3991	END DO
3992	END DO
3993	END DO
3994
3995
3996	DO j = 1, na
3997	DO k = 1, na
3998	DO i = 1, ncum
3999	IF (k>=icb(i) .AND. k<=inb(i)) THEN
4000	eplaMm(i, j) = eplamm(i, j) + &
4001	ep(i, j)clw(i, j)ment(i, j, k)*(1.-sigij(i,j,k))
4002	END IF
4003	END DO
4004	END DO
4005	END DO
4006
4007	DO j = 1, na
4008	DO k = 1, j - 1
4009	DO i = 1, ncum
4010	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
4011	epmlmMm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
4012	END IF
4013	END DO
4014	END DO
4015	END DO
4016
4017	! matrices pour calculer la tendance des concentrations dans cvltr.F90
4018	DO j = 1, na
4019	DO k = 1, na
4020	DO i = 1, ncum
4021	da(i, j) = da(i, j) + (1.-sigij(i,k,j))*ment(i, k, j)
4022	phi(i, j, k) = sigij(i, k, j)*ment(i, k, j)
4023	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sigij(i,k,j))
4024	IF (k<=j) THEN
4025	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, k, j)(1.-ep(i,k))*(1.-sigij(i,k,j))
4026	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
4027	END IF
4028	END DO
4029	END DO
4030	END DO
4031
4032	RETURN
4033	END SUBROUTINE cv3_tracer
4034	!AC! et !RomP <<<
4035
4036	SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
4037	iflag, &
4038	precip, sig, w0, &
4039	ft, fq, fu, fv, ftra, &
4040	Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
4041	iflag1, &
4042	precip1, sig1, w01, &
4043	ft1, fq1, fu1, fv1, ftra1, &
4044	Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1)
4045	IMPLICIT NONE
4046
4047	include "cv3param.h"
4048
4049	!inputs:
4050	INTEGER len, ncum, nd, ntra, nloc
4051	INTEGER idcum(nloc)
4052	INTEGER iflag(nloc)
4053	REAL precip(nloc)
4054	REAL sig(nloc, nd), w0(nloc, nd)
4055	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
4056	REAL ftra(nloc, nd, ntra)
4057	REAL ma(nloc, nd)
4058	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
4059	REAL qcondc(nloc, nd)
4060	REAL wd(nloc), cape(nloc)
4061
4062	!outputs:
4063	INTEGER iflag1(len)
4064	REAL precip1(len)
4065	REAL sig1(len, nd), w01(len, nd)
4066	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
4067	REAL ftra1(len, nd, ntra)
4068	REAL ma1(len, nd)
4069	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
4070	REAL qcondc1(nloc, nd)
4071	REAL wd1(nloc), cape1(nloc)
4072
4073	!local variables:
4074	INTEGER i, k, j
4075
4076	DO i = 1, ncum
4077	precip1(idcum(i)) = precip(i)
4078	iflag1(idcum(i)) = iflag(i)
4079	wd1(idcum(i)) = wd(i)
4080	cape1(idcum(i)) = cape(i)
4081	END DO
4082
4083	DO k = 1, nl
4084	DO i = 1, ncum
4085	sig1(idcum(i), k) = sig(i, k)
4086	w01(idcum(i), k) = w0(i, k)
4087	ft1(idcum(i), k) = ft(i, k)
4088	fq1(idcum(i), k) = fq(i, k)
4089	fu1(idcum(i), k) = fu(i, k)
4090	fv1(idcum(i), k) = fv(i, k)
4091	ma1(idcum(i), k) = ma(i, k)
4092	upwd1(idcum(i), k) = upwd(i, k)
4093	dnwd1(idcum(i), k) = dnwd(i, k)
4094	dnwd01(idcum(i), k) = dnwd0(i, k)
4095	qcondc1(idcum(i), k) = qcondc(i, k)
4096	END DO
4097	END DO
4098
4099	DO i = 1, ncum
4100	sig1(idcum(i), nd) = sig(i, nd)
4101	END DO
4102
4103
4104	!AC! do 2100 j=1,ntra
4105	!AC!c oct3 do 2110 k=1,nl
4106	!AC! do 2110 k=1,nd ! oct3
4107	!AC! do 2120 i=1,ncum
4108	!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
4109	!AC! 2120 continue
4110	!AC! 2110 continue
4111	!AC! 2100 continue
4112	!
4113	RETURN
4114	END SUBROUTINE cv3_uncompress

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

Download in other formats:

Original Format