module cv3p_mixing_mod implicit none contains ! ************************************************************** ! * ! CV3P_MIXING : compute mixed draught properties and, * ! within a scaling factor, mixed draught * ! mass fluxes. * ! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15* ! modified by : * ! ************************************************************** SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, & ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, & unk, vnk, hp, tv, tvp, ep, clw, sig, & Ment, Qent, hent, uent, vent, nent, & Sigij, elij, supmax, Ments, Qents, traent) !inputs: INTEGER, INTENT(IN) :: ncum, nd, na INTEGER, INTENT(IN) :: ntra, nloc INTEGER, DIMENSION(nloc), INTENT(IN) :: icb, inb, nk REAL, DIMENSION(nloc, nd), INTENT(IN) :: sig REAL, DIMENSION(nloc), INTENT(IN) :: unk, vnk REAL, DIMENSION(nloc, nd), INTENT(IN) :: qta REAL, DIMENSION(nloc, nd + 1), INTENT(IN) :: ph REAL, DIMENSION(nloc, nd), INTENT(IN) :: t, rr, rs REAL, DIMENSION(nloc, nd), INTENT(IN) :: u, v REAL, DIMENSION(nloc, nd, ntra), INTENT(IN) :: tra ! input of convect3 REAL, DIMENSION(nloc, na), INTENT(IN) :: lv REAL, DIMENSION(nloc, na), INTENT(IN) :: lf REAL, DIMENSION(nloc, na), INTENT(IN) :: frac !ice fraction in condensate REAL, DIMENSION(nloc, na), INTENT(IN) :: h !liquid water static energy of environment REAL, DIMENSION(nloc, na), INTENT(IN) :: hp !liquid water static energy of air shed from adiab. asc. REAL, DIMENSION(nloc, na), INTENT(IN) :: tv, tvp REAL, DIMENSION(nloc, na), INTENT(IN) :: ep, clw !outputs: REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Ment, Qent REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: uent, vent REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Sigij, elij REAL, DIMENSION(nloc, na), INTENT(OUT) :: supmax ! Highest mixing fraction of mixed ! updraughts with the sign of (h-hp) REAL, DIMENSION(nloc, nd, nd, ntra), INTENT(OUT) :: traent REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: Ments, Qents REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: hent INTEGER, DIMENSION(nloc, nd), INTENT(OUT) :: nent logical, parameter :: use_old = .false. if (use_old) then call cv3p_mixing_old(nloc, ncum, nd, na, ntra, icb, nk, inb, & ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, & unk, vnk, hp, tv, tvp, ep, clw, sig, & Ment, Qent, hent, uent, vent, nent, & Sigij, elij, supmax, Ments, Qents, traent) else call cv3p_mixing_new(nloc, ncum, nd, na, icb, inb, & t, rr, rs, u, v, h, lv, lf, frac, qta, & unk, vnk, hp, ep, clw, & Ment, Qent, hent, uent, vent, nent, & Sigij, elij, supmax) endif END SUBROUTINE ! modified by : Arnaud Durocher (04/2020) : changed loops order SUBROUTINE cv3p_mixing_new(nloc, ncum, nd, na, icb, inb, & t, rr, rs, u, v, h, lv, lf, frac, qta, & unk, vnk, hp, ep, clw, & Ment, Qent, hent, uent, vent, nent, & Sigij, elij, supmax) USE print_control_mod, ONLY: prt_level USE add_phys_tend_mod, ONLY: fl_cor_ebil IMPLICIT NONE include "cvthermo.h" include "cv3param.h" include "YOMCST2.h" include "cvflag.h" !inputs: INTEGER, INTENT(IN) :: ncum, nd, na INTEGER, INTENT(IN) :: nloc INTEGER, DIMENSION(nloc), INTENT(IN) :: icb, inb REAL, DIMENSION(nloc), INTENT(IN) :: unk, vnk REAL, DIMENSION(nloc, nd), INTENT(IN) :: qta REAL, DIMENSION(nloc, nd), INTENT(IN) :: t, rr, rs REAL, DIMENSION(nloc, nd), INTENT(IN) :: u, v REAL, DIMENSION(nloc, na), INTENT(IN) :: lv REAL, DIMENSION(nloc, na), INTENT(IN) :: lf REAL, DIMENSION(nloc, na), INTENT(IN) :: frac !ice fraction in condensate REAL, DIMENSION(nloc, na), INTENT(IN) :: h !liquid water static energy of environment REAL, DIMENSION(nloc, na), INTENT(IN) :: hp !liquid water static energy of air shed from adiab. asc. REAL, DIMENSION(nloc, na), INTENT(IN) :: ep, clw !outputs: REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Ment, Qent REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: uent, vent REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Sigij, elij REAL, DIMENSION(nloc, na), INTENT(OUT) :: supmax ! Highest mixing fraction of mixed ! updraughts with the sign of (h-hp) REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: hent INTEGER, DIMENSION(nloc, nd), INTENT(OUT) :: nent !local variables: INTEGER i, j, k, il REAL, DIMENSION(nloc, nd) :: ASij REAL, DIMENSION(nloc, nd, nd) :: Sij real :: Sjmin(nd), Sjmax(nd), Qmixmax(nd), Qmixmin(nd), Rmixmax(nd), Rmixmin(nd) INTEGER, SAVE :: igout = 1 !$omp THREADPRIVATE(igout) ! -- Mixing probability distribution functions REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F Qcoef1(F) = tanh(F/gammas) Qcoef2(F) = (tanh(F/gammas) + gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas))) QFF(F) = max(min(F, 1.), 0.) QFFf(F) = min(QFF(F), scut) Qmix1(F) = (tanh((QFF(F) - Fmax)/gammas) + Qcoef1max)/Qcoef2max Rmix1(F) = (gammas*log(cosh((QFF(F) - Fmax)/gammas)) + QFF(F)*Qcoef1max)/Qcoef2max Qmix2(F) = -log(1.-QFFf(F))/scut Rmix2(F) = (QFFf(F) + (1.-QFF(F))*log(1.-QFFf(F)))/scut Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F) Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F) INTEGER, SAVE :: ifrst = 0 !$omp THREADPRIVATE(ifrst) ! ===================================================================== ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS ! ===================================================================== ! -- Initialize mixing PDF coefficients IF (ifrst == 0) THEN ifrst = 1 Qcoef1max = Qcoef1(Fmax) Qcoef2max = Qcoef2(Fmax) print *, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', & fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max END IF nent(:ncum, :nl) = 0 elij(:ncum, :nl, :nl) = 0 hent(:ncum, :nl, :nl) = 0 DO j = 1, nl DO k = 1, nl DO i = 1, ncum Qent(i, k, j) = rr(i, j) uent(i, k, j) = u(i, j) vent(i, k, j) = v(i, j) END DO END DO END DO Ment(1:ncum, 1:nd, 1:nd) = 0.0 Sij(1:ncum, 1:nd, 1:nd) = 0.0 ! ===================================================================== ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING ! --- FRACTION (Sij) ! ===================================================================== IF (ok_entrain) THEN ! NOTE : (04/2020 AD) this loop order gave good results in cv3a_tracer ! if( icb(il) <= i <= inb(il) ) is probably more expensive than reading non-contiguously DO il = 1, ncum DO i = icb(il), inb(il) DO j = (icb(il) - 1), inb(il) block real :: rti, bf2, anum, denom, dei, altem, cwat, stemp rti = qta(il, i - 1) - ep(il, i)*clw(il, i) bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) IF (cvflag_ice .AND. t(il, j) <= 263.15) THEN bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* & lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) END IF anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j)) denom = h(il, i) - hp(il, i) + (cpd - cpv)*(rr(il, i) - rti)*t(il, j) dei = denom IF (abs(dei) < 0.01) dei = 0.01 Sij(il, i, j) = anum/dei Sij(il, i, i) = 1.0 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) altem = altem/bf2 cwat = clw(il, j)*(1.-ep(il, j)) stemp = Sij(il, i, j) IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN IF (cvflag_ice) THEN anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2) denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti) ELSE anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2) denom = denom + lv(il, j)*(rr(il, i) - rti) END IF IF (abs(denom) < 0.01) denom = 0.01 Sij(il, i, j) = anum/denom altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) altem = altem - (bf2 - 1.)*cwat END IF IF (Sij(il, i, j) > 0.0) THEN Ment(il, i, j) = 1. elij(il, i, j) = altem elij(il, i, j) = amax1(0.0, elij(il, i, j)) nent(il, i) = nent(il, i) + 1 END IF Sij(il, i, j) = amax1(0.0, Sij(il, i, j)) Sij(il, i, j) = amin1(1.0, Sij(il, i, j)) end block END DO IF (nent(il, i) == 0) THEN Ment(il, i, i) = 1. Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) IF (fl_cor_ebil .GE. 2) THEN hent(il, i, i) = hp(il, i) ENDIF elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) Sij(il, i, i) = 0.0 END IF END DO END DO ELSE DO i = minorig + 1, nl DO il = 1, ncum IF ((i >= icb(il)) .AND. (i <= inb(il))) THEN Ment(il, i, i) = 1. Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) IF (fl_cor_ebil .GE. 2) THEN hent(il, i, i) = hp(il, i) ENDIF elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) Sij(il, i, i) = 0.0 END IF END DO END DO ENDIF DO i = minorig + 1, nl DO j = minorig, nl IF (prt_level >= 10) THEN print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout) !IF (nent(igout, i) .gt. 0) THEN ! print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout)) !ENDIF ENDIF END DO END DO Sigij(1:ncum, 1:nd, 1:nd) = 0.0 DO il = 1, ncum DO i = icb(il), inb(il) DO j = (icb(il) - 1), inb(il) Sigij(il, i, j) = Sij(il, i, j) END DO END DO END DO ! ===================================================================== ! --- NORMALIZE ENTRAINED AIR MASS FLUXES ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING ! ===================================================================== DO il = 1, ncum DO i = icb(il), inb(il) !Loop on origin level "i" block real :: signhpmh, Scrit block real :: rti, anum, denom, alt, Sx signhpmh = sign(1., hp(il, i) - h(il, i)) rti = qta(il, i - 1) - ep(il, i)*clw(il, i) IF (cvflag_ice) THEN anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* & (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i)) denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* & (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti) ELSE anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + & (cpv - cpd)*t(il, i)*(rti - rr(il, i)) denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + & (cpd - cpv)*t(il, i)*(rr(il, i) - rti) END IF IF (abs(denom) < 0.01) denom = 0.01 Scrit = min(anum/denom, 1.) alt = rti - rs(il, i) + Scrit*(rr(il, i) - rti) if (alt <= 0) then Scrit = 1.0 else ! Get max of Sij Sx = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh) ! Find new critical value Scrit ! such that : Sij > Scrit => mixed draught will detrain at J mixed draught will detrain at J>I Scrit = max(0., min(Scrit, Sx*max(0., -signhpmh)) + Scrit*max(0., signhpmh)) endif end block ASij(il, i) = 0.0 ! Compute all Sjmax(j) and Sjmin(j) block real :: smin, smax, sup, Sbef, Smid smin = 1 smax = 0.0 sup = 0. ! upper S-value reached by descending draughts ! Glitchy : why so complicated? if (i < inb(il)) then Sbef = Sij(il, i, inb(il)) else Sbef = max(0., signhpmh) endif DO j = (icb(il) - 1), i - 1 IF (Sij(il, i, j) > 0.0) THEN Smid = max(Sij(il, i, j), Scrit) Sjmax(j) = Smid Sjmin(j) = Smid IF (Smid > smax .AND. Sij(il, i, j + 1) > Smid) THEN smax = Smid Sjmax(j) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) Sjmax(j) = max(Sjmax(j), Scrit) Sjmin(j) = min((Sbef + Sij(il, i, j))/2., Sij(il, i, j)) Sjmin(j) = max(Sjmin(j), Scrit) Sbef = Sij(il, i, j) END IF IF (abs(Sjmin(j) - Sjmax(j)) > 1.E-10) & sup = max(Sjmin(j), Sjmax(j), sup) ENDIF END DO j = i IF (Sij(il, i, j) > 0.0) THEN Smid = 1. Sjmin(j) = max((Sij(il, i, j - 1) + Smid)/2., Scrit)*max(0., -signhpmh) + & min((Sij(il, i, j + 1) + Smid)/2., Scrit)*max(0., signhpmh) Sjmin(j) = max(Sjmin(j), sup) Sjmax(j) = 1. ! preparation des variables Scrit, Smin et Sbef pour la partie j>i Scrit = min(Sjmin(j), Sjmax(j), Scrit) Sbef = max(0., signhpmh) supmax(il, i) = sign(Scrit, -signhpmh) ENDIF DO j = i + 1, inb(il) IF (Sij(il, i, j) > 0.0) THEN Smid = min(Sij(il, i, j), Scrit) Sjmax(j) = Smid Sjmin(j) = Smid IF (Smid < smin .AND. Sij(il, i, j + 1) < Smid) THEN smin = min(smin, Smid) Sjmax(j) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit) Sjmin(j) = max((Sbef + Sij(il, i, j))/2., Sij(il, i, j)) Sjmin(j) = min(Sjmin(j), Scrit) Sbef = Sij(il, i, j) END IF ENDIF END DO end block DO j = (icb(il) - 1), inb(il) Qmixmax(j) = Qmix(Sjmax(j)) Qmixmin(j) = Qmix(Sjmin(j)) Rmixmax(j) = Rmix(Sjmax(j)) Rmixmin(j) = Rmix(Sjmin(j)) END DO DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" IF (Sij(il, i, j) > 0.0) THEN block real :: sqmrmax, sqmrmin, cpm, awat, Tm, rti rti = qta(il, i - 1) - ep(il, i)*clw(il, i) sqmrmax = Sjmax(j)*Qmixmax(j) - Rmixmax(j) sqmrmin = Sjmin(j)*Qmixmin(j) - Rmixmin(j) Ment(il, i, j) = abs(Qmixmax(j) - Qmixmin(j))*Ment(il, i, j) IF (abs(Qmixmax(j) - Qmixmin(j)) > 1.E-10) THEN Sigij(il, i, j) = (sqmrmax - sqmrmin)/(Qmixmax(j) - Qmixmin(j)) ELSE Sigij(il, i, j) = 0. END IF ! Compute Qent, uent, vent according to the true mixing fraction Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i) uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i) vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i) ! Compute liquid water static energy of mixed draughts hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i) ! Heat capacity of mixed draught cpm = cpd + Qent(il, i, j)*(cpv - cpd) elij(il, i, j) = Qent(il, i, j) - rs(il, j) elij(il, i, j) = elij(il, i, j) + & (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN elij(il, i, j) = elij(il, i, j)/ & (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ & (cpm*rrv*t(il, j)*t(il, j))) ELSE elij(il, i, j) = elij(il, i, j)/ & (1.+lv(il, j)*lv(il, j)*rs(il, j)/ & (cpm*rrv*t(il, j)*t(il, j))) ENDIF elij(il, i, j) = max(elij(il, i, j), 0.) elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j)) IF (j > i) THEN awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j) awat = amax1(awat, 0.0) ELSE awat = 0. END IF ! Mixed draught temperature at level j Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j)) IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat ELSE hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat ENDIF ! ASij is the integral of P(F) over the relevant F interval ASij(il, i) = ASij(il, i) + abs(Qmixmax(j)*(1.-Sjmax(j)) + Rmixmax(j) - & Qmixmin(j)*(1.-Sjmin(j)) - Rmixmin(j)) ! If I=J (detrainement and entrainement at the same level), then only the ! adiabatic ascent part of the mixture is considered IF (i == j) THEN rti = qta(il, i - 1) - ep(il, i)*clw(il, i) Ment(il, i, i) = abs(Qmixmax(j)*(1.-Sjmax(j)) + Rmixmax(j) - & Qmixmin(j)*(1.-Sjmin(j)) - Rmixmin(j)) Qent(il, i, i) = rti uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) hent(il, i, i) = hp(il, i) elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) Sigij(il, i, i) = 0. END IF end block END IF END DO ! Loop j = (icb(il) - 1), inb(il) end block END DO ! Loop i = icb(il), inb(il) END DO ! Loop il = 1, ncum DO il = 1, ncum DO i = icb(il), inb(il) !Loop on origin level "i" block real :: csum, ASij_inv csum = 0 DO j = (icb(il) - 1), inb(il) ASij(il, i) = amax1(1.0E-16, ASij(il, i)) ASij_inv = 1.0/ASij(il, i) ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs IF (ASij_inv > 100.) ASij_inv = 0. Ment(il, i, j) = Ment(il, i, j)*ASij_inv csum = csum + Ment(il, i, j) END DO IF (csum < 1.) THEN nent(il, i) = 0 Ment(il, i, i) = 1. Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) IF (fl_cor_ebil .GE. 2) THEN hent(il, i, i) = hp(il, i) Sigij(il, i, i) = 0.0 ELSE Sij(il, i, i) = 0.0 ENDIF ENDIF end block END DO ! Loop il = 1, ncum END DO ! End loop on origin level "i" END SUBROUTINE ! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15 SUBROUTINE cv3p_mixing_old(nloc, ncum, nd, na, ntra, icb, nk, inb, & ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, & unk, vnk, hp, tv, tvp, ep, clw, sig, & Ment, Qent, hent, uent, vent, nent, & Sigij, elij, supmax, Ments, Qents, traent) USE print_control_mod, ONLY: prt_level USE add_phys_tend_mod, ONLY: fl_cor_ebil IMPLICIT NONE include "cvthermo.h" include "cv3param.h" include "YOMCST2.h" include "cvflag.h" !inputs: INTEGER, INTENT(IN) :: ncum, nd, na INTEGER, INTENT(IN) :: ntra, nloc INTEGER, DIMENSION(nloc), INTENT(IN) :: icb, inb, nk REAL, DIMENSION(nloc, nd), INTENT(IN) :: sig REAL, DIMENSION(nloc), INTENT(IN) :: unk, vnk REAL, DIMENSION(nloc, nd), INTENT(IN) :: qta REAL, DIMENSION(nloc, nd + 1), INTENT(IN) :: ph REAL, DIMENSION(nloc, nd), INTENT(IN) :: t, rr, rs REAL, DIMENSION(nloc, nd), INTENT(IN) :: u, v REAL, DIMENSION(nloc, nd, ntra), INTENT(IN) :: tra ! input of convect3 REAL, DIMENSION(nloc, na), INTENT(IN) :: lv REAL, DIMENSION(nloc, na), INTENT(IN) :: lf REAL, DIMENSION(nloc, na), INTENT(IN) :: frac !ice fraction in condensate REAL, DIMENSION(nloc, na), INTENT(IN) :: h !liquid water static energy of environment REAL, DIMENSION(nloc, na), INTENT(IN) :: hp !liquid water static energy of air shed from adiab. asc. REAL, DIMENSION(nloc, na), INTENT(IN) :: tv, tvp REAL, DIMENSION(nloc, na), INTENT(IN) :: ep, clw !outputs: REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Ment, Qent REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: uent, vent REAL, DIMENSION(nloc, na, na), INTENT(OUT) :: Sigij, elij REAL, DIMENSION(nloc, na), INTENT(OUT) :: supmax ! Highest mixing fraction of mixed ! updraughts with the sign of (h-hp) REAL, DIMENSION(nloc, nd, nd, ntra), INTENT(OUT) :: traent REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: Ments, Qents REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: hent INTEGER, DIMENSION(nloc, nd), INTENT(OUT) :: nent !local variables: INTEGER i, j, k, il, im, jm INTEGER num1, num2 REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp REAL :: alt, delp, delm REAL, DIMENSION(nloc) :: Qmixmax, Rmixmax, sqmrmax REAL, DIMENSION(nloc) :: Qmixmin, Rmixmin, sqmrmin REAL, DIMENSION(nloc) :: signhpmh REAL, DIMENSION(nloc) :: Sx REAL :: Scrit2 REAL, DIMENSION(nloc) :: Smid, Sjmin, Sjmax REAL, DIMENSION(nloc) :: Sbef, sup, smin REAL, DIMENSION(nloc) :: ASij, ASij_inv, smax, Scrit REAL, DIMENSION(nloc, nd, nd) :: Sij REAL, DIMENSION(nloc, nd) :: csum REAL :: awat REAL :: cpm !Mixed draught heat capacity REAL :: Tm !Mixed draught temperature LOGICAL, DIMENSION(nloc) :: lwork REAL amxupcrit, df, ff INTEGER nstep INTEGER, SAVE :: igout = 1 !$omp THREADPRIVATE(igout) ! -- Mixing probability distribution functions REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F Qcoef1(F) = tanh(F/gammas) Qcoef2(F) = (tanh(F/gammas) + gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas))) QFF(F) = max(min(F, 1.), 0.) QFFf(F) = min(QFF(F), scut) Qmix1(F) = (tanh((QFF(F) - Fmax)/gammas) + Qcoef1max)/Qcoef2max Rmix1(F) = (gammas*log(cosh((QFF(F) - Fmax)/gammas)) + QFF(F)*Qcoef1max)/Qcoef2max Qmix2(F) = -log(1.-QFFf(F))/scut Rmix2(F) = (QFFf(F) + (1.-QFF(F))*log(1.-QFFf(F)))/scut Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F) Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F) INTEGER, SAVE :: ifrst DATA ifrst/0/ !$omp THREADPRIVATE(ifrst) ! ===================================================================== ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS ! ===================================================================== ! -- Initialize mixing PDF coefficients IF (ifrst == 0) THEN ifrst = 1 Qcoef1max = Qcoef1(Fmax) Qcoef2max = Qcoef2(Fmax) !jyg ! END IF ! ori do 360 i=1,ncum*nlp DO j = 1, nl DO i = 1, ncum nent(i, j) = 0 ! in convect3, m is computed in cv3_closure ! ori m(i,1)=0.0 END DO END DO ! ori do 400 k=1,nlp ! ori do 390 j=1,nlp DO j = 1, nl DO k = 1, nl DO i = 1, ncum Qent(i, k, j) = rr(i, j) uent(i, k, j) = u(i, j) vent(i, k, j) = v(i, j) elij(i, k, j) = 0.0 hent(i, k, j) = 0.0 !AC ! Ment(i,k,j)=0.0 !AC ! Sij(i,k,j)=0.0 END DO END DO END DO !AC ! Ment(1:ncum, 1:nd, 1:nd) = 0.0 Sij(1:ncum, 1:nd, 1:nd) = 0.0 !AC ! !ym Sigij(1:ncum, 1:nd, 1:nd) = 0.0 !ym !jyg ! DO k = 1, ntra !jyg ! DO j = 1, nd ! instead nlp !jyg ! DO i = 1, nd ! instead nlp !jyg ! DO il = 1, ncum !jyg ! traent(il, i, j, k) = tra(il, j, k) !jyg ! END DO !jyg ! END DO !jyg ! END DO !jyg ! END DO ! ===================================================================== ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING ! --- FRACTION (Sij) ! ===================================================================== DO i = minorig + 1, nl IF (ok_entrain) THEN DO j = minorig, nl DO il = 1, ncum IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (j >= (icb(il) - 1)) & .AND. (j <= inb(il))) THEN ! ! rti = qnk(il) - ep(il, i)*clw(il, i) rti = qta(il, i - 1) - ep(il, i)*clw(il, i) bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) !jyg(from aj)< IF (cvflag_ice) THEN ! print*,cvflag_ice,'cvflag_ice dans do 700' IF (t(il, j) <= 263.15) THEN bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* & lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd) END IF END IF !>jyg anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j)) denom = h(il, i) - hp(il, i) + (cpd - cpv)*(rr(il, i) - rti)*t(il, j) dei = denom IF (abs(dei) < 0.01) dei = 0.01 Sij(il, i, j) = anum/dei Sij(il, i, i) = 1.0 altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) altem = altem/bf2 cwat = clw(il, j)*(1.-ep(il, j)) stemp = Sij(il, i, j) IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN !jyg(from aj)< IF (cvflag_ice) THEN anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2) denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti) ELSE anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2) denom = denom + lv(il, j)*(rr(il, i) - rti) END IF !>jyg IF (abs(denom) < 0.01) denom = 0.01 Sij(il, i, j) = anum/denom altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j) altem = altem - (bf2 - 1.)*cwat END IF IF (Sij(il, i, j) > 0.0) THEN ! ! ! Ment(il,i,j)=m(il,i) Ment(il, i, j) = 1. elij(il, i, j) = altem elij(il, i, j) = amax1(0.0, elij(il, i, j)) nent(il, i) = nent(il, i) + 1 END IF Sij(il, i, j) = amax1(0.0, Sij(il, i, j)) Sij(il, i, j) = amin1(1.0, Sij(il, i, j)) ELSE IF (j > i) THEN IF (prt_level >= 10) THEN print *, 'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il, i, j) ENDIF END IF ! new END DO END DO ELSE ! (ok_entrain) DO il = 1, ncum nent(il, i) = 0 ENDDO ENDIF ! (ok_entrain) !jygdebug< IF (prt_level >= 10) THEN print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout) IF (nent(igout, i) .gt. 0) THEN print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout)) ENDIF ENDIF !>jygdebug ! *** if no air can entrain at level i assume that updraft detrains *** ! *** at that level and calculate detrained air flux and properties *** ! @ do 170 i=icb(il),inb(il) DO il = 1, ncum IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (nent(il, i) == 0)) THEN ! @ if(nent(il,i).eq.0)then ! ! ! Ment(il,i,i)=m(il,i) Ment(il, i, i) = 1. ! ! Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i) Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) IF (fl_cor_ebil .GE. 2) THEN hent(il, i, i) = hp(il, i) ENDIF elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) Sij(il, i, i) = 0.0 END IF END DO END DO ! i = minorig + 1, nl !jyg ! DO j = 1, ntra !jyg ! DO i = minorig + 1, nl !jyg ! DO il = 1, ncum !jyg ! IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN !jyg ! traent(il, i, i, j) = tra(il, nk(il), j) !jyg ! END IF !jyg ! END DO !jyg ! END DO !jyg ! END DO DO j = minorig, nl DO i = minorig, nl DO il = 1, ncum IF ((j >= (icb(il) - 1)) .AND. (j <= inb(il)) .AND. & (i >= icb(il)) .AND. (i <= inb(il))) THEN Sigij(il, i, j) = Sij(il, i, j) END IF END DO END DO END DO ! @ enddo ! @170 continue ! ===================================================================== ! --- NORMALIZE ENTRAINED AIR MASS FLUXES ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING ! ===================================================================== CALL zilch(csum, nloc*nd) DO il = 1, ncum lwork(il) = .FALSE. END DO ! --------------------------------------------------------------- DO i = minorig + 1, nl !Loop on origin level "i" ! --------------------------------------------------------------- num1 = 0 DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il)) num1 = num1 + 1 END DO IF (num1 <= 0) GO TO 789 !JYG1 Find maximum of SIJ for J>I, if any. Sx(:) = 0. DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il)) THEN signhpmh(il) = sign(1., hp(il, i) - h(il, i)) Sbef(il) = max(0., signhpmh(il)) END IF END DO DO j = i + 1, nl DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. j <= inb(il)) THEN IF (Sbef(il) < Sij(il, i, j)) THEN Sx(il) = max(Sij(il, i, j), Sx(il)) END IF Sbef(il) = Sij(il, i, j) END IF END DO END DO DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il)) THEN lwork(il) = (nent(il, i) /= 0) ! ! rti = qnk(il) - ep(il, i)*clw(il, i) rti = qta(il, i - 1) - ep(il, i)*clw(il, i) !jyg< IF (cvflag_ice) THEN anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* & (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i)) denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* & (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti) ELSE anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + & (cpv - cpd)*t(il, i)*(rti - rr(il, i)) denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + & (cpd - cpv)*t(il, i)*(rr(il, i) - rti) END IF !>jyg IF (abs(denom) < 0.01) denom = 0.01 Scrit(il) = min(anum/denom, 1.) alt = rti - rs(il, i) + Scrit(il)*(rr(il, i) - rti) !JYG1 Find new critical value Scrit2 ! such that : Sij > Scrit2 => mixed draught will detrain at J mixed draught will detrain at J>I Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + & Scrit(il)*max(0., signhpmh(il)) Scrit(il) = Scrit2 !JYG Correction pour la nouvelle logique; la correction pour ALT ! est un peu au hazard IF (Scrit(il) <= 0.0) Scrit(il) = 0.0 IF (alt <= 0.0) Scrit(il) = 1.0 smax(il) = 0.0 ASij(il) = 0.0 sup(il) = 0. ! upper S-value reached by descending draughts END IF END DO ! --------------------------------------------------------------- DO j = minorig, nl !Loop on destination level "j" ! --------------------------------------------------------------- num2 = 0 DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. & j >= (icb(il) - 1) .AND. j <= inb(il) .AND. & lwork(il)) num2 = num2 + 1 END DO IF (num2 <= 0) GO TO 175 ! ----------------------------------------------- IF (j > i) THEN ! ----------------------------------------------- DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. & j >= (icb(il) - 1) .AND. j <= inb(il) .AND. & lwork(il)) THEN IF (Sij(il, i, j) > 0.0) THEN Smid(il) = min(Sij(il, i, j), Scrit(il)) Sjmax(il) = Smid(il) Sjmin(il) = Smid(il) IF (Smid(il) < smin(il) .AND. Sij(il, i, j + 1) < Smid(il)) THEN smin(il) = Smid(il) Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il)) Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) Sjmin(il) = min(Sjmin(il), Scrit(il)) Sbef(il) = Sij(il, i, j) END IF END IF END IF END DO ! ----------------------------------------------- ELSE IF (j == i) THEN ! ----------------------------------------------- DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. & j >= (icb(il) - 1) .AND. j <= inb(il) .AND. & lwork(il)) THEN IF (Sij(il, i, j) > 0.0) THEN Smid(il) = 1. Sjmin(il) = max((Sij(il, i, j - 1) + Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + & min((Sij(il, i, j + 1) + Smid(il))/2., Scrit(il))*max(0., signhpmh(il)) Sjmin(il) = max(Sjmin(il), sup(il)) Sjmax(il) = 1. ! - preparation des variables Scrit, Smin et Sbef pour la partie j>i Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il)) smin(il) = 1. Sbef(il) = max(0., signhpmh(il)) supmax(il, i) = sign(Scrit(il), -signhpmh(il)) END IF END IF END DO ! ----------------------------------------------- ELSE IF (j < i) THEN ! ----------------------------------------------- DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. & j >= (icb(il) - 1) .AND. j <= inb(il) .AND. & lwork(il)) THEN IF (Sij(il, i, j) > 0.0) THEN Smid(il) = max(Sij(il, i, j), Scrit(il)) Sjmax(il) = Smid(il) Sjmin(il) = Smid(il) IF (Smid(il) > smax(il) .AND. Sij(il, i, j + 1) > Smid(il)) THEN smax(il) = Smid(il) Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) Sjmax(il) = max(Sjmax(il), Scrit(il)) Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j)) Sjmin(il) = max(Sjmin(il), Scrit(il)) Sbef(il) = Sij(il, i, j) END IF IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) & sup(il) = max(Sjmin(il), Sjmax(il), sup(il)) END IF END IF END DO ! ----------------------------------------------- END IF ! ----------------------------------------------- DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. & j >= (icb(il) - 1) .AND. j <= inb(il) .AND. & lwork(il)) THEN IF (Sij(il, i, j) > 0.0) THEN ! ! rti = qnk(il) - ep(il, i)*clw(il, i) rti = qta(il, i - 1) - ep(il, i)*clw(il, i) Qmixmax(il) = Qmix(Sjmax(il)) Qmixmin(il) = Qmix(Sjmin(il)) Rmixmax(il) = Rmix(Sjmax(il)) Rmixmin(il) = Rmix(Sjmin(il)) sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il)) sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il)) Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j) ! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j) IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il)) ELSE Sigij(il, i, j) = 0. END IF ! -- Compute Qent, uent, vent according to the true mixing fraction Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i) uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i) vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i) ! -- Compute liquid water static energy of mixed draughts ! IF (j .GT. i) THEN ! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j) ! awat=amax1(awat,0.0) ! ELSE ! awat = 0. ! ENDIF ! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i) ! : + Sigij(il,i,j)*H(il,i) ! : + (LV(il,j)+(cpd-cpv)*t(il,j))*awat !IM 301008 beg hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i) !jyg< ! elij(il, i, j) = Qent(il, i, j) - rs(il, j) ! elij(il, i, j) = elij(il, i, j) + & ! ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / & ! ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j))) ! elij(il, i, j) = elij(il, i, j) / & ! (1.+lv(il,j)*lv(il,j)*rs(il,j) / & ! ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j))) ! ! Computation of condensate amount Elij, taking into account the ice fraction frac ! Warning : the same saturation humidity rs is used over both liquid water and ice; this ! should be corrected. ! ! Heat capacity of mixed draught cpm = cpd + Qent(il, i, j)*(cpv - cpd) ! IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN elij(il, i, j) = Qent(il, i, j) - rs(il, j) elij(il, i, j) = elij(il, i, j) + & (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) elij(il, i, j) = elij(il, i, j)/ & (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ & (cpm*rrv*t(il, j)*t(il, j))) ELSE elij(il, i, j) = Qent(il, i, j) - rs(il, j) elij(il, i, j) = elij(il, i, j) + & (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* & rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j)) elij(il, i, j) = elij(il, i, j)/ & (1.+lv(il, j)*lv(il, j)*rs(il, j)/ & (cpm*rrv*t(il, j)*t(il, j))) ENDIF !>jyg elij(il, i, j) = max(elij(il, i, j), 0.) elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j)) IF (j > i) THEN awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j) awat = amax1(awat, 0.0) ELSE awat = 0. END IF ! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)* ! : t(il,j)) !jyg< ! hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat ! Mixed draught temperature at level j IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j)) hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat ELSE Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j)) hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat ENDIF !>jyg !IM 301008 end ! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ', ! : i,j,hent(il,i,j),Sigij(il,i,j) ! -- ASij is the integral of P(F) over the relevant F interval ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) END IF END IF END DO !jyg ! DO k = 1, ntra !jyg ! DO il = 1, ncum !jyg ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. & !jyg ! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. & !jyg ! lwork(il)) THEN !jyg ! IF (Sij(il,i,j)>0.0) THEN !jyg ! traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + & !jyg ! (1.-Sigij(il,i,j))*tra(il, nk(il), k) !jyg ! END IF !jyg ! END IF !jyg ! END DO !jyg ! END DO ! -- If I=J (detrainement and entrainement at the same level), then only the ! -- adiabatic ascent part of the mixture is considered IF (i == j) THEN DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. & j >= (icb(il) - 1) .AND. j <= inb(il) .AND. & lwork(il)) THEN IF (Sij(il, i, j) > 0.0) THEN ! ! rti = qnk(il) - ep(il, i)*clw(il, i) rti = qta(il, i - 1) - ep(il, i)*clw(il, i) ! ! ! Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il)) Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - & Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) Qent(il, i, i) = rti uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) hent(il, i, i) = hp(il, i) elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) Sigij(il, i, i) = 0. END IF END IF END DO !jyg ! DO k = 1, ntra !jyg ! DO il = 1, ncum !jyg ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. & !jyg ! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. & !jyg ! lwork(il)) THEN !jyg ! IF (Sij(il,i,j)>0.0) THEN !jyg ! traent(il, i, i, k) = tra(il, nk(il), k) !jyg ! END IF !jyg ! END IF !jyg ! END DO !jyg ! END DO END IF ! --------------------------------------------------------------- 175 END DO ! End loop on destination level "j" ! --------------------------------------------------------------- DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il)) THEN ASij(il) = amax1(1.0E-16, ASij(il)) !jyg+lluis< ! ! ASij(il) = 1.0/ASij(il) ASij_inv(il) = 1.0/ASij(il) ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs IF (ASij_inv(il) > 100.) ASij_inv(il) = 0. !>jyg+lluis csum(il, i) = 0.0 END IF END DO DO j = minorig, nl DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. & j >= (icb(il) - 1) .AND. j <= inb(il)) THEN !jyg Ment(il, i, j) = Ment(il, i, j)*ASij(il) Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il) END IF END DO END DO DO j = minorig, nl DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. & j >= (icb(il) - 1) .AND. j <= inb(il)) THEN csum(il, i) = csum(il, i) + Ment(il, i, j) END IF END DO END DO DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. csum(il, i) < 1.) THEN ! cc : .and. csum(il,i).lt.m(il,i) ) then nent(il, i) = 0 ! cc Ment(il,i,i)=m(il,i) Ment(il, i, i) = 1. ! ! Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i) Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i) uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) elij(il, i, i) = clw(il, i)*(1.-ep(il, i)) IF (fl_cor_ebil .GE. 2) THEN hent(il, i, i) = hp(il, i) Sigij(il, i, i) = 0.0 ELSE Sij(il, i, i) = 0.0 ENDIF END IF END DO ! il !jyg ! DO j = 1, ntra !jyg ! DO il = 1, ncum !jyg ! IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN !jyg ! ! cc : .and. csum(il,i).lt.m(il,i) ) then !jyg ! traent(il, i, i, j) = tra(il, nk(il), j) !jyg ! END IF !jyg ! END DO !jyg ! END DO ! --------------------------------------------------------------- 789 END DO ! End loop on origin level "i" ! --------------------------------------------------------------- RETURN END SUBROUTINE end module