Context Navigation

cv3_routines.F90 @ 2039

Last change on this file since 2039 was 2039, checked in by fhourdin, 11 years ago

Inclusion d'un appel supplementaire a geopot dans leapfrog dans

les versions dyn3d et dyn3dpar pour garantir la convergence
des 3 versions dyn3d/dyn3dpar/dyn3dmem.
La convergence fonctionne avec un compilateur gfortran 4.6.3 / openmpi
distribué sous ubuntu avec la nouvelle physique (NPv3.2) et le guidage.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Une seule exception : si on guide avec dyn3dpar et openMP.

Correction du guidage dans dyn3dmem.

un print supprime dans cv3_toutines.F90

Rm : On ne devrait perdre la convergence numérique avec les précédentes
versions que pour la nouvelle physique puisque le géopotentiel
n'est utilisé dans la physique que par les theermiques et sisvat.

=====================================================================

Added call to geopot in leapfrog in dyn3d and dyn3dpar in order to

insure numerical convergence with dyn3dmem.
The convergence dyn3d/dyn3dpar/dyn3dmem is satistied with
gfortran 4.6.3 / openmpi, new physics (NPv3.2) and nudging.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Only one exception : with nudging & dynd3dpar & openMP.

Bug fixing for nudging in dyn3dmem.

print supressed in cv3_routines.F90

Rm : Numerical convergence with previous releases should be lost for new
physics only since the geopotential is used only by thermals and
sisvat in the physics.

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

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

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

Download in other formats:

Original Format