! $Id: cv_routines.f90 5894 2025-11-28 16:34:54Z snguyen $ MODULE cv_routines_mod PRIVATE PUBLIC cv_param, cv_prelim, cv_feed, cv_undilute1, cv_trigger, cv_compress, cv_undilute2, & cv_closure, cv_mixing, cv_yield, cv_unsat, cv_uncompress CONTAINS SUBROUTINE cv_param(nd) USE lmdz_cv_ini, ONLY : alpha,betad,coeffr,coeffs,cu,damp,delta,dtmax,elcrit,entp,minorig,nl,nlm,nlp,noff,omtrain,omtsnow,sigd,sigs,tlcrit IMPLICIT NONE ! ------------------------------------------------------------ ! Set parameters for convectL ! (includes microphysical parameters and parameters that ! control the rate of approach to quasi-equilibrium) ! ------------------------------------------------------------ ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *** ! *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *** ! *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *** ! *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *** ! *** BETWEEN 0 C AND TLCRIT) *** ! *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *** ! *** FORMULATION *** ! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** ! *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** ! *** OF CLOUD *** ! *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *** ! *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *** ! *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** ! *** OF RAIN *** ! *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** ! *** OF SNOW *** ! *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *** ! *** TRANSPORT *** ! *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *** ! *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *** ! *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *** ! *** APPROACH TO QUASI-EQUILIBRIUM *** ! *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *** ! *** (DAMP MUST BE LESS THAN 1) *** INTEGER nd CHARACTER (LEN=20) :: modname = 'cv_routines' CHARACTER (LEN=80) :: abort_message ! noff: integer limit for convection (nd-noff) ! minorig: First level of convection noff = 2 minorig = 2 nl = nd - noff nlp = nl + 1 nlm = nl - 1 elcrit = 0.0011 tlcrit = -55.0 entp = 1.5 sigs = 0.12 sigd = 0.05 omtrain = 50.0 omtsnow = 5.5 coeffr = 1.0 coeffs = 0.8 dtmax = 0.9 cu = 0.70 betad = 10.0 damp = 0.1 alpha = 0.2 delta = 0.01 ! cld RETURN END SUBROUTINE cv_param SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm) USE lmdz_cv_ini, ONLY : cl,clmcpv,cpd,cpv,epsim1,hrd,lv0,nlp,t0 IMPLICIT NONE ! ===================================================================== ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY ! ===================================================================== ! inputs: INTEGER len, nd, ndp1 REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1) ! outputs: REAL lv(len, nd), cpn(len, nd), tv(len, nd) REAL gz(len, nd), h(len, nd), hm(len, nd) ! local variables: INTEGER k, i REAL cpx(len, nd) DO k = 1, nlp DO i = 1, len lv(i, k) = lv0 - clmcpv*(t(i,k)-t0) cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k) cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k) tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1) ENDDO ENDDO ! gz = phi at the full levels (same as p). DO i = 1, len gz(i, 1) = 0.0 ENDDO DO k = 2, nlp DO i = 1, len gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, & k) ENDDO ENDDO ! h = phi + cpT (dry static energy). ! hm = phi + cp(T-Tbase)+Lq DO k = 1, nlp DO i = 1, len h(i, k) = gz(i, k) + cpn(i, k)*t(i, k) hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k) ENDDO ENDDO RETURN END SUBROUTINE cv_prelim SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, & qnk, gznk, plcl) USE lmdz_cv_ini, ONLY : minorig,nl,nlm,nlp IMPLICIT NONE ! ================================================================ ! Purpose: CONVECTIVE FEED ! ================================================================ ! inputs: INTEGER len, nd REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd) REAL hm(len, nd), gz(len, nd) ! outputs: INTEGER iflag(len), nk(len), icb(len), icbmax REAL tnk(len), qnk(len), gznk(len), plcl(len) ! local variables: INTEGER i, k INTEGER ihmin(len) REAL work(len) REAL pnk(len), qsnk(len), rh(len), chi(len) ! ------------------------------------------------------------------- ! --- Find level of minimum moist static energy ! --- If level of minimum moist static energy coincides with ! --- or is lower than minimum allowable parcel origin level, ! --- set iflag to 6. ! ------------------------------------------------------------------- DO i = 1, len work(i) = 1.0E12 ihmin(i) = nl ENDDO DO k = 2, nlp DO i = 1, len IF ((hm(i,k)work(i)) .AND. (k<=ihmin(i))) THEN work(i) = hm(i, k) nk(i) = k ENDIF ENDDO ENDDO ! ------------------------------------------------------------------- ! --- Check whether parcel level temperature and specific humidity ! --- are reasonable ! ------------------------------------------------------------------- DO i = 1, len IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< & 400.0)) .AND. (iflag(i)==0)) iflag(i) = 7 ENDDO ! ------------------------------------------------------------------- ! --- Calculate lifted condensation level of air at parcel origin level ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980) ! ------------------------------------------------------------------- DO i = 1, len tnk(i) = t(i, nk(i)) qnk(i) = q(i, nk(i)) gznk(i) = gz(i, nk(i)) pnk(i) = p(i, nk(i)) qsnk(i) = qs(i, nk(i)) rh(i) = qnk(i)/qsnk(i) rh(i) = min(1.0, rh(i)) chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i)) plcl(i) = pnk(i)*(rh(i)**chi(i)) IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i & ) = 8 ENDDO ! ------------------------------------------------------------------- ! --- Calculate first level above lcl (=icb) ! ------------------------------------------------------------------- DO i = 1, len icb(i) = nlm ENDDO DO k = minorig, nl DO i = 1, len IF ((k>=(nk(i)+1)) .AND. (p(i,k)=nlm) .AND. (iflag(i)==0)) iflag(i) = 9 ENDDO ! Compute icbmax. !ym do not do that, independance between column !ym icbmax = 2 !ym DO i = 1, len !ym icbmax = max(icbmax, icb(i)) !ym END DO RETURN END SUBROUTINE cv_feed SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, & clw) USE lmdz_cv_ini, ONLY : cl,clmcpv,cpd,cpv,eps,epsi,lv0,minorig,rrv,t0,nl IMPLICIT NONE ! inputs: INTEGER len, nd INTEGER nk(len), icb(len), icbmax REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd) REAL p(len, nd) ! outputs: REAL tp(len, nd), tvp(len, nd), clw(len, nd) ! local variables: INTEGER i, k REAL tg, qg, alv, s, ahg, tc, denom, es, rg REAL ah0(len), cpp(len) REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len) ! ------------------------------------------------------------------- ! --- Calculates the lifted parcel virtual temperature at nk, ! --- the actual temperature, and the adiabatic ! --- liquid water content. The procedure is to solve the equation. ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. ! ------------------------------------------------------------------- DO i = 1, len tnk(i) = t(i, nk(i)) qnk(i) = q(i, nk(i)) gznk(i) = gz(i, nk(i)) ticb(i) = t(i, icb(i)) gzicb(i) = gz(i, icb(i)) END DO ! *** Calculate certain parcel quantities, including static energy *** DO i = 1, len ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- & 273.15)) + gznk(i) cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv END DO ! *** Calculate lifted parcel quantities below cloud base *** !ym bad dependance between column => icbmax computed in cv_feed !ym DO k = minorig, icbmax - 1 DO k = minorig, nd DO i = 1, len IF (k <= MAX(2,icb(i))-1) THEN tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i) tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi) ENDIF ENDDO ENDDO ! *** Find lifted parcel quantities above cloud base *** DO i = 1, len tg = ticb(i) qg = qs(i, icb(i)) alv = lv0 - clmcpv*(ticb(i)-t0) ! First iteration. s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i)) s = 1./s ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i) tg = tg + s*(ah0(i)-ahg) tg = max(tg, 35.0) tc = tg - t0 denom = 243.5 + tc IF (tc>=0.0) THEN es = 6.112*exp(17.67*tc/denom) ELSE es = exp(23.33086-6111.72784/tg+0.15215*log(tg)) ENDIF qg = eps*es/(p(i,icb(i))-es*(1.-eps)) ! Second iteration. s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i)) s = 1./s ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i) tg = tg + s*(ah0(i)-ahg) tg = max(tg, 35.0) tc = tg - t0 denom = 243.5 + tc IF (tc>=0.0) THEN es = 6.112*exp(17.67*tc/denom) ELSE es = exp(23.33086-6111.72784/tg+0.15215*log(tg)) ENDIF qg = eps*es/(p(i,icb(i))-es*(1.-eps)) alv = lv0 - clmcpv*(ticb(i)-273.15) tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd clw(i, icb(i)) = qnk(i) - qg clw(i, icb(i)) = max(0.0, clw(i,icb(i))) rg = qg/(1.-qnk(i)) tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi) ENDDO !ym bad dependance between column => ibmax computed in cv_feed !ym DO k = minorig, icbmax DO k = minorig, nd DO i = 1, len IF (k <= MAX(2,icb(i))) THEN tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i) ENDIF ENDDO ENDDO RETURN END SUBROUTINE cv_undilute1 SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag) USE lmdz_cv_ini, ONLY : dtmax IMPLICIT NONE ! ------------------------------------------------------------------- ! --- Test for instability. ! --- If there was no convection at last time step and parcel ! --- is stable at icb, then set iflag to 4. ! ------------------------------------------------------------------- ! inputs: INTEGER len, nd, icb(len) REAL cbmf(len), tv(len, nd), tvp(len, nd) ! outputs: INTEGER iflag(len) ! also an input ! local variables: INTEGER i DO i = 1, len IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, & icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4 END DO RETURN END SUBROUTINE cv_trigger SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, compress, nk1, icb1, cbmf1, plcl1, & tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, & tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, & v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph) USE lmdz_cv_ini, ONLY : nl USE print_control_mod, ONLY: lunout IMPLICIT NONE ! inputs: INTEGER len, ncum, nd, nloc INTEGER iflag1(len), nk1(len), icb1(len) LOGICAL compress REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len) REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd) REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd) REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd) REAL tvp1(len, nd), clw1(len, nd) ! outputs: INTEGER iflag(nloc), nk(nloc), icb(nloc) REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc) REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd) REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd) REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd) REAL tvp(nloc, nd), clw(nloc, nd) REAL dph(nloc, nd) ! local variables: INTEGER i, k, nn CHARACTER (LEN=20) :: modname = 'cv_compress' CHARACTER (LEN=80) :: abort_message IF (compress) THEN DO k = 1, nl + 1 nn = 0 DO i = 1, len IF (iflag1(i)==0) THEN nn = nn + 1 t(nn, k) = t1(i, k) q(nn, k) = q1(i, k) qs(nn, k) = qs1(i, k) u(nn, k) = u1(i, k) v(nn, k) = v1(i, k) gz(nn, k) = gz1(i, k) h(nn, k) = h1(i, k) lv(nn, k) = lv1(i, k) cpn(nn, k) = cpn1(i, k) p(nn, k) = p1(i, k) ph(nn, k) = ph1(i, k) tv(nn, k) = tv1(i, k) tp(nn, k) = tp1(i, k) tvp(nn, k) = tvp1(i, k) clw(nn, k) = clw1(i, k) ENDIF ENDDO ENDDO IF (nn/=ncum) THEN WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum abort_message = '' CALL abort_physic(modname, abort_message, 1) ENDIF nn = 0 DO i = 1, len IF (iflag1(i)==0) THEN nn = nn + 1 cbmf(nn) = cbmf1(i) plcl(nn) = plcl1(i) tnk(nn) = tnk1(i) qnk(nn) = qnk1(i) gznk(nn) = gznk1(i) nk(nn) = nk1(i) icb(nn) = icb1(i) iflag(nn) = iflag1(i) ENDIF ENDDO ELSE !compress t(:, 1:nl+1) = t1(:, 1:nl+1) q(:, 1:nl+1) = q1(:, 1:nl+1) qs(:, 1:nl+1) = qs1(:, 1:nl+1) u(:, 1:nl+1) = u1(:, 1:nl+1) v(:, 1:nl+1) = v1(:, 1:nl+1) gz(:, 1:nl+1) = gz1(:, 1:nl+1) h(:, 1:nl+1) = h1(:, 1:nl+1) lv(:, 1:nl+1) = lv1(:, 1:nl+1) cpn(:, 1:nl+1) = cpn1(:, 1:nl+1) p(:, 1:nl+1) = p1(:, 1:nl+1) ph(:, 1:nl+1) = ph1(:, 1:nl+1) tv(:, 1:nl+1) = tv1(:, 1:nl+1) tp(:, 1:nl+1) = tp1(:, 1:nl+1) tvp(:, 1:nl+1) = tvp1(:, 1:nl+1) clw(:, 1:nl+1) = clw1(:, 1:nl+1) cbmf(:) = cbmf1(:) plcl(:) = plcl1(:) tnk(:) = tnk1(:) qnk(:) = qnk1(:) gznk(:) = gznk1(:) nk(:) = nk1(:) icb(:) = icb1(:) iflag(:) = iflag1(:) ENDIF DO k = 1, nl DO i = 1, ncum dph(i, k) = ph(i, k) - ph(i, k+1) ENDDO ENDDO RETURN END SUBROUTINE cv_compress SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, & gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac) USE lmdz_cv_ini, ONLY : nl,cl,clmcpv,cpd,cpv,elcrit,eps,epsi,lv0,minorig,nlp,rrv,sigs,t0,tlcrit IMPLICIT NONE ! --------------------------------------------------------------------- ! Purpose: ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES ! & ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD ! & ! FIND THE LEVEL OF NEUTRAL BUOYANCY ! --------------------------------------------------------------------- ! inputs: INTEGER ncum, nd, nloc INTEGER icb(nloc), nk(nloc) REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd) REAL p(nloc, nd), dph(nloc, nd) REAL tnk(nloc), qnk(nloc), gznk(nloc) REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd) ! outputs: INTEGER inb(nloc), inb1(nloc) REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd) REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd) REAL frac(nloc) ! local variables: INTEGER i, k REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit REAL by, defrac REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc) LOGICAL lcape(nloc) ! ===================================================================== ! --- SOME INITIALIZATIONS ! ===================================================================== DO k = 1, nl DO i = 1, ncum ep(i, k) = 0.0 sigp(i, k) = sigs ENDDO ENDDO ! ===================================================================== ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES ! ===================================================================== ! --- The procedure is to solve the equation. ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. ! *** Calculate certain parcel quantities, including static energy *** DO i = 1, ncum ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- & t0)) + gznk(i) ENDDO ! *** Find lifted parcel quantities above cloud base *** DO k = minorig + 1, nl DO i = 1, ncum IF (k>=(icb(i)+1)) THEN tg = t(i, k) qg = qs(i, k) alv = lv0 - clmcpv*(t(i,k)-t0) ! First iteration. s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k)) s = 1./s ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k) tg = tg + s*(ah0(i)-ahg) tg = max(tg, 35.0) tc = tg - t0 denom = 243.5 + tc IF (tc>=0.0) THEN es = 6.112*exp(17.67*tc/denom) ELSE es = exp(23.33086-6111.72784/tg+0.15215*log(tg)) ENDIF qg = eps*es/(p(i,k)-es*(1.-eps)) ! Second iteration. s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k)) s = 1./s ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k) tg = tg + s*(ah0(i)-ahg) tg = max(tg, 35.0) tc = tg - t0 denom = 243.5 + tc IF (tc>=0.0) THEN es = 6.112*exp(17.67*tc/denom) ELSE es = exp(23.33086-6111.72784/tg+0.15215*log(tg)) ENDIF qg = eps*es/(p(i,k)-es*(1.-eps)) alv = lv0 - clmcpv*(t(i,k)-t0) ! print*,'cpd dans convect2 ',cpd ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd' ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd ! if (.not.cpd.gt.1000.) then ! print*,'CPD=',cpd ! stop ! endif clw(i, k) = qnk(i) - qg clw(i, k) = max(0.0, clw(i,k)) rg = qg/(1.-qnk(i)) tvp(i, k) = tp(i, k)*(1.+rg*epsi) ENDIF ENDDO ENDDO ! ===================================================================== ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) ! ===================================================================== DO k = minorig + 1, nl DO i = 1, ncum IF (k>=(nk(i)+1)) THEN tca = tp(i, k) - t0 IF (tca>=0.0) THEN elacrit = elcrit ELSE elacrit = elcrit*(1.0-tca/tlcrit) ENDIF elacrit = max(elacrit, 0.0) ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8) ep(i, k) = max(ep(i,k), 0.0) ep(i, k) = min(ep(i,k), 1.0) sigp(i, k) = sigs ENDIF ENDDO ENDDO ! ===================================================================== ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL ! --- VIRTUAL TEMPERATURE ! ===================================================================== DO k = minorig + 1, nl DO i = 1, ncum IF (k>=(icb(i)+1)) THEN tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k)) ! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)' ! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k) ENDIF ENDDO ENDDO DO i = 1, ncum tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd ENDDO ! ===================================================================== ! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S ! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY ! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB) ! ===================================================================== DO i = 1, ncum cape(i) = 0.0 capem(i) = 0.0 inb(i) = icb(i) + 1 inb1(i) = inb(i) ENDDO ! Originial Code ! do 530 k=minorig+1,nl-1 ! do 520 i=1,ncum ! if(k.ge.(icb(i)+1))then ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) ! cape(i)=cape(i)+by ! if(by.ge.0.0)inb1(i)=k+1 ! if(cape(i).gt.0.0)then ! inb(i)=k+1 ! capem(i)=cape(i) ! endif ! endif ! 520 continue ! 530 continue ! do 540 i=1,ncum ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl) ! cape(i)=capem(i)+byp ! defrac=capem(i)-cape(i) ! defrac=max(defrac,0.001) ! frac(i)=-cape(i)/defrac ! frac(i)=min(frac(i),1.0) ! frac(i)=max(frac(i),0.0) ! 540 continue ! K Emanuel fix ! call zilch(byp,ncum) ! do 530 k=minorig+1,nl-1 ! do 520 i=1,ncum ! if(k.ge.(icb(i)+1))then ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) ! cape(i)=cape(i)+by ! if(by.ge.0.0)inb1(i)=k+1 ! if(cape(i).gt.0.0)then ! inb(i)=k+1 ! capem(i)=cape(i) ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) ! endif ! endif ! 520 continue ! 530 continue ! do 540 i=1,ncum ! inb(i)=max(inb(i),inb1(i)) ! cape(i)=capem(i)+byp(i) ! defrac=capem(i)-cape(i) ! defrac=max(defrac,0.001) ! frac(i)=-cape(i)/defrac ! frac(i)=min(frac(i),1.0) ! frac(i)=max(frac(i),0.0) ! 540 continue ! J Teixeira fix byp(1:ncum) = 0 DO i = 1, ncum lcape(i) = .TRUE. ENDDO DO k = minorig + 1, nl - 1 DO i = 1, ncum IF (cape(i)<0.0) lcape(i) = .FALSE. IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k) byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1) cape(i) = cape(i) + by IF (by>=0.0) inb1(i) = k + 1 IF (cape(i)>0.0) THEN inb(i) = k + 1 capem(i) = cape(i) ENDIF ENDIF ENDDO ENDDO DO i = 1, ncum cape(i) = capem(i) + byp(i) defrac = capem(i) - cape(i) defrac = max(defrac, 0.001) frac(i) = -cape(i)/defrac frac(i) = min(frac(i), 1.0) frac(i) = max(frac(i), 0.0) ENDDO ! ===================================================================== ! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL ! ===================================================================== ! initialization: !ym very bad !ym DO i = 1, ncum*nlp !ym hp(i, 1) = h(i, 1) !ym END DO DO k=1,nlp DO i=1,ncum hp(i, k) = h(i, k) ENDDO ENDDO DO k = minorig + 1, nl DO i = 1, ncum IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k & ) ENDIF ENDDO ENDDO RETURN END SUBROUTINE cv_undilute2 SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, & cpn, iflag, cbmf) USE lmdz_cv_ini, ONLY : alpha,damp,dtmax,minorig,rrd IMPLICIT NONE ! inputs: INTEGER ncum, nd, nloc INTEGER nk(nloc), icb(nloc) REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd) REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent... REAL plcl(nloc), cpn(nloc, nd) ! outputs: INTEGER iflag(nloc) REAL cbmf(nloc) ! also an input ! local variables: INTEGER i, k, icbmax REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc) REAL work(nloc) ! ------------------------------------------------------------------- ! Compute icbmax. ! ------------------------------------------------------------------- !ym independance between columns !ym icbmax = 2 !ym DO i = 1, ncum !ym icbmax = max(icbmax, icb(i)) !ym END DO ! ===================================================================== ! --- CALCULATE CLOUD BASE MASS FLUX ! ===================================================================== ! tvpplcl = parcel temperature lifted adiabatically from level ! icb-1 to the LCL. ! tvaplcl = virtual temperature at the LCL. DO i = 1, ncum dtpbl(i) = 0.0 tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( & i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1)) tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i & ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1)) ENDDO ! ------------------------------------------------------------------- ! --- Interpolate difference between lifted parcel and ! --- environmental temperatures to lifted condensation level ! ------------------------------------------------------------------- ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1). !ym independance betwwen column !ym DO k = minorig, icbmax DO k = minorig, nd DO i = 1, ncum IF (k<=MAX(2,icb(i))) THEN IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k) ENDIF ENDIF ENDDO ENDDO DO i = 1, ncum dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i))) dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i) ENDDO ! ------------------------------------------------------------------- ! --- Adjust cloud base mass flux ! ------------------------------------------------------------------- DO i = 1, ncum work(i) = cbmf(i) cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i)) IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN iflag(i) = 3 ENDIF ENDDO RETURN END SUBROUTINE cv_closure SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, & h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, & sij, elij) USE lmdz_cv_ini, ONLY : cpd,cpv,entp,minorig,nl,nlp,rrv IMPLICIT NONE ! inputs: INTEGER ncum, nd, nloc INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc) REAL cbmf(nloc), qnk(nloc) REAL ph(nloc, nd+1) REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd) REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd) REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd) ! outputs: INTEGER nent(nloc, nd) REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd) REAL uent(nloc, nd, nd), vent(nloc, nd, nd) REAL sij(nloc, nd, nd), elij(nloc, nd, nd) ! local variables: INTEGER i, j, k, ij INTEGER num1, num2 REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp REAL alt, qp1, smid, sjmin, sjmax, delp, delm REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc) REAL bsum(nloc, nd) LOGICAL lwork(nloc) ! ===================================================================== ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS ! ===================================================================== !ym very bad !ym DO i = 1, ncum*nlp !ym nent(i, 1) = 0 !ym m(i, 1) = 0.0 !ym END DO DO k = 1, nlp DO i = 1, ncum nent(i, k) = 0 m(i, k) = 0.0 ENDDO ENDDO DO k = 1, nlp DO j = 1, nlp DO i = 1, ncum qent(i, k, j) = q(i, j) uent(i, k, j) = u(i, j) vent(i, k, j) = v(i, j) elij(i, k, j) = 0.0 ment(i, k, j) = 0.0 sij(i, k, j) = 0.0 ENDDO ENDDO ENDDO ! ------------------------------------------------------------------- ! --- Calculate rates of mixing, m(i) ! ------------------------------------------------------------------- work(1:ncum) = 0. DO j = minorig + 1, nl DO i = 1, ncum IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN k = min(j, inb1(i)) dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + & entp*0.04*(ph(i,k)-ph(i,k+1)) work(i) = work(i) + dbo m(i, j) = cbmf(i)*dbo ENDIF ENDDO ENDDO DO k = minorig + 1, nl DO i = 1, ncum IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN m(i, k) = m(i, k)/work(i) ENDIF ENDDO ENDDO ! ===================================================================== ! --- 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 + 1, nl DO ij = 1, ncum IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= & inb(ij))) THEN qti = qnk(ij) - ep(ij, i)*clw(ij, i) bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd) anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j)) denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j) dei = denom IF (abs(dei)<0.01) dei = 0.01 sij(ij, i, j) = anum/dei sij(ij, i, i) = 1.0 altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j) altem = altem/bf2 cwat = clw(ij, j)*(1.-ep(ij,j)) stemp = sij(ij, i, j) IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2) denom = denom + lv(ij, j)*(q(ij,i)-qti) IF (abs(denom)<0.01) denom = 0.01 sij(ij, i, j) = anum/denom altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j) altem = altem - (bf2-1.)*cwat ENDIF IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + & (1.-sij(ij,i,j))*u(ij, nk(ij)) vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + & (1.-sij(ij,i,j))*v(ij, nk(ij)) elij(ij, i, j) = altem elij(ij, i, j) = max(0.0, elij(ij,i,j)) ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j)) nent(ij, i) = nent(ij, i) + 1 ENDIF sij(ij, i, j) = max(0.0, sij(ij,i,j)) sij(ij, i, j) = min(1.0, sij(ij,i,j)) ENDIF ENDDO ENDDO ! *** If no air can entrain at level i assume that updraft detrains ! *** ! *** at that level and calculate detrained air flux and properties ! *** DO ij = 1, ncum IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN ment(ij, i, i) = m(ij, i) qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i) uent(ij, i, i) = u(ij, nk(ij)) vent(ij, i, i) = v(ij, nk(ij)) elij(ij, i, i) = clw(ij, i) sij(ij, i, i) = 1.0 ENDIF ENDDO ENDDO DO i = 1, ncum sij(i, inb(i), inb(i)) = 1.0 ENDDO ! ===================================================================== ! --- NORMALIZE ENTRAINED AIR MASS FLUXES ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING ! ===================================================================== bsum(1:ncum,1:nlp) = 0. DO ij = 1, ncum lwork(ij) = .FALSE. ENDDO DO i = minorig + 1, nl !789 ENDDO num1 = 0 DO ij = 1, ncum IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1 ENDDO !ym IF (num1<=0) GO TO 789 IF (num1<=0) CYCLE DO ij = 1, ncum IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN lwork(ij) = (nent(ij,i)/=0) qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i) anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i)) denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1) IF (abs(denom)<0.01) denom = 0.01 scrit(ij) = anum/denom alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1) IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0 asij(ij) = 0.0 smin(ij) = 1.0 ENDIF ENDDO DO j = minorig, nl ! 783 ENDDO num2 = 0 DO ij = 1, ncum IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( & ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1 END DO !ym IF (num2<=0) GO TO 783 IF (num2<=0) CYCLE DO ij = 1, ncum IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( & ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN IF (j>i) THEN smid = min(sij(ij,i,j), scrit(ij)) sjmax = smid sjmin = smid IF (smid1) sjmin = sij(ij, i, j-1) sjmin = max(sjmin, scrit(ij)) ENDIF delp = abs(sjmax-smid) delm = abs(sjmin-smid) asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1)) ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1)) ENDIF ENDIF ENDDO 783 ENDDO DO ij = 1, ncum IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN asij(ij) = max(1.0E-21, asij(ij)) asij(ij) = 1.0/asij(ij) bsum(ij, i) = 0.0 ENDIF ENDDO DO j = minorig, nl + 1 DO ij = 1, ncum IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( & ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN ment(ij, i, j) = ment(ij, i, j)*asij(ij) bsum(ij, i) = bsum(ij, i) + ment(ij, i, j) ENDIF ENDDO ENDDO DO ij = 1, ncum IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, & i)<1.0E-18) .AND. lwork(ij)) THEN nent(ij, i) = 0 ment(ij, i, i) = m(ij, i) qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i) uent(ij, i, i) = u(ij, nk(ij)) vent(ij, i, i) = v(ij, nk(ij)) elij(ij, i, i) = clw(ij, i) sij(ij, i, i) = 1.0 ENDIF ENDDO 789 ENDDO RETURN END SUBROUTINE cv_mixing SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, & ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap) USE lmdz_cv_ini, ONLY : cl,coeffr,coeffs,cpd,g,ginv,nl,omtrain,omtsnow,sigd IMPLICIT NONE ! inputs: INTEGER ncum, nd, nloc INTEGER inb(nloc) REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd) REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd) REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd) REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd) REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd) ! outputs: INTEGER iflag(nloc) ! also an input REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd) REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd) ! local variables: INTEGER i, j, k, ij, num1 INTEGER jtt(nloc) REAL awat, coeff, qsm, afac, sigt, b6, c6, revap REAL dhdp, fac, qstm, rat REAL wdtrain(nloc) LOGICAL lwork(nloc) ! ===================================================================== ! --- PRECIPITATING DOWNDRAFT CALCULATION ! ===================================================================== ! Initializations: DO i = 1, ncum DO k = 1, nl + 1 wt(i, k) = omtsnow mp(i, k) = 0.0 evap(i, k) = 0.0 water(i, k) = 0.0 ENDDO ENDDO DO i = 1, ncum qp(i, 1) = q(i, 1) up(i, 1) = u(i, 1) vp(i, 1) = v(i, 1) ENDDO DO k = 2, nl + 1 DO i = 1, ncum qp(i, k) = q(i, k-1) up(i, k) = u(i, k-1) vp(i, k) = v(i, k-1) ENDDO ENDDO ! *** Check whether ep(inb)=0, if so, skip precipitating *** ! *** downdraft calculation *** ! *** Integrate liquid water equation to find condensed water *** ! *** and condensed water flux *** DO i = 1, ncum jtt(i) = 2 IF (ep(i,inb(i))<=0.0001) iflag(i) = 2 IF (iflag(i)==0) THEN lwork(i) = .TRUE. ELSE lwork(i) = .FALSE. ENDIF ENDDO ! *** Begin downdraft loop *** wdtrain(1:ncum) = 0. DO i = nl + 1, 1, -1 ! 899 ENDDO num1 = 0 DO ij = 1, ncum IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1 END DO !ym IF (num1<=0) GO TO 899 IF (num1<=0) CYCLE ! *** Calculate detrained precipitation *** DO ij = 1, ncum IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i) ENDIF ENDDO IF (i>1) THEN DO j = 1, i - 1 DO ij = 1, ncum IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i) awat = max(0.0, awat) wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i) ENDIF ENDDO ENDDO ENDIF ! *** Find rain water and evaporation using provisional *** ! *** estimates of qp(i)and qp(i-1) *** ! *** Value of terminal velocity and coeffecient of evaporation for snow ! *** DO ij = 1, ncum IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN coeff = coeffs wt(ij, i) = omtsnow ! *** Value of terminal velocity and coeffecient of evaporation for ! rain *** IF (t(ij,i)>273.0) THEN coeff = coeffr wt(ij, i) = omtrain ENDIF qsm = 0.5*(q(ij,i)+qp(ij,i+1)) afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0E4+2.0E3*ph(ij,i)*qs(ij,i)) afac = max(afac, 0.0) sigt = sigp(ij, i) sigt = max(0.0, sigt) sigt = min(1.0, sigt) b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i) c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i) revap = 0.5*(-b6+sqrt(b6*b6+4.*c6)) evap(ij, i) = sigt*afac*revap water(ij, i) = revap*revap ! *** Calculate precipitating downdraft mass flux under *** ! *** hydrostatic approximation *** IF (i>1) THEN dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i)) dhdp = max(dhdp, 10.0) mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp mp(ij, i) = max(mp(ij,i), 0.0) ! *** Add small amount of inertia to downdraft *** fac = 20.0/(ph(ij,i-1)-ph(ij,i)) mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac) ! *** Force mp to decrease linearly to zero ! *** ! *** between about 950 mb and the surface ! *** IF (p(ij,i)>(0.949*p(ij,1))) THEN jtt(ij) = max(jtt(ij), i) mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ & (p(ij,1)-p(ij,jtt(ij))) ENDIF ENDIF ! *** Find mixing ratio of precipitating downdraft *** IF (i/=inb(ij)) THEN IF (i==1) THEN qstm = qs(ij, 1) ELSE qstm = qs(ij, i-1) ENDIF IF (mp(ij,i)>mp(ij,i+1)) THEN rat = mp(ij, i+1)/mp(ij, i) qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + & 100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i)) up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat) vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat) ELSE IF (mp(ij,i+1)>0.0) THEN qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, & i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, & i)))/(lv(ij,i)+t(ij,i)*(cl-cpd)) up(ij, i) = up(ij, i+1) vp(ij, i) = vp(ij, i+1) ENDIF ENDIF qp(ij, i) = min(qp(ij,i), qstm) qp(ij, i) = max(qp(ij,i), 0.0) ENDIF ENDIF ENDDO 899 ENDDO RETURN END SUBROUTINE cv_unsat SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, & ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, & ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, & precip, cbmf, ft, fq, fu, fv, ma, qcondc) USE lmdz_cv_ini, ONLY : g,lv0,nl,rrd,sigd,betad,cl,cpd,cpv,cu,delta IMPLICIT NONE ! inputs INTEGER ncum, nd, nloc INTEGER nk(nloc), icb(nloc), inb(nloc) INTEGER nent(nloc, nd) REAL delt REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd) REAL gz(nloc, nd) REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd) REAL hp(nloc, nd), lv(nloc, nd) REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc) REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd) REAL up(nloc, nd), vp(nloc, nd) REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd) REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd) REAL uent(nloc, nd, nd), vent(nloc, nd, nd) REAL tv(nloc, nd), tvp(nloc, nd) ! outputs INTEGER iflag(nloc) ! also an input REAL cbmf(nloc) ! also an input REAL wd(nloc), tprime(nloc), qprime(nloc) REAL precip(nloc) REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd) REAL ma(nloc, nd) REAL qcondc(nloc, nd) ! local variables INTEGER i, j, ij, k, num1 REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti REAL work(nloc), am(nloc), amp1(nloc), ad(nloc) REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd) REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld ! -- initializations: delti = 1.0/delt DO i = 1, ncum precip(i) = 0.0 wd(i) = 0.0 tprime(i) = 0.0 qprime(i) = 0.0 DO k = 1, nl + 1 ft(i, k) = 0.0 fu(i, k) = 0.0 fv(i, k) = 0.0 fq(i, k) = 0.0 lvcp(i, k) = lv(i, k)/cpn(i, k) qcondc(i, k) = 0.0 ! cld qcond(i, k) = 0.0 ! cld nqcond(i, k) = 0.0 ! cld ENDDO ENDDO ! *** Calculate surface precipitation in mm/day *** DO i = 1, ncum IF (iflag(i)<=1) THEN ! c precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000. ! c & /(rowl*g) ! c precip(i)=precip(i)*delt/86400. precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g ENDIF ENDDO ! *** Calculate downdraft velocity scale and surface temperature and *** ! *** water vapor fluctuations *** DO i = 1, ncum wd(i) = betad*abs(mp(i,icb(i)))*0.01*rrd*t(i, icb(i))/(sigd*p(i,icb(i))) qprime(i) = 0.5*(qp(i,1)-q(i,1)) tprime(i) = lv0*qprime(i)/cpd ENDDO ! *** Calculate tendencies of lowest level potential temperature *** ! *** and mixing ratio *** DO i = 1, ncum work(i) = 0.01/(ph(i,1)-ph(i,2)) am(i) = 0.0 ENDDO DO k = 2, nl DO i = 1, ncum IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN am(i) = am(i) + m(i, k) ENDIF ENDDO ENDDO DO i = 1, ncum IF ((g*work(i)*am(i))>=delti) iflag(i) = 1 ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, & 1))/cpn(i,1)) ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1) ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* & work(i)/cpn(i, 1) fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + & sigd*evap(i, 1) fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i) fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, & 2)-u(i,1))) fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, & 2)-v(i,1))) ENDDO DO j = 2, nl DO i = 1, ncum IF (j<=inb(i)) THEN fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1)) fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1)) fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1)) ENDIF ENDDO ENDDO ! *** Calculate tendencies of potential temperature and mixing ratio *** ! *** at levels above the lowest level *** ! *** First find the net saturated updraft and downdraft mass fluxes *** ! *** through each level *** DO i = 2, nl + 1 ! 1500 ENDDO num1 = 0 DO ij = 1, ncum IF (i<=inb(ij)) num1 = num1 + 1 ENDDO !ym IF (num1<=0) GO TO 1500 IF (num1<=0) CYCLE amp1(1:ncum)=0. ad(1:ncum)=0. DO k = i + 1, nl + 1 DO ij = 1, ncum IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN amp1(ij) = amp1(ij) + m(ij, k) ENDIF ENDDO ENDDO DO k = 1, i DO j = i + 1, nl + 1 DO ij = 1, ncum IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN amp1(ij) = amp1(ij) + ment(ij, k, j) ENDIF ENDDO ENDDO ENDDO DO k = 1, i - 1 DO j = i, nl + 1 DO ij = 1, ncum IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN ad(ij) = ad(ij) + ment(ij, j, k) ENDIF ENDDO ENDDO ENDDO DO ij = 1, ncum IF (i<=inb(ij)) THEN dpinv = 0.01/(ph(ij,i)-ph(ij,i+1)) cpinv = 1.0/cpn(ij, i) ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, & i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, & i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i) ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij & ,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( & ij,i+1)-t(ij,i))*dpinv*cpinv fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, & i))-ad(ij)*(q(ij,i)-q(ij,i-1))) fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, & i))-ad(ij)*(u(ij,i)-u(ij,i-1))) fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, & i))-ad(ij)*(v(ij,i)-v(ij,i-1))) ENDIF ENDDO DO k = 1, i - 1 DO ij = 1, ncum IF (i<=inb(ij)) THEN dpinv = 0.01/(ph(ij,i)-ph(ij,i+1)) awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i) awat = max(awat, 0.0) fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q & (ij,i)) fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i & )) fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i & )) ! (saturated updrafts resulting from mixing) ! cld qcond(ij, i) = qcond(ij, i) + (elij(ij,k,i)-awat) ! cld nqcond(ij, i) = nqcond(ij, i) + 1. ! cld ENDIF ENDDO ENDDO DO k = i, nl + 1 DO ij = 1, ncum IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN dpinv = 0.01/(ph(ij,i)-ph(ij,i+1)) fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i & )) fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i & )) fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i & )) ENDIF ENDDO ENDDO DO ij = 1, ncum IF (i<=inb(ij)) THEN dpinv = 0.01/(ph(ij,i)-ph(ij,i+1)) fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, & i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, & i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, & i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv ! (saturated downdrafts resulting from mixing) ! cld DO k = i + 1, inb(ij) ! cld qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld nqcond(ij, i) = nqcond(ij, i) + 1. ! cld ENDDO ! cld ! (particular case: no detraining level is found) ! cld IF (nent(ij,i)==0) THEN ! cld qcond(ij, i) = qcond(ij, i) + (1.-ep(ij,i))*clw(ij, i) ! cld nqcond(ij, i) = nqcond(ij, i) + 1. ! cld ENDIF ! cld IF (nqcond(ij,i)/=0.) THEN ! cld qcond(ij, i) = qcond(ij, i)/nqcond(ij, i) ! cld ENDIF ! cld ENDIF ENDDO 1500 ENDDO ! *** Adjust tendencies at top of convection layer to reflect *** ! *** actual position of the level zero cape *** DO ij = 1, ncum fqold = fq(ij, inb(ij)) fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij)) fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, & inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, & inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1) ftold = ft(ij, inb(ij)) ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij)) ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, & inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, & inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1) fuold = fu(ij, inb(ij)) fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij)) fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, & inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij)))) fvold = fv(ij, inb(ij)) fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij)) fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, & inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij)))) ENDDO ! *** Very slightly adjust tendencies to force exact *** ! *** enthalpy, momentum and tracer conservation *** DO ij = 1, ncum ents(ij) = 0.0 uav(ij) = 0.0 vav(ij) = 0.0 DO i = 1, inb(ij) ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- & ph(ij,i+1)) uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1)) vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1)) ENDDO ENDDO DO ij = 1, ncum ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) ENDDO DO ij = 1, ncum DO i = 1, inb(ij) ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i) fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij)) fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij)) ENDDO ENDDO DO k = 1, nl + 1 DO i = 1, ncum IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10 ENDDO ENDDO DO i = 1, ncum IF (iflag(i)>2) THEN precip(i) = 0.0 cbmf(i) = 0.0 ENDIF ENDDO DO k = 1, nl DO i = 1, ncum IF (iflag(i)>2) THEN ft(i, k) = 0.0 fq(i, k) = 0.0 fu(i, k) = 0.0 fv(i, k) = 0.0 qcondc(i, k) = 0.0 ! cld ENDIF ENDDO ENDDO DO k = 1, nl + 1 DO i = 1, ncum ma(i, k) = 0. ENDDO ENDDO DO k = nl, 1, -1 DO i = 1, ncum ma(i, k) = ma(i, k+1) + m(i, k) ENDDO ENDDO ! *** diagnose the in-cloud mixing ratio *** ! cld ! *** of condensed water *** ! cld ! ! cld DO ij = 1, ncum ! cld DO i = 1, nd ! cld mac(ij, i) = 0.0 ! cld wa(ij, i) = 0.0 ! cld siga(ij, i) = 0.0 ! cld ENDDO ! cld DO i = nk(ij), inb(ij) ! cld DO k = i + 1, inb(ij) + 1 ! cld mac(ij, i) = mac(ij, i) + m(ij, k) ! cld ENDDO ! cld ENDDO ! cld DO i = icb(ij), inb(ij) - 1 ! cld ax(ij, i) = 0. ! cld DO j = icb(ij), i ! cld ax(ij, i) = ax(ij, i) + rrd*(tvp(ij,j)-tv(ij,j)) & ! cld *(ph(ij,j)-ph(ij,j+1))/p(ij, j) ! cld ENDDO ! cld IF (ax(ij,i)>0.0) THEN ! cld wa(ij, i) = sqrt(2.*ax(ij,i)) ! cld ENDIF ! cld ENDDO ! cld DO i = 1, nl ! cld IF (wa(ij,i)>0.0) & ! cld siga(ij, i) = mac(ij, i)/wa(ij, i) & ! cld *rrd*tvp(ij, i)/p(ij, i)/100./delta ! cld siga(ij, i) = min(siga(ij,i), 1.0) ! cld qcondc(ij, i) = siga(ij, i)*clw(ij, i)*(1.-ep(ij,i)) & ! cld +(1.-siga(ij,i))*qcond(ij, i) ! cld ENDDO ! cld ENDDO ! cld RETURN END SUBROUTINE cv_yield SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, is_convect, compress, iflag, precip, cbmf, ft, & fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, & qcondc1) USE lmdz_cv_ini, ONLY : nl IMPLICIT NONE ! inputs: INTEGER len, ncum, nd, nloc INTEGER idcum(nloc) LOGICAL is_convect(nloc) LOGICAL compress INTEGER iflag(nloc) REAL precip(nloc), cbmf(nloc) REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd) REAL ma(nloc, nd) REAL qcondc(nloc, nd) !cld ! outputs: INTEGER iflag1(len) REAL precip1(len), cbmf1(len) REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd) REAL ma1(len, nd) REAL qcondc1(len, nd) !cld ! local variables: INTEGER i, k IF (compress) THEN DO i = 1, ncum precip1(idcum(i)) = precip(i) cbmf1(idcum(i)) = cbmf(i) iflag1(idcum(i)) = iflag(i) ENDDO DO k = 1, nl DO i = 1, ncum ft1(idcum(i), k) = ft(i, k) fq1(idcum(i), k) = fq(i, k) fu1(idcum(i), k) = fu(i, k) fv1(idcum(i), k) = fv(i, k) ma1(idcum(i), k) = ma(i, k) qcondc1(idcum(i), k) = qcondc(i, k) ENDDO ENDDO ELSE DO i = 1, len IF (is_convect(i)) THEN precip1(i) = precip(i) cbmf1(i) = cbmf(i) iflag1(i) = iflag(i) ENDIF ENDDO DO k = 1, nl DO i = 1, ncum IF (is_convect(i)) THEN ft1(i, k) = ft(i, k) fq1(i, k) = fq(i, k) fu1(i, k) = fu(i, k) fv1(i, k) = fv(i, k) ma1(i, k) = ma(i, k) qcondc1(i, k) = qcondc(i, k) ENDIF ENDDO ENDDO ENDIF RETURN END SUBROUTINE cv_uncompress END MODULE cv_routines_mod