Context Navigation

cv3_routines.F90 @ 2094

Last change on this file since 2094 was 2078, checked in by idelkadi, 10 years ago
Rajout des directives OPENMP pour la lecture des cles et parametres de convection et wakes dans les fichier conv_param.data, ep_param.data et wake_param.data
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 2078 2014-07-04 10:44:45Z lguez $
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	REAL qi(len, nl)
997	REAL t(len, nl), clw(len, nl)
998	REAL fracg
999	INTEGER nl, len, k, i
1000
1001	DO k = 3, nl
1002	DO i = 1, len
1003	IF (t(i,k)>263.15) THEN
1004	qi(i, k) = 0.
1005	ELSE
1006	IF (t(i,k)<243.15) THEN
1007	qi(i, k) = clw(i, k)
1008	ELSE
1009	fracg = (263.15-t(i,k))/20
1010	qi(i, k) = clw(i, k)*fracg
1011	END IF
1012	END IF
1013	! print*,t(i,k),qi(i,k),'temp,testglace'
1014	END DO
1015	END DO
1016
1017	RETURN
1018
1019	END SUBROUTINE icefrac
1020
1021	SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, &
1022	tnk, qnk, gznk, hnk, t, q, qs, gz, &
1023	p, h, tv, lv, lf, pbase, buoybase, plcl, &
1024	inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
1025	IMPLICIT NONE
1026
1027	! ---------------------------------------------------------------------
1028	! Purpose:
1029	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1030	! &
1031	! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
1032	! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
1033	! &
1034	! FIND THE LEVEL OF NEUTRAL BUOYANCY
1035
1036	! Main differences convect3/convect4:
1037	! - icbs (input) is the first level above LCL (may differ from icb)
1038	! - many minor differences in the iterations
1039	! - condensed water not removed from tvp in convect3
1040	! - vertical profile of buoyancy computed here (use of buoybase)
1041	! - the determination of inb is different
1042	! - no inb1, only inb in output
1043	! ---------------------------------------------------------------------
1044
1045	include "cvthermo.h"
1046	include "cv3param.h"
1047	include "conema3.h"
1048	include "cvflag.h"
1049
1050	!inputs:
1051	INTEGER ncum, nd, nloc, j
1052	INTEGER icb(nloc), icbs(nloc), nk(nloc)
1053	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
1054	REAL p(nloc, nd)
1055	REAL tnk(nloc), qnk(nloc), gznk(nloc)
1056	REAL hnk(nloc)
1057	REAL lv(nloc, nd), lf(nloc, nd), tv(nloc, nd), h(nloc, nd)
1058	REAL pbase(nloc), buoybase(nloc), plcl(nloc)
1059
1060	!outputs:
1061	INTEGER inb(nloc)
1062	REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
1063	REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
1064	REAL buoy(nloc, nd)
1065
1066	!local variables:
1067	INTEGER i, k
1068	REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
1069	REAL als
1070	REAL qsat_new, snew, qi(nloc, nd)
1071	REAL by, defrac, pden, tbis
1072	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
1073	REAL frac(nloc, nd)
1074	LOGICAL lcape(nloc)
1075	INTEGER iposit(nloc)
1076	REAL fracg
1077
1078	! =====================================================================
1079	! --- SOME INITIALIZATIONS
1080	! =====================================================================
1081
1082	DO k = 1, nl
1083	DO i = 1, ncum
1084	ep(i, k) = 0.0
1085	sigp(i, k) = spfac
1086	qi(i, k) = 0.
1087	END DO
1088	END DO
1089
1090	! =====================================================================
1091	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1092	! =====================================================================
1093
1094	! --- The procedure is to solve the equation.
1095	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
1096
1097	! * Calculate certain parcel quantities, including static energy *
1098
1099
1100	DO i = 1, ncum
1101	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))*tnk(i)+ &
1102	! debug qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
1103	qnk(i)(lv0-clmcpv(tnk(i)-273.15)) + gznk(i)
1104	END DO
1105
1106
1107	! * Find lifted parcel quantities above cloud base *
1108
1109
1110	DO k = minorig + 1, nl
1111	DO i = 1, ncum
1112	! ori if(k.ge.(icb(i)+1))then
1113	IF (k>=(icbs(i)+1)) THEN ! convect3
1114	tg = t(i, k)
1115	qg = qs(i, k)
1116	! debug alv=lv0-clmcpv*(t(i,k)-t0)
1117	alv = lv0 - clmcpv*(t(i,k)-273.15)
1118
1119	! First iteration.
1120
1121	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1122	s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
1123	alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
1124	s = 1./s
1125	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1126	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1127	tg = tg + s*(ah0(i)-ahg)
1128	! ori tg=max(tg,35.0)
1129	! debug tc=tg-t0
1130	tc = tg - 273.15
1131	denom = 243.5 + tc
1132	denom = max(denom, 1.0) ! convect3
1133	! ori if(tc.ge.0.0)then
1134	es = 6.112exp(17.67tc/denom)
1135	! ori else
1136	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1137	! ori endif
1138	qg = epses/(p(i,k)-es(1.-eps))
1139
1140	! Second iteration.
1141
1142	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1143	! ori s=1./s
1144	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1145	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1146	tg = tg + s*(ah0(i)-ahg)
1147	! ori tg=max(tg,35.0)
1148	! debug tc=tg-t0
1149	tc = tg - 273.15
1150	denom = 243.5 + tc
1151	denom = max(denom, 1.0) ! convect3
1152	! ori if(tc.ge.0.0)then
1153	es = 6.112exp(17.67tc/denom)
1154	! ori else
1155	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1156	! ori endif
1157	qg = epses/(p(i,k)-es(1.-eps))
1158
1159	! debug alv=lv0-clmcpv*(t(i,k)-t0)
1160	alv = lv0 - clmcpv*(t(i,k)-273.15)
1161	! print*,'cpd dans convect2 ',cpd
1162	! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
1163	! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
1164
1165	! ori c approximation here:
1166	! ori tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
1167
1168	! convect3: no approximation:
1169	IF (cvflag_ice) THEN
1170	tp(i, k) = max(0., (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i)))
1171	ELSE
1172	tp(i, k) = (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
1173	END IF
1174
1175	clw(i, k) = qnk(i) - qg
1176	clw(i, k) = max(0.0, clw(i,k))
1177	rg = qg/(1.-qnk(i))
1178	! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
1179	! convect3: (qg utilise au lieu du vrai mixing ratio rg):
1180	tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
1181	IF (cvflag_ice) THEN
1182	IF (clw(i,k)<1.E-11) THEN
1183	tp(i, k) = tv(i, k)
1184	tvp(i, k) = tv(i, k)
1185	END IF
1186	END IF
1187	END IF
1188
1189	IF (cvflag_ice) THEN
1190	!CR:attention boucle en klon dans Icefrac
1191	! Call Icefrac(t,clw,qi,nl,nloc)
1192	IF (t(i,k)>263.15) THEN
1193	qi(i, k) = 0.
1194	ELSE
1195	IF (t(i,k)<243.15) THEN
1196	qi(i, k) = clw(i, k)
1197	ELSE
1198	fracg = (263.15-t(i,k))/20
1199	qi(i, k) = clw(i, k)*fracg
1200	END IF
1201	END IF
1202	!CR: fin test
1203	IF (t(i,k)<263.15) THEN
1204	!CR: on commente les calculs d'Arnaud car division par zero
1205	! nouveau calcul propose par JYG
1206	! alv=lv0-clmcpv*(t(i,k)-273.15)
1207	! alf=lf0-clmci*(t(i,k)-273.15)
1208	! tg=tp(i,k)
1209	! tc=tp(i,k)-273.15
1210	! denom=243.5+tc
1211	! do j=1,3
1212	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1213	! il faudra que esi vienne en argument de la convection
1214	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1215	! tbis=t(i,k)+(tp(i,k)-tg)
1216	! esi=exp(23.33086-(6111.72784/tbis) + &
1217	! 0.15215*log(tbis))
1218	! qsat_new=epsesi/(p(i,k)-esi(1.-eps))
1219	! snew=cpd(1.-qnk(i))+clqnk(i)+alvalvqsat_new/ &
1220	! (rrvtbistbis)
1221	! snew=1./snew
1222	! print*,esi,qsat_new,snew,'esi,qsat,snew'
1223	! tp(i,k)=tg+(alfqi(i,k)+alvqg(1.-(esi/es)))snew
1224	! print*,k,tp(i,k),qnk(i),'avec glace'
1225	! print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
1226	! enddo
1227
1228	alv = lv0 - clmcpv*(t(i,k)-273.15)
1229	alf = lf0 + clmci*(t(i,k)-273.15)
1230	als = alf + alv
1231	tg = tp(i, k)
1232	tp(i, k) = t(i, k)
1233	DO j = 1, 3
1234	esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
1235	qsat_new = epsesi/(p(i,k)-esi(1.-eps))
1236	snew = cpd(1.-qnk(i)) + clqnk(i) + alvalsqsat_new/ &
1237	(rrvtp(i,k)tp(i,k))
1238	snew = 1./snew
1239	! c print*,esi,qsat_new,snew,'esi,qsat,snew'
1240	tp(i, k) = tp(i, k) + &
1241	((cpd(1.-qnk(i))+clqnk(i))*(tg-tp(i,k)) + &
1242	alv(qg-qsat_new)+alfqi(i,k))*snew
1243	! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
1244	! 'k,tp,q,qt,qi avec glace'
1245	END DO
1246
1247	!CR:reprise du code AJ
1248	clw(i, k) = qnk(i) - qsat_new
1249	clw(i, k) = max(0.0, clw(i,k))
1250	tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
1251	! print*,tvp(i,k),'tvp'
1252	END IF
1253	IF (clw(i,k)<1.E-11) THEN
1254	tp(i, k) = tv(i, k)
1255	tvp(i, k) = tv(i, k)
1256	END IF
1257	END IF ! (cvflag_ice)
1258
1259	END DO
1260	END DO
1261
1262	! =====================================================================
1263	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
1264	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
1265	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
1266	! =====================================================================
1267
1268	IF (flag_epkeorig/=1) THEN
1269	DO k = 1, nl ! convect3
1270	DO i = 1, ncum
1271	pden = ptcrit - pbcrit
1272	ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1273	ep(i, k) = max(ep(i,k), 0.0)
1274	ep(i, k) = min(ep(i,k), epmax)
1275	sigp(i, k) = spfac
1276	END DO
1277	END DO
1278	ELSE
1279	DO k = 1, nl
1280	DO i = 1, ncum
1281	IF (k>=(nk(i)+1)) THEN
1282	tca = tp(i, k) - t0
1283	IF (tca>=0.0) THEN
1284	elacrit = elcrit
1285	ELSE
1286	elacrit = elcrit*(1.0-tca/tlcrit)
1287	END IF
1288	elacrit = max(elacrit, 0.0)
1289	ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
1290	ep(i, k) = max(ep(i,k), 0.0)
1291	ep(i, k) = min(ep(i,k), epmax)
1292	sigp(i, k) = spfac
1293	END IF
1294	END DO
1295	END DO
1296	END IF
1297	! =====================================================================
1298	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1299	! --- VIRTUAL TEMPERATURE
1300	! =====================================================================
1301
1302	! dans convect3, tvp est calcule en une seule fois, et sans retirer
1303	! l'eau condensee (~> reversible CAPE)
1304
1305	! ori do 340 k=minorig+1,nl
1306	! ori do 330 i=1,ncum
1307	! ori if(k.ge.(icb(i)+1))then
1308	! ori tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
1309	! oric print*,'i,k,tvp(i,k),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	! ori endif
1312	! ori 330 continue
1313	! ori 340 continue
1314
1315	! ori do 350 i=1,ncum
1316	! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1317	! ori 350 continue
1318
1319	DO i = 1, ncum ! convect3
1320	tp(i, nlp) = tp(i, nl) ! convect3
1321	END DO ! convect3
1322
1323	! =====================================================================
1324	! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1325	! =====================================================================
1326
1327	! -- this is for convect3 only:
1328
1329	! first estimate of buoyancy:
1330
1331	DO i = 1, ncum
1332	DO k = 1, nl
1333	buoy(i, k) = tvp(i, k) - tv(i, k)
1334	END DO
1335	END DO
1336
1337	! set buoyancy=buoybase for all levels below base
1338	! for safety, set buoy(icb)=buoybase
1339
1340	DO i = 1, ncum
1341	DO k = 1, nl
1342	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1343	buoy(i, k) = buoybase(i)
1344	END IF
1345	END DO
1346	! buoy(icb(i),k)=buoybase(i)
1347	buoy(i, icb(i)) = buoybase(i)
1348	END DO
1349
1350	! -- end convect3
1351
1352	! =====================================================================
1353	! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1354	! --- LEVEL OF NEUTRAL BUOYANCY
1355	! =====================================================================
1356
1357	! -- this is for convect3 only:
1358
1359	DO i = 1, ncum
1360	inb(i) = nl - 1
1361	iposit(i) = nl
1362	END DO
1363
1364
1365	! -- iposit(i) = first level, above icb, with positive buoyancy
1366	DO k = 1, nl - 1
1367	DO i = 1, ncum
1368	IF (k>=icb(i) .AND. buoy(i,k)>0.) THEN
1369	iposit(i) = min(iposit(i), k)
1370	END IF
1371	END DO
1372	END DO
1373
1374	DO i = 1, ncum
1375	IF (iposit(i)==nl) THEN
1376	iposit(i) = icb(i)
1377	END IF
1378	END DO
1379
1380	DO k = 1, nl - 1
1381	DO i = 1, ncum
1382	IF ((k>=iposit(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1383	inb(i) = min(inb(i), k)
1384	END IF
1385	END DO
1386	END DO
1387
1388	! -- end convect3
1389
1390	! ori do 510 i=1,ncum
1391	! ori cape(i)=0.0
1392	! ori capem(i)=0.0
1393	! ori inb(i)=icb(i)+1
1394	! ori inb1(i)=inb(i)
1395	! ori 510 continue
1396
1397	! Originial Code
1398
1399	! do 530 k=minorig+1,nl-1
1400	! do 520 i=1,ncum
1401	! if(k.ge.(icb(i)+1))then
1402	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1403	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1404	! cape(i)=cape(i)+by
1405	! if(by.ge.0.0)inb1(i)=k+1
1406	! if(cape(i).gt.0.0)then
1407	! inb(i)=k+1
1408	! capem(i)=cape(i)
1409	! endif
1410	! endif
1411	!520 continue
1412	!530 continue
1413	! do 540 i=1,ncum
1414	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1415	! cape(i)=capem(i)+byp
1416	! defrac=capem(i)-cape(i)
1417	! defrac=max(defrac,0.001)
1418	! frac(i)=-cape(i)/defrac
1419	! frac(i)=min(frac(i),1.0)
1420	! frac(i)=max(frac(i),0.0)
1421	!540 continue
1422
1423	! K Emanuel fix
1424
1425	! call zilch(byp,ncum)
1426	! do 530 k=minorig+1,nl-1
1427	! do 520 i=1,ncum
1428	! if(k.ge.(icb(i)+1))then
1429	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1430	! cape(i)=cape(i)+by
1431	! if(by.ge.0.0)inb1(i)=k+1
1432	! if(cape(i).gt.0.0)then
1433	! inb(i)=k+1
1434	! capem(i)=cape(i)
1435	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1436	! endif
1437	! endif
1438	!520 continue
1439	!530 continue
1440	! do 540 i=1,ncum
1441	! inb(i)=max(inb(i),inb1(i))
1442	! cape(i)=capem(i)+byp(i)
1443	! defrac=capem(i)-cape(i)
1444	! defrac=max(defrac,0.001)
1445	! frac(i)=-cape(i)/defrac
1446	! frac(i)=min(frac(i),1.0)
1447	! frac(i)=max(frac(i),0.0)
1448	!540 continue
1449
1450	! J Teixeira fix
1451
1452	! ori call zilch(byp,ncum)
1453	! ori do 515 i=1,ncum
1454	! ori lcape(i)=.true.
1455	! ori 515 continue
1456	! ori do 530 k=minorig+1,nl-1
1457	! ori do 520 i=1,ncum
1458	! ori if(cape(i).lt.0.0)lcape(i)=.false.
1459	! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1460	! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1461	! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1462	! ori cape(i)=cape(i)+by
1463	! ori if(by.ge.0.0)inb1(i)=k+1
1464	! ori if(cape(i).gt.0.0)then
1465	! ori inb(i)=k+1
1466	! ori capem(i)=cape(i)
1467	! ori endif
1468	! ori endif
1469	! ori 520 continue
1470	! ori 530 continue
1471	! ori do 540 i=1,ncum
1472	! ori cape(i)=capem(i)+byp(i)
1473	! ori defrac=capem(i)-cape(i)
1474	! ori defrac=max(defrac,0.001)
1475	! ori frac(i)=-cape(i)/defrac
1476	! ori frac(i)=min(frac(i),1.0)
1477	! ori frac(i)=max(frac(i),0.0)
1478	! ori 540 continue
1479
1480	! =====================================================================
1481	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1482	! =====================================================================
1483
1484	DO k = 1, nd
1485	DO i = 1, ncum
1486	hp(i, k) = h(i, k)
1487	END DO
1488	END DO
1489
1490	DO k = minorig + 1, nl
1491	DO i = 1, ncum
1492	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1493
1494	IF (cvflag_ice) THEN
1495	frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
1496	frac(i, k) = min(max(frac(i,k),0.0), 1.0)
1497	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac(i,k)lf(i,k))* &
1498	ep(i, k)*clw(i, k)
1499
1500	ELSE
1501	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k)
1502	END IF
1503
1504	END IF
1505	END DO
1506	END DO
1507
1508	RETURN
1509	END SUBROUTINE cv3_undilute2
1510
1511	SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
1512	pbase, p, ph, tv, buoy, &
1513	sig, w0, cape, m, iflag)
1514	IMPLICIT NONE
1515
1516	! ===================================================================
1517	! --- CLOSURE OF CONVECT3
1518	!
1519	! vectorization: S. Bony
1520	! ===================================================================
1521
1522	include "cvthermo.h"
1523	include "cv3param.h"
1524
1525	!input:
1526	INTEGER ncum, nd, nloc
1527	INTEGER icb(nloc), inb(nloc)
1528	REAL pbase(nloc)
1529	REAL p(nloc, nd), ph(nloc, nd+1)
1530	REAL tv(nloc, nd), buoy(nloc, nd)
1531
1532	!input/output:
1533	REAL sig(nloc, nd), w0(nloc, nd)
1534	INTEGER iflag(nloc)
1535
1536	!output:
1537	REAL cape(nloc)
1538	REAL m(nloc, nd)
1539
1540	!local variables:
1541	INTEGER i, j, k, icbmax
1542	REAL deltap, fac, w, amu
1543	REAL dtmin(nloc, nd), sigold(nloc, nd)
1544	REAL cbmflast(nloc)
1545
1546
1547	! -------------------------------------------------------
1548	! -- Initialization
1549	! -------------------------------------------------------
1550
1551	DO k = 1, nl
1552	DO i = 1, ncum
1553	m(i, k) = 0.0
1554	END DO
1555	END DO
1556
1557	! -------------------------------------------------------
1558	! -- Reset sig(i) and w0(i) for i>inb and i<icb
1559	! -------------------------------------------------------
1560
1561	! update sig and w0 above LNB:
1562
1563	DO k = 1, nl - 1
1564	DO i = 1, ncum
1565	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
1566	sig(i, k) = beta*sig(i, k) + &
1567	2.alphabuoy(i, inb(i))*abs(buoy(i,inb(i)))
1568	sig(i, k) = amax1(sig(i,k), 0.0)
1569	w0(i, k) = beta*w0(i, k)
1570	END IF
1571	END DO
1572	END DO
1573
1574	! compute icbmax:
1575
1576	icbmax = 2
1577	DO i = 1, ncum
1578	icbmax = max(icbmax, icb(i))
1579	END DO
1580
1581	! update sig and w0 below cloud base:
1582
1583	DO k = 1, icbmax
1584	DO i = 1, ncum
1585	IF (k<=icb(i)) THEN
1586	sig(i, k) = beta*sig(i, k) - &
1587	2.alphabuoy(i, icb(i))*buoy(i, icb(i))
1588	sig(i, k) = max(sig(i,k), 0.0)
1589	w0(i, k) = beta*w0(i, k)
1590	END IF
1591	END DO
1592	END DO
1593
1594	!! if(inb.lt.(nl-1))then
1595	!! do 85 i=inb+1,nl-1
1596	!! sig(i)=betasig(i)+2.alphabuoy(inb)
1597	!! 1 abs(buoy(inb))
1598	!! sig(i)=max(sig(i),0.0)
1599	!! w0(i)=beta*w0(i)
1600	!! 85 continue
1601	!! end if
1602
1603	!! do 87 i=1,icb
1604	!! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
1605	!! sig(i)=max(sig(i),0.0)
1606	!! w0(i)=beta*w0(i)
1607	!! 87 continue
1608
1609	! -------------------------------------------------------------
1610	! -- Reset fractional areas of updrafts and w0 at initial time
1611	! -- and after 10 time steps of no convection
1612	! -------------------------------------------------------------
1613
1614	DO k = 1, nl - 1
1615	DO i = 1, ncum
1616	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
1617	sig(i, k) = 0.0
1618	w0(i, k) = 0.0
1619	END IF
1620	END DO
1621	END DO
1622
1623	! -------------------------------------------------------------
1624	! -- Calculate convective available potential energy (cape),
1625	! -- vertical velocity (w), fractional area covered by
1626	! -- undilute updraft (sig), and updraft mass flux (m)
1627	! -------------------------------------------------------------
1628
1629	DO i = 1, ncum
1630	cape(i) = 0.0
1631	END DO
1632
1633	! compute dtmin (minimum buoyancy between ICB and given level k):
1634
1635	DO i = 1, ncum
1636	DO k = 1, nl
1637	dtmin(i, k) = 100.0
1638	END DO
1639	END DO
1640
1641	DO i = 1, ncum
1642	DO k = 1, nl
1643	DO j = minorig, nl
1644	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
1645	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
1646	END IF
1647	END DO
1648	END DO
1649	END DO
1650
1651	! the interval on which cape is computed starts at pbase :
1652
1653	DO k = 1, nl
1654	DO i = 1, ncum
1655
1656	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
1657
1658	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
1659	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
1660	cape(i) = amax1(0.0, cape(i))
1661	sigold(i, k) = sig(i, k)
1662
1663	! dtmin(i,k)=100.0
1664	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
1665	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1666	! 97 continue
1667
1668	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
1669	sig(i, k) = max(sig(i,k), 0.0)
1670	sig(i, k) = amin1(sig(i,k), 0.01)
1671	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
1672	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
1673	amu = 0.5(sig(i,k)+sigold(i,k))w
1674	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
1675	w0(i, k) = w
1676	END IF
1677
1678	END DO
1679	END DO
1680
1681	DO i = 1, ncum
1682	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
1683	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))
1684	sig(i, icb(i)) = sig(i, icb(i)+1)
1685	sig(i, icb(i)-1) = sig(i, icb(i))
1686	END DO
1687
1688	! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
1689	! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
1690	! ccc the final mass flux (cbmflast) is greater than the target mass flux
1691	! ccc (cbmf) ??).
1692	! cc
1693	! c do i = 1,ncum
1694	! c cbmflast(i) = 0.
1695	! c enddo
1696	! cc
1697	! c do k= 1,nl
1698	! c do i = 1,ncum
1699	! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
1700	! c cbmflast(i) = cbmflast(i)+M(i,k)
1701	! c ENDIF
1702	! c enddo
1703	! c enddo
1704	! cc
1705	! c do i = 1,ncum
1706	! c IF (cbmflast(i) .lt. 1.e-6) THEN
1707	! c iflag(i) = 3
1708	! c ENDIF
1709	! c enddo
1710	! cc
1711	! c do k= 1,nl
1712	! c do i = 1,ncum
1713	! c IF (iflag(i) .ge. 3) THEN
1714	! c M(i,k) = 0.
1715	! c sig(i,k) = 0.
1716	! c w0(i,k) = 0.
1717	! c ENDIF
1718	! c enddo
1719	! c enddo
1720	! cc
1721	!! cape=0.0
1722	!! do 98 i=icb+1,inb
1723	!! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1724	!! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
1725	!! dcape=rrdbuoy(i-1)deltap/p(i-1)
1726	!! dlnp=deltap/p(i-1)
1727	!! cape=max(0.0,cape)
1728	!! sigold=sig(i)
1729
1730	!! dtmin=100.0
1731	!! do 97 j=icb,i-1
1732	!! dtmin=amin1(dtmin,buoy(j))
1733	!! 97 continue
1734
1735	!! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
1736	!! sig(i)=max(sig(i),0.0)
1737	!! sig(i)=amin1(sig(i),0.01)
1738	!! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1739	!! w=(1.-beta)facsqrt(cape)+beta*w0(i)
1740	!! amu=0.5(sig(i)+sigold)w
1741	!! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
1742	!! w0(i)=w
1743	!! 98 continue
1744	!! w0(icb)=0.5*w0(icb+1)
1745	!! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1746	!! sig(icb)=sig(icb+1)
1747	!! sig(icb-1)=sig(icb)
1748
1749	RETURN
1750	END SUBROUTINE cv3_closure
1751
1752	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
1753	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
1754	unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
1755	ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
1756	IMPLICIT NONE
1757
1758	! ---------------------------------------------------------------------
1759	! a faire:
1760	! - vectorisation de la partie normalisation des flux (do 789...)
1761	! ---------------------------------------------------------------------
1762
1763	include "cvthermo.h"
1764	include "cv3param.h"
1765	include "cvflag.h"
1766
1767	!inputs:
1768	INTEGER ncum, nd, na, ntra, nloc
1769	INTEGER icb(nloc), inb(nloc), nk(nloc)
1770	REAL sig(nloc, nd)
1771	REAL qnk(nloc), unk(nloc), vnk(nloc)
1772	REAL ph(nloc, nd+1)
1773	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1774	REAL u(nloc, nd), v(nloc, nd)
1775	REAL tra(nloc, nd, ntra) ! input of convect3
1776	REAL lv(nloc, na), h(nloc, na), hp(nloc, na)
1777	REAL lf(nloc, na), frac(nloc, na)
1778	REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
1779	REAL m(nloc, na) ! input of convect3
1780
1781	!outputs:
1782	REAL ment(nloc, na, na), qent(nloc, na, na)
1783	REAL uent(nloc, na, na), vent(nloc, na, na)
1784	REAL sij(nloc, na, na), elij(nloc, na, na)
1785	REAL traent(nloc, nd, nd, ntra)
1786	REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
1787	REAL sigij(nloc, nd, nd)
1788	INTEGER nent(nloc, nd)
1789
1790	!local variables:
1791	INTEGER i, j, k, il, im, jm
1792	INTEGER num1, num2
1793	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1794	REAL alt, smid, sjmin, sjmax, delp, delm
1795	REAL asij(nloc), smax(nloc), scrit(nloc)
1796	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
1797	REAL wgh
1798	REAL zm(nloc, na)
1799	LOGICAL lwork(nloc)
1800
1801	! =====================================================================
1802	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1803	! =====================================================================
1804
1805	! ori do 360 i=1,ncum*nlp
1806	DO j = 1, nl
1807	DO i = 1, ncum
1808	nent(i, j) = 0
1809	! in convect3, m is computed in cv3_closure
1810	! ori m(i,1)=0.0
1811	END DO
1812	END DO
1813
1814	! ori do 400 k=1,nlp
1815	! ori do 390 j=1,nlp
1816	DO j = 1, nl
1817	DO k = 1, nl
1818	DO i = 1, ncum
1819	qent(i, k, j) = rr(i, j)
1820	uent(i, k, j) = u(i, j)
1821	vent(i, k, j) = v(i, j)
1822	elij(i, k, j) = 0.0
1823	!ym ment(i,k,j)=0.0
1824	!ym sij(i,k,j)=0.0
1825	END DO
1826	END DO
1827	END DO
1828
1829	!ym
1830	ment(1:ncum, 1:nd, 1:nd) = 0.0
1831	sij(1:ncum, 1:nd, 1:nd) = 0.0
1832
1833	!AC! do k=1,ntra
1834	!AC! do j=1,nd ! instead nlp
1835	!AC! do i=1,nd ! instead nlp
1836	!AC! do il=1,ncum
1837	!AC! traent(il,i,j,k)=tra(il,j,k)
1838	!AC! enddo
1839	!AC! enddo
1840	!AC! enddo
1841	!AC! enddo
1842	zm(:, :) = 0.
1843
1844	! =====================================================================
1845	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1846	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1847	! --- FRACTION (sij)
1848	! =====================================================================
1849
1850	DO i = minorig + 1, nl
1851
1852	DO j = minorig, nl
1853	DO il = 1, ncum
1854	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
1855
1856	rti = qnk(il) - ep(il, i)*clw(il, i)
1857	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
1858
1859
1860	IF (cvflag_ice) THEN
1861	! print*,cvflag_ice,'cvflag_ice dans do 700'
1862	IF (t(il,j)<=263.15) THEN
1863	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
1864	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
1865	END IF
1866	END IF
1867
1868	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
1869	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
1870	dei = denom
1871	IF (abs(dei)<0.01) dei = 0.01
1872	sij(il, i, j) = anum/dei
1873	sij(il, i, i) = 1.0
1874	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
1875	altem = altem/bf2
1876	cwat = clw(il, j)*(1.-ep(il,j))
1877	stemp = sij(il, i, j)
1878	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
1879
1880	IF (cvflag_ice) THEN
1881	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
1882	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
1883	ELSE
1884	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
1885	denom = denom + lv(il, j)*(rr(il,i)-rti)
1886	END IF
1887
1888	IF (abs(denom)<0.01) denom = 0.01
1889	sij(il, i, j) = anum/denom
1890	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
1891	altem = altem - (bf2-1.)*cwat
1892	END IF
1893	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
1894	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
1895	uent(il, i, j) = sij(il, i, j)u(il, i) + (1.-sij(il,i,j))unk(il)
1896	vent(il, i, j) = sij(il, i, j)v(il, i) + (1.-sij(il,i,j))vnk(il)
1897	!!!! do k=1,ntra
1898	!!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1899	!!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1900	!!!! end do
1901	elij(il, i, j) = altem
1902	elij(il, i, j) = max(0.0, elij(il,i,j))
1903	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
1904	nent(il, i) = nent(il, i) + 1
1905	END IF
1906	sij(il, i, j) = max(0.0, sij(il,i,j))
1907	sij(il, i, j) = amin1(1.0, sij(il,i,j))
1908	END IF ! new
1909	END DO
1910	END DO
1911
1912	!AC! do k=1,ntra
1913	!AC! do j=minorig,nl
1914	!AC! do il=1,ncum
1915	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1916	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1917	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1918	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1919	!AC! endif
1920	!AC! enddo
1921	!AC! enddo
1922	!AC! enddo
1923
1924
1925	! * if no air can entrain at level i assume that updraft detrains *
1926	! * at that level and calculate detrained air flux and properties *
1927
1928
1929	! @ do 170 i=icb(il),inb(il)
1930
1931	DO il = 1, ncum
1932	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
1933	! @ if(nent(il,i).eq.0)then
1934	ment(il, i, i) = m(il, i)
1935	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
1936	uent(il, i, i) = unk(il)
1937	vent(il, i, i) = vnk(il)
1938	elij(il, i, i) = clw(il, i)
1939	! MAF sij(il,i,i)=1.0
1940	sij(il, i, i) = 0.0
1941	END IF
1942	END DO
1943	END DO
1944
1945	!AC! do j=1,ntra
1946	!AC! do i=minorig+1,nl
1947	!AC! do il=1,ncum
1948	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1949	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
1950	!AC! endif
1951	!AC! enddo
1952	!AC! enddo
1953	!AC! enddo
1954
1955	DO j = minorig, nl
1956	DO i = minorig, nl
1957	DO il = 1, ncum
1958	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
1959	sigij(il, i, j) = sij(il, i, j)
1960	END IF
1961	END DO
1962	END DO
1963	END DO
1964	! @ enddo
1965
1966	! @170 continue
1967
1968	! =====================================================================
1969	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
1970	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
1971	! =====================================================================
1972
1973	CALL zilch(asum, nloc*nd)
1974	CALL zilch(csum, nloc*nd)
1975	CALL zilch(csum, nloc*nd)
1976
1977	DO il = 1, ncum
1978	lwork(il) = .FALSE.
1979	END DO
1980
1981	DO i = minorig + 1, nl
1982
1983	num1 = 0
1984	DO il = 1, ncum
1985	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
1986	END DO
1987	IF (num1<=0) GO TO 789
1988
1989
1990	DO il = 1, ncum
1991	IF (i>=icb(il) .AND. i<=inb(il)) THEN
1992	lwork(il) = (nent(il,i)/=0)
1993	qp = qnk(il) - ep(il, i)*clw(il, i)
1994
1995	IF (cvflag_ice) THEN
1996
1997	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
1998	(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
1999	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
2000	(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
2001	ELSE
2002
2003	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
2004	(cpv-cpd)t(il, i)(qp-rr(il,i))
2005	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
2006	(cpd-cpv)t(il, i)(rr(il,i)-qp)
2007	END IF
2008
2009	IF (abs(denom)<0.01) denom = 0.01
2010	scrit(il) = anum/denom
2011	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
2012	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
2013	smax(il) = 0.0
2014	asij(il) = 0.0
2015	END IF
2016	END DO
2017
2018	DO j = nl, minorig, -1
2019
2020	num2 = 0
2021	DO il = 1, ncum
2022	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2023	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2024	lwork(il)) num2 = num2 + 1
2025	END DO
2026	IF (num2<=0) GO TO 175
2027
2028	DO il = 1, ncum
2029	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2030	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2031	lwork(il)) THEN
2032
2033	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
2034	wgh = 1.0
2035	IF (j>i) THEN
2036	sjmax = max(sij(il,i,j+1), smax(il))
2037	sjmax = amin1(sjmax, scrit(il))
2038	smax(il) = max(sij(il,i,j), smax(il))
2039	sjmin = max(sij(il,i,j-1), smax(il))
2040	sjmin = amin1(sjmin, scrit(il))
2041	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
2042	smid = amin1(sij(il,i,j), scrit(il))
2043	ELSE
2044	sjmax = max(sij(il,i,j+1), scrit(il))
2045	smid = max(sij(il,i,j), scrit(il))
2046	sjmin = 0.0
2047	IF (j>1) sjmin = sij(il, i, j-1)
2048	sjmin = max(sjmin, scrit(il))
2049	END IF
2050	delp = abs(sjmax-smid)
2051	delm = abs(sjmin-smid)
2052	asij(il) = asij(il) + wgh*(delp+delm)
2053	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
2054	END IF
2055	END IF
2056	END DO
2057
2058	175 END DO
2059
2060	DO il = 1, ncum
2061	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2062	asij(il) = max(1.0E-16, asij(il))
2063	asij(il) = 1.0/asij(il)
2064	asum(il, i) = 0.0
2065	bsum(il, i) = 0.0
2066	csum(il, i) = 0.0
2067	END IF
2068	END DO
2069
2070	DO j = minorig, nl
2071	DO il = 1, ncum
2072	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2073	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2074	ment(il, i, j) = ment(il, i, j)*asij(il)
2075	END IF
2076	END DO
2077	END DO
2078
2079	DO j = minorig, nl
2080	DO il = 1, ncum
2081	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2082	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2083	asum(il, i) = asum(il, i) + ment(il, i, j)
2084	ment(il, i, j) = ment(il, i, j)*sig(il, j)
2085	bsum(il, i) = bsum(il, i) + ment(il, i, j)
2086	END IF
2087	END DO
2088	END DO
2089
2090	DO il = 1, ncum
2091	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2092	bsum(il, i) = max(bsum(il,i), 1.0E-16)
2093	bsum(il, i) = 1.0/bsum(il, i)
2094	END IF
2095	END DO
2096
2097	DO j = minorig, nl
2098	DO il = 1, ncum
2099	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2100	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2101	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
2102	END IF
2103	END DO
2104	END DO
2105
2106	DO j = minorig, nl
2107	DO il = 1, ncum
2108	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2109	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2110	csum(il, i) = csum(il, i) + ment(il, i, j)
2111	END IF
2112	END DO
2113	END DO
2114
2115	DO il = 1, ncum
2116	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2117	csum(il,i)<m(il,i)) THEN
2118	nent(il, i) = 0
2119	ment(il, i, i) = m(il, i)
2120	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2121	uent(il, i, i) = unk(il)
2122	vent(il, i, i) = vnk(il)
2123	elij(il, i, i) = clw(il, i)
2124	! MAF sij(il,i,i)=1.0
2125	sij(il, i, i) = 0.0
2126	END IF
2127	END DO ! il
2128
2129	!AC! do j=1,ntra
2130	!AC! do il=1,ncum
2131	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
2132	!AC! : .and. csum(il,i).lt.m(il,i) ) then
2133	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2134	!AC! endif
2135	!AC! enddo
2136	!AC! enddo
2137	789 END DO
2138
2139	! MAF: renormalisation de MENT
2140	CALL zilch(zm, nloc*na)
2141	DO jm = 1, nd
2142	DO im = 1, nd
2143	DO il = 1, ncum
2144	zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
2145	END DO
2146	END DO
2147	END DO
2148
2149	DO jm = 1, nd
2150	DO im = 1, nd
2151	DO il = 1, ncum
2152	IF (zm(il,im)/=0.) THEN
2153	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
2154	END IF
2155	END DO
2156	END DO
2157	END DO
2158
2159	DO jm = 1, nd
2160	DO im = 1, nd
2161	DO il = 1, ncum
2162	qents(il, im, jm) = qent(il, im, jm)
2163	ments(il, im, jm) = ment(il, im, jm)
2164	END DO
2165	END DO
2166	END DO
2167
2168	RETURN
2169	END SUBROUTINE cv3_mixing
2170
2171	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
2172	t, rr, rs, gz, u, v, tra, p, ph, &
2173	th, tv, lv, lf, cpn, ep, sigp, clw, &
2174	m, ment, elij, delt, plcl, coef_clos, &
2175	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
2176	faci, b, sigd, &
2177	wdtrainA, wdtrainM) ! RomP
2178	IMPLICIT NONE
2179
2180
2181	include "cvthermo.h"
2182	include "cv3param.h"
2183	include "cvflag.h"
2184
2185	!inputs:
2186	INTEGER ncum, nd, na, ntra, nloc
2187	INTEGER icb(nloc), inb(nloc)
2188	REAL delt, plcl(nloc)
2189	REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd), gz(nloc, na)
2190	REAL u(nloc, nd), v(nloc, nd)
2191	REAL tra(nloc, nd, ntra)
2192	REAL p(nloc, nd), ph(nloc, nd+1)
2193	REAL ep(nloc, na), sigp(nloc, na), clw(nloc, na)
2194	REAL th(nloc, na), tv(nloc, na), lv(nloc, na), cpn(nloc, na)
2195	REAL lf(nloc, na)
2196	REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
2197	REAL coef_clos(nloc)
2198
2199	!input/output
2200	INTEGER iflag(nloc)
2201
2202	!outputs:
2203	REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
2204	REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
2205	REAL ice(nloc, na), fondue(nloc, na), faci(nloc, na)
2206	REAL trap(nloc, na, ntra)
2207	REAL b(nloc, na), sigd(nloc)
2208	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
2209	! de l ascendance adiabatique et des flux melanges Pa et Pm.
2210	! Distinction des wdtrain
2211	! Pa = wdtrainA Pm = wdtrainM
2212	REAL wdtrainA(nloc, na), wdtrainM(nloc, na)
2213
2214	!local variables
2215	INTEGER i, j, k, il, num1, ndp1
2216	REAL tinv, delti, coef
2217	REAL awat, afac, afac1, afac2, bfac
2218	REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
2219	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
2220	REAL ampmax, thaw
2221	REAL tevap(nloc)
2222	REAL lvcp(nloc, na), lfcp(nloc, na)
2223	REAL h(nloc, na), hm(nloc, na)
2224	REAL frac(nloc, na)
2225	REAL fraci(nloc, na), prec(nloc, na)
2226	REAL wdtrain(nloc)
2227	LOGICAL lwork(nloc), mplus(nloc)
2228
2229
2230	! ------------------------------------------------------
2231
2232	delti = 1./delt
2233	tinv = 1./3.
2234
2235	mp(:, :) = 0.
2236
2237	DO i = 1, nl
2238	DO il = 1, ncum
2239	mp(il, i) = 0.0
2240	rp(il, i) = rr(il, i)
2241	up(il, i) = u(il, i)
2242	vp(il, i) = v(il, i)
2243	wt(il, i) = 0.001
2244	water(il, i) = 0.0
2245	frac(il, i) = 0.0
2246	faci(il, i) = 0.0
2247	fraci(il, i) = 0.0
2248	ice(il, i) = 0.0
2249	prec(il, i) = 0.0
2250	fondue(il, i) = 0.0
2251	evap(il, i) = 0.0
2252	b(il, i) = 0.0
2253	lvcp(il, i) = lv(il, i)/cpn(il, i)
2254	lfcp(il, i) = lf(il, i)/cpn(il, i)
2255	END DO
2256	END DO
2257	!AC! do k=1,ntra
2258	!AC! do i=1,nd
2259	!AC! do il=1,ncum
2260	!AC! trap(il,i,k)=tra(il,i,k)
2261	!AC! enddo
2262	!AC! enddo
2263	!AC! enddo
2264	!! RomP >>>
2265	DO i = 1, nd
2266	DO il = 1, ncum
2267	wdtrainA(il, i) = 0.0
2268	wdtrainM(il, i) = 0.0
2269	END DO
2270	END DO
2271	!! RomP <<<
2272
2273	! * check whether ep(inb)=0, if so, skip precipitating *
2274	! * downdraft calculation *
2275
2276
2277	DO il = 1, ncum
2278	!! lwork(il)=.TRUE.
2279	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
2280	lwork(il) = ep(il, inb(il)) >= 0.0001
2281	END DO
2282
2283	! *** Set the fractionnal area sigd of precipitating downdraughts
2284	DO il = 1, ncum
2285	sigd(il) = sigdz*coef_clos(il)
2286	END DO
2287
2288
2289	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2290	!
2291	! * begin downdraft loop *
2292	!
2293	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2294
2295	DO i = nl + 1, 1, -1
2296
2297	num1 = 0
2298	DO il = 1, ncum
2299	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
2300	END DO
2301	IF (num1<=0) GO TO 400
2302
2303	CALL zilch(wdtrain, ncum)
2304
2305
2306	! * integrate liquid water equation to find condensed water *
2307	! * and condensed water flux *
2308	!
2309	!
2310	! * calculate detrained precipitation *
2311
2312	DO il = 1, ncum
2313	IF (i<=inb(il) .AND. lwork(il)) THEN
2314	IF (cvflag_grav) THEN
2315	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
2316	wdtrainA(il, i) = wdtrain(il)/grav ! Pa RomP
2317	ELSE
2318	wdtrain(il) = 10.0ep(il, i)m(il, i)*clw(il, i)
2319	wdtrainA(il, i) = wdtrain(il)/10. ! Pa RomP
2320	END IF
2321	END IF
2322	END DO
2323
2324	IF (i>1) THEN
2325	DO j = 1, i - 1
2326	DO il = 1, ncum
2327	IF (i<=inb(il) .AND. lwork(il)) THEN
2328	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
2329	awat = max(awat, 0.0)
2330	IF (cvflag_grav) THEN
2331	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
2332	wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP
2333	ELSE
2334	wdtrain(il) = wdtrain(il) + 10.0awatment(il, j, i)
2335	wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i) ! Pm RomP
2336	END IF
2337	END IF
2338	END DO
2339	END DO
2340	END IF
2341
2342
2343	! * find rain water and evaporation using provisional *
2344	! * estimates of rp(i)and rp(i-1) *
2345
2346
2347	DO il = 1, ncum
2348	IF (i<=inb(il) .AND. lwork(il)) THEN
2349
2350	wt(il, i) = 45.0
2351
2352	IF (cvflag_ice) THEN
2353	frac(il, inb(il)) = 1. - (t(il,inb(il))-243.15)/(263.15-243.15)
2354	frac(il, inb(il)) = min(max(frac(il,inb(il)),0.), 1.)
2355	fraci(il, inb(il)) = frac(il, inb(il))
2356	ELSE
2357	CONTINUE
2358	END IF
2359
2360	IF (i<inb(il)) THEN
2361
2362	IF (cvflag_ice) THEN
2363	thaw = (t(il,i)-273.15)/(275.15-273.15)
2364	thaw = min(max(thaw,0.0), 1.0)
2365	frac(il, i) = frac(il, i)*(1.-thaw)
2366	ELSE
2367	CONTINUE
2368	END IF
2369
2370	rp(il, i) = rp(il, i+1) + &
2371	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
2372	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
2373	END IF
2374	fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
2375	fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
2376	rp(il, i) = max(rp(il,i), 0.0)
2377	rp(il, i) = amin1(rp(il,i), rs(il,i))
2378	rp(il, inb(il)) = rr(il, inb(il))
2379
2380	IF (i==1) THEN
2381	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
2382	IF (cvflag_ice) THEN
2383	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
2384	END IF
2385	ELSE
2386	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)
2387	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
2388	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
2389	rp(il, i-1) = max(rp(il,i-1), 0.0)
2390	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
2391	afac2 = p(il, i-1)(rs(il,i-1)-rp(il,i-1))/(1.0E4+2000.0p(il,i-1)*rs(il,i-1))
2392	afac = 0.5*(afac1+afac2)
2393	END IF
2394	IF (i==inb(il)) afac = 0.0
2395	afac = max(afac, 0.0)
2396	bfac = 1./(sigd(il)*wt(il,i))
2397
2398	!JYG1
2399	! cc sigt=1.0
2400	! cc if(i.ge.icb)sigt=sigp(i)
2401	! prise en compte de la variation progressive de sigt dans
2402	! les couches icb et icb-1:
2403	! pour plcl<ph(i+1), pr1=0 & pr2=1
2404	! pour plcl>ph(i), pr1=1 & pr2=0
2405	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2406	! sur le nuage, et pr2 est la proportion sous la base du
2407	! nuage.
2408	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2409	pr1 = max(0., min(1.,pr1))
2410	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2411	pr2 = max(0., min(1.,pr2))
2412	sigt = sigp(il, i)*pr1 + pr2
2413	!JYG2
2414
2415	!JYG----
2416	! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2417	! c6 = water(il,i+1) + wdtrain(il)*bfac
2418	! c6 = prec(il,i+1) + wdtrain(il)*bfac
2419	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
2420	! evap(il,i)=sigtafacrevap
2421	! water(il,i)=revap*revap
2422	! prec(il,i)=revap*revap
2423	!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
2424	!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
2425	!!---end jyg---
2426
2427	! --------retour à la formulation originale d'Emanuel.
2428	IF (cvflag_ice) THEN
2429
2430	! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2431	! c6=prec(il,i+1)+bfac*wdtrain(il) &
2432	! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2433	! if(c6.gt.0.0)then
2434	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
2435
2436	!JAM Attention: evap=sigt*E
2437	! Modification: evap devient l'évaporation en milieu de couche
2438	! car nécessaire dans cv3_yield
2439	! Du coup, il faut modifier pas mal d'équations...
2440	! et l'expression de afac qui devient afac1
2441	! revap=sqrt((prec(i+1)+prec(i))/2)
2442
2443	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
2444	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
2445	! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
2446	! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
2447	! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
2448	IF (c6>b6*b6+1.E-20) THEN
2449	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
2450	ELSE
2451	revap = (-b6+sqrt(b6b6+4.c6))/2.
2452	END IF
2453	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
2454	! print*,prec(il,i),'neige'
2455
2456	!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
2457	! c evap(il,i)=sigtafacrevap
2458	! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
2459	! Ici,l'evaporation evap est simplement calculee par l'equation de
2460	! conservation.
2461	! prec(il,i)=revap*revap
2462	! else
2463	!JYG---- Correction : si c6 <= 0, water(il,i)=0.
2464	! prec(il,i)=0.
2465	! endif
2466
2467	!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
2468	! moins [tt ce qui sort de la couche i]
2469	! print *, 'evap avec ice'
2470	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
2471	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2472
2473	d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
2474	e6 = bfac*wdtrain(il)
2475	f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
2476
2477	thaw = (t(il,i)-273.15)/(275.15-273.15)
2478	thaw = min(max(thaw,0.0), 1.0)
2479	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
2480	water(il, i) = max(water(il,i), 0.)
2481	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
2482	ice(il, i) = max(ice(il,i), 0.)
2483	fondue(il, i) = ice(il, i)*thaw
2484	water(il, i) = water(il, i) + fondue(il, i)
2485	ice(il, i) = ice(il, i) - fondue(il, i)
2486
2487	IF (water(il,i)+ice(il,i)<1.E-30) THEN
2488	faci(il, i) = 0.
2489	ELSE
2490	faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
2491	END IF
2492
2493	! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
2494	! water(il,i)=max(water(il,i),0.)
2495	! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
2496	! ice(il,i)=max(ice(il,i),0.)
2497	! fondue(il,i)=ice(il,i)*thaw
2498	! water(il,i)=water(il,i)+fondue(il,i)
2499	! ice(il,i)=ice(il,i)-fondue(il,i)
2500
2501	! if((water(il,i)+ice(il,i)).lt.1.e-30)then
2502	! faci(il,i)=0.
2503	! else
2504	! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
2505	! endif
2506
2507	ELSE
2508	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2509	c6 = water(il, i+1) + bfac*wdtrain(il) - &
2510	50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
2511	IF (c6>0.0) THEN
2512	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
2513	water(il, i) = revap*revap
2514	ELSE
2515	water(il, i) = 0.
2516	END IF
2517	! print *, 'evap sans ice'
2518	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
2519	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2520
2521	END IF
2522	END IF !(i.le.inb(il) .and. lwork(il))
2523	END DO
2524	! ----------------------------------------------------------------
2525
2526	! cc
2527	! * calculate precipitating downdraft mass flux under *
2528	! * hydrostatic approximation *
2529
2530	DO il = 1, ncum
2531	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
2532
2533	tevap(il) = max(0.0, evap(il,i))
2534	delth = max(0.001, (th(il,i)-th(il,i-1)))
2535	IF (cvflag_ice) THEN
2536	IF (cvflag_grav) THEN
2537	mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
2538	(p(il,i-1)-p(il,i))/delth + &
2539	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
2540	(p(il,i-1)-p(il,i))/delth + &
2541	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
2542	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
2543	ELSE
2544	mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
2545	(p(il,i-1)-p(il,i))/delth + &
2546	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
2547	(p(il,i-1)-p(il,i))/delth + &
2548	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
2549	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
2550
2551	END IF
2552	ELSE
2553	IF (cvflag_grav) THEN
2554	mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
2555	(p(il,i-1)-p(il,i))/delth
2556	ELSE
2557	mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
2558	(p(il,i-1)-p(il,i))/delth
2559	END IF
2560
2561	END IF
2562
2563	END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2564	END DO
2565	! ----------------------------------------------------------------
2566
2567	! * if hydrostatic assumption fails, *
2568	! * solve cubic difference equation for downdraft theta *
2569	! * and mass flux from two simultaneous differential eqns *
2570
2571	DO il = 1, ncum
2572	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
2573
2574	amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
2575	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2576	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2577
2578	IF (amp2>(0.1*amfac)) THEN
2579	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
2580	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
2581	(lvcp(il,i)sigd(il)th(il,i))
2582	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
2583
2584	IF (cvflag_ice) THEN
2585	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2586	50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
2587	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
2588	ELSE
2589
2590	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2591	50.(p(il,i-1)-p(il,i))xf*tevap(il)
2592	END IF
2593
2594	fac2 = 1.0
2595	IF (bf<0.0) fac2 = -1.0
2596	bf = abs(bf)
2597	ur = 0.25bfbf - afafaftinvtinv*tinv
2598	IF (ur>=0.0) THEN
2599	sru = sqrt(ur)
2600	fac = 1.0
2601	IF ((0.5*bf-sru)<0.0) fac = -1.0
2602	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
2603	fac(abs(0.5bf-sru))**tinv
2604	ELSE
2605	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
2606	IF (fac2<0.0) d = 3.14159 - d
2607	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
2608	END IF
2609	mp(il, i) = max(0.0, mp(il,i))
2610
2611	IF (cvflag_ice) THEN
2612	IF (cvflag_grav) THEN
2613	!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2614	! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
2615	! Et il faut bien revoir les facteurs 100.
2616	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i)) &
2617	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
2618	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
2619	(ph(il,i)-ph(il,i+1))) / &
2620	(mp(il,i)+sigd(il)*0.1) - &
2621	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2622	(lvcp(il,i)sigd(il)th(il,i))
2623	ELSE
2624	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
2625	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
2626	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
2627	(ph(il,i)-ph(il,i+1))) / &
2628	(mp(il,i)+sigd(il)*0.1) - &
2629	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2630	(lvcp(il,i)sigd(il)th(il,i))
2631	END IF
2632	ELSE
2633	IF (cvflag_grav) THEN
2634	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
2635	(mp(il,i)+sigd(il)*0.1) - &
2636	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2637	(lvcp(il,i)sigd(il)th(il,i))
2638	ELSE
2639	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
2640	(mp(il,i)+sigd(il)*0.1) - &
2641	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2642	(lvcp(il,i)sigd(il)th(il,i))
2643	END IF
2644	END IF
2645	b(il, i-1) = max(b(il,i-1), 0.0)
2646
2647	END IF !(amp2.gt.(0.1*amfac))
2648
2649	! * limit magnitude of mp(i) to meet cfl condition *
2650
2651	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
2652	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
2653	ampmax = min(ampmax, amp2)
2654	mp(il, i) = min(mp(il,i), ampmax)
2655
2656	! * force mp to decrease linearly to zero *
2657	! * between cloud base and the surface *
2658
2659
2660	! c if(p(il,i).gt.p(il,icb(il)))then
2661	! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2662	! c endif
2663	IF (ph(il,i)>0.9*plcl(il)) THEN
2664	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
2665	END IF
2666
2667	END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2668	END DO
2669	! ----------------------------------------------------------------
2670
2671	! * find mixing ratio of precipitating downdraft *
2672
2673	DO il = 1, ncum
2674	IF (i<inb(il) .AND. lwork(il)) THEN
2675	mplus(il) = mp(il, i) > mp(il, i+1)
2676	END IF ! (i.lt.inb(il) .and. lwork(il))
2677	END DO
2678
2679	DO il = 1, ncum
2680	IF (i<inb(il) .AND. lwork(il)) THEN
2681
2682	rp(il, i) = rr(il, i)
2683
2684	IF (mplus(il)) THEN
2685
2686	IF (cvflag_grav) THEN
2687	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
2688	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2689	ELSE
2690	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
2691	5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2692	END IF
2693	rp(il, i) = rp(il, i)/mp(il, i)
2694	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
2695	up(il, i) = up(il, i)/mp(il, i)
2696	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
2697	vp(il, i) = vp(il, i)/mp(il, i)
2698
2699	ELSE ! if (mplus(il))
2700
2701	IF (mp(il,i+1)>1.0E-16) THEN
2702	IF (cvflag_grav) THEN
2703	rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
2704	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2705	ELSE
2706	rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
2707	(evap(il,i+1)+evap(il,i))/mp(il, i+1)
2708	END IF
2709	up(il, i) = up(il, i+1)
2710	vp(il, i) = vp(il, i+1)
2711	END IF ! (mp(il,i+1).gt.1.0e-16)
2712	END IF ! (mplus(il)) else if (.not.mplus(il))
2713
2714	rp(il, i) = amin1(rp(il,i), rs(il,i))
2715	rp(il, i) = max(rp(il,i), 0.0)
2716
2717	END IF ! (i.lt.inb(il) .and. lwork(il))
2718	END DO
2719	! ----------------------------------------------------------------
2720
2721	! * find tracer concentrations in precipitating downdraft *
2722
2723	!AC! do j=1,ntra
2724	!AC! do il = 1,ncum
2725	!AC! if (i.lt.inb(il) .and. lwork(il)) then
2726	!AC!c
2727	!AC! if(mplus(il))then
2728	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2729	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2730	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2731	!AC! else ! if (mplus(il))
2732	!AC! if(mp(il,i+1).gt.1.0e-16)then
2733	!AC! trap(il,i,j)=trap(il,i+1,j)
2734	!AC! endif
2735	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
2736	!AC!c
2737	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
2738	!AC! enddo
2739	!AC! end do
2740
2741	400 END DO
2742	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2743
2744	! * end of downdraft loop *
2745
2746	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2747
2748
2749	RETURN
2750	END SUBROUTINE cv3_unsat
2751
2752	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
2753	icb, inb, delt, &
2754	t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
2755	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
2756	ep, clw, m, tp, mp, rp, up, vp, trap, &
2757	wt, water, ice, evap, fondue, faci, b, sigd, &
2758	ment, qent, hent, iflag_mix, uent, vent, &
2759	nent, elij, traent, sig, &
2760	tv, tvp, wghti, &
2761	iflag, precip, Vprecip, ft, fr, fu, fv, ftra, &
2762	cbmf, upwd, dnwd, dnwd0, ma, mip, &
2763	tls, tps, qcondc, wd, &
2764	ftd, fqd)
2765
2766	IMPLICIT NONE
2767
2768	include "cvthermo.h"
2769	include "cv3param.h"
2770	include "cvflag.h"
2771	include "conema3.h"
2772
2773	!inputs:
2774	INTEGER iflag_mix
2775	INTEGER ncum, nd, na, ntra, nloc
2776	LOGICAL ok_conserv_q
2777	INTEGER icb(nloc), inb(nloc)
2778	REAL delt
2779	REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
2780	REAL t_wake(nloc, nd), rr_wake(nloc, nd)
2781	REAL s_wake(nloc)
2782	REAL tra(nloc, nd, ntra), sig(nloc, nd)
2783	REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
2784	REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
2785	REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
2786	REAL lf(nloc, na)
2787	REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
2788	REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
2789	REAL water(nloc, na), evap(nloc, na), b(nloc, na), sigd(nloc)
2790	REAL fondue(nloc, na), faci(nloc, na), ice(nloc, na)
2791	REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
2792	REAL hent(nloc, na, na)
2793	!IM bug real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2794	REAL vent(nloc, na, na), elij(nloc, na, na)
2795	INTEGER nent(nloc, nd)
2796	REAL traent(nloc, na, na, ntra)
2797	REAL tv(nloc, nd), tvp(nloc, nd), wghti(nloc, nd)
2798	!
2799	!input/output:
2800	INTEGER iflag(nloc)
2801	!
2802	!outputs:
2803	REAL precip(nloc)
2804	REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
2805	REAL ftd(nloc, nd), fqd(nloc, nd)
2806	REAL ftra(nloc, nd, ntra)
2807	REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
2808	REAL dnwd0(nloc, nd), mip(nloc, nd)
2809	REAL Vprecip(nloc, nd+1)
2810	REAL tls(nloc, nd), tps(nloc, nd)
2811	REAL qcondc(nloc, nd) ! cld
2812	REAL wd(nloc) ! gust
2813	REAL cbmf(nloc)
2814	!
2815	!local variables:
2816	INTEGER i, k, il, n, j, num1
2817	REAL rat, delti
2818	REAL ax, bx, cx, dx, ex
2819	REAL cpinv, rdcp, dpinv
2820	REAL awat(nloc)
2821	REAL lvcp(nloc, na), lfcp(nloc, na), mke(nloc, na)
2822	REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
2823	!! real up1(nloc), dn1(nloc)
2824	REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
2825	REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2826	REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
2827	REAL th_wake(nloc, nd)
2828	REAL alpha_qpos(nloc), alpha_qpos1(nloc)
2829	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
2830	REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
2831
2832	REAL sumdq !jyg
2833	!
2834	! -------------------------------------------------------------
2835
2836	! initialization:
2837
2838	delti = 1.0/delt
2839	! print*,'cv3_yield initialisation delt', delt
2840	!
2841	DO il = 1, ncum
2842	precip(il) = 0.0
2843	Vprecip(il, nd+1) = 0.0
2844	wd(il) = 0.0 ! gust
2845	END DO
2846
2847	DO i = 1, nd
2848	DO il = 1, ncum
2849	Vprecip(il, i) = 0.0
2850	ft(il, i) = 0.0
2851	fr(il, i) = 0.0
2852	fu(il, i) = 0.0
2853	fv(il, i) = 0.0
2854	upwd(il, i) = 0.0
2855	dnwd(il, i) = 0.0
2856	dnwd0(il, i) = 0.0
2857	mip(il, i) = 0.0
2858	ftd(il, i) = 0.0
2859	fqd(il, i) = 0.0
2860	qcondc(il, i) = 0.0 ! cld
2861	qcond(il, i) = 0.0 ! cld
2862	nqcond(il, i) = 0.0 ! cld
2863	END DO
2864	END DO
2865	! print*,'cv3_yield initialisation 2'
2866	!AC! do j=1,ntra
2867	!AC! do i=1,nd
2868	!AC! do il=1,ncum
2869	!AC! ftra(il,i,j)=0.0
2870	!AC! enddo
2871	!AC! enddo
2872	!AC! enddo
2873	! print*,'cv3_yield initialisation 3'
2874	DO i = 1, nl
2875	DO il = 1, ncum
2876	lvcp(il, i) = lv(il, i)/cpn(il, i)
2877	lfcp(il, i) = lf(il, i)/cpn(il, i)
2878	END DO
2879	END DO
2880
2881
2882
2883	! * calculate surface precipitation in mm/day *
2884
2885	DO il = 1, ncum
2886	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
2887	IF (cvflag_ice) THEN
2888	precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
2889	86400.1000./(rowl*grav)
2890	ELSE
2891	precip(il) = wt(il, 1)sigd(il)water(il, 1) &
2892	86400.1000./(rowl*grav)
2893	END IF
2894	END IF
2895	END DO
2896	! print*,'cv3_yield apres calcul precip'
2897
2898
2899	! === calculate vertical profile of precipitation in kg/m2/s ===
2900
2901	DO i = 1, nl
2902	DO il = 1, ncum
2903	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
2904	IF (cvflag_ice) THEN
2905	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
2906	ELSE
2907	Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
2908	END IF
2909	END IF
2910	END DO
2911	END DO
2912
2913
2914	! * Calculate downdraft velocity scale *
2915	! * NE PAS UTILISER POUR L'INSTANT *
2916
2917	!! do il=1,ncum
2918	!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
2919	!! /(sigd(il)*p(il,icb(il)))
2920	!! enddo
2921
2922
2923	! * calculate tendencies of lowest level potential temperature *
2924	! * and mixing ratio *
2925
2926	DO il = 1, ncum
2927	work(il) = 1.0/(ph(il,1)-ph(il,2))
2928	cbmf(il) = 0.0
2929	END DO
2930
2931	DO k = 2, nl
2932	DO il = 1, ncum
2933	IF (k>=icb(il)) THEN
2934	cbmf(il) = cbmf(il) + m(il, k)
2935	END IF
2936	END DO
2937	END DO
2938
2939	! print*,'cv3_yield avant ft'
2940	! am is the part of cbmf taken from the first level
2941	DO il = 1, ncum
2942	am(il) = cbmf(il)*wghti(il, 1)
2943	END DO
2944
2945	DO il = 1, ncum
2946	IF (iflag(il)<=1) THEN
2947	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
2948	!JYG Correction pour conserver l'eau
2949	! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
2950	IF (cvflag_ice) THEN
2951	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
2952	lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
2953	lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
2954	(100.*(ph(il,1)-ph(il,2))) !precip
2955	ELSE
2956	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
2957	END IF
2958
2959	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
2960
2961	IF (cvflag_ice) THEN
2962	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
2963	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
2964	0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
2965	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
2966	ELSE
2967	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
2968	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
2969	END IF
2970
2971	ftd(il, 1) = ft(il, 1) ! fin precip
2972
2973	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
2974	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
2975	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
2976	END IF ! iflag
2977	END DO
2978
2979
2980	DO j = 2, nl
2981	IF (iflag_mix>0) THEN
2982	DO il = 1, ncum
2983	! FH WARNING a modifier :
2984	cpinv = 0.
2985	! cpinv=1.0/cpn(il,1)
2986	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
2987	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
2988	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
2989	END IF ! j
2990	END DO
2991	END IF
2992	END DO
2993	! fin sature
2994
2995
2996	DO il = 1, ncum
2997	IF (iflag(il)<=1) THEN
2998	!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
2999	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
3000	sigd(il)*evap(il, 1)
3001	!!! sigd(il)0.5(evap(il,1)+evap(il,2))
3002
3003	fqd(il, 1) = fr(il, 1) !precip
3004
3005	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
3006
3007	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
3008	am(il)*(u(il,2)-u(il,1)))
3009	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
3010	am(il)*(v(il,2)-v(il,1)))
3011	END IF ! iflag
3012	END DO ! il
3013
3014
3015	!AC! do j=1,ntra
3016	!AC! do il=1,ncum
3017	!AC! if (iflag(il) .le. 1) then
3018	!AC! if (cvflag_grav) then
3019	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
3020	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3021	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3022	!AC! else
3023	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
3024	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3025	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3026	!AC! endif
3027	!AC! endif ! iflag
3028	!AC! enddo
3029	!AC! enddo
3030
3031	DO j = 2, nl
3032	DO il = 1, ncum
3033	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3034	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3035	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
3036	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
3037	END IF ! j
3038	END DO
3039	END DO
3040
3041	!AC! do k=1,ntra
3042	!AC! do j=2,nl
3043	!AC! do il=1,ncum
3044	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
3045	!AC!
3046	!AC! if (cvflag_grav) then
3047	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
3048	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3049	!AC! else
3050	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
3051	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3052	!AC! endif
3053	!AC!
3054	!AC! endif
3055	!AC! enddo
3056	!AC! enddo
3057	!AC! enddo
3058	! print*,'cv3_yield apres ft'
3059
3060	! * calculate tendencies of potential temperature and mixing ratio *
3061	! * at levels above the lowest level *
3062
3063	! * first find the net saturated updraft and downdraft mass fluxes *
3064	! * through each level *
3065
3066
3067	DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
3068
3069	num1 = 0
3070	DO il = 1, ncum
3071	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
3072	END DO
3073	IF (num1<=0) GO TO 500
3074
3075	CALL zilch(amp1, ncum)
3076	CALL zilch(ad, ncum)
3077
3078	DO k = 1, nl + 1
3079	DO il = 1, ncum
3080	IF (i>=icb(il)) THEN
3081	IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
3082	amp1(il) = amp1(il) + m(il, k)
3083	END IF
3084	ELSE
3085	! AMP1 is the part of cbmf taken from layers I and lower
3086	IF (k<=i) THEN
3087	amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
3088	END IF
3089	END IF
3090	END DO
3091	END DO
3092
3093	DO k = 1, i
3094	DO j = i + 1, nl + 1
3095	DO il = 1, ncum
3096	IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
3097	amp1(il) = amp1(il) + ment(il, k, j)
3098	END IF
3099	END DO
3100	END DO
3101	END DO
3102
3103	DO k = 1, i - 1
3104	DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
3105	DO il = 1, ncum
3106	IF (i<=inb(il) .AND. j<=inb(il)) THEN
3107	ad(il) = ad(il) + ment(il, j, k)
3108	END IF
3109	END DO
3110	END DO
3111	END DO
3112
3113	DO il = 1, ncum
3114	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3115	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3116	cpinv = 1.0/cpn(il, i)
3117
3118	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
3119	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
3120
3121	! precip
3122	! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
3123	IF (cvflag_ice) THEN
3124	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
3125	sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
3126	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
3127	ELSE
3128	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
3129	END IF
3130
3131	rat = cpn(il, i-1)*cpinv
3132
3133	ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
3134	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
3135	IF (cvflag_ice) THEN
3136	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
3137	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
3138	0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
3139	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
3140	ELSE
3141	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
3142	(t_wake(il,i+1)-t_wake(il,i))dpinv &
3143	cpinv
3144	END IF
3145
3146	ftd(il, i) = ft(il, i)
3147	! fin precip
3148
3149	! sature
3150	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
3151	(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
3152	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
3153
3154
3155	IF (iflag_mix==0) THEN
3156	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
3157	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
3158	END IF
3159
3160
3161
3162	! sb: on ne fait pas encore la correction permettant de mieux
3163	! conserver l'eau:
3164	!JYG: correction permettant de mieux conserver l'eau:
3165	! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
3166	fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
3167	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
3168	fqd(il, i) = fr(il, i) ! precip
3169
3170	fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
3171	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
3172	fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
3173	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
3174
3175
3176	fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
3177	ad(il)*(rr(il,i)-rr(il,i-1)))
3178	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
3179	ad(il)*(u(il,i)-u(il,i-1)))
3180	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
3181	ad(il)*(v(il,i)-v(il,i-1)))
3182
3183	END IF ! i
3184	END DO
3185
3186	!AC! do k=1,ntra
3187	!AC! do il=1,ncum
3188	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3189	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3190	!AC! cpinv=1.0/cpn(il,i)
3191	!AC! if (cvflag_grav) then
3192	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
3193	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
3194	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
3195	!AC! else
3196	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
3197	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
3198	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
3199	!AC! endif
3200	!AC! endif
3201	!AC! enddo
3202	!AC! enddo
3203
3204	DO k = 1, i - 1
3205
3206	DO il = 1, ncum
3207	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
3208	awat(il) = max(awat(il), 0.0)
3209	END DO
3210
3211	IF (iflag_mix/=0) THEN
3212	DO il = 1, ncum
3213	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3214	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3215	cpinv = 1.0/cpn(il, i)
3216	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
3217	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
3218	!
3219	!
3220	END IF ! i
3221	END DO
3222	END IF
3223
3224	DO il = 1, ncum
3225	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3226	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3227	cpinv = 1.0/cpn(il, i)
3228	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
3229	(qent(il,k,i)-awat(il)-rr(il,i))
3230	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
3231	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
3232
3233	! (saturated updrafts resulting from mixing) ! cld
3234	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
3235	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3236	END IF ! i
3237	END DO
3238	END DO
3239
3240	!AC! do j=1,ntra
3241	!AC! do k=1,i-1
3242	!AC! do il=1,ncum
3243	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3244	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3245	!AC! cpinv=1.0/cpn(il,i)
3246	!AC! if (cvflag_grav) then
3247	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3248	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3249	!AC! else
3250	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3251	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3252	!AC! endif
3253	!AC! endif
3254	!AC! enddo
3255	!AC! enddo
3256	!AC! enddo
3257
3258	DO k = i, nl + 1
3259
3260	IF (iflag_mix/=0) THEN
3261	DO il = 1, ncum
3262	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
3263	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3264	cpinv = 1.0/cpn(il, i)
3265	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
3266	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
3267
3268
3269	END IF ! i
3270	END DO
3271	END IF
3272
3273	DO il = 1, ncum
3274	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
3275	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3276	cpinv = 1.0/cpn(il, i)
3277
3278	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
3279	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
3280	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
3281	END IF ! i and k
3282	END DO
3283	END DO
3284
3285	!AC! do j=1,ntra
3286	!AC! do k=i,nl+1
3287	!AC! do il=1,ncum
3288	!AC! if (i.le.inb(il) .and. k.le.inb(il)
3289	!AC! $ .and. iflag(il) .le. 1) then
3290	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3291	!AC! cpinv=1.0/cpn(il,i)
3292	!AC! if (cvflag_grav) then
3293	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3294	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3295	!AC! else
3296	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3297	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3298	!AC! endif
3299	!AC! endif ! i and k
3300	!AC! enddo
3301	!AC! enddo
3302	!AC! enddo
3303
3304	! sb: interface with the cloud parameterization: ! cld
3305
3306	DO k = i + 1, nl
3307	DO il = 1, ncum
3308	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
3309	! (saturated downdrafts resulting from mixing) ! cld
3310	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
3311	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3312	END IF ! cld
3313	END DO ! cld
3314	END DO ! cld
3315
3316	! (particular case: no detraining level is found) ! cld
3317	DO il = 1, ncum ! cld
3318	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
3319	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
3320	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3321	END IF ! cld
3322	END DO ! cld
3323
3324	DO il = 1, ncum ! cld
3325	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
3326	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
3327	END IF ! cld
3328	END DO
3329
3330	!AC! do j=1,ntra
3331	!AC! do il=1,ncum
3332	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3333	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3334	!AC! cpinv=1.0/cpn(il,i)
3335	!AC!
3336	!AC! if (cvflag_grav) then
3337	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
3338	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3339	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3340	!AC! else
3341	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
3342	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3343	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3344	!AC! endif
3345	!AC! endif ! i
3346	!AC! enddo
3347	!AC! enddo
3348
3349
3350	500 END DO
3351
3352	!JYG<
3353	!Conservation de l'eau
3354	! sumdq = 0.
3355	! DO k = 1, nl
3356	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3357	! END DO
3358	! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3359	!JYG>
3360	! * move the detrainment at level inb down to level inb-1 *
3361	! * in such a way as to preserve the vertically *
3362	! * integrated enthalpy and water tendencies *
3363
3364	! Correction bug le 18-03-09
3365	DO il = 1, ncum
3366	IF (iflag(il)<=1) THEN
3367	ax = 0.01gravment(il, inb(il), inb(il))* &
3368	(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
3369	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
3370	ft(il, inb(il)) = ft(il, inb(il)) - ax
3371	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3372	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
3373
3374	bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
3375	(ph(il,inb(il))-ph(il,inb(il)+1))
3376	fr(il, inb(il)) = fr(il, inb(il)) - bx
3377	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3378	(ph(il,inb(il)-1)-ph(il,inb(il)))
3379
3380	cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
3381	(ph(il,inb(il))-ph(il,inb(il)+1))
3382	fu(il, inb(il)) = fu(il, inb(il)) - cx
3383	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3384	(ph(il,inb(il)-1)-ph(il,inb(il)))
3385
3386	dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
3387	(ph(il,inb(il))-ph(il,inb(il)+1))
3388	fv(il, inb(il)) = fv(il, inb(il)) - dx
3389	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3390	(ph(il,inb(il)-1)-ph(il,inb(il)))
3391	END IF !iflag
3392	END DO
3393
3394	!JYG<
3395	!Conservation de l'eau
3396	! sumdq = 0.
3397	! DO k = 1, nl
3398	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3399	! END DO
3400	! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3401	!JYG>
3402
3403	!AC! do j=1,ntra
3404	!AC! do il=1,ncum
3405	!AC! IF (iflag(il) .le. 1) THEN
3406	!AC! IF (cvflag_grav) then
3407	!AC! ex=0.01gravment(il,inb(il),inb(il))
3408	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3409	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3410	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3411	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3412	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3413	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3414	!AC! else
3415	!AC! ex=0.1*ment(il,inb(il),inb(il))
3416	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3417	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3418	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3419	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3420	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3421	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3422	!AC! ENDIF !cvflag grav
3423	!AC! ENDIF !iflag
3424	!AC! enddo
3425	!AC! enddo
3426
3427
3428	! * homogenize tendencies below cloud base *
3429
3430
3431	DO il = 1, ncum
3432	asum(il) = 0.0
3433	bsum(il) = 0.0
3434	csum(il) = 0.0
3435	dsum(il) = 0.0
3436	esum(il) = 0.0
3437	fsum(il) = 0.0
3438	gsum(il) = 0.0
3439	hsum(il) = 0.0
3440	END DO
3441
3442	!do i=1,nl
3443	!do il=1,ncum
3444	!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
3445	!enddo
3446	!enddo
3447
3448	DO i = 1, nl
3449	DO il = 1, ncum
3450	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
3451	!jyg Saturated part : use T profile
3452	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
3453	!jyg<20140311
3454	!Correction pour conserver l eau
3455	IF (ok_conserv_q) THEN
3456	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
3457	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
3458
3459	ELSE
3460	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
3461	(ph(il,i)-ph(il,i+1))
3462	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
3463	(ph(il,i)-ph(il,i+1))
3464	ENDIF ! (ok_conserv_q)
3465	!jyg>
3466	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
3467	!jyg Unsaturated part : use T_wake profile
3468	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
3469	!jyg<20140311
3470	!Correction pour conserver l eau
3471	IF (ok_conserv_q) THEN
3472	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
3473	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
3474	ELSE
3475	fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
3476	(ph(il,i)-ph(il,i+1))
3477	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
3478	(ph(il,i)-ph(il,i+1))
3479	ENDIF ! (ok_conserv_q)
3480	!jyg>
3481	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
3482	END IF
3483	END DO
3484	END DO
3485
3486	!!!! do 700 i=1,icb(il)-1
3487	DO i = 1, nl
3488	DO il = 1, ncum
3489	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
3490	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
3491	fqd(il, i) = fsum(il)/gsum(il)
3492	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
3493	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
3494	END IF
3495	END DO
3496	END DO
3497
3498	!jyg<
3499	!Conservation de l'eau
3500	!! sumdq = 0.
3501	!! DO k = 1, nl
3502	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3503	!! END DO
3504	!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3505	!jyg>
3506
3507
3508	! *** Check that moisture stays positive. If not, scale tendencies
3509	! in order to ensure moisture positivity
3510	DO il = 1, ncum
3511	alpha_qpos(il) = 1.
3512	IF (iflag(il)<=1) THEN
3513	IF (fr(il,1)<=0.) THEN
3514	alpha_qpos(il) = max(alpha_qpos(il), (-deltfr(il,1))/(s_wake(il)rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
3515	END IF
3516	END IF
3517	END DO
3518	DO i = 2, nl
3519	DO il = 1, ncum
3520	IF (iflag(il)<=1) THEN
3521	IF (fr(il,i)<=0.) THEN
3522	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
3523	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
3524	END IF
3525	END IF
3526	END DO
3527	END DO
3528	DO il = 1, ncum
3529	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
3530	alpha_qpos(il) = alpha_qpos(il)*1.1
3531	END IF
3532	END DO
3533	DO il = 1, ncum
3534	IF (iflag(il)<=1) THEN
3535	sigd(il) = sigd(il)/alpha_qpos(il)
3536	precip(il) = precip(il)/alpha_qpos(il)
3537	END IF
3538	END DO
3539	DO i = 1, nl
3540	DO il = 1, ncum
3541	IF (iflag(il)<=1) THEN
3542	fr(il, i) = fr(il, i)/alpha_qpos(il)
3543	ft(il, i) = ft(il, i)/alpha_qpos(il)
3544	fqd(il, i) = fqd(il, i)/alpha_qpos(il)
3545	ftd(il, i) = ftd(il, i)/alpha_qpos(il)
3546	fu(il, i) = fu(il, i)/alpha_qpos(il)
3547	fv(il, i) = fv(il, i)/alpha_qpos(il)
3548	m(il, i) = m(il, i)/alpha_qpos(il)
3549	mp(il, i) = mp(il, i)/alpha_qpos(il)
3550	Vprecip(il, i) = vprecip(il, i)/alpha_qpos(il)
3551	END IF
3552	END DO
3553	END DO
3554	DO i = 1, nl
3555	DO j = 1, nl
3556	DO il = 1, ncum
3557	IF (iflag(il)<=1) THEN
3558	ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
3559	END IF
3560	END DO
3561	END DO
3562	END DO
3563
3564	!AC! DO j = 1,ntra
3565	!AC! DO i = 1,nl
3566	!AC! DO il = 1,ncum
3567	!AC! IF (iflag(il) .le. 1) THEN
3568	!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
3569	!AC! ENDIF
3570	!AC! ENDDO
3571	!AC! ENDDO
3572	!AC! ENDDO
3573
3574
3575	! * reset counter and return *
3576
3577	DO il = 1, ncum
3578	sig(il, nd) = 2.0
3579	END DO
3580
3581
3582	DO i = 1, nd
3583	DO il = 1, ncum
3584	upwd(il, i) = 0.0
3585	dnwd(il, i) = 0.0
3586	END DO
3587	END DO
3588
3589	DO i = 1, nl
3590	DO il = 1, ncum
3591	dnwd0(il, i) = -mp(il, i)
3592	END DO
3593	END DO
3594	DO i = nl + 1, nd
3595	DO il = 1, ncum
3596	dnwd0(il, i) = 0.
3597	END DO
3598	END DO
3599
3600
3601	DO i = 1, nl
3602	DO il = 1, ncum
3603	IF (i>=icb(il) .AND. i<=inb(il)) THEN
3604	upwd(il, i) = 0.0
3605	dnwd(il, i) = 0.0
3606	END IF
3607	END DO
3608	END DO
3609
3610	DO i = 1, nl
3611	DO k = 1, nl
3612	DO il = 1, ncum
3613	up1(il, k, i) = 0.0
3614	dn1(il, k, i) = 0.0
3615	END DO
3616	END DO
3617	END DO
3618
3619	DO i = 1, nl
3620	DO k = i, nl
3621	DO n = 1, i - 1
3622	DO il = 1, ncum
3623	IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
3624	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
3625	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
3626	END IF
3627	END DO
3628	END DO
3629	END DO
3630	END DO
3631
3632	DO i = 1, nl
3633	DO k = 1, nl
3634	DO il = 1, ncum
3635	IF (i>=icb(il)) THEN
3636	IF (k>=i .AND. k<=(inb(il))) THEN
3637	upwd(il, i) = upwd(il, i) + m(il, k)
3638	END IF
3639	ELSE
3640	IF (k<i) THEN
3641	upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
3642	END IF
3643	END IF
3644	! c print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
3645	END DO
3646	END DO
3647	END DO
3648
3649	DO i = 2, nl
3650	DO k = i, nl
3651	DO il = 1, ncum
3652	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
3653	IF (i<=inb(il) .AND. k<=inb(il)) THEN
3654	upwd(il, i) = upwd(il, i) + up1(il, k, i)
3655	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
3656	END IF
3657	! c print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
3658	END DO
3659	END DO
3660	END DO
3661
3662
3663	!!!! DO il=1,ncum
3664	!!!! do i=icb(il),inb(il)
3665	!!!!
3666	!!!! upwd(il,i)=0.0
3667	!!!! dnwd(il,i)=0.0
3668	!!!! do k=i,inb(il)
3669	!!!! up1=0.0
3670	!!!! dn1=0.0
3671	!!!! do n=1,i-1
3672	!!!! up1=up1+ment(il,n,k)
3673	!!!! dn1=dn1-ment(il,k,n)
3674	!!!! enddo
3675	!!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
3676	!!!! dnwd(il,i)=dnwd(il,i)+dn1
3677	!!!! enddo
3678	!!!! enddo
3679	!!!!
3680	!!!! ENDDO
3681
3682	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3683	! determination de la variation de flux ascendant entre
3684	! deux niveau non dilue mip
3685	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3686
3687	DO i = 1, nl
3688	DO il = 1, ncum
3689	mip(il, i) = m(il, i)
3690	END DO
3691	END DO
3692
3693	DO i = nl + 1, nd
3694	DO il = 1, ncum
3695	mip(il, i) = 0.
3696	END DO
3697	END DO
3698
3699	DO i = 1, nd
3700	DO il = 1, ncum
3701	ma(il, i) = 0
3702	END DO
3703	END DO
3704
3705	DO i = 1, nl
3706	DO j = i, nl
3707	DO il = 1, ncum
3708	ma(il, i) = ma(il, i) + m(il, j)
3709	END DO
3710	END DO
3711	END DO
3712
3713	DO i = nl + 1, nd
3714	DO il = 1, ncum
3715	ma(il, i) = 0.
3716	END DO
3717	END DO
3718
3719	DO i = 1, nl
3720	DO il = 1, ncum
3721	IF (i<=(icb(il)-1)) THEN
3722	ma(il, i) = 0
3723	END IF
3724	END DO
3725	END DO
3726
3727	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3728	! icb represente de niveau ou se trouve la
3729	! base du nuage , et inb le top du nuage
3730	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3731
3732	DO i = 1, nd
3733	DO il = 1, ncum
3734	mke(il, i) = upwd(il, i) + dnwd(il, i)
3735	END DO
3736	END DO
3737
3738	DO i = 1, nd
3739	DO il = 1, ncum
3740	rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd(1.-rr(il,i))+rr(il,i)cpv)
3741	tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp
3742	tps(il, i) = tp(il, i)
3743	END DO
3744	END DO
3745
3746
3747	! * diagnose the in-cloud mixing ratio * ! cld
3748	! * of condensed water * ! cld
3749	!! cld
3750
3751	DO i = 1, nd ! cld
3752	DO il = 1, ncum ! cld
3753	mac(il, i) = 0.0 ! cld
3754	wa(il, i) = 0.0 ! cld
3755	siga(il, i) = 0.0 ! cld
3756	sax(il, i) = 0.0 ! cld
3757	END DO ! cld
3758	END DO ! cld
3759
3760	DO i = minorig, nl ! cld
3761	DO k = i + 1, nl + 1 ! cld
3762	DO il = 1, ncum ! cld
3763	IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
3764	mac(il, i) = mac(il, i) + m(il, k) ! cld
3765	END IF ! cld
3766	END DO ! cld
3767	END DO ! cld
3768	END DO ! cld
3769
3770	DO i = 1, nl ! cld
3771	DO j = 1, i ! cld
3772	DO il = 1, ncum ! cld
3773	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
3774	.AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
3775	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
3776	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
3777	END IF ! cld
3778	END DO ! cld
3779	END DO ! cld
3780	END DO ! cld
3781
3782	DO i = 1, nl ! cld
3783	DO il = 1, ncum ! cld
3784	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
3785	.AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
3786	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
3787	END IF ! cld
3788	END DO ! cld
3789	END DO ! cld
3790
3791	DO i = 1, nl ! cld
3792	DO il = 1, ncum ! cld
3793	IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
3794	siga(il, i) = mac(il, i)/wa(il, i) & ! cld
3795	rrdtvp(il, i)/p(il, i)/100./delta ! cld
3796	siga(il, i) = min(siga(il,i), 1.0) ! cld
3797	! IM cf. FH
3798	IF (iflag_clw==0) THEN ! cld
3799	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
3800	+(1.-siga(il,i))*qcond(il, i) ! cld
3801	ELSE IF (iflag_clw==1) THEN ! cld
3802	qcondc(il, i) = qcond(il, i) ! cld
3803	END IF ! cld
3804
3805	END DO ! cld
3806	END DO
3807	! print*,'cv3_yield fin'
3808
3809	RETURN
3810	END SUBROUTINE cv3_yield
3811
3812	!AC! et !RomP >>>
3813	SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
3814	ment, sigij, da, phi, phi2, d1a, dam, &
3815	ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
3816	icb, inb)
3817	IMPLICIT NONE
3818
3819	include "cv3param.h"
3820
3821	!inputs:
3822	INTEGER ncum, nd, na, nloc, len
3823	REAL ment(nloc, na, na), sigij(nloc, na, na)
3824	REAL clw(nloc, nd), elij(nloc, na, na)
3825	REAL ep(nloc, na)
3826	INTEGER icb(nloc), inb(nloc)
3827	REAL Vprecip(nloc, nd+1)
3828	!ouputs:
3829	REAL da(nloc, na), phi(nloc, na, na)
3830	REAL phi2(nloc, na, na)
3831	REAL d1a(nloc, na), dam(nloc, na)
3832	REAL epmlmMm(nloc, na, na), eplaMm(nloc, na)
3833	! variables pour tracer dans precip de l'AA et des mel
3834	!local variables:
3835	INTEGER i, j, k
3836	REAL epm(nloc, na, na)
3837
3838	! variables d'Emanuel : du second indice au troisieme
3839	! ---> tab(i,k,j) -> de l origine k a l arrivee j
3840	! ment, sigij, elij
3841	! variables personnelles : du troisieme au second indice
3842	! ---> tab(i,j,k) -> de k a j
3843	! phi, phi2
3844
3845	! initialisations
3846
3847	da(:, :) = 0.
3848	d1a(:, :) = 0.
3849	dam(:, :) = 0.
3850	epm(:, :, :) = 0.
3851	eplaMm(:, :) = 0.
3852	epmlmMm(:, :, :) = 0.
3853	phi(:, :, :) = 0.
3854	phi2(:, :, :) = 0.
3855
3856	! fraction deau condensee dans les melanges convertie en precip : epm
3857	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3858	DO j = 1, na
3859	DO k = 1, na
3860	DO i = 1, ncum
3861	IF (k>=icb(i) .AND. k<=inb(i) .AND. &
3862	!!jyg j.ge.k.and.j.le.inb(i)) then
3863	!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3864	j>k .AND. j<=inb(i)) THEN
3865	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
3866	!!
3867	epm(i, j, k) = max(epm(i,j,k), 0.0)
3868	END IF
3869	END DO
3870	END DO
3871	END DO
3872
3873
3874	DO j = 1, na
3875	DO k = 1, na
3876	DO i = 1, ncum
3877	IF (k>=icb(i) .AND. k<=inb(i)) THEN
3878	eplaMm(i, j) = eplamm(i, j) + &
3879	ep(i, j)clw(i, j)ment(i, j, k)*(1.-sigij(i,j,k))
3880	END IF
3881	END DO
3882	END DO
3883	END DO
3884
3885	DO j = 1, na
3886	DO k = 1, j - 1
3887	DO i = 1, ncum
3888	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3889	epmlmMm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
3890	END IF
3891	END DO
3892	END DO
3893	END DO
3894
3895	! matrices pour calculer la tendance des concentrations dans cvltr.F90
3896	DO j = 1, na
3897	DO k = 1, na
3898	DO i = 1, ncum
3899	da(i, j) = da(i, j) + (1.-sigij(i,k,j))*ment(i, k, j)
3900	phi(i, j, k) = sigij(i, k, j)*ment(i, k, j)
3901	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sigij(i,k,j))
3902	IF (k<=j) THEN
3903	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, k, j)(1.-ep(i,k))*(1.-sigij(i,k,j))
3904	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
3905	END IF
3906	END DO
3907	END DO
3908	END DO
3909
3910	RETURN
3911	END SUBROUTINE cv3_tracer
3912	!AC! et !RomP <<<
3913
3914	SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
3915	iflag, &
3916	precip, sig, w0, &
3917	ft, fq, fu, fv, ftra, &
3918	Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
3919	iflag1, &
3920	precip1, sig1, w01, &
3921	ft1, fq1, fu1, fv1, ftra1, &
3922	Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1)
3923	IMPLICIT NONE
3924
3925	include "cv3param.h"
3926
3927	!inputs:
3928	INTEGER len, ncum, nd, ntra, nloc
3929	INTEGER idcum(nloc)
3930	INTEGER iflag(nloc)
3931	REAL precip(nloc)
3932	REAL sig(nloc, nd), w0(nloc, nd)
3933	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
3934	REAL ftra(nloc, nd, ntra)
3935	REAL ma(nloc, nd)
3936	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
3937	REAL qcondc(nloc, nd)
3938	REAL wd(nloc), cape(nloc)
3939
3940	!outputs:
3941	INTEGER iflag1(len)
3942	REAL precip1(len)
3943	REAL sig1(len, nd), w01(len, nd)
3944	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
3945	REAL ftra1(len, nd, ntra)
3946	REAL ma1(len, nd)
3947	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
3948	REAL qcondc1(nloc, nd)
3949	REAL wd1(nloc), cape1(nloc)
3950
3951	!local variables:
3952	INTEGER i, k, j
3953
3954	DO i = 1, ncum
3955	precip1(idcum(i)) = precip(i)
3956	iflag1(idcum(i)) = iflag(i)
3957	wd1(idcum(i)) = wd(i)
3958	cape1(idcum(i)) = cape(i)
3959	END DO
3960
3961	DO k = 1, nl
3962	DO i = 1, ncum
3963	sig1(idcum(i), k) = sig(i, k)
3964	w01(idcum(i), k) = w0(i, k)
3965	ft1(idcum(i), k) = ft(i, k)
3966	fq1(idcum(i), k) = fq(i, k)
3967	fu1(idcum(i), k) = fu(i, k)
3968	fv1(idcum(i), k) = fv(i, k)
3969	ma1(idcum(i), k) = ma(i, k)
3970	upwd1(idcum(i), k) = upwd(i, k)
3971	dnwd1(idcum(i), k) = dnwd(i, k)
3972	dnwd01(idcum(i), k) = dnwd0(i, k)
3973	qcondc1(idcum(i), k) = qcondc(i, k)
3974	END DO
3975	END DO
3976
3977	DO i = 1, ncum
3978	sig1(idcum(i), nd) = sig(i, nd)
3979	END DO
3980
3981
3982	!AC! do 2100 j=1,ntra
3983	!AC!c oct3 do 2110 k=1,nl
3984	!AC! do 2110 k=1,nd ! oct3
3985	!AC! do 2120 i=1,ncum
3986	!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
3987	!AC! 2120 continue
3988	!AC! 2110 continue
3989	!AC! 2100 continue
3990	!
3991	RETURN
3992	END SUBROUTINE cv3_uncompress

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

Download in other formats:

Original Format