Context Navigation

cv3p_mixing.F90 @ 4491

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

Rev	Line
[3927]	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
[4143]	23	USE infotrac_phy, ONLY: ntraciso=>ntiso
[3927]	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), &
[4491]	453	ridicule,deltalim_snow,'cv3p_mixing 1993')
[3927]	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
[4033]	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)
[3927]	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)
[4143]	1303	endif !if (iso_eau>0) then
[3927]	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')
[4143]	1355	! if (iso_verif_positif_nostop(xtclw(itZonIso( &
[3927]	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
[4143]	1363	! endif !if (iso_verif_positif_nostop(xtclw(itZonIso(
[3927]	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