Context Navigation

cv3_routines.F90 @ 2198

Last change on this file since 2198 was 2110, checked in by lguez, 10 years ago

Abort if surface pressure becomes negative. The call to abort_gcm
already was in the sequential version but not in dyn3dpar nor in
dyn3dmem. In dyn3dmem, there is a variable checksum_all, which is
always true. The call to MPI_ALLREDUCE which should update
checksum_all is commented out. I do not know why (performance?).

Non-ASCII characters in comments are not always rendered properly and
they risk being lost. See revision [1740].

Bug fix in procedure convect2: the dimension len must be declared
before the array idcum which has this dimension. Bug fix in procedure
icefrac: the dimensions nl and len must be declared before the arrays
which have these dimensions.

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

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

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

Download in other formats:

Original Format