Context Navigation

cv3_routines.f90 @ 5289

Last change on this file since 5289 was 5285, checked in by abarral, 8 months ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
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: 180.0 KB

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

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

Download in other formats:

Original Format