Context Navigation

cv3p_mixing.f90 @ 3777

Last change on this file since 3777 was 3717, checked in by adurocher, 5 years ago

cv3p_mixing_new : VLAs removed from block construct (intel 19 bug)

Intel compiler 19 crashes (internal compiler error) when VLA (variable length arrays) are used within block constructs

File size: 47.4 KB

Rev	Line
[3709]	1	module cv3p_mixing_mod
	2	implicit none
	3
	4	contains
	5
	6	! **************************************************************
	7	! *
	8	! CV3P_MIXING : compute mixed draught properties and, *
	9	! within a scaling factor, mixed draught *
	10	! mass fluxes. *
	11	! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
	12	! modified by : *
	13	! **************************************************************
	14	SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
	15	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
	16	unk, vnk, hp, tv, tvp, ep, clw, sig, &
	17	Ment, Qent, hent, uent, vent, nent, &
	18	Sigij, elij, supmax, Ments, Qents, traent)
	19
	20	!inputs:
	21	INTEGER, INTENT(IN) :: ncum, nd, na
	22	INTEGER, INTENT(IN) :: ntra, nloc
	23	INTEGER, DIMENSION(nloc), INTENT(IN) :: icb, inb, nk
	24	REAL, DIMENSION(nloc, nd), INTENT(IN) :: sig
	25	REAL, DIMENSION(nloc), INTENT(IN) :: unk, vnk
	26	REAL, DIMENSION(nloc, nd), INTENT(IN) :: qta
	27	REAL, DIMENSION(nloc, nd + 1), INTENT(IN) :: ph
	28	REAL, DIMENSION(nloc, nd), INTENT(IN) :: t, rr, rs
	29	REAL, DIMENSION(nloc, nd), INTENT(IN) :: u, v
	30	REAL, DIMENSION(nloc, nd, ntra), INTENT(IN) :: tra ! input of convect3
	31	REAL, DIMENSION(nloc, na), INTENT(IN) :: lv
	32	REAL, DIMENSION(nloc, na), INTENT(IN) :: lf
	33	REAL, DIMENSION(nloc, na), INTENT(IN) :: frac !ice fraction in condensate
	34	REAL, DIMENSION(nloc, na), INTENT(IN) :: h !liquid water static energy of environment
	35	REAL, DIMENSION(nloc, na), INTENT(IN) :: hp !liquid water static energy of air shed from adiab. asc.
	36	REAL, DIMENSION(nloc, na), INTENT(IN) :: tv, tvp
	37	REAL, DIMENSION(nloc, na), INTENT(IN) :: ep, clw
	38
	39	!outputs:
	40	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Ment, Qent
	41	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: uent, vent
	42	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Sigij, elij
	43	REAL, DIMENSION(nloc, na), INTENT(OUT) :: supmax ! Highest mixing fraction of mixed
	44	! updraughts with the sign of (h-hp)
	45	REAL, DIMENSION(nloc, nd, nd, ntra), INTENT(OUT) :: traent
	46	REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: Ments, Qents
	47	REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: hent
	48	INTEGER, DIMENSION(nloc, nd), INTENT(OUT) :: nent
	49
	50	logical, parameter :: use_old = .false.
	51
	52	if (use_old) then
	53	call cv3p_mixing_old(nloc, ncum, nd, na, ntra, icb, nk, inb, &
	54	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
	55	unk, vnk, hp, tv, tvp, ep, clw, sig, &
	56	Ment, Qent, hent, uent, vent, nent, &
	57	Sigij, elij, supmax, Ments, Qents, traent)
	58	else
	59	call cv3p_mixing_new(nloc, ncum, nd, na, icb, inb, &
	60	t, rr, rs, u, v, h, lv, lf, frac, qta, &
	61	unk, vnk, hp, ep, clw, &
	62	Ment, Qent, hent, uent, vent, nent, &
	63	Sigij, elij, supmax)
	64	endif
	65
	66	END SUBROUTINE
	67
	68	! modified by : Arnaud Durocher (04/2020) : changed loops order
	69	SUBROUTINE cv3p_mixing_new(nloc, ncum, nd, na, icb, inb, &
	70	t, rr, rs, u, v, h, lv, lf, frac, qta, &
	71	unk, vnk, hp, ep, clw, &
	72	Ment, Qent, hent, uent, vent, nent, &
	73	Sigij, elij, supmax)
	74
	75	USE print_control_mod, ONLY: prt_level
	76	USE add_phys_tend_mod, ONLY: fl_cor_ebil
	77
	78	IMPLICIT NONE
	79
	80	include "cvthermo.h"
	81	include "cv3param.h"
	82	include "YOMCST2.h"
	83	include "cvflag.h"
	84
	85	!inputs:
	86	INTEGER, INTENT(IN) :: ncum, nd, na
	87	INTEGER, INTENT(IN) :: nloc
	88	INTEGER, DIMENSION(nloc), INTENT(IN) :: icb, inb
	89	REAL, DIMENSION(nloc), INTENT(IN) :: unk, vnk
	90	REAL, DIMENSION(nloc, nd), INTENT(IN) :: qta
	91	REAL, DIMENSION(nloc, nd), INTENT(IN) :: t, rr, rs
	92	REAL, DIMENSION(nloc, nd), INTENT(IN) :: u, v
	93	REAL, DIMENSION(nloc, na), INTENT(IN) :: lv
	94	REAL, DIMENSION(nloc, na), INTENT(IN) :: lf
	95	REAL, DIMENSION(nloc, na), INTENT(IN) :: frac !ice fraction in condensate
	96	REAL, DIMENSION(nloc, na), INTENT(IN) :: h !liquid water static energy of environment
	97	REAL, DIMENSION(nloc, na), INTENT(IN) :: hp !liquid water static energy of air shed from adiab. asc.
	98	REAL, DIMENSION(nloc, na), INTENT(IN) :: ep, clw
	99
	100	!outputs:
	101	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Ment, Qent
	102	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: uent, vent
	103	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Sigij, elij
	104	REAL, DIMENSION(nloc, na), INTENT(OUT) :: supmax ! Highest mixing fraction of mixed
	105	! updraughts with the sign of (h-hp)
	106	REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: hent
	107	INTEGER, DIMENSION(nloc, nd), INTENT(OUT) :: nent
[3710]	108
	109	!local variables:
[3716]	110	INTEGER i, j, k, il
[3715]	111	REAL, DIMENSION(nloc, nd) :: ASij
	112	REAL, DIMENSION(nloc, nd, nd) :: Sij
