Context Navigation

cv3_routines.F90 @ 2393

Last change on this file since 2393 was 2393, checked in by jyg, 9 years ago
Add various intializations of arrays in lmdz1d.F90 and in the convection scheme. Add output variables for boundary layer splitting.
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: 134.2 KB

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

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

Download in other formats:

Original Format