Context Navigation

cv3p_mixing.F90 @ 5276

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

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

Download in other formats:

Original Format