Index: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.F90 =================================================================== --- LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.F90 (revision 3708) +++ (revision ) @@ -1,701 +1,0 @@ -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) -! ************************************************************** -! * -! 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 : * -! ************************************************************** - - USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level - USE ioipsl_getin_p_mod, ONLY: getin_p - 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)=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)=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=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 cv3p_mixing - 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 3709) @@ -0,0 +1,1256 @@ +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, im, jm + REAL Smid + 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) :: 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 = 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 + 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) + ! ===================================================================== + + 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 = 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) THEN + 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 + 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)) + 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) + + 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 + + ! *** 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 + 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 ! i = minorig + 1, nl + + 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 + + ! ===================================================================== + ! --- 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) ) + + lwork(il) = (nent(il, i) /= 0) + 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) = 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 + END IF + END DO + + ! --------------------------------------------------------------- + DO j = minorig, nl !Loop on destination level "j" + ! --------------------------------------------------------------- + ! ----------------------------------------------- + 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 = 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 + 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) + 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)) + + 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" + ! --------------------------------------------------------------- + + 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 + + ! --------------------------------------------------------------- + 789 END DO ! End loop on origin level "i" + ! --------------------------------------------------------------- + RETURN + 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 Index: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cva_driver.F90 =================================================================== --- LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cva_driver.F90 (revision 3708) +++ LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cva_driver.F90 (revision 3709) @@ -42,4 +42,5 @@ USE add_phys_tend_mod, ONLY: fl_cor_ebil USE cv3a_compress_mod + USE cv3p_mixing_mod IMPLICIT NONE