SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, & ph, t, rr, rs, u, v, tra, h, lv, qnk, & 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 IMPLICIT NONE include "cvthermo.h" include "cv3param.h" include "YOMCST2.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) :: qnk, unk, vnk 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) :: 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, qp 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 !jyg REAL, DIMENSION (nloc) :: ASij, smax, Scrit REAL, DIMENSION (nloc) :: ASij, ASij_inv, smax, Scrit REAL, DIMENSION (nloc, nd, nd) :: Sij REAL, DIMENSION (nloc, nd) :: csum REAL :: awat 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! !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 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) bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd) 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 anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2) denom = denom + lv(il, j)*(rr(il,i)-rti) 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)) END IF ! new END DO END DO !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) uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) elij(il, i, i) = clw(il, i)*(1.-ep(il,i)) Sij(il, i, i) = 0.0 END IF END DO END DO !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) 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) 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) 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) 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))) 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)) hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat !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) !!! 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) uent(il, i, i) = unk(il) vent(il, i, i) = vnk(il) elij(il, i, i) = clw(il, i)*(1.-ep(il,i)) Sij(il, i, i) = 0.0 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