Index: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 =================================================================== --- LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 (revision 3714) +++ LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 (revision 3715) @@ -112,15 +112,6 @@ 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, smin - REAL :: Scrit2 - REAL, DIMENSION(nloc) :: Sjmin, Sjmax - REAL, DIMENSION(nloc) :: Sbef, sup - REAL, DIMENSION(nloc, nd) :: ASij - REAL, DIMENSION(nloc) :: ASij_inv, smax, Scrit - REAL, DIMENSION(nloc, nd, nd) :: Sij - REAL, DIMENSION(nloc, nd) :: csum + REAL, DIMENSION(nloc, nd) :: ASij + REAL, DIMENSION(nloc, nd, nd) :: Sij REAL :: awat REAL :: cpm !Mixed draught heat capacity @@ -179,5 +170,4 @@ Ment(1:ncum, 1:nd, 1:nd) = 0.0 Sij(1:ncum, 1:nd, 1:nd) = 0.0 - Sigij(1:ncum, 1:nd, 1:nd) = 0.0 ! ===================================================================== @@ -287,171 +277,174 @@ ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING ! ===================================================================== - DO il = 1, ncum DO i = icb(il), inb(il) !Loop on origin level "i" - signhpmh(il) = sign(1., hp(il, i) - h(il, i)) - - Sx(il) = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il)) - - 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(il) = min(anum/denom, 1.) - alt = rti - rs(il, i) + Scrit(il)*(rr(il, i) - rti) - - ! 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 - - ! 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, i) = 0.0 - sup(il) = 0. ! upper S-value reached by descending draughts - - if (i < inb(il)) then - Sbef(il) = Sij(il, i, inb(il)) - else - Sbef(il) = max(0., signhpmh(il)) - endif - - DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" - IF (Sij(il, i, j) > 0.0) THEN - ! ----------------------------------------------- - IF (j > i) THEN - Smid = min(Sij(il, i, j), Scrit(il)) - Sjmax(il) = Smid - Sjmin(il) = Smid - IF (Smid < smin(il) .AND. Sij(il, i, j + 1) < Smid) THEN - smin(il) = min(smin(il), Smid) - 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 - ELSE IF (j == i) THEN - Smid = 1. - Sjmin(il) = max((Sij(il, i, j - 1) + Smid)/2., Scrit(il))*max(0., -signhpmh(il)) + & - min((Sij(il, i, j + 1) + Smid)/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)) - ELSE IF (j < i) THEN - Smid = max(Sij(il, i, j), Scrit(il)) - Sjmax(il) = Smid - Sjmin(il) = Smid - IF (Smid > smax(il) .AND. Sij(il, i, j + 1) > Smid) THEN - smax(il) = Smid - 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)) + block + real :: signhpmh, Sx, Scrit + real :: smax, sup, Sbef, smin + + 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) + + ! 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 = min(Scrit, Sx*max(0., -signhpmh)) + Scrit*max(0., signhpmh) + + ! Correction pour la nouvelle logique; la correction pour ALT + ! est un peu au hazard + IF (Scrit <= 0.0) Scrit = 0.0 + IF (alt <= 0.0) Scrit = 1.0 + + smax = 0.0 + ASij(il, i) = 0.0 + sup = 0. ! upper S-value reached by descending draughts + + ! Glitchy : why? + if (i < inb(il)) then + Sbef = Sij(il, i, inb(il)) + else + Sbef = max(0., signhpmh) + endif + + DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" + IF (Sij(il, i, j) > 0.0) THEN + block + real :: Sjmax, Sjmin, Qmixmax, Qmixmin, Rmixmax, Rmixmin, sqmrmax, sqmrmin + ! ----------------------------------------------- + IF (j > i) THEN + Smid = min(Sij(il, i, j), Scrit) + Sjmax = Smid + Sjmin = Smid + IF (Smid < smin .AND. Sij(il, i, j + 1) < Smid) THEN + smin = min(smin, Smid) + Sjmax = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit) + Sjmin = max((Sbef + Sij(il, i, j))/2., Sij(il, i, j)) + Sjmin = min(Sjmin, Scrit) + Sbef = Sij(il, i, j) + END IF + ELSE IF (j == i) THEN + Smid = 1. + Sjmin = 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 = max(Sjmin, sup) + Sjmax = 1. + ! preparation des variables Scrit, Smin et Sbef pour la partie j>i + Scrit = min(Sjmin, Sjmax, Scrit) + + smin = 1. + Sbef = max(0., signhpmh) + supmax(il, i) = sign(Scrit, -signhpmh) + ELSE IF (j < i) THEN + Smid = max(Sij(il, i, j), Scrit) + Sjmax = Smid + Sjmin = Smid + IF (Smid > smax .AND. Sij(il, i, j + 1) > Smid) THEN + smax = Smid + Sjmax = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j)) + Sjmax = max(Sjmax, Scrit) + Sjmin = min((Sbef + Sij(il, i, j))/2., Sij(il, i, j)) + Sjmin = max(Sjmin, Scrit) + Sbef = Sij(il, i, j) + END IF + IF (abs(Sjmin - Sjmax) > 1.E-10) & + sup = max(Sjmin, Sjmax, sup) + END IF + ! ----------------------------------------------- + + rti = qta(il, i - 1) - ep(il, i)*clw(il, i) + Qmixmax = Qmix(Sjmax) + Qmixmin = Qmix(Sjmin) + Rmixmax = Rmix(Sjmax) + Rmixmin = Rmix(Sjmin) + sqmrmax = Sjmax*Qmix(Sjmax) - Rmix(Sjmax) + sqmrmin = Sjmin*Qmix(Sjmin) - Rmix(Sjmin) + + Ment(il, i, j) = abs(Qmixmax - Qmixmin)*Ment(il, i, j) + IF (abs(Qmixmax - Qmixmin) > 1.E-10) THEN + Sigij(il, i, j) = (sqmrmax - sqmrmin)/(Qmixmax - Qmixmin) + 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*(1.-Sjmax) + Rmixmax - & + Qmixmin*(1.-Sjmin) - Rmixmin) + + ! 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*(1.-Sjmax) + Rmixmax - & + Qmixmin*(1.-Sjmin) - Rmixmin) + 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 - ! ----------------------------------------------- - - 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) - 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 - 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) - 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 - 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(il)*(1.-Sjmax(il)) + Rmixmax(il) - & - Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il)) - - ! 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(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 ! Loop j = (icb(il) - 1), inb(il) + 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 @@ -459,27 +452,30 @@ DO il = 1, ncum DO i = icb(il), inb(il) !Loop on origin level "i" - csum(il, i) = 0 - DO j = (icb(il) - 1), inb(il) - ASij(il, i) = amax1(1.0E-16, ASij(il, i)) - ASij_inv(il) = 1.0/ASij(il, i) - ! 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. - Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il) - csum(il, i) = csum(il, i) + Ment(il, i, j) - END DO - IF (csum(il, i) < 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 + 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 - ENDIF + end block END DO ! Loop il = 1, ncum END DO ! End loop on origin level "i"