Context Navigation

cv3p_mixing.F90 @ 5020

Last change on this file since 5020 was 4491, checked in by crisi, 20 months ago
Bug corrections in LMDZiso, especially for water tagging
Property svn:keywords set to `Id`
File size: 50.8 KB

Line
1	SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
2	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
3	unk, vnk, hp, tv, tvp, ep, clw, sig, &
4	Ment, Qent, hent, uent, vent, nent, &
5	Sigij, elij, supmax, Ments, Qents, traent &
6	#ifdef ISO
7	& ,xt,xtta,xtclw,xtent,xtelij &
8	#endif
9	& )
10	! **************************************************************
11	! *
12	! CV3P_MIXING : compute mixed draught properties and, *
13	! within a scaling factor, mixed draught *
14	! mass fluxes. *
15	! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
16	! modified by : *
17	! **************************************************************
18
19	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
20	USE ioipsl_getin_p_mod, ONLY: getin_p
21	USE add_phys_tend_mod, ONLY: fl_cor_ebil
22	#ifdef ISO
23	USE infotrac_phy, ONLY: ntraciso=>ntiso
24	USE isotopes_mod, ONLY: pxtmelt,pxtice
25	USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall
26	#ifdef ISOTRAC
27	USE infotrac_phy, ONLY: niso
28	USE isotrac_mod, ONLY: option_cond,izone_cond, &
29	option_tmin,option_traceurs,seuil_tag_tmin, &
30	index_zone,index_iso
31	USE isotrac_routines_mod, ONLY: iso_recolorise_condensation
32	USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall_trac
33	#endif
34	#ifdef ISOVERIF
35	USE isotopes_verif_mod
36	!, ONLY: iso_verif_positif_nostop, iso_verif_egalite_nostop, &
37	! errmax,errmaxrel,deltalim
38	USE isotopes_mod, ONLY: ridicule,iso_eau,iso_hdo
39	#endif
40	#endif
41	IMPLICIT NONE
42
43	include "cvthermo.h"
44	include "cv3param.h"
45	include "YOMCST2.h"
46	include "cvflag.h"
47
48	!inputs:
49	INTEGER, INTENT (IN) :: ncum, nd, na
50	INTEGER, INTENT (IN) :: ntra, nloc
51	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
52	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
53	REAL, DIMENSION (nloc), INTENT (IN) :: unk, vnk
54	REAL, DIMENSION (nloc, nd), INTENT (IN) :: qta
55	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
56	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
57	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
58	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
59	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv
60	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
61	REAL, DIMENSION (nloc, na), INTENT (IN) :: frac !ice fraction in condensate
62	REAL, DIMENSION (nloc, na), INTENT (IN) :: h !liquid water static energy of environment
63	REAL, DIMENSION (nloc, na), INTENT (IN) :: hp !liquid water static energy of air shed from adiab. asc.
64	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp
65	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, clw
66	#ifdef ISO
67	REAL, DIMENSION (ntraciso,nloc, na), INTENT (IN) :: xtclw
68	REAL, DIMENSION (ntraciso,nloc, nd), INTENT (IN) :: xt
69	REAL, DIMENSION (ntraciso,nloc, nd), INTENT (IN) :: xtta
70	#endif
71
72	!outputs:
73	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Ment, Qent
74	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
75	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Sigij, elij
76	REAL, DIMENSION (nloc, na), INTENT (OUT) :: supmax ! Highest mixing fraction of mixed
77	! updraughts with the sign of (h-hp)
78	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
79	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: Ments, Qents
80	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: hent
81	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
82	#ifdef ISO
83	REAL, DIMENSION (ntraciso,nloc, na, na), INTENT (OUT) :: xtent
84	REAL, DIMENSION (ntraciso,nloc, na, na), INTENT (OUT) :: xtelij
85	#endif
86
87
88	!local variables:
89	INTEGER i, j, k, il, im, jm
90	INTEGER num1, num2
91	REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp
92	REAL :: alt, delp, delm
93	REAL, DIMENSION (nloc) :: Qmixmax, Rmixmax, sqmrmax
94	REAL, DIMENSION (nloc) :: Qmixmin, Rmixmin, sqmrmin
95	REAL, DIMENSION (nloc) :: signhpmh
96	REAL, DIMENSION (nloc) :: Sx
97	REAL :: Scrit2
98	REAL, DIMENSION (nloc) :: Smid, Sjmin, Sjmax
99	REAL, DIMENSION (nloc) :: Sbef, sup, smin
100	REAL, DIMENSION (nloc) :: ASij, ASij_inv, smax, Scrit
101	REAL, DIMENSION (nloc, nd, nd) :: Sij
102	REAL, DIMENSION (nloc, nd) :: csum
103	REAL :: awat
104	REAL :: cpm !Mixed draught heat capacity
105	REAL :: Tm !Mixed draught temperature
106	LOGICAL, DIMENSION (nloc) :: lwork
107
108	REAL amxupcrit, df, ff
109	INTEGER nstep
110
111	INTEGER,SAVE :: igout=1
112	!$OMP THREADPRIVATE(igout)
113
114	! -- Mixing probability distribution functions
115
116	REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
117
118	Qcoef1(F) = tanh(F/gammas)
119	Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
120	QFF(F) = max(min(F,1.), 0.)
121	QFFf(F) = min(QFF(F), scut)
122	Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
123	Rmix1(F) = (gammaslog(cosh((QFF(F)-Fmax)/gammas))+QFF(F)Qcoef1max)/Qcoef2max
124	Qmix2(F) = -log(1.-QFFf(F))/scut
125	Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
126	Qmix(F) = qqa1Qmix1(F) + qqa2Qmix2(F)
127	Rmix(F) = qqa1Rmix1(F) + qqa2Rmix2(F)
128
129	INTEGER, SAVE :: ifrst
130	DATA ifrst/0/
131	!$OMP THREADPRIVATE(ifrst)
132
133	#ifdef ISO
134	integer ixt
135	real xtrti(ntraciso)
136	real xtres(ntraciso)
137	real zfice(nloc),zxtliq(ntraciso,nloc),zxtice(ntraciso,nloc)
138	real qent_diag(nloc,nd,nd)
139	real xtent_diag(ntraciso,nloc,nd,nd)
140	#endif
141
142
143	! =====================================================================
144	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
145	! =====================================================================
146
147	! -- Initialize mixing PDF coefficients
148	IF (ifrst==0) THEN
149	ifrst = 1
150	Qcoef1max = Qcoef1(Fmax)
151	Qcoef2max = Qcoef2(Fmax)
152	!<jyg
153	print*, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
154	fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
155	!>jyg
156	!
157	END IF
158
159
160	! ori do 360 i=1,ncum*nlp
161	DO j = 1, nl
162	DO i = 1, ncum
163	nent(i, j) = 0
164	! in convect3, m is computed in cv3_closure
165	! ori m(i,1)=0.0
166	END DO
167	END DO
168
169	#ifdef ISOVERIF
170	Qent(:,:,:)=0.0
171	xtent(:,:,:,:)=0.0
172	#endif
173
174	! ori do 400 k=1,nlp
175	! ori do 390 j=1,nlp
176	DO j = 1, nl
177	DO k = 1, nl
178	DO i = 1, ncum
179	Qent(i, k, j) = rr(i, j)
180	uent(i, k, j) = u(i, j)
181	vent(i, k, j) = v(i, j)
182	elij(i, k, j) = 0.0
183	hent(i, k, j) = 0.0
184	!AC! Ment(i,k,j)=0.0
185	!AC! Sij(i,k,j)=0.0
186	#ifdef ISO
187	xtelij(:,i, k, j) = 0.0
188	xtent(:,i, k, j) = xt(:,i, j)
189	#endif
190	END DO
191	END DO
192	END DO
193	#ifdef ISO
194	! on initialise mieux par sécurité:
195	elij(:, :, :)=0.0
196	xtelij(:,:, :, :)=0.0
197	#endif
198
199	!AC!
200	Ment(1:ncum, 1:nd, 1:nd) = 0.0
201	Sij(1:ncum, 1:nd, 1:nd) = 0.0
202	!AC!
203	!ym
204	Sigij(1:ncum, 1:nd, 1:nd) = 0.0
205	!ym
206
207	!jyg! DO k = 1, ntra
208	!jyg! DO j = 1, nd ! instead nlp
209	!jyg! DO i = 1, nd ! instead nlp
210	!jyg! DO il = 1, ncum
211	!jyg! traent(il, i, j, k) = tra(il, j, k)
212	!jyg! END DO
213	!jyg! END DO
214	!jyg! END DO
215	!jyg! END DO
216
217	! =====================================================================
218	! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
219	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
220	! --- FRACTION (Sij)
221	! =====================================================================
222
223	DO i = minorig + 1, nl
224
225	IF (ok_entrain) THEN
226	DO j = minorig, nl
227	DO il = 1, ncum
228	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
229	.AND. (j<=inb(il))) THEN
230
231	!! rti = qnk(il) - ep(il, i)*clw(il, i)
232	rti = qta(il,i-1) - ep(il, i)*clw(il, i)
233	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
234	!jyg(from aj)<
235	IF (cvflag_ice) THEN
236	! print*,cvflag_ice,'cvflag_ice dans do 700'
237	IF (t(il,j)<=263.15) THEN
238	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
239	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
240	END IF
241	END IF ! IF (cvflag_ice) THEN
242	!>jyg
243	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
244	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
245	dei = denom
246	IF (abs(dei)<0.01) dei = 0.01
247	Sij(il, i, j) = anum/dei
248	Sij(il, i, i) = 1.0
249	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il,i,j))rti - rs(il, j)
250	altem = altem/bf2
251	cwat = clw(il, j)*(1.-ep(il,j))
252	stemp = Sij(il, i, j)
253	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
254	!jyg(from aj)<
255	IF (cvflag_ice) THEN
256	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
257	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
258	ELSE
259	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
260	denom = denom + lv(il, j)*(rr(il,i)-rti)
261	END IF
262	!>jyg
263	IF (abs(denom)<0.01) denom = 0.01
264	Sij(il, i, j) = anum/denom
265	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il,i,j))rti - rs(il, j)
266	altem = altem - (bf2-1.)*cwat
267	END IF
268	IF (Sij(il,i,j)>0.0) THEN
269	!!! Ment(il,i,j)=m(il,i)
270	Ment(il, i, j) = 1.
271	elij(il, i, j) = altem
272	elij(il, i, j) = amax1(0.0, elij(il,i,j))
273	nent(il, i) = nent(il, i) + 1
274	END IF
275
276	Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
277	Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
278	ELSE IF (j > i) THEN
279	IF (prt_level >= 10) THEN
280	print *,'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il,i,j)
281	ENDIF
282	END IF ! new ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
283	#ifdef ISO
284	if(sij(il,i,j).gt.0.0)then
285	qent_diag(il,i,j)=sij(il,i,j)*rr(il,i) &
286	+(1.-sij(il,i,j))*rti
287	! question: pourquoi qent n'est-il pas recalculé ici?
288	endif !if(sij(il,i,j).gt.0.0)then
289	#endif
290	END DO ! il
291
292	#ifdef ISO
293	do il=1,ncum
294	zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
295	zfice(il) = MIN(MAX(zfice(il),0.0),1.0)
296
297	if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &
298	(j.ge.(icb(il)-1)).and.(j.le.inb(il))) then
299	if (sij(il,i,j).gt.0.0) then
300	do ixt=1,ntraciso
301	xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
302	xtent_diag(ixt,il,i,j)=sij(il,i,j)*xt(ixt,il,i) &
303	+(1.-sij(il,i,j))*xtrti(ixt)
304	enddo ! do ixt=1,ntraciso
305	#ifdef ISOVERIF
306	if (iso_eau.gt.0) then
307	call iso_verif_egalite(xtta(iso_eau,il,i-1), &
308	qta(il,i-1),'cv3p_mixing 256a')
309	call iso_verif_egalite(xt(iso_eau,il,i), &
310	rr(il,i),'cv3p_mixing 256b')
311	call iso_verif_egalite(xtclw(iso_eau,il,i), &
312	clw(il,i),'cv3p_mixing 256c')
313	call iso_verif_egalite(xtent_diag(iso_eau,il,i,j), &
314	qent_diag(il,i,j),'cv3p_mixing 256d')
315	endif ! if (iso_eau.gt.0) then
316	#ifdef ISOTRAC
317	call iso_verif_traceur(xtrti,'cv3p_mixing 1968')
318	#endif
319	#endif
320	if (qent_diag(il,i,j).lt.elij(il,i,j)) then
321	do ixt=1,ntraciso
322	xtent_diag(ixt,il,i,j)= xtent_diag(ixt,il,i,j) &
323	/ qent_diag(il,i,j)*elij(il,i,j)
324	enddo ! do ixt=1,ntraciso
325	qent_diag(il,i,j)=elij(il,i,j)
326	endif !if (qent_diag(il,i,j).lt.elij(il,i,j)) then
327
328	#ifdef ISOVERIFIDRIS
329	call iso_verif_noNaN(qent_diag(il,i,j),'cv3p_routines 261')
330	#endif
331	#ifdef ISOVERIF
332	if (iso_eau.gt.0) then
333	if (iso_verif_egalite_nostop(xtent_diag(iso_eau,il,i,j), &
334	qent_diag(il,i,j),'cv3p_mixing 258').eq.1) then
335	write(,) 'il,i,j=',il,i,j
336	write(,) 'qent_diag(il,i,j),elij(il,i,j)=', &
337	qent_diag(il,i,j),elij(il,i,j)
338	write(,) 'sij(il,i,j)=',sij(il,i,j)
339	write(,) 'xt=',xt(iso_eau,il,i),rr(il,i)
340	write(,) 'xtclw=',xtclw(iso_eau,il,i),clw(il,i)
341	write(,) 'ep(il,i)=',ep(il,i)
342	stop
343	endif
344	endif
345	call iso_verif_positif(qent_diag(il,i,j)-elij(il,i,j), &
346	'cv3p_routines 262')
347	#endif
348	else !if (sij(il,i,j).gt.0.0) then
349	! ajout le 4 mai 2012: initialiser au cas où
350	elij(il,i,j)=0.0
351	qent_diag(il,i,j)=rr(il,i)
352	do ixt=1,ntraciso
353	xtent_diag(ixt,il,i,j)=xt(ixt,il,i)
354	enddo !do ixt=1,ntraciso
355	endif !if(sij(il,i,j).gt.0.0)then
356	else !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
357	! ajout le 4 mai 2012: initialiser au cas où
358	elij(il,i,j)=0.0
359	qent_diag(il,i,j)=rr(il,i)
360	do ixt=1,ntraciso
361	xtent_diag(ixt,il,i,j)=xt(ixt,il,i)
362	enddo !do ixt=1,ntraciso
363	endif !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
364	#ifdef ISOVERIF
365	#ifdef VERIFNEGATIF
366	call iso_verif_positif(qent_diag(il,i,j), &
367	'cv3p_mixing 303: qt<0')
368	#endif
369	call iso_verif_positif(qent_diag(il,i,j)-elij(il,i,j), &
370	'cv3p_mixing 305: qt<elij')
371	#endif
372	enddo !do il=1,ncum
373
374	!write(,) 'cv3p_mixing 265: call condiso'
375	call condiso_liq_ice_vectall(xtent_diag(1,1,i,j), &
376	qent_diag(1,i,j),elij(1,i,j), &
377	t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
378	#ifdef ISOTRAC
379	call condiso_liq_ice_vectall_trac(xtent_diag(1,1,i,j), &
380	qent_diag(1,i,j),elij(1,i,j), &
381	t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
382	#ifdef ISOVERIF
383	do il=1,ncum
384	call iso_verif_traceur(xt(1,il,i),'cv3p_mixing 1967')
385	call iso_verif_traceur(xtent(1,il,i,j),'cv3p_mixing 1969')
386	enddo !do il=1,ncum
387	#endif
388	#endif
389	do il=1,ncum
390	do ixt = 1, ntraciso
391	xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
392	enddo !do ixt = 1, ntraciso
393	enddo !do il=1,ncum
394
395	#ifdef ISOTRAC
396	! write(,) 'cv3p_mixing tmp 1987,option_traceurs=',
397	! : option_traceurs
398	if (option_tmin.ge.1) then
399	do il=1,ncum
400	! colorier la vapeur résiduelle selon température de
401	! condensation, et le condensat en un tag spécifique
402	if ((elij(il,i,j).gt.0.0).and. &
403	(qent_diag(il,i,j).gt.0.0)) then
404	if (option_traceurs.eq.17) then
405	call iso_recolorise_condensation(qent_diag(il,i,j), &
406	elij(il,i,j), &
407	xtent_diag(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
408	0.0,xtres, &
409	seuil_tag_tmin)
410	else !if (option_traceurs.eq.17) then
411	! write(,) 'cv3 2002: il,i,j =',il,i,j
412	call iso_recolorise_condensation(qent_diag(il,i,j), &
413	elij(il,i,j), &
414	xtent_diag(1,il,i,j),xtelij(1,il,i,j),rs(il,j), &
415	0.0,xtres,seuil_tag_tmin)
416	endif !if (option_traceurs.eq.17) then
417	do ixt=1+niso,ntraciso
418	xtent_diag(ixt,il,i,j)=xtres(ixt)
419	enddo
420	endif !if (cond.gt.0.0) then
421	enddo !do il=1,ncum
422	#ifdef ISOVERIF
423	do il=1,ncum
424	call iso_verif_traceur(xtent(1,il,i,j),'cv3p_mixing 1996')
425	call iso_verif_traceur(xtelij(1,il,i,j),'cv3p_mixing 1997')
426	call iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
427	'cv3p_mixing 2042')
428	enddo !do il=1,ncum
429	#endif
430	endif !if (option_tmin.ge.1) then
431	#endif
432
433	! fractionation:
434	#ifdef ISOVERIF
435	do il=1,ncum
436	if ((i.ge.icb(il)).and.(i.le.inb(il)).and. &
437	(j.ge.(icb(il)-1)).and.(j.le.inb(il))) then
438	if (sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.scut) then
439	if (iso_eau.gt.0) then
440	call iso_verif_egalite_choix(xtent(iso_eau,il,i,j), &
441	qent(il,i,j),'cv3p_mixing 1889',errmax,errmaxrel)
442	call iso_verif_egalite_choix(xtelij(iso_eau,il,i,j), &
443	elij(il,i,j),'cv3p_mixing 1890',errmax,errmaxrel)
444	endif
445	if (iso_hdo.gt.0) then
446	call iso_verif_aberrant_choix(xt(iso_HDO,il,i),rr(il,i), &
447	ridicule,deltalim,'cv3p_mixing 1997')
448	call iso_verif_aberrant_choix( &
449	xtent(iso_HDO,il,i,j),qent(il,i,j), &
450	ridicule,deltalim,'cv3p_mixing 1931')
451	call iso_verif_aberrant_choix( &
452	xtelij(iso_HDO,il,i,j),elij(il,i,j), &
453	ridicule,deltalim_snow,'cv3p_mixing 1993')
454	endif !if (iso_hdo.gt.0) then
455	#ifdef ISOTRAC
456	! write(,) 'cv3p_mixing tmp 2039 il=',il
457	call iso_verif_traceur(xtent(1,il,i,j), &
458	'cv3p_mixing 2031')
459	call iso_verif_traceur(xtelij(1,il,i,j), &
460	'cv3p_mixing 2033')
461	#endif
462
463	endif !if(sij(il,i,j).gt.0.0)then
464	endif !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
465	enddo !do il=1,ncum
466	#endif
467	! write(,) 'cv3p_routine tmp 1984: cond=',elij(il,i,j)
468	#endif
469
470	END DO ! j
471	ELSE ! (ok_entrain)
472	DO il = 1,ncum
473	nent(il,i) = 0
474	ENDDO
475	ENDIF ! (ok_entrain)
476
477	!jygdebug<
478	IF (prt_level >= 10) THEN
479	print *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
480	IF (nent(igout,i) .gt. 0) THEN
481	print *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
482	ENDIF
483	ENDIF
484	!>jygdebug
485
486	! * if no air can entrain at level i assume that updraft detrains *
487	! * at that level and calculate detrained air flux and properties *
488
489
490	! @ do 170 i=icb(il),inb(il)
491
492	DO il = 1, ncum
493	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
494	! @ if(nent(il,i).eq.0)then
495	!!! Ment(il,i,i)=m(il,i)
496	Ment(il, i, i) = 1.
497	!! Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
498	Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
499	uent(il, i, i) = unk(il)
500	vent(il, i, i) = vnk(il)
501	IF (fl_cor_ebil .GE. 2) THEN
502	hent(il, i, i) = hp(il,i)
503	ENDIF
504	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
505	Sij(il, i, i) = 0.0
506
507
508	#ifdef ISO
509	do ixt = 1, ntraciso
510	xtent(ixt,il,i,i)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
511	xtelij(ixt,il,i,i)=xtclw(ixt,il,i) ! rq: ne sera pas utilise ensuite
512	enddo ! do ixt = 1, ntraciso
513	#ifdef ISOVERIF
514	if (iso_eau.gt.0) then
515	call iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
516	qent(il,i,i),'cv3p_mixing 414',errmax,errmaxrel)
517	endif
518	#endif
519	#ifdef ISOTRAC
520	if (option_tmin.ge.1) then
521	! colorier la vapeur résiduelle selon température de
522	! condensation, et le condensat en un tag spécifique
523	! write(,) 'cv3 tmp 2095 il,i,j,xtent(:,il,i,j)=',
524	! : il,i,j,xtent(:,il,i,j)
525	if ((elij(il,i,i).gt.0.0).and.(qent(il,i,i).gt.0.0)) then
526	if (option_traceurs.eq.17) then
527	call iso_recolorise_condensation(qent(il,i,i), &
528	elij(il,i,i), &
529	xt(1,il,nk(il)),xtclw(1,il,i),t(il,i),ep(il,i), &
530	xtres, &
531	seuil_tag_tmin)
532	else !if (option_traceurs.eq.17) then
533	call iso_recolorise_condensation(qent(il,i,i), &
534	elij(il,i,i), &
535	xt(1,il,nk(il)),xtclw(1,il,i),rs(il,i),ep(il,i), &
536	xtres, &
537	seuil_tag_tmin)
538	endif !if (option_traceurs.eq.17) then
539	do ixt=1+niso,ntraciso
540	xtent(ixt,il,i,i)=xtres(ixt)
541	enddo
542	#ifdef ISOVERIF
543	do ixt=1,niso
544	call iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), &
545	'cv3p_mixing 2102',errmax,errmaxrel)
546	call iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
547	'cv3p_mixing 2154')
548	enddo
549	#endif
550	endif !if (cond.gt.0.0) then
551
552	#ifdef ISOVERIF
553	call iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
554	qent(il,i,i),'cv3p_mixing 2103',errmax,errmaxrel)
555	call iso_verif_traceur(xtent(1,il,i,i),'cv3p_mixing 2095')
556	call iso_verif_traceur(xtelij(1,il,i,i),'cv3p_mixing 2096')
557	#endif
558	endif !if (option_tmin.ge.1) then
559	#endif
560	! endif ISOTRAC
561	#endif
562	! endif ISO
563
564	END IF ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
565	END DO ! il
566	END DO ! i = minorig + 1, nl
567
568	!jyg! DO j = 1, ntra
569	!jyg! DO i = minorig + 1, nl
570	!jyg! DO il = 1, ncum
571	!jyg! IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
572	!jyg! traent(il, i, i, j) = tra(il, nk(il), j)
573	!jyg! END IF
574	!jyg! END DO
575	!jyg! END DO
576	!jyg! END DO
577
578	DO j = minorig, nl
579	DO i = minorig, nl
580	DO il = 1, ncum
581	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
582	(i>=icb(il)) .AND. (i<=inb(il))) THEN
583	Sigij(il, i, j) = Sij(il, i, j)
584	END IF
585	END DO
586	END DO
587	END DO
588	! @ enddo
589
590	! @170 continue
591
592	! =====================================================================
593	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
594	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
595	! =====================================================================
596
597	CALL zilch(csum, nloc*nd)
598
599	DO il = 1, ncum
600	lwork(il) = .FALSE.
601	END DO
602
603	! ---------------------------------------------------------------
604	DO i = minorig + 1, nl !Loop on origin level "i"
605	! ---------------------------------------------------------------
606
607	num1 = 0
608	DO il = 1, ncum
609	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
610	END DO
611	IF (num1<=0) GO TO 789
612
613
614	!JYG1 Find maximum of SIJ for J>I, if any.
615
616	Sx(:) = 0.
617
618	DO il = 1, ncum
619	IF (i>=icb(il) .AND. i<=inb(il)) THEN
620	signhpmh(il) = sign(1., hp(il,i)-h(il,i))
621	Sbef(il) = max(0., signhpmh(il))
622	END IF
623	END DO
624
625	DO j = i + 1, nl
626	DO il = 1, ncum
627	IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
628	IF (Sbef(il)<Sij(il,i,j)) THEN
629	Sx(il) = max(Sij(il,i,j), Sx(il))
630	END IF
631	Sbef(il) = Sij(il, i, j)
632	END IF
633	END DO
634	END DO
635
636
637	DO il = 1, ncum
638	IF (i>=icb(il) .AND. i<=inb(il)) THEN
639	lwork(il) = (nent(il,i)/=0)
640	!! rti = qnk(il) - ep(il, i)*clw(il, i)
641	rti = qta(il,i-1) - ep(il, i)*clw(il, i)
642	!jyg<
643	IF (cvflag_ice) THEN
644
645	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
646	(rti-rs(il,i)) + (cpv-cpd)t(il, i)(rti-rr(il,i))
647	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
648	(rr(il,i)-rti) + (cpd-cpv)t(il, i)(rr(il,i)-rti)
649	ELSE
650
651	anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
652	(cpv-cpd)t(il, i)(rti-rr(il,i))
653	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
654	(cpd-cpv)t(il, i)(rr(il,i)-rti)
655	END IF
656	!>jyg
657	IF (abs(denom)<0.01) denom = 0.01
658	Scrit(il) = min(anum/denom, 1.)
659	alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
660
661	!JYG1 Find new critical value Scrit2
662	! such that : Sij > Scrit2 => mixed draught will detrain at J<I
663	! Sij < Scrit2 => mixed draught will detrain at J>I
664
665	Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
666	Scrit(il)*max(0., signhpmh(il))
667
668	Scrit(il) = Scrit2
669
670	!JYG Correction pour la nouvelle logique; la correction pour ALT
671	! est un peu au hazard
672	IF (Scrit(il)<=0.0) Scrit(il) = 0.0
673	IF (alt<=0.0) Scrit(il) = 1.0
674
675	smax(il) = 0.0
676	ASij(il) = 0.0
677	sup(il) = 0. ! upper S-value reached by descending draughts
678	END IF
679	END DO
680
681	! ---------------------------------------------------------------
682	DO j = minorig, nl !Loop on destination level "j"
683	! ---------------------------------------------------------------
684
685	num2 = 0
686	DO il = 1, ncum
687	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
688	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
689	lwork(il)) num2 = num2 + 1
690	END DO
691	IF (num2<=0) GO TO 175
692
693	! -----------------------------------------------
694	IF (j>i) THEN
695	! -----------------------------------------------
696	DO il = 1, ncum
697	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
698	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
699	lwork(il)) THEN
700	IF (Sij(il,i,j)>0.0) THEN
701	Smid(il) = min(Sij(il,i,j), Scrit(il))
702	Sjmax(il) = Smid(il)
703	Sjmin(il) = Smid(il)
704	IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
705	smin(il) = Smid(il)
706	Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
707	Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
708	Sjmin(il) = min(Sjmin(il), Scrit(il))
709	Sbef(il) = Sij(il, i, j)
710	END IF
711	END IF
712	END IF
713	END DO
714	! -----------------------------------------------
715	ELSE IF (j==i) THEN
716	! -----------------------------------------------
717	DO il = 1, ncum
718	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
719	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
720	lwork(il)) THEN
721	IF (Sij(il,i,j)>0.0) THEN
722	Smid(il) = 1.
723	Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
724	min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
725	Sjmin(il) = max(Sjmin(il), sup(il))
726	Sjmax(il) = 1.
727
728	! - preparation des variables Scrit, Smin et Sbef pour la partie j>i
729	Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
730
731	smin(il) = 1.
732	Sbef(il) = max(0., signhpmh(il))
733	supmax(il, i) = sign(Scrit(il), -signhpmh(il))
734	END IF
735	END IF
736	END DO
737	! -----------------------------------------------
738	ELSE IF (j<i) THEN
739	! -----------------------------------------------
740	DO il = 1, ncum
741	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
742	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
743	lwork(il)) THEN
744	IF (Sij(il,i,j)>0.0) THEN
745	Smid(il) = max(Sij(il,i,j), Scrit(il))
746	Sjmax(il) = Smid(il)
747	Sjmin(il) = Smid(il)
748	IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
749	smax(il) = Smid(il)
750	Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
751	Sjmax(il) = max(Sjmax(il), Scrit(il))
752	Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
753	Sjmin(il) = max(Sjmin(il), Scrit(il))
754	Sbef(il) = Sij(il, i, j)
755	END IF
756	IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
757	sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
758	END IF
759	END IF
760	END DO
761	! -----------------------------------------------
762	END IF
763	! -----------------------------------------------
764
765
766	DO il = 1, ncum
767	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
768	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
769	lwork(il)) THEN
770	IF (Sij(il,i,j)>0.0) THEN
771	!! rti = qnk(il) - ep(il, i)*clw(il, i)
772	rti = qta(il,i-1) - ep(il, i)*clw(il, i)
773	Qmixmax(il) = Qmix(Sjmax(il))
774	Qmixmin(il) = Qmix(Sjmin(il))
775	Rmixmax(il) = Rmix(Sjmax(il))
776	Rmixmin(il) = Rmix(Sjmin(il))
777	sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
778	sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
779
780	Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
781
782	! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
783	IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
784	Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
785	ELSE
786	Sigij(il, i, j) = 0.
787	END IF
788
789	! -- Compute Qent, uent, vent according to the true mixing fraction
790	Qent(il, i, j) = (1.-Sigij(il,i,j))rti + Sigij(il, i, j)rr(il, i)
791	uent(il, i, j) = (1.-Sigij(il,i,j))unk(il) + Sigij(il, i, j)u(il, i)
792	vent(il, i, j) = (1.-Sigij(il,i,j))vnk(il) + Sigij(il, i, j)v(il, i)
793
794	! -- Compute liquid water static energy of mixed draughts
795	! IF (j .GT. i) THEN
796	! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
797	! awat=amax1(awat,0.0)
798	! ELSE
799	! awat = 0.
800	! ENDIF
801	! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
802	! : + Sigij(il,i,j)*H(il,i)
803	! : + (LV(il,j)+(cpd-cpv)t(il,j))awat
804	!IM 301008 beg
805	hent(il, i, j) = (1.-Sigij(il,i,j))hp(il, i) + Sigij(il, i, j)h(il, i)
806
807	!jyg<
808	! elij(il, i, j) = Qent(il, i, j) - rs(il, j)
809	! elij(il, i, j) = elij(il, i, j) + &
810	! ((h(il,j)-hent(il,i,j))rs(il,j)lv(il,j) / &
811	! ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
812	! elij(il, i, j) = elij(il, i, j) / &
813	! (1.+lv(il,j)lv(il,j)rs(il,j) / &
814	! ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
815	!
816	! Computation of condensate amount Elij, taking into account the ice fraction frac
817	! Warning : the same saturation humidity rs is used over both liquid water and ice; this
818	! should be corrected.
819	!
820	! Heat capacity of mixed draught
821	cpm = cpd+Qent(il,i,j)*(cpv-cpd)
822	!
823	IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
824	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
825	elij(il, i, j) = elij(il, i, j) + &
826	(h(il,j)-hent(il,i,j)+(cpv-cpd)(Qent(il,i,j)-rr(il,j))t(il,j))* &
827	rs(il,j)lv(il,j) / (cpmrrvt(il,j)t(il,j))
828	elij(il, i, j) = elij(il, i, j) / &
829	(1.+(lv(il,j)+frac(il,j)lf(il,j))lv(il,j)*rs(il,j) / &
830	(cpmrrvt(il,j)*t(il,j)))
831	ELSE
832	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
833	elij(il, i, j) = elij(il, i, j) + &
834	(h(il,j)-hent(il,i,j)+(cpv-cpd)(Qent(il,i,j)-rr(il,j))t(il,j))* &
835	rs(il,j)lv(il,j) / (cpmrrvt(il,j)t(il,j))
836	elij(il, i, j) = elij(il, i, j) / &
837	(1.+lv(il,j)lv(il,j)rs(il,j) / &
838	(cpmrrvt(il,j)*t(il,j)))
839	ENDIF
840	!>jyg
841	elij(il, i, j) = max(elij(il,i,j), 0.)
842
843	elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
844
845	IF (j>i) THEN
846	awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
847	awat = amax1(awat, 0.0)
848	ELSE
849	awat = 0.
850	END IF
851
852	! print ,h(il,j)-hent(il,i,j),LV(il,j)rs(il,j)/(cpdrrvt(il,j)*
853	! : t(il,j))
854
855	!jyg<
856	! hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)t(il,j))awat
857	! Mixed draught temperature at level j
858	IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
859	Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))rrvt(il,j)t(il,j)/(lv(il,j)rs(il,j))
860	hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)lf(il,j)+(cpd-cpv)Tm)*awat
861	ELSE
862	Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))rrvt(il,j)t(il,j)/(lv(il,j)rs(il,j))
863	hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)Tm)awat
864	ENDIF
865	!>jyg
866
867	!IM 301008 end
868
869	! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
870	! : i,j,hent(il,i,j),Sigij(il,i,j)
871
872	! -- ASij is the integral of P(F) over the relevant F interval
873	ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
874	Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
875
876	END IF
877	END IF
878	END DO ! il
879
880
881
882	#ifdef ISO
883	do il=1,ncum
884	zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
885	zfice(il) = MIN(MAX(zfice(il),0.0),1.0)
886	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
887	j.ge.(icb(il)-1) .and. j.le.inb(il) &
888	.and. lwork(il) ) then
889	if(sij(il,i,j).gt.0.0)then
890
891	do ixt=1,ntraciso
892	xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
893	xtent(ixt,il,i,j) = (1.-Sigij(il,i,j))*xtrti(ixt) &
894	+ Sigij(il,i,j)*xt(ixt,il,i)
895	enddo !do ixt = 1, ntraciso
896	#ifdef ISOVERIF
897	#ifdef ISOTRAC
898	call iso_verif_traceur(xtrti,'cv3p_routines 714')
899	call iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 714')
900	#endif
901	if (iso_eau.gt.0) then
902	call iso_verif_egalite(xtent(iso_eau,il,i,j), &
903	qent(il,i,j),'cv3p_routines 718')
904	endif
905	#endif
906	endif !if(sij(il,i,j).gt.0.0)then
907	endif !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
908	enddo !do il=1,ncum
909
910	#ifdef ISOVERIFIDRIS
911	do il=1,ncum
912	call iso_verif_noNaN(qent(il,i,j),'cv3p_routines 771')
913	enddo !do il=1,ncum
914	#endif
915	#ifdef ISOVERIF
916	if (iso_eau.gt.0) then
917	do il=1,ncum
918	call iso_verif_egalite(xtent(iso_eau,il,i,j), &
919	qent(il,i,j),'cv3p_routines 744')
920	enddo !do il=1,ncum
921	endif
922	do il=1,ncum
923	call iso_verif_positif(qent(il,i,j)-elij(il,i,j), &
924	'cv3p_routines 749')
925	enddo !do il=1,ncum
926	!write(,) 'cv3p_mixing 788: call condiso_liq_ice_vectall'
927	#endif
928
929	call condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), &
930	elij(1,i,j), &
931	t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
932	#ifdef ISOTRAC
933	call condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), &
934	elij(1,i,j), &
935	t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
936	#endif
937	do il=1,ncum
938	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
939	j.ge.(icb(il)-1) .and. j.le.inb(il) &
940	.and. lwork(il) ) then
941	if(sij(il,i,j).gt.0.0)then
942	do ixt = 1, ntraciso
943	xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
944	enddo !do ixt = 1, ntraciso
945	endif !if(sij(il,i,j).gt.0.0)then
946	endif !if ( i.ge.icb(il) .and. i.le.inb(il) .and.
947	enddo !do il=1,ncum
948
949	#ifdef ISOVERIF
950	if (iso_eau.gt.0) then
951	do il=1,ncum
952	call iso_verif_egalite(xtent(iso_eau,il,i,j), &
953	qent(il,i,j),'cv3p_routines 744')
954	call iso_verif_egalite(xtelij(iso_eau,il,i,j), &
955	elij(il,i,j),'cv3p_routines 746')
956	enddo !do il=1,ncum
957	endif
958	#endif
959	!#ifdef ISOVERIF
960
961	#ifdef ISOTRAC
962	! write(,) 'cv3_routines tmp 1987,option_traceurs=',
963	! : option_traceurs
964	if (option_tmin.ge.1) then
965	do il=1,ncum
966	! write(,) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),',
967	! : 'tcond(il),rs(il,j)=',
968	! : il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j)
969	! colorier la vapeur résiduelle selon température de
970	! condensation, et le condensat en un tag spécifique
971	if ((elij(il,i,j).gt.0.0).and.(qent(il,i,j).gt.0.0)) then
972	if (option_traceurs.eq.17) then
973	call iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
974	xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
975	0.0,xtres, &
976	seuil_tag_tmin)
977	else !if (option_traceurs.eq.17) then
978	! write(,) 'cv3 2002: il,i,j =',il,i,j
979	call iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
980	xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, &
981	seuil_tag_tmin)
982	endif !if (option_traceurs.eq.17) then
983	do ixt=1+niso,ntraciso
984	xtent(ixt,il,i,j)=xtres(ixt)
985	enddo
986	endif !if (cond.gt.0.0) then
987	enddo !do il=1,ncum
988	#ifdef ISOVERIF
989	do il=1,ncum
990	call iso_verif_traceur(xtent(1,il,i,j),'cv3_routines 1996')
991	call iso_verif_traceur(xtelij(1,il,i,j),'cv3_routines 1997')
992	call iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
993	'cv3_routines 2042')
994	enddo !do il=1,ncum
995	#endif
996	endif !if (option_tmin.ge.1) then
997	#endif
998	! endif ISOTRAC
999
1000	#endif
1001	!!#ifdef ISO
1002
1003	!jyg! DO k = 1, ntra
1004	!jyg! DO il = 1, ncum
1005	!jyg! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
1006	!jyg! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
1007	!jyg! lwork(il)) THEN
1008	!jyg! IF (Sij(il,i,j)>0.0) THEN
1009	!jyg! traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
1010	!jyg! (1.-Sigij(il,i,j))*tra(il, nk(il), k)
1011	!jyg! END IF
1012	!jyg! END IF
1013	!jyg! END DO
1014	!jyg! END DO
1015
1016	! -- If I=J (detrainement and entrainement at the same level), then only the
1017	! -- adiabatic ascent part of the mixture is considered
1018	IF (i==j) THEN
1019	DO il = 1, ncum
1020	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
1021	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
1022	lwork(il)) THEN
1023	IF (Sij(il,i,j)>0.0) THEN
1024	!! rti = qnk(il) - ep(il, i)*clw(il, i)
1025	rti = qta(il,i-1) - ep(il, i)*clw(il, i)
1026	!!! Ment(il,i,i) = m(il,i)abs(Qmixmax(il)(1.-Sjmax(il))
1027	Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
1028	Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
1029	Qent(il, i, i) = rti
1030	uent(il, i, i) = unk(il)
1031	vent(il, i, i) = vnk(il)
1032	hent(il, i, i) = hp(il, i)
1033	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
1034	Sigij(il, i, i) = 0.
1035	END IF
1036	END IF
1037	END DO
1038
1039
1040	#ifdef ISO
1041	IF (I .EQ. J) THEN
1042	do il=1,ncum
1043	zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
1044	zfice(il) = MIN(MAX(zfice(il),0.0),1.0)
1045	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
1046	j.ge.(icb(il)-1) .and. j.le.inb(il) &
1047	.and. lwork(il) ) then
1048	if(sij(il,i,j).gt.0.0)then
1049	do ixt=1,ntraciso
1050	xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
1051	xtent(ixt,il,i,j) = xtrti(ixt)
1052	enddo !do ixt = 1, ntraciso
1053	#ifdef ISOVERIFIDRIS
1054	call iso_verif_noNaN(qent(il,i,j),'cv3p_routines 923')
1055	#endif
1056	#ifdef ISOVERIF
1057	#ifdef ISOTRAC
1058	call iso_verif_traceur(xtrti,'cv3p_routines 714')
1059	call iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 714')
1060	#endif
1061	if (iso_eau.gt.0) then
1062	call iso_verif_egalite(xtent(iso_eau,il,i,j), &
1063	qent(il,i,j),'cv3p_routines 718')
1064	endif
1065	call iso_verif_positif(qent(il,i,j)-elij(il,i,j), &
1066	'cv3p_routines 887')
1067	#endif
1068	endif !if(sij(il,i,j).gt.0.0)then
1069	endif !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1070	enddo !do il=1,ncum
1071
1072	!write(,) 'cv3p_mixing 808: call condiso'
1073	call condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), &
1074	elij(1,i,j), &
1075	t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
1076
1077	do il=1,ncum
1078	do ixt = 1, ntraciso
1079	xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
1080	enddo !do ixt = 1, ntraciso
1081	#ifdef ISOVERIF
1082	if (iso_eau.gt.0) then
1083	call iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1084	elij(il,i,j),'cv3p_routines 1074')
1085	endif ! if (iso_eau.gt.0) then
1086	#endif
1087	enddo !do il=1,ncum
1088
1089	#ifdef ISOTRAC
1090	call condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), &
1091	elij(1,i,j), &
1092	t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
1093	#ifdef ISOVERIF
1094	do il=1,ncum
1095	call iso_verif_traceur(xt(1,il,i),'cv3p_routines 1967')
1096	call iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 1969')
1097	enddo !do il=1,ncum
1098	#endif
1099
1100	if (option_tmin.ge.1) then
1101	do il=1,ncum
1102	! write(,) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),',
1103	! : 'tcond(il),rs(il,j)=',
1104	! : il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j)
1105	! colorier la vapeur résiduelle selon température de
1106	! condensation, et le condensat en un tag spécifique
1107	if ((elij(il,i,j).gt.0.0).and.(qent(il,i,j).gt.0.0)) then
1108	if (option_traceurs.eq.17) then
1109	call iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
1110	xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
1111	0.0,xtres, &
1112	seuil_tag_tmin)
1113	else !if (option_traceurs.eq.17) then
1114	! write(,) 'cv3 2002: il,i,j =',il,i,j
1115	call iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
1116	xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, &
1117	seuil_tag_tmin)
1118	endif !if (option_traceurs.eq.17) then
1119	do ixt=1+niso,ntraciso
1120	xtent(ixt,il,i,j)=xtres(ixt)
1121	enddo
1122	endif !if (cond.gt.0.0) then
1123	enddo !do il=1,ncum
1124	#ifdef ISOVERIF
1125	do il=1,ncum
1126	call iso_verif_traceur(xtent(1,il,i,j),'cv3_routines 1996')
1127	call iso_verif_traceur(xtelij(1,il,i,j),'cv3_routines 1997')
1128	call iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
1129	'cv3_routines 2042')
1130	enddo !do il=1,ncum
1131	#endif
1132	endif !if (option_tmin.ge.1) then
1133	!#ifdef ISOVERIF
1134	#endif
1135	!#ifdef ISOTRAC
1136	#ifdef ISOVERIF
1137	if (iso_eau.gt.0) then
1138	do il=1,ncum
1139	call iso_verif_egalite(xtent(iso_eau,il,i,j), &
1140	qent(il,i,j),'cv3p_routines 883')
1141	call iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1142	elij(il,i,j),'cv3p_routines 885')
1143	enddo !do il=1,ncum
1144	endif !if (iso_eau.gt.0) then
1145	#endif
1146	!#ifdef ISOVERIF
1147	endif !IF (I .EQ. J) THEN
1148	#endif
1149	!#ifdef ISO
1150
1151	!jyg! DO k = 1, ntra
1152	!jyg! DO il = 1, ncum
1153	!jyg! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
1154	!jyg! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
1155	!jyg! lwork(il)) THEN
1156	!jyg! IF (Sij(il,i,j)>0.0) THEN
1157	!jyg! traent(il, i, i, k) = tra(il, nk(il), k)
1158	!jyg! END IF
1159	!jyg! END IF
1160	!jyg! END DO
1161	!jyg! END DO
1162
1163	END IF
1164
1165	! ---------------------------------------------------------------
1166	175 END DO ! End loop on destination level "j"
1167	! ---------------------------------------------------------------
1168
1169	DO il = 1, ncum
1170	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1171	ASij(il) = amax1(1.0E-16, ASij(il))
1172	!jyg+lluis<
1173	!! ASij(il) = 1.0/ASij(il)
1174	ASij_inv(il) = 1.0/ASij(il)
1175	! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
1176	IF (ASij_inv(il) > 100.) ASij_inv(il) = 0.
1177	!>jyg+lluis
1178	csum(il, i) = 0.0
1179	END IF
1180	END DO
1181
1182	DO j = minorig, nl
1183	DO il = 1, ncum
1184	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1185	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
1186	!jyg Ment(il, i, j) = Ment(il, i, j)*ASij(il)
1187	Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
1188	END IF
1189	END DO
1190	END DO
1191
1192	DO j = minorig, nl
1193	DO il = 1, ncum
1194	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1195	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
1196	csum(il, i) = csum(il, i) + Ment(il, i, j)
1197	END IF
1198	END DO
1199	END DO
1200
1201	DO il = 1, ncum
1202	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
1203	! cc : .and. csum(il,i).lt.m(il,i) ) then
1204	nent(il, i) = 0
1205	! cc Ment(il,i,i)=m(il,i)
1206	Ment(il, i, i) = 1.
1207	!! Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
1208	Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
1209	uent(il, i, i) = unk(il)
1210	vent(il, i, i) = vnk(il)
1211	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
1212
1213	#ifdef ISO
1214	do ixt = 1, ntraciso
1215	xtent(ixt,il,i,i)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
1216	xtelij(ixt,il,i,i)=xtclw(ixt,il,i)*(1.-ep(il,i))
1217	enddo
1218	#ifdef ISOVERIF
1219	if (iso_eau.gt.0) then
1220	call iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1221	elij(il,i,j),'cv3p_routines 1204')
1222	endif ! if (iso_eau.gt.0) then
1223	#endif
1224	#endif
1225
1226	#ifdef ISOTRAC
1227	if (option_tmin.ge.1) then
1228	! colorier la vapeur résiduelle selon température de
1229	! condensation, et le condensat en un tag spécifique
1230	! write(,) 'cv3 tmp 2314 il,i,j,xtent(:,il,i,j)=',
1231	! : il,i,j,xtent(:,il,i,j)
1232	if ((elij(il,i,i).gt.0.0).and.(qent(il,i,i).gt.0.0)) then
1233	if (option_traceurs.eq.17) then
1234	call iso_recolorise_condensation(qent(il,i,i), &
1235	elij(il,i,i), &
1236	xt(1,il,1),xtclw(1,il,i),t(il,i),ep(il,i), &
1237	xtres, &
1238	seuil_tag_tmin)
1239	else !if (option_traceurs.eq.17) then
1240	call iso_recolorise_condensation(qent(il,i,i), &
1241	elij(il,i,i), &
1242	xt(1,il,1),xtclw(1,il,i),rs(il,i),ep(il,i), &
1243	xtres, &
1244	seuil_tag_tmin)
1245	endif ! if (option_traceurs.eq.17) then
1246	do ixt=1+niso,ntraciso
1247	xtent(ixt,il,i,i)=xtres(ixt)
1248	enddo
1249	#ifdef ISOVERIF
1250	do ixt=1,niso
1251	call iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), &
1252	'cv3p_mixing 2318',errmax,errmaxrel)
1253	call iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
1254	'cv3p_mixing 2383')
1255	enddo
1256	#endif
1257	endif !if (cond.gt.0.0) then
1258	#ifdef ISOVERIF
1259	call iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
1260	qent(il,i,i),'cv3p_mixing 2321',errmax,errmaxrel)
1261	call iso_verif_egalite_choix(xtelij(iso_eau,il,i,i), &
1262	elij(il,i,i),'cv3p_mixing 2321b',errmax,errmaxrel)
1263	call iso_verif_traceur(xtent(1,il,i,i),'cv3p_mixing 2322')
1264	call iso_verif_traceur(xtelij(1,il,i,i),'cv3p_mixing 2323')
1265	#endif
1266	endif !if (option_tmin.ge.1) then
1267	#endif
1268	IF (fl_cor_ebil .GE. 2) THEN
1269	hent(il, i, i) = hp(il,i)
1270	Sigij(il, i, i) = 0.0
1271	ELSE
1272	Sij(il, i, i) = 0.0
1273	ENDIF
1274	END IF
1275	END DO ! il
1276
1277	!jyg! DO j = 1, ntra
1278	!jyg! DO il = 1, ncum
1279	!jyg! IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
1280	!jyg!! cc : .and. csum(il,i).lt.m(il,i) ) then
1281	!jyg! traent(il, i, i, j) = tra(il, nk(il), j)
1282	!jyg! END IF
1283	!jyg! END DO
1284	!jyg! END DO
1285
1286	! ---------------------------------------------------------------
1287	789 END DO ! End loop on origin level "i"
1288	! ---------------------------------------------------------------
1289
1290	#ifdef ISO
1291	#ifdef ISOVERIF
1292	! write(,) 'cv3p_mixing 2540: ', &
1293	! 'verif finale en sortant de cv3p_mixing'
1294	! write(,) 'qent,xtent(1,1,1)=',qent(1,1,1),xtent(iso_eau,1,1,1)
1295	do im = 1, nd
1296	do jm = 1, nd
1297	do il = 1, ncum
1298	if (iso_eau.gt.0) then
1299	call iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), &
1300	elij(il,im,jm),'cv3p_mixing 2110',errmax,errmaxrel)
1301	call iso_verif_egalite_choix(xtent(iso_eau,il,im,jm), &
1302	qent(il,im,jm),'cv3p_mixing 2112',errmax,errmaxrel)
1303	endif !if (iso_eau>0) then
1304	#ifdef ISOTRAC
1305	call iso_verif_traceur_justmass(xtelij(1,il,im,jm), &
1306	'cv3p_routine 2250')
1307	#endif
1308	enddo !do il = 1, nloc
1309	enddo !do jm = 1, klev
1310	enddo !do im = 1, klev
1311	#endif
1312	#endif
1313
1314	#ifdef ISO
1315	#ifdef ISOTRAC
1316	! seulement à la fin on taggue le condensat
1317	if (option_cond.ge.1) then
1318	do im = 1, nd
1319	do jm = 1, nd
1320	do il = 1, ncum
1321	! colorier le condensat en un tag spécifique
1322	do ixt=niso+1,ntraciso
1323	if (index_zone(ixt).eq.izone_cond) then
1324	xtelij(ixt,il,im,jm)=xtelij(index_iso(ixt),il,im,jm)
1325	else !if (index_zone(ixt).eq.izone_cond) then
1326	xtelij(ixt,il,im,jm)=0.0
1327	endif !if (index_zone(ixt).eq.izone_cond) then
1328	enddo !do ixt=1,ntraciso
1329	#ifdef ISOVERIF
1330	call iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), &
1331	elij(il,im,jm),'cv3p_mixing 2408',errmax,errmaxrel)
1332	call iso_verif_traceur(xtelij(1,il,im,jm), &
1333	'cv3p_mixing 358')
1334	#endif
1335	enddo !do il = 1, ncum
1336	enddo !do jm = 1, nd
1337	enddo !do im = 1, nd
1338	! 17 juin 2020: on supprime ci-dessous car normalement xtclw a déjà été
1339	! recolorié ds cv3_undilute
1340	! do im = 1, nd
1341	! do il = 1, ncum
1342	! ! colorier le condensat en un tag spécifique
1343	! do ixt=niso+1,ntraciso
1344	! if (index_zone(ixt).eq.izone_cond) then
1345	! xtclw(ixt,il,im)=xtclw(index_iso(ixt),il,im)
1346	! else !if (index_zone(ixt).eq.izone_cond) then
1347	! xtclw(ixt,il,im)=0.0
1348	! endif !if (index_zone(ixt).eq.izone_cond) then
1349	! enddo !do ixt=1,ntraciso
1350	!#ifdef ISOVERIF
1351	! call iso_verif_egalite_choix(xtclw(iso_eau,il,im), &
1352	! clw(il,im),'cv3p_mixing 2427',errmax,errmaxrel)
1353	! call iso_verif_traceur(xtclw(1,il,im), &
1354	! 'cv3p_mixing 358')
1355	! if (iso_verif_positif_nostop(xtclw(itZonIso( &
1356	! izone_cond,iso_eau),i,k)-xtclw(iso_eau,i,k) &
1357	! ,'cv3p_mixing 909').eq.1) then
1358	! write(,) 'i,k=',i,k
1359	! write(,) 'xtclw=',xtclw(:,i,k)
1360	! write(,) 'niso,ntraciso,index_zone,izone_cond=', &
1361	! niso,ntraciso,index_zone,izone_cond
1362	! stop
1363	! endif !if (iso_verif_positif_nostop(xtclw(itZonIso(
1364	!#endif
1365	! enddo !do il = 1, ncum
1366	! enddo !do im = 1, nd
1367	! write(,) 'xtclw(:,1,2)=',xtclw(:,1,2)
1368	endif !if (option_tmin.eq.1) then
1369	#endif
1370	#endif
1371
1372
1373	RETURN
1374	END SUBROUTINE cv3p_mixing
1375

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

Download in other formats:

Original Format