Context Navigation

cv3_routines.F90 @ 2500

Last change on this file since 2500 was 2490, checked in by jyg, 8 years ago
Bug fix in cv3_unsat (in cv3_routines.F90) in order to guarantee water conservation.
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: 140.9 KB

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

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

Download in other formats:

Original Format