Context Navigation

cv3_routines.F90 @ 2852

Last change on this file since 2852 was 2818, checked in by Laurent Fairhead, 7 years ago
Bug corrected to get same restarts when using ok_conserv_q=y and different parallel domain sectionning
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: 149.7 KB

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

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

Download in other formats:

Original Format