Context Navigation

cv3p_mixing.F90 @ 5275

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

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

Download in other formats:

Original Format