Context Navigation

cv3_routines.F90 @ 2415

Last change on this file since 2415 was 2398, checked in by jyg, 9 years ago

Introduction of two new flags (ok_conv_stop
[Def=F], ok_intermittent [Def=F]) and one new
parameter (tau_stop [Def=15000]) in
conv_param.data:

. ok_conv_stop=T => convective mass fluxes are

set to zero if there is no trigger for a time
lapse longer than tau_stop.

. ok_intermittent=T => intermittent convection is

represented; the change concerns the upper
bound on the cloud base mass flux in
cv3p2_closure.F90.

The upper bound on Alp in physiq.F90 should

also be changed: this is still to be done.

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

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

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

Download in other formats:

Original Format