Index: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 =================================================================== --- LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 (revision 3710) +++ LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 (revision 3711) @@ -131,7 +131,7 @@ INTEGER, SAVE :: igout = 1 - !$omp THREADPRIVATE(igout) - - ! -- Mixing probability distribution functions +!$omp THREADPRIVATE(igout) + +! -- Mixing probability distribution functions REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F @@ -149,10 +149,10 @@ INTEGER, SAVE :: ifrst = 0 - !$omp THREADPRIVATE(ifrst) - - ! ===================================================================== - ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS - ! ===================================================================== - ! -- Initialize mixing PDF coefficients +!$omp THREADPRIVATE(ifrst) + +! ===================================================================== +! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS +! ===================================================================== +! -- Initialize mixing PDF coefficients IF (ifrst == 0) THEN ifrst = 1 @@ -181,9 +181,9 @@ Sigij(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) - ! ===================================================================== +! ===================================================================== +! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING +! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING +! --- FRACTION (Sij) +! ===================================================================== DO i = minorig + 1, nl @@ -253,6 +253,6 @@ ENDIF - ! *** if no air can entrain at level i assume that updraft detrains *** - ! *** at that level and calculate detrained air flux and properties *** +! *** if no air can entrain at level i assume that updraft detrains *** +! *** at that level and calculate detrained air flux and properties *** DO il = 1, ncum @@ -282,10 +282,8 @@ END DO - ! ===================================================================== - ! --- NORMALIZE ENTRAINED AIR MASS FLUXES - ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING - ! ===================================================================== - - CALL zilch(csum, nloc*nd) +! ===================================================================== +! --- NORMALIZE ENTRAINED AIR MASS FLUXES +! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING +! ===================================================================== DO il = 1, ncum @@ -293,7 +291,7 @@ END DO - ! --------------------------------------------------------------- +! --------------------------------------------------------------- DO i = minorig + 1, nl !Loop on origin level "i" - ! --------------------------------------------------------------- +! --------------------------------------------------------------- DO il = 1, ncum @@ -341,185 +339,156 @@ Sbef(il) = max(0., signhpmh(il)) endif - END IF - END DO - - DO il = 1, ncum - DO j = minorig, nl !Loop on destination level "j" - IF (i >= icb(il) .AND. i <= inb(il) .AND. & - j >= (icb(il) - 1) .AND. j <= inb(il) .AND. & - lwork(il) .AND. 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) + + DO j = (icb(il) - 1), inb(il) !Loop on destination level "j" + IF (lwork(il) .AND. 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)) 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) + ! ----------------------------------------------- + + 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 - IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) & - sup(il) = max(Sjmin(il), Sjmax(il), sup(il)) + + ! 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) = ASij(il) + 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 - ! ----------------------------------------------- - - 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 - 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 - - ! 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)) - - ! 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 + END DO ! Loop j = (icb(il) - 1), inb(il) + + IF (lwork(il)) THEN + csum(il, i) = 0 + DO j = (icb(il) - 1), inb(il) + ASij(il) = amax1(1.0E-16, 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. + 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) - 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 ! End loop on destination level "j" - ! --------------------------------------------------------------- - END DO - - DO il = 1, ncum - IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il)) THEN - ASij(il) = amax1(1.0E-16, 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. - 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 - 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 - 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 - END IF - END DO ! il - - ! --------------------------------------------------------------- - END DO ! End loop on origin level "i" - ! --------------------------------------------------------------- - RETURN + 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 IF + END IF ! i >= icb(il) .AND. i <= inb(il) + END DO ! Loop il = 1, ncum + END DO ! End loop on origin level "i" END SUBROUTINE