Context Navigation

cv3_routines.F90 @ 5105

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

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

Download in other formats:

Original Format