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 : * ! ************************************************************** IMPLICIT NONE include "cvthermo.h" include "cv3param.h" include "YOMCST2.h" ! inputs: INTEGER ncum, nd, na, ntra, nloc INTEGER icb(nloc), inb(nloc), nk(nloc) REAL sig(nloc, nd) REAL qnk(nloc), unk(nloc), vnk(nloc) REAL ph(nloc, nd+1) REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd) REAL u(nloc, nd), v(nloc, nd) REAL tra(nloc, nd, ntra) ! input of convect3 REAL lv(nloc, na) REAL h(nloc, na) !liquid water static energy of environment REAL hp(nloc, na) !liquid water static energy of air shed from adiab. asc. REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na) ! outputs: REAL ment(nloc, na, na), qent(nloc, na, na) REAL uent(nloc, na, na), vent(nloc, na, na) REAL sigij(nloc, na, na), elij(nloc, na, na) REAL supmax(nloc, na) ! Highest mixing fraction of mixed updraughts ! with the sign of (h-hp) REAL traent(nloc, nd, nd, ntra) REAL ments(nloc, nd, nd), qents(nloc, nd, nd) REAL hent(nloc, nd, nd) INTEGER nent(nloc, nd) ! 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 qmixmax(nloc), rmixmax(nloc), sqmrmax(nloc) REAL qmixmin(nloc), rmixmin(nloc), sqmrmin(nloc) REAL signhpmh(nloc) REAL sx(nloc), scrit2 REAL smid(nloc), sjmin(nloc), sjmax(nloc) REAL sbef(nloc), sup(nloc), smin(nloc) REAL asij(nloc), smax(nloc), scrit(nloc) REAL sij(nloc, nd, nd) REAL csum(nloc, nd) REAL awat LOGICAL lwork(nloc) REAL amxupcrit, df, ff INTEGER nstep ! -- 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) 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! DO k = 1, ntra DO j = 1, nd ! instead nlp DO i = 1, nd ! instead nlp DO il = 1, ncum traent(il, i, j, k) = tra(il, j, k) END DO END DO END DO 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 ! cc 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 ! *** 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 ! 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 END DO DO j = 1, ntra DO i = minorig + 1, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN traent(il, i, i, j) = tra(il, nk(il), j) END IF END DO END DO 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) qp = qnk(il) - ep(il, i)*clw(il, i) anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + & (cpv-cpd)*t(il, i)*(qp-rr(il,i)) denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + & (cpd-cpv)*t(il, i)*(rr(il,i)-qp) IF (abs(denom)<0.01) denom = 0.01 scrit(il) = min(anum/denom, 1.) alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp) ! 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 DO k = 1, ntra 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 traent(il, i, j, k) = sigij(il, i, j)*tra(il, i, k) + & (1.-sigij(il,i,j))*tra(il, nk(il), k) END IF END IF END DO 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) ! cc 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 DO k = 1, ntra 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 traent(il, i, i, k) = tra(il, nk(il), k) END IF END IF END DO END DO END IF 175 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(il) = 1.0/asij(il) 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(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 DO j = 1, ntra 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 traent(il, i, i, j) = tra(il, nk(il), j) END IF END DO END DO 789 END DO RETURN END SUBROUTINE cv3p_mixing