Context Navigation

cv3_routines.F90 @ 5139

Last change on this file since 5139 was 5139, checked in by abarral, 3 months ago
Put nuage.h, flux_arp.h, compbl.h into modules Move unused phylmd/ini_hist* to obsolete
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: 179.8 KB

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

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

Download in other formats:

Original Format