Context Navigation

cv3_routines.F90 @ 2343

Last change on this file since 2343 was 2327, checked in by jyg, 9 years ago
Changes in output levels of some variables and Fixes of units for Ale and Alp variables. Cleaning of variable declarations in cv3_yield.
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: 132.1 KB

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

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

Download in other formats:

Original Format