Context Navigation

cv3_routines.F90 @ 2021

Last change on this file since 2021 was 2007, checked in by fhourdin, 10 years ago

Modification pour garantir la conservation des espèces traces.

Modification for tracer conservation

Jean-Yves Grandpeix

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

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

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

Download in other formats:

Original Format