[3717]	113	real :: Sjmin(nd), Sjmax(nd), Qmixmax(nd), Qmixmin(nd), Rmixmax(nd), Rmixmin(nd)
[3710]	114
[3716]	115	INTEGER, SAVE :: igout = 1
[3712]	116	!$omp THREADPRIVATE(igout)
[3710]	117
[3712]	118	! -- Mixing probability distribution functions
[3710]	119
[3709]	120	REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
[3710]	121
[3709]	122	Qcoef1(F) = tanh(F/gammas)
	123	Qcoef2(F) = (tanh(F/gammas) + gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
	124	QFF(F) = max(min(F, 1.), 0.)
	125	QFFf(F) = min(QFF(F), scut)
	126	Qmix1(F) = (tanh((QFF(F) - Fmax)/gammas) + Qcoef1max)/Qcoef2max
	127	Rmix1(F) = (gammaslog(cosh((QFF(F) - Fmax)/gammas)) + QFF(F)Qcoef1max)/Qcoef2max
	128	Qmix2(F) = -log(1.-QFFf(F))/scut
	129	Rmix2(F) = (QFFf(F) + (1.-QFF(F))*log(1.-QFFf(F)))/scut
	130	Qmix(F) = qqa1Qmix1(F) + qqa2Qmix2(F)
	131	Rmix(F) = qqa1Rmix1(F) + qqa2Rmix2(F)
[3710]	132
[3709]	133	INTEGER, SAVE :: ifrst = 0
[3712]	134	!$omp THREADPRIVATE(ifrst)
[3710]	135
[3712]	136	! =====================================================================
	137	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	138	! =====================================================================
	139	! -- Initialize mixing PDF coefficients
[3709]	140	IF (ifrst == 0) THEN
	141	ifrst = 1
	142	Qcoef1max = Qcoef1(Fmax)
	143	Qcoef2max = Qcoef2(Fmax)
	144	print *, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
	145	fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
	146	END IF
[3710]	147
[3709]	148	nent(:ncum, :nl) = 0
	149	elij(:ncum, :nl, :nl) = 0
	150	hent(:ncum, :nl, :nl) = 0
[3710]	151
[3709]	152	DO j = 1, nl
	153	DO k = 1, nl
	154	DO i = 1, ncum
	155	Qent(i, k, j) = rr(i, j)
	156	uent(i, k, j) = u(i, j)
	157	vent(i, k, j) = v(i, j)
	158	END DO
	159	END DO
	160	END DO
[3710]	161
[3709]	162	Ment(1:ncum, 1:nd, 1:nd) = 0.0
	163	Sij(1:ncum, 1:nd, 1:nd) = 0.0
[3710]	164
[3712]	165	! =====================================================================
	166	! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
	167	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	168	! --- FRACTION (Sij)
	169	! =====================================================================
[3710]	170
[3714]	171	IF (ok_entrain) THEN
	172	! NOTE : (04/2020 AD) this loop order gave good results in cv3a_tracer
	173	! if( icb(il) <= i <= inb(il) ) is probably more expensive than reading non-contiguously
	174	DO il = 1, ncum
	175	DO i = icb(il), inb(il)
	176	DO j = (icb(il) - 1), inb(il)
[3716]	177	block
	178	real :: rti, bf2, anum, denom, dei, altem, cwat, stemp
	179
	180	rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
	181	bf2 = 1.+lv(il, j)lv(il, j)rs(il, j)/(rrvt(il, j)t(il, j)*cpd)
	182	IF (cvflag_ice .AND. t(il, j) <= 263.15) THEN
	183	bf2 = 1.+(lf(il, j) + lv(il, j))(lv(il, j) + frac(il, j) &
	184	lf(il, j))rs(il, j)/(rrvt(il, j)t(il, j)cpd)
	185	END IF
	186	anum = h(il, j) - hp(il, i) + (cpv - cpd)t(il, j)(rti - rr(il, j))
[3714]	187	denom = h(il, i) - hp(il, i) + (cpd - cpv)(rr(il, i) - rti)t(il, j)
	188	dei = denom
	189	IF (abs(dei) < 0.01) dei = 0.01
	190	Sij(il, i, j) = anum/dei
	191	Sij(il, i, i) = 1.0
[3713]	192	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il, i, j))rti - rs(il, j)
[3714]	193	altem = altem/bf2
	194	cwat = clw(il, j)*(1.-ep(il, j))
[3716]	195	stemp = Sij(il, i, j)
	196	IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN
	197	IF (cvflag_ice) THEN
	198	anum = anum - (lv(il, j) + frac(il, j)lf(il, j))(rti - rs(il, j) - cwat*bf2)
	199	denom = denom + (lv(il, j) + frac(il, j)lf(il, j))(rr(il, i) - rti)
	200	ELSE
	201	anum = anum - lv(il, j)(rti - rs(il, j) - cwatbf2)
	202	denom = denom + lv(il, j)*(rr(il, i) - rti)
	203	END IF
	204	IF (abs(denom) < 0.01) denom = 0.01
	205	Sij(il, i, j) = anum/denom
	206	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il, i, j))rti - rs(il, j)
	207	altem = altem - (bf2 - 1.)*cwat
[3714]	208	END IF
[3716]	209	IF (Sij(il, i, j) > 0.0) THEN
	210	Ment(il, i, j) = 1.
	211	elij(il, i, j) = altem
	212	elij(il, i, j) = amax1(0.0, elij(il, i, j))
	213	nent(il, i) = nent(il, i) + 1
	214	END IF
[3714]	215
[3716]	216	Sij(il, i, j) = amax1(0.0, Sij(il, i, j))
	217	Sij(il, i, j) = amin1(1.0, Sij(il, i, j))
	218	end block
[3714]	219	END DO
	220	IF (nent(il, i) == 0) THEN
	221	Ment(il, i, i) = 1.
	222	Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
	223	uent(il, i, i) = unk(il)
	224	vent(il, i, i) = vnk(il)
	225	IF (fl_cor_ebil .GE. 2) THEN
	226	hent(il, i, i) = hp(il, i)
	227	ENDIF
	228	elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
	229	Sij(il, i, i) = 0.0
[3713]	230	END IF
[3714]	231	END DO
	232	END DO
	233	ELSE
	234	DO i = minorig + 1, nl
	235	DO il = 1, ncum
	236	IF ((i >= icb(il)) .AND. (i <= inb(il))) THEN
	237	Ment(il, i, i) = 1.
	238	Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
	239	uent(il, i, i) = unk(il)
	240	vent(il, i, i) = vnk(il)
	241	IF (fl_cor_ebil .GE. 2) THEN
	242	hent(il, i, i) = hp(il, i)
	243	ENDIF
	244	elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
	245	Sij(il, i, i) = 0.0
[3713]	246	END IF
[3709]	247	END DO
[3713]	248	END DO
[3714]	249	ENDIF
	250
[3713]	251	DO i = minorig + 1, nl
[3714]	252	DO j = minorig, nl
	253	IF (prt_level >= 10) THEN
	254	print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout)
	255	!IF (nent(igout, i) .gt. 0) THEN
	256	! print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout))
	257	!ENDIF
	258	ENDIF
[3709]	259	END DO
[3713]	260	END DO
[3710]	261
[3714]	262	Sigij(1:ncum, 1:nd, 1:nd) = 0.0
	263	DO il = 1, ncum
	264	DO i = icb(il), inb(il)
	265	DO j = (icb(il) - 1), inb(il)
	266	Sigij(il, i, j) = Sij(il, i, j)
	267	END DO
[3709]	268	END DO
	269	END DO
[3710]	270
[3712]	271	! =====================================================================
	272	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	273	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	274	! =====================================================================
[3709]	275	DO il = 1, ncum
[3712]	276	DO i = icb(il), inb(il) !Loop on origin level "i"
[3715]	277	block
[3716]	278	real :: signhpmh, Scrit
[3710]	279
[3716]	280	block
	281	real :: rti, anum, denom, alt, Sx
[3710]	282
[3716]	283	signhpmh = sign(1., hp(il, i) - h(il, i))
[3710]	284
[3716]	285	rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
	286	IF (cvflag_ice) THEN
	287	anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)lf(il, i)) &
	288	(rti - rs(il, i)) + (cpv - cpd)t(il, i)(rti - rr(il, i))
	289	denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)lf(il, i)) &
	290	(rr(il, i) - rti) + (cpd - cpv)t(il, i)(rr(il, i) - rti)
	291	ELSE
	292	anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + &
	293	(cpv - cpd)t(il, i)(rti - rr(il, i))
	294	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + &
	295	(cpd - cpv)t(il, i)(rr(il, i) - rti)
	296	END IF
	297	IF (abs(denom) < 0.01) denom = 0.01
	298	Scrit = min(anum/denom, 1.)
	299	alt = rti - rs(il, i) + Scrit*(rr(il, i) - rti)
[3710]	300
[3716]	301	if (alt <= 0) then
	302	Scrit = 1.0
	303	else
	304	! Get max of Sij
	305	Sx = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh)
	306	! Find new critical value Scrit
	307	! such that : Sij > Scrit => mixed draught will detrain at J<I
	308	! Sij < Scrit => mixed draught will detrain at J>I
	309	Scrit = max(0., min(Scrit, Sxmax(0., -signhpmh)) + Scritmax(0., signhpmh))
	310	endif
[3710]	311
[3716]	312	end block
	313
[3715]	314	ASij(il, i) = 0.0
[3710]	315
[3716]	316	! Compute all Sjmax(j) and Sjmin(j)
	317	block
	318	real :: smin, smax, sup, Sbef, Smid
	319	smin = 1
	320	smax = 0.0
	321	sup = 0. ! upper S-value reached by descending draughts
	322	! Glitchy : why so complicated?
	323	if (i < inb(il)) then
	324	Sbef = Sij(il, i, inb(il))
	325	else
	326	Sbef = max(0., signhpmh)
	327	endif
[3710]	328
[3716]	329	DO j = (icb(il) - 1), i - 1
	330	IF (Sij(il, i, j) > 0.0) THEN
	331	Smid = max(Sij(il, i, j), Scrit)
	332	Sjmax(j) = Smid
	333	Sjmin(j) = Smid
	334	IF (Smid > smax .AND. Sij(il, i, j + 1) > Smid) THEN
	335	smax = Smid
	336	Sjmax(j) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j))
	337	Sjmax(j) = max(Sjmax(j), Scrit)
	338	Sjmin(j) = min((Sbef + Sij(il, i, j))/2., Sij(il, i, j))
	339	Sjmin(j) = max(Sjmin(j), Scrit)
	340	Sbef = Sij(il, i, j)
	341	END IF
	342	IF (abs(Sjmin(j) - Sjmax(j)) > 1.E-10) &
	343	sup = max(Sjmin(j), Sjmax(j), sup)
	344	ENDIF
	345	END DO
	346	j = i
	347	IF (Sij(il, i, j) > 0.0) THEN
	348	Smid = 1.
	349	Sjmin(j) = max((Sij(il, i, j - 1) + Smid)/2., Scrit)*max(0., -signhpmh) + &
	350	min((Sij(il, i, j + 1) + Smid)/2., Scrit)*max(0., signhpmh)
	351	Sjmin(j) = max(Sjmin(j), sup)
	352	Sjmax(j) = 1.
	353	! preparation des variables Scrit, Smin et Sbef pour la partie j>i
	354	Scrit = min(Sjmin(j), Sjmax(j), Scrit)
	355	Sbef = max(0., signhpmh)
	356	supmax(il, i) = sign(Scrit, -signhpmh)
	357	ENDIF
	358	DO j = i + 1, inb(il)
	359	IF (Sij(il, i, j) > 0.0) THEN
	360	Smid = min(Sij(il, i, j), Scrit)
	361	Sjmax(j) = Smid
	362	Sjmin(j) = Smid
	363	IF (Smid < smin .AND. Sij(il, i, j + 1) < Smid) THEN
	364	smin = min(smin, Smid)
	365	Sjmax(j) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit)
	366	Sjmin(j) = max((Sbef + Sij(il, i, j))/2., Sij(il, i, j))
	367	Sjmin(j) = min(Sjmin(j), Scrit)
	368	Sbef = Sij(il, i, j)
	369	END IF
	370	ENDIF
	371	END DO
	372	end block
	373
	374	DO j = (icb(il) - 1), inb(il)
	375	Qmixmax(j) = Qmix(Sjmax(j))
	376	Qmixmin(j) = Qmix(Sjmin(j))
	377	Rmixmax(j) = Rmix(Sjmax(j))
	378	Rmixmin(j) = Rmix(Sjmin(j))
	379	END DO
	380
[3715]	381	DO j = (icb(il) - 1), inb(il) !Loop on destination level "j"
	382	IF (Sij(il, i, j) > 0.0) THEN
	383	block
[3716]	384	real :: sqmrmax, sqmrmin, cpm, awat, Tm, rti
[3710]	385
[3715]	386	rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
[3710]	387
[3716]	388	sqmrmax = Sjmax(j)*Qmixmax(j) - Rmixmax(j)
	389	sqmrmin = Sjmin(j)*Qmixmin(j) - Rmixmin(j)
	390
	391	Ment(il, i, j) = abs(Qmixmax(j) - Qmixmin(j))*Ment(il, i, j)
	392	IF (abs(Qmixmax(j) - Qmixmin(j)) > 1.E-10) THEN
	393	Sigij(il, i, j) = (sqmrmax - sqmrmin)/(Qmixmax(j) - Qmixmin(j))
[3715]	394	ELSE
	395	Sigij(il, i, j) = 0.
	396	END IF
[3710]	397
[3715]	398	! Compute Qent, uent, vent according to the true mixing fraction
	399	Qent(il, i, j) = (1.-Sigij(il, i, j))rti + Sigij(il, i, j)rr(il, i)
	400	uent(il, i, j) = (1.-Sigij(il, i, j))unk(il) + Sigij(il, i, j)u(il, i)
	401	vent(il, i, j) = (1.-Sigij(il, i, j))vnk(il) + Sigij(il, i, j)v(il, i)
[3710]	402
[3715]	403	! Compute liquid water static energy of mixed draughts
	404	hent(il, i, j) = (1.-Sigij(il, i, j))hp(il, i) + Sigij(il, i, j)h(il, i)
	405	! Heat capacity of mixed draught
	406	cpm = cpd + Qent(il, i, j)*(cpv - cpd)
[3710]	407
[3715]	408	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
	409	elij(il, i, j) = elij(il, i, j) + &
	410	(h(il, j) - hent(il, i, j) + (cpv - cpd)(Qent(il, i, j) - rr(il, j))t(il, j))* &
	411	rs(il, j)lv(il, j)/(cpmrrvt(il, j)t(il, j))
	412	IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
	413	elij(il, i, j) = elij(il, i, j)/ &
	414	(1.+(lv(il, j) + frac(il, j)lf(il, j))lv(il, j)*rs(il, j)/ &
	415	(cpmrrvt(il, j)*t(il, j)))
	416	ELSE
[3710]	417
[3715]	418	elij(il, i, j) = elij(il, i, j)/ &
	419	(1.+lv(il, j)lv(il, j)rs(il, j)/ &
	420	(cpmrrvt(il, j)*t(il, j)))
	421	ENDIF
	422	elij(il, i, j) = max(elij(il, i, j), 0.)
	423	elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j))
[3710]	424
[3715]	425	IF (j > i) THEN
	426	awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j)
	427	awat = amax1(awat, 0.0)
	428	ELSE
	429	awat = 0.
	430	END IF
[3710]	431
[3715]	432	! Mixed draught temperature at level j
	433	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))
	434	IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
	435	hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)lf(il, j) + (cpd - cpv)Tm)*awat
	436	ELSE
	437	hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)Tm)awat
	438	ENDIF
[3710]	439
[3715]	440	! ASij is the integral of P(F) over the relevant F interval
[3716]	441	ASij(il, i) = ASij(il, i) + abs(Qmixmax(j)*(1.-Sjmax(j)) + Rmixmax(j) - &
	442	Qmixmin(j)*(1.-Sjmin(j)) - Rmixmin(j))
[3715]	443
	444	! If I=J (detrainement and entrainement at the same level), then only the
	445	! adiabatic ascent part of the mixture is considered
	446	IF (i == j) THEN
	447	rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
[3716]	448	Ment(il, i, i) = abs(Qmixmax(j)*(1.-Sjmax(j)) + Rmixmax(j) - &
	449	Qmixmin(j)*(1.-Sjmin(j)) - Rmixmin(j))
[3715]	450	Qent(il, i, i) = rti
	451	uent(il, i, i) = unk(il)
	452	vent(il, i, i) = vnk(il)
	453	hent(il, i, i) = hp(il, i)
	454	elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
	455	Sigij(il, i, i) = 0.
	456	END IF
	457	end block
[3712]	458	END IF
[3715]	459	END DO ! Loop j = (icb(il) - 1), inb(il)
	460	end block
[3712]	461	END DO ! Loop i = icb(il), inb(il)
	462	END DO ! Loop il = 1, ncum
	463
	464	DO il = 1, ncum
	465	DO i = icb(il), inb(il) !Loop on origin level "i"
[3715]	466	block
	467	real :: csum, ASij_inv
	468	csum = 0
	469	DO j = (icb(il) - 1), inb(il)
	470	ASij(il, i) = amax1(1.0E-16, ASij(il, i))
	471	ASij_inv = 1.0/ASij(il, i)
	472	! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
	473	IF (ASij_inv > 100.) ASij_inv = 0.
	474	Ment(il, i, j) = Ment(il, i, j)*ASij_inv
	475	csum = csum + Ment(il, i, j)
	476	END DO
	477	IF (csum < 1.) THEN
	478	nent(il, i) = 0
	479	Ment(il, i, i) = 1.
	480	Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
	481	uent(il, i, i) = unk(il)
	482	vent(il, i, i) = vnk(il)
	483	elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
	484	IF (fl_cor_ebil .GE. 2) THEN
	485	hent(il, i, i) = hp(il, i)
	486	Sigij(il, i, i) = 0.0
	487	ELSE
	488	Sij(il, i, i) = 0.0
	489	ENDIF
[3712]	490	ENDIF
[3715]	491	end block
[3711]	492	END DO ! Loop il = 1, ncum
	493	END DO ! End loop on origin level "i"
[3709]	494	END SUBROUTINE
	495
	496	! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15
	497	SUBROUTINE cv3p_mixing_old(nloc, ncum, nd, na, ntra, icb, nk, inb, &
	498	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
	499	unk, vnk, hp, tv, tvp, ep, clw, sig, &
	500	Ment, Qent, hent, uent, vent, nent, &
	501	Sigij, elij, supmax, Ments, Qents, traent)
	502
	503	USE print_control_mod, ONLY: prt_level
	504	USE add_phys_tend_mod, ONLY: fl_cor_ebil
	505
	506	IMPLICIT NONE
	507
	508	include "cvthermo.h"
	509	include "cv3param.h"
	510	include "YOMCST2.h"
	511	include "cvflag.h"
	512
	513	!inputs:
	514	INTEGER, INTENT(IN) :: ncum, nd, na
	515	INTEGER, INTENT(IN) :: ntra, nloc
	516	INTEGER, DIMENSION(nloc), INTENT(IN) :: icb, inb, nk
	517	REAL, DIMENSION(nloc, nd), INTENT(IN) :: sig
	518	REAL, DIMENSION(nloc), INTENT(IN) :: unk, vnk
	519	REAL, DIMENSION(nloc, nd), INTENT(IN) :: qta
	520	REAL, DIMENSION(nloc, nd + 1), INTENT(IN) :: ph
	521	REAL, DIMENSION(nloc, nd), INTENT(IN) :: t, rr, rs
	522	REAL, DIMENSION(nloc, nd), INTENT(IN) :: u, v
	523	REAL, DIMENSION(nloc, nd, ntra), INTENT(IN) :: tra ! input of convect3
	524	REAL, DIMENSION(nloc, na), INTENT(IN) :: lv
	525	REAL, DIMENSION(nloc, na), INTENT(IN) :: lf
	526	REAL, DIMENSION(nloc, na), INTENT(IN) :: frac !ice fraction in condensate
	527	REAL, DIMENSION(nloc, na), INTENT(IN) :: h !liquid water static energy of environment
	528	REAL, DIMENSION(nloc, na), INTENT(IN) :: hp !liquid water static energy of air shed from adiab. asc.
	529	REAL, DIMENSION(nloc, na), INTENT(IN) :: tv, tvp
	530	REAL, DIMENSION(nloc, na), INTENT(IN) :: ep, clw
	531
	532	!outputs:
	533	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Ment, Qent
	534	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: uent, vent
	535	REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Sigij, elij
	536	REAL, DIMENSION(nloc, na), INTENT(OUT) :: supmax ! Highest mixing fraction of mixed
	537	! updraughts with the sign of (h-hp)
	538	REAL, DIMENSION(nloc, nd, nd, ntra), INTENT(OUT) :: traent
	539	REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: Ments, Qents
	540	REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: hent
	541	INTEGER, DIMENSION(nloc, nd), INTENT(OUT) :: nent
	542
	543	!local variables:
	544	INTEGER i, j, k, il, im, jm
	545	INTEGER num1, num2
	546	REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp
	547	REAL :: alt, delp, delm
	548	REAL, DIMENSION(nloc) :: Qmixmax, Rmixmax, sqmrmax
	549	REAL, DIMENSION(nloc) :: Qmixmin, Rmixmin, sqmrmin
	550	REAL, DIMENSION(nloc) :: signhpmh
	551	REAL, DIMENSION(nloc) :: Sx
	552	REAL :: Scrit2
	553	REAL, DIMENSION(nloc) :: Smid, Sjmin, Sjmax
	554	REAL, DIMENSION(nloc) :: Sbef, sup, smin
	555	REAL, DIMENSION(nloc) :: ASij, ASij_inv, smax, Scrit
	556	REAL, DIMENSION(nloc, nd, nd) :: Sij
	557	REAL, DIMENSION(nloc, nd) :: csum
	558	REAL :: awat
	559	REAL :: cpm !Mixed draught heat capacity
	560	REAL :: Tm !Mixed draught temperature
	561	LOGICAL, DIMENSION(nloc) :: lwork
	562
	563	REAL amxupcrit, df, ff
	564	INTEGER nstep
	565
	566	INTEGER, SAVE :: igout = 1
	567	!$omp THREADPRIVATE(igout)
	568
	569	! -- Mixing probability distribution functions
	570
	571	REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
	572
	573	Qcoef1(F) = tanh(F/gammas)
	574	Qcoef2(F) = (tanh(F/gammas) + gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
	575	QFF(F) = max(min(F, 1.), 0.)
	576	QFFf(F) = min(QFF(F), scut)
	577	Qmix1(F) = (tanh((QFF(F) - Fmax)/gammas) + Qcoef1max)/Qcoef2max
	578	Rmix1(F) = (gammaslog(cosh((QFF(F) - Fmax)/gammas)) + QFF(F)Qcoef1max)/Qcoef2max
	579	Qmix2(F) = -log(1.-QFFf(F))/scut
	580	Rmix2(F) = (QFFf(F) + (1.-QFF(F))*log(1.-QFFf(F)))/scut
	581	Qmix(F) = qqa1Qmix1(F) + qqa2Qmix2(F)
	582	Rmix(F) = qqa1Rmix1(F) + qqa2Rmix2(F)
	583
	584	INTEGER, SAVE :: ifrst
	585	DATA ifrst/0/
	586	!$omp THREADPRIVATE(ifrst)
	587
	588	! =====================================================================
	589	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	590	! =====================================================================
	591
	592	! -- Initialize mixing PDF coefficients
	593	IF (ifrst == 0) THEN
	594	ifrst = 1
	595	Qcoef1max = Qcoef1(Fmax)
	596	Qcoef2max = Qcoef2(Fmax)
	597	!<jyg
	598	print *, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
	599	fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
	600	!>jyg
	601	!
	602	END IF
	603
	604	! ori do 360 i=1,ncum*nlp
	605	DO j = 1, nl
	606	DO i = 1, ncum
	607	nent(i, j) = 0
	608	! in convect3, m is computed in cv3_closure
	609	! ori m(i,1)=0.0
	610	END DO
	611	END DO
	612
	613	! ori do 400 k=1,nlp
	614	! ori do 390 j=1,nlp
	615	DO j = 1, nl
	616	DO k = 1, nl
	617	DO i = 1, ncum
	618	Qent(i, k, j) = rr(i, j)
	619	uent(i, k, j) = u(i, j)
	620	vent(i, k, j) = v(i, j)
	621	elij(i, k, j) = 0.0
	622	hent(i, k, j) = 0.0
	623	!AC ! Ment(i,k,j)=0.0
	624	!AC ! Sij(i,k,j)=0.0
	625	END DO
	626	END DO
	627	END DO
	628
	629	!AC !
	630	Ment(1:ncum, 1:nd, 1:nd) = 0.0
	631	Sij(1:ncum, 1:nd, 1:nd) = 0.0
	632	!AC !
	633	!ym
	634	Sigij(1:ncum, 1:nd, 1:nd) = 0.0
	635	!ym
	636
	637	!jyg ! DO k = 1, ntra
	638	!jyg ! DO j = 1, nd ! instead nlp
	639	!jyg ! DO i = 1, nd ! instead nlp
	640	!jyg ! DO il = 1, ncum
	641	!jyg ! traent(il, i, j, k) = tra(il, j, k)
	642	!jyg ! END DO
	643	!jyg ! END DO
	644	!jyg ! END DO
	645	!jyg ! END DO
	646
	647	! =====================================================================
	648	! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
	649	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	650	! --- FRACTION (Sij)
	651	! =====================================================================
	652
	653	DO i = minorig + 1, nl
	654
	655	IF (ok_entrain) THEN
	656	DO j = minorig, nl
	657	DO il = 1, ncum
	658	IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (j >= (icb(il) - 1)) &
	659	.AND. (j <= inb(il))) THEN
	660
	661	! ! rti = qnk(il) - ep(il, i)*clw(il, i)
	662	rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
	663	bf2 = 1.+lv(il, j)lv(il, j)rs(il, j)/(rrvt(il, j)t(il, j)*cpd)
	664	!jyg(from aj)<
	665	IF (cvflag_ice) THEN
	666	! print*,cvflag_ice,'cvflag_ice dans do 700'
	667	IF (t(il, j) <= 263.15) THEN
	668	bf2 = 1.+(lf(il, j) + lv(il, j))(lv(il, j) + frac(il, j) &
	669	lf(il, j))rs(il, j)/(rrvt(il, j)t(il, j)cpd)
	670	END IF
	671	END IF
	672	!>jyg
	673	anum = h(il, j) - hp(il, i) + (cpv - cpd)t(il, j)(rti - rr(il, j))
	674	denom = h(il, i) - hp(il, i) + (cpd - cpv)(rr(il, i) - rti)t(il, j)
	675	dei = denom
	676	IF (abs(dei) < 0.01) dei = 0.01
	677	Sij(il, i, j) = anum/dei
	678	Sij(il, i, i) = 1.0
	679	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il, i, j))rti - rs(il, j)
	680	altem = altem/bf2
	681	cwat = clw(il, j)*(1.-ep(il, j))
	682	stemp = Sij(il, i, j)
	683	IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN
	684	!jyg(from aj)<
	685	IF (cvflag_ice) THEN
	686	anum = anum - (lv(il, j) + frac(il, j)lf(il, j))(rti - rs(il, j) - cwat*bf2)
	687	denom = denom + (lv(il, j) + frac(il, j)lf(il, j))(rr(il, i) - rti)
	688	ELSE
	689	anum = anum - lv(il, j)(rti - rs(il, j) - cwatbf2)
	690	denom = denom + lv(il, j)*(rr(il, i) - rti)
	691	END IF
	692	!>jyg
	693	IF (abs(denom) < 0.01) denom = 0.01
	694	Sij(il, i, j) = anum/denom
	695	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il, i, j))rti - rs(il, j)
	696	altem = altem - (bf2 - 1.)*cwat
	697	END IF
	698	IF (Sij(il, i, j) > 0.0) THEN
	699	! ! ! Ment(il,i,j)=m(il,i)
	700	Ment(il, i, j) = 1.
	701	elij(il, i, j) = altem
	702	elij(il, i, j) = amax1(0.0, elij(il, i, j))
	703	nent(il, i) = nent(il, i) + 1
	704	END IF
	705
	706	Sij(il, i, j) = amax1(0.0, Sij(il, i, j))
	707	Sij(il, i, j) = amin1(1.0, Sij(il, i, j))
	708	ELSE IF (j > i) THEN
	709	IF (prt_level >= 10) THEN
	710	print *, 'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il, i, j)
	711	ENDIF
	712	END IF ! new
	713	END DO
	714	END DO
	715	ELSE ! (ok_entrain)
	716	DO il = 1, ncum
	717	nent(il, i) = 0
	718	ENDDO
	719	ENDIF ! (ok_entrain)
	720
	721	!jygdebug<
	722	IF (prt_level >= 10) THEN
	723	print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout)
	724	IF (nent(igout, i) .gt. 0) THEN
	725	print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout))
	726	ENDIF
	727	ENDIF
	728	!>jygdebug
	729
	730	! * if no air can entrain at level i assume that updraft detrains *
	731	! * at that level and calculate detrained air flux and properties *
	732
	733	! @ do 170 i=icb(il),inb(il)
	734
	735	DO il = 1, ncum
	736	IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (nent(il, i) == 0)) THEN
	737	! @ if(nent(il,i).eq.0)then
	738	! ! ! Ment(il,i,i)=m(il,i)
	739	Ment(il, i, i) = 1.
	740	! ! Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	741	Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
	742	uent(il, i, i) = unk(il)
	743	vent(il, i, i) = vnk(il)
	744	IF (fl_cor_ebil .GE. 2) THEN
	745	hent(il, i, i) = hp(il, i)
	746	ENDIF
	747	elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
	748	Sij(il, i, i) = 0.0
	749	END IF
	750	END DO
	751	END DO ! i = minorig + 1, nl
	752
	753	!jyg ! DO j = 1, ntra
	754	!jyg ! DO i = minorig + 1, nl
	755	!jyg ! DO il = 1, ncum
	756	!jyg ! IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
	757	!jyg ! traent(il, i, i, j) = tra(il, nk(il), j)
	758	!jyg ! END IF
	759	!jyg ! END DO
	760	!jyg ! END DO
	761	!jyg ! END DO
	762
	763	DO j = minorig, nl
	764	DO i = minorig, nl
	765	DO il = 1, ncum
	766	IF ((j >= (icb(il) - 1)) .AND. (j <= inb(il)) .AND. &
	767	(i >= icb(il)) .AND. (i <= inb(il))) THEN
	768	Sigij(il, i, j) = Sij(il, i, j)
	769	END IF
	770	END DO
	771	END DO
	772	END DO
	773	! @ enddo
	774
	775	! @170 continue
	776
	777	! =====================================================================
	778	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
	779	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	780	! =====================================================================
	781
	782	CALL zilch(csum, nloc*nd)
	783
	784	DO il = 1, ncum
	785	lwork(il) = .FALSE.
	786	END DO
	787
	788	! ---------------------------------------------------------------
	789	DO i = minorig + 1, nl !Loop on origin level "i"
	790	! ---------------------------------------------------------------
	791
	792	num1 = 0
	793	DO il = 1, ncum
	794	IF (i >= icb(il) .AND. i <= inb(il)) num1 = num1 + 1
	795	END DO
	796	IF (num1 <= 0) GO TO 789
	797
	798	!JYG1 Find maximum of SIJ for J>I, if any.
	799
	800	Sx(:) = 0.
	801
	802	DO il = 1, ncum
	803	IF (i >= icb(il) .AND. i <= inb(il)) THEN
	804	signhpmh(il) = sign(1., hp(il, i) - h(il, i))
	805	Sbef(il) = max(0., signhpmh(il))
	806	END IF
	807	END DO
	808
	809	DO j = i + 1, nl
	810	DO il = 1, ncum
	811	IF (i >= icb(il) .AND. i <= inb(il) .AND. j <= inb(il)) THEN
	812	IF (Sbef(il) < Sij(il, i, j)) THEN
	813	Sx(il) = max(Sij(il, i, j), Sx(il))
	814	END IF
	815	Sbef(il) = Sij(il, i, j)
	816	END IF
	817	END DO
	818	END DO
	819
	820	DO il = 1, ncum
	821	IF (i >= icb(il) .AND. i <= inb(il)) THEN
	822	lwork(il) = (nent(il, i) /= 0)
	823	! ! rti = qnk(il) - ep(il, i)*clw(il, i)
	824	rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
	825	!jyg<
	826	IF (cvflag_ice) THEN
	827
	828	anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)lf(il, i)) &
	829	(rti - rs(il, i)) + (cpv - cpd)t(il, i)(rti - rr(il, i))
	830	denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)lf(il, i)) &
	831	(rr(il, i) - rti) + (cpd - cpv)t(il, i)(rr(il, i) - rti)
	832	ELSE
	833
	834	anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + &
	835	(cpv - cpd)t(il, i)(rti - rr(il, i))
	836	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + &
	837	(cpd - cpv)t(il, i)(rr(il, i) - rti)
	838	END IF
	839	!>jyg
	840	IF (abs(denom) < 0.01) denom = 0.01
	841	Scrit(il) = min(anum/denom, 1.)
	842	alt = rti - rs(il, i) + Scrit(il)*(rr(il, i) - rti)
	843
	844	!JYG1 Find new critical value Scrit2
	845	! such that : Sij > Scrit2 => mixed draught will detrain at J<I
	846	! Sij < Scrit2 => mixed draught will detrain at J>I
	847
	848	Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
	849	Scrit(il)*max(0., signhpmh(il))
	850
	851	Scrit(il) = Scrit2
	852
	853	!JYG Correction pour la nouvelle logique; la correction pour ALT
	854	! est un peu au hazard
	855	IF (Scrit(il) <= 0.0) Scrit(il) = 0.0
	856	IF (alt <= 0.0) Scrit(il) = 1.0
	857
	858	smax(il) = 0.0
	859	ASij(il) = 0.0
	860	sup(il) = 0. ! upper S-value reached by descending draughts
	861	END IF
	862	END DO
	863
	864	! ---------------------------------------------------------------
	865	DO j = minorig, nl !Loop on destination level "j"
	866	! ---------------------------------------------------------------
	867
	868	num2 = 0
	869	DO il = 1, ncum
	870	IF (i >= icb(il) .AND. i <= inb(il) .AND. &
	871	j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
	872	lwork(il)) num2 = num2 + 1
	873	END DO
	874	IF (num2 <= 0) GO TO 175
	875
	876	! -----------------------------------------------
	877	IF (j > i) THEN
	878	! -----------------------------------------------
	879	DO il = 1, ncum
	880	IF (i >= icb(il) .AND. i <= inb(il) .AND. &
	881	j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
	882	lwork(il)) THEN
	883	IF (Sij(il, i, j) > 0.0) THEN
	884	Smid(il) = min(Sij(il, i, j), Scrit(il))
	885	Sjmax(il) = Smid(il)
	886	Sjmin(il) = Smid(il)
	887	IF (Smid(il) < smin(il) .AND. Sij(il, i, j + 1) < Smid(il)) THEN
	888	smin(il) = Smid(il)
	889	Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il))
	890	Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
	891	Sjmin(il) = min(Sjmin(il), Scrit(il))
	892	Sbef(il) = Sij(il, i, j)
	893	END IF
	894	END IF
	895	END IF
	896	END DO
	897	! -----------------------------------------------
	898	ELSE IF (j == i) THEN
	899	! -----------------------------------------------
	900	DO il = 1, ncum
	901	IF (i >= icb(il) .AND. i <= inb(il) .AND. &
	902	j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
	903	lwork(il)) THEN
	904	IF (Sij(il, i, j) > 0.0) THEN
	905	Smid(il) = 1.
	906	Sjmin(il) = max((Sij(il, i, j - 1) + Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
	907	min((Sij(il, i, j + 1) + Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
	908	Sjmin(il) = max(Sjmin(il), sup(il))
	909	Sjmax(il) = 1.
	910
	911	! - preparation des variables Scrit, Smin et Sbef pour la partie j>i
	912	Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
	913
	914	smin(il) = 1.
	915	Sbef(il) = max(0., signhpmh(il))
	916	supmax(il, i) = sign(Scrit(il), -signhpmh(il))
	917	END IF
	918	END IF
	919	END DO
	920	! -----------------------------------------------
	921	ELSE IF (j < i) THEN
	922	! -----------------------------------------------
	923	DO il = 1, ncum
	924	IF (i >= icb(il) .AND. i <= inb(il) .AND. &
	925	j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
	926	lwork(il)) THEN
	927	IF (Sij(il, i, j) > 0.0) THEN
	928	Smid(il) = max(Sij(il, i, j), Scrit(il))
	929	Sjmax(il) = Smid(il)
	930	Sjmin(il) = Smid(il)
	931	IF (Smid(il) > smax(il) .AND. Sij(il, i, j + 1) > Smid(il)) THEN
	932	smax(il) = Smid(il)
	933	Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j))
	934	Sjmax(il) = max(Sjmax(il), Scrit(il))
	935	Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
	936	Sjmin(il) = max(Sjmin(il), Scrit(il))
	937	Sbef(il) = Sij(il, i, j)
	938	END IF
	939	IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) &
	940	sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
	941	END IF
	942	END IF
	943	END DO
	944	! -----------------------------------------------
	945	END IF
	946	! -----------------------------------------------
	947
	948	DO il = 1, ncum
	949	IF (i >= icb(il) .AND. i <= inb(il) .AND. &
	950	j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
	951	lwork(il)) THEN
	952	IF (Sij(il, i, j) > 0.0) THEN
	953	! ! rti = qnk(il) - ep(il, i)*clw(il, i)
	954	rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
	955	Qmixmax(il) = Qmix(Sjmax(il))
	956	Qmixmin(il) = Qmix(Sjmin(il))
	957	Rmixmax(il) = Rmix(Sjmax(il))
	958	Rmixmin(il) = Rmix(Sjmin(il))
	959	sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
	960	sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
	961
	962	Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j)
	963
	964	! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
	965	IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN
	966	Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il))
	967	ELSE
	968	Sigij(il, i, j) = 0.
	969	END IF
	970
	971	! -- Compute Qent, uent, vent according to the true mixing fraction
	972	Qent(il, i, j) = (1.-Sigij(il, i, j))rti + Sigij(il, i, j)rr(il, i)
	973	uent(il, i, j) = (1.-Sigij(il, i, j))unk(il) + Sigij(il, i, j)u(il, i)
	974	vent(il, i, j) = (1.-Sigij(il, i, j))vnk(il) + Sigij(il, i, j)v(il, i)
	975
	976	! -- Compute liquid water static energy of mixed draughts
	977	! IF (j .GT. i) THEN
	978	! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
	979	! awat=amax1(awat,0.0)
	980	! ELSE
	981	! awat = 0.
	982	! ENDIF
	983	! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
	984	! : + Sigij(il,i,j)*H(il,i)
	985	! : + (LV(il,j)+(cpd-cpv)t(il,j))awat
	986	!IM 301008 beg
	987	hent(il, i, j) = (1.-Sigij(il, i, j))hp(il, i) + Sigij(il, i, j)h(il, i)
	988
	989	!jyg<
	990	! elij(il, i, j) = Qent(il, i, j) - rs(il, j)
	991	! elij(il, i, j) = elij(il, i, j) + &
	992	! ((h(il,j)-hent(il,i,j))rs(il,j)lv(il,j) / &
	993	! ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
	994	! elij(il, i, j) = elij(il, i, j) / &
	995	! (1.+lv(il,j)lv(il,j)rs(il,j) / &
	996	! ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
	997	!
	998	! Computation of condensate amount Elij, taking into account the ice fraction frac
	999	! Warning : the same saturation humidity rs is used over both liquid water and ice; this
	1000	! should be corrected.
	1001	!
	1002	! Heat capacity of mixed draught
	1003	cpm = cpd + Qent(il, i, j)*(cpv - cpd)
	1004	!
	1005	IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
	1006	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
	1007	elij(il, i, j) = elij(il, i, j) + &
	1008	(h(il, j) - hent(il, i, j) + (cpv - cpd)(Qent(il, i, j) - rr(il, j))t(il, j))* &
	1009	rs(il, j)lv(il, j)/(cpmrrvt(il, j)t(il, j))
	1010	elij(il, i, j) = elij(il, i, j)/ &
	1011	(1.+(lv(il, j) + frac(il, j)lf(il, j))lv(il, j)*rs(il, j)/ &
	1012	(cpmrrvt(il, j)*t(il, j)))
	1013	ELSE
	1014	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
	1015	elij(il, i, j) = elij(il, i, j) + &
	1016	(h(il, j) - hent(il, i, j) + (cpv - cpd)(Qent(il, i, j) - rr(il, j))t(il, j))* &
	1017	rs(il, j)lv(il, j)/(cpmrrvt(il, j)t(il, j))
	1018	elij(il, i, j) = elij(il, i, j)/ &
	1019	(1.+lv(il, j)lv(il, j)rs(il, j)/ &
	1020	(cpmrrvt(il, j)*t(il, j)))
	1021	ENDIF
	1022	!>jyg
	1023	elij(il, i, j) = max(elij(il, i, j), 0.)
	1024
	1025	elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j))
	1026
	1027	IF (j > i) THEN
	1028	awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j)
	1029	awat = amax1(awat, 0.0)
	1030	ELSE
	1031	awat = 0.
	1032	END IF
	1033
	1034	! print ,h(il,j)-hent(il,i,j),LV(il,j)rs(il,j)/(cpdrrvt(il,j)*
	1035	! : t(il,j))
	1036
	1037	!jyg<
	1038	! hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)t(il,j))awat
	1039	! Mixed draught temperature at level j
	1040	IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
	1041	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))
	1042	hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)lf(il, j) + (cpd - cpv)Tm)*awat
	1043	ELSE
	1044	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))
	1045	hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)Tm)awat
	1046	ENDIF
	1047	!>jyg
	1048
	1049	!IM 301008 end
	1050
	1051	! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
	1052	! : i,j,hent(il,i,j),Sigij(il,i,j)
	1053
	1054	! -- ASij is the integral of P(F) over the relevant F interval
	1055	ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
	1056	Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
	1057
	1058	END IF
	1059	END IF
	1060	END DO
	1061	!jyg ! DO k = 1, ntra
	1062	!jyg ! DO il = 1, ncum
	1063	!jyg ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
	1064	!jyg ! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
	1065	!jyg ! lwork(il)) THEN
	1066	!jyg ! IF (Sij(il,i,j)>0.0) THEN
	1067	!jyg ! traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
	1068	!jyg ! (1.-Sigij(il,i,j))*tra(il, nk(il), k)
	1069	!jyg ! END IF
	1070	!jyg ! END IF
	1071	!jyg ! END DO
	1072	!jyg ! END DO
	1073
	1074	! -- If I=J (detrainement and entrainement at the same level), then only the
	1075	! -- adiabatic ascent part of the mixture is considered
	1076	IF (i == j) THEN
	1077	DO il = 1, ncum
	1078	IF (i >= icb(il) .AND. i <= inb(il) .AND. &
	1079	j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
	1080	lwork(il)) THEN
	1081	IF (Sij(il, i, j) > 0.0) THEN
	1082	! ! rti = qnk(il) - ep(il, i)*clw(il, i)
	1083	rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
	1084	! ! ! Ment(il,i,i) = m(il,i)abs(Qmixmax(il)(1.-Sjmax(il))
	1085	Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
	1086	Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
	1087	Qent(il, i, i) = rti
	1088	uent(il, i, i) = unk(il)
	1089	vent(il, i, i) = vnk(il)
	1090	hent(il, i, i) = hp(il, i)
	1091	elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
	1092	Sigij(il, i, i) = 0.
	1093	END IF
	1094	END IF
	1095	END DO
	1096	!jyg ! DO k = 1, ntra
	1097	!jyg ! DO il = 1, ncum
	1098	!jyg ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
	1099	!jyg ! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
	1100	!jyg ! lwork(il)) THEN
	1101	!jyg ! IF (Sij(il,i,j)>0.0) THEN
	1102	!jyg ! traent(il, i, i, k) = tra(il, nk(il), k)
	1103	!jyg ! END IF
	1104	!jyg ! END IF
	1105	!jyg ! END DO
	1106	!jyg ! END DO
	1107
	1108	END IF
	1109
	1110	! ---------------------------------------------------------------
	1111	175 END DO ! End loop on destination level "j"
	1112	! ---------------------------------------------------------------
	1113
	1114	DO il = 1, ncum
	1115	IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il)) THEN
	1116	ASij(il) = amax1(1.0E-16, ASij(il))
	1117	!jyg+lluis<
	1118	! ! ASij(il) = 1.0/ASij(il)
	1119	ASij_inv(il) = 1.0/ASij(il)
	1120	! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
	1121	IF (ASij_inv(il) > 100.) ASij_inv(il) = 0.
	1122	!>jyg+lluis
	1123	csum(il, i) = 0.0
	1124	END IF
	1125	END DO
	1126
	1127	DO j = minorig, nl
	1128	DO il = 1, ncum
	1129	IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. &
	1130	j >= (icb(il) - 1) .AND. j <= inb(il)) THEN
	1131	!jyg Ment(il, i, j) = Ment(il, i, j)*ASij(il)
	1132	Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
	1133	END IF
	1134	END DO
	1135	END DO
	1136
	1137	DO j = minorig, nl
	1138	DO il = 1, ncum
	1139	IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. &
	1140	j >= (icb(il) - 1) .AND. j <= inb(il)) THEN
	1141	csum(il, i) = csum(il, i) + Ment(il, i, j)
	1142	END IF
	1143	END DO
	1144	END DO
	1145
	1146	DO il = 1, ncum
	1147	IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. csum(il, i) < 1.) THEN
	1148	! cc : .and. csum(il,i).lt.m(il,i) ) then
	1149	nent(il, i) = 0
	1150	! cc Ment(il,i,i)=m(il,i)
	1151	Ment(il, i, i) = 1.
	1152	! ! Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
	1153	Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
	1154	uent(il, i, i) = unk(il)
	1155	vent(il, i, i) = vnk(il)
	1156	elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
	1157	IF (fl_cor_ebil .GE. 2) THEN
	1158	hent(il, i, i) = hp(il, i)
	1159	Sigij(il, i, i) = 0.0
	1160	ELSE
	1161	Sij(il, i, i) = 0.0
	1162	ENDIF
	1163	END IF
	1164	END DO ! il
	1165
	1166	!jyg ! DO j = 1, ntra
	1167	!jyg ! DO il = 1, ncum
	1168	!jyg ! IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
	1169	!jyg ! ! cc : .and. csum(il,i).lt.m(il,i) ) then
	1170	!jyg ! traent(il, i, i, j) = tra(il, nk(il), j)
	1171	!jyg ! END IF
	1172	!jyg ! END DO
	1173	!jyg ! END DO
	1174
	1175	! ---------------------------------------------------------------
	1176	789 END DO ! End loop on origin level "i"
	1177	! ---------------------------------------------------------------
	1178
	1179	RETURN
	1180	END SUBROUTINE
	1181	end module

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

Download in other formats:

Original Format