Index: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 =================================================================== --- LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 (revision 3709) +++ LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 (revision 3710) @@ -106,6 +106,6 @@ REAL, DIMENSION(nloc, nd, nd), INTENT(OUT) :: hent INTEGER, DIMENSION(nloc, nd), INTENT(OUT) :: nent - - !local variables: + +!local variables: INTEGER i, j, k, il, im, jm REAL Smid @@ -126,15 +126,15 @@ 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 - + !$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))) @@ -147,12 +147,12 @@ 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 + !$omp THREADPRIVATE(ifrst) + + ! ===================================================================== + ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS + ! ===================================================================== + ! -- Initialize mixing PDF coefficients IF (ifrst == 0) THEN ifrst = 1 @@ -162,9 +162,9 @@ 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 @@ -176,15 +176,15 @@ END DO END DO - + 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 - - ! ===================================================================== - ! --- 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 IF (ok_entrain) THEN @@ -230,5 +230,5 @@ 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)) @@ -245,5 +245,5 @@ ENDDO ENDIF ! (ok_entrain) - + IF (prt_level >= 10) THEN print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout) @@ -252,8 +252,8 @@ ENDIF 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 IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (nent(il, i) == 0)) THEN @@ -270,5 +270,5 @@ END DO END DO ! i = minorig + 1, nl - + DO j = minorig, nl DO i = minorig, nl @@ -281,26 +281,26 @@ END DO END DO - - ! ===================================================================== - ! --- NORMALIZE ENTRAINED AIR MASS FLUXES - ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING - ! ===================================================================== - + + ! ===================================================================== + ! --- 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" - ! --------------------------------------------------------------- - + ! --------------------------------------------------------------- + DO il = 1, ncum IF (i >= icb(il) .AND. i <= inb(il)) THEN signhpmh(il) = sign(1., hp(il, i) - h(il, i)) - - Sx(il) = max( maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il) ) - + + Sx(il) = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il)) + lwork(il) = (nent(il, i) /= 0) rti = qta(il, i - 1) - ep(il, i)*clw(il, i) @@ -319,5 +319,5 @@ 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 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 = 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 - 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 = 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)) - 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 = 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 - END IF - END DO - ! ----------------------------------------------- - END IF - ! ----------------------------------------------- - - DO il = 1, ncum + + 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)) THEN - IF (Sij(il, i, j) > 0.0) THEN + 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 + ! ----------------------------------------------- + + 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) - 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)) - + 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 - - ! -- 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 = 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 - END IF - - ! --------------------------------------------------------------- - 175 END DO ! End loop on destination level "j" - ! --------------------------------------------------------------- - + 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 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 @@ -536,5 +491,5 @@ END DO END DO - + DO j = minorig, nl DO il = 1, ncum @@ -545,5 +500,5 @@ 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 @@ -562,8 +517,8 @@ END IF END DO ! il - - ! --------------------------------------------------------------- - 789 END DO ! End loop on origin level "i" - ! --------------------------------------------------------------- + + ! --------------------------------------------------------------- + END DO ! End loop on origin level "i" + ! --------------------------------------------------------------- RETURN END SUBROUTINE