subroutine temf2d(j,ux,vx,tx,thx,qvx,qcx,qix,p2d,p2di,pi2d,rho, & rubltenx,rvbltenx,rthbltenx, & rqvbltenx,rqcbltenx,rqibltenx, & g,cp,rcp,r_d,r_v,cpv, & z2d, & xlv,psfcpa, & znt,zsrf,ust,zol,hol,hpbl,psim,psih, & xland,hfx,qfx,tsk,qsfc,gz1oz0,wspd,br, & dt,dtmin,kpbl1d, & svp1,svp2,svp3,svpt0,ep1,ep2,karman,eomeg,stbolt, & kh_temfx,km_temfx, & u10,v10,t2, & te_temfx,shf_temfx,qf_temfx,uw_temfx,vw_temfx, & wupd_temfx,mf_temfx,thup_temfx,qtup_temfx,qlup_temfx, & cf3d_temfx,cfm_temfx, & hd_temfx,lcl_temfx,hct_temfx,exch_temfx, & flhc,flqc, & fCor, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte & ) !------------------------------------------------------------------- implicit none !------------------------------------------------------------------- ! ! This is the Total Energy - Mass Flux (TEMF) PBL scheme. ! Initial implementation 2010 by Wayne Angevine, CIRES/NOAA ESRL. ! References: ! Angevine et al., 2010, MWR ! Angevine, 2005, JAM ! Mauritsen et al., 2007, JAS ! !------------------------------------------------------------------- ! integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, j ! real, intent(in ) :: dt,dtmin,g,cp,rcp,r_d,r_v,cpv,xlv ! real, intent(in ) :: svp1,svp2,svp3,svpt0 real, intent(in ) :: ep1,ep2,karman,eomeg,stbolt ! real, dimension( ims:ime, kms:kme ), & intent(in) :: z2d ! real, dimension( ims:ime, kms:kme ) , & intent(in ) :: ux, vx real, dimension( ims:ime, kms:kme ) , & intent(inout) :: te_temfx real, dimension( ims:ime, kms:kme ) , & intent( out) :: shf_temfx, qf_temfx, uw_temfx, vw_temfx , & wupd_temfx, mf_temfx,thup_temfx, & qtup_temfx, qlup_temfx, cf3d_temfx real, dimension( ims:ime ) , & intent( out) :: hd_temfx, lcl_temfx, hct_temfx, cfm_temfx real, dimension( ims:ime ) , & intent(in ) :: fCor real, dimension( ims:ime ) , & intent(inout) :: flhc, flqc, exch_temfx real, dimension( ims:ime, kms:kme ) , & intent(in ) :: tx, thx, qvx, qcx, qix, pi2d, rho real, dimension( ims:ime, kms:kme ) , & intent(in ) :: p2di, p2d ! real, dimension( ims:ime, kms:kme ) , & intent(inout) :: rubltenx, rvbltenx, rthbltenx, & rqvbltenx, rqcbltenx, rqibltenx ! real, dimension( ims:ime ) , & intent(inout) :: hol, ust, hpbl, znt real, dimension( ims:ime ) , & intent(in ) :: xland, zsrf real, dimension( ims:ime ) , & intent(inout) :: hfx, qfx ! real, dimension( ims:ime ), intent(inout) :: wspd real, dimension( ims:ime ), intent(in ) :: br ! real, dimension( ims:ime ), intent(in ) :: psim, psih real, dimension( ims:ime ), intent(in ) :: gz1oz0 ! real, dimension( ims:ime ), intent(in ) :: psfcpa real, dimension( ims:ime ), intent(in ) :: tsk, qsfc real, dimension( ims:ime ), intent(inout) :: zol integer, dimension( ims:ime ), intent(out ) :: kpbl1d real, dimension( ims:ime, kms:kme ) , & intent(inout) :: kh_temfx, km_temfx ! real, dimension( ims:ime ) , & intent(inout) :: u10, v10, t2 ! ! !----------------------------------------------------------- ! Local variables ! ! TE model constants logical, parameter :: MFopt = .true. ! Use mass flux or not real, parameter :: visc_temf = 1.57e-4 ! WA TEST bigger minimum K real, parameter :: conduc_temf = 1.57e-4 / 0.733 real, parameter :: Pr_temf = 0.733 real, parameter :: TEmin = 1e-3 real, parameter :: ftau0 = 0.17 real, parameter :: fth0 = 0.145 real, parameter :: critRi = 0.25 real, parameter :: Cf = 0.185 real, parameter :: CN = 2.0 real, parameter :: Ceps = 0.070 real, parameter :: Cgamma = Ceps real, parameter :: Cphi = Ceps real, parameter :: PrT0 = Cphi/Ceps * ftau0**2 / 2. / fth0**2 ! EDMF constants ! WA TEST 10/16/17 reduce initial M to reduce need for flux limiter real, parameter :: CM = 0.03 ! Proportionality constant for subcloud MF (orig) ! real, parameter :: CM = 0.02 ! Proportionality constant for subcloud MF real, parameter :: Cdelt = 0.006 ! Prefactor for detrainment rate real, parameter :: Cw = 0.5 ! Prefactor for surface wUPD real, parameter :: Cc = 3.0 ! Prefactor for convective length scale real, parameter :: lasymp = 200.0 ! Asymptotic length scale WA 11/20/09 real, parameter :: hmax = 4000.0 ! Max hd,hct WA 11/20/09 ! integer :: i, k, kt ! Loop variable integer, dimension( its:ite) :: h0idx real, dimension( its:ite) :: h0 real, dimension( its:ite) :: wstr, ang, wm real, dimension( its:ite) :: hd,lcl,hct,ht real, dimension( its:ite) :: convection_TKE_surface_src, sfcFTE real, dimension( its:ite) :: sfcTHVF real, dimension( its:ite) :: z0t integer, dimension( its:ite) :: hdidx,lclidx,hctidx,htidx integer, dimension( its:ite) :: hmax_idx integer, dimension( its:ite) :: tval real, dimension( its:ite, kts:kte) :: thetal, qt real, dimension( its:ite, kts:kte) :: u_temf, v_temf real, dimension( its:ite, kts:kte) :: rv, rl, rt real, dimension( its:ite, kts:kte) :: chi_poisson, gam real, dimension( its:ite, kts:kte) :: dthdz, dqtdz, dudz, dvdz real, dimension( its:ite, kts:kte) :: lepsmin real, dimension( its:ite, kts:kte) :: thetav real, dimension( its:ite, kts:kte) :: MFCth, MFCq, MFCu, MFCv real, dimension( its:ite, kts:kte) :: MFCql, MFCthv, MFCTE real, dimension( its:ite, kts:kte) :: epsmf, deltmf, dMdz real, dimension( its:ite, kts:kte) :: UUPD, VUPD real, dimension( its:ite, kts:kte) :: thetavUPD, qlUPD, TEUPD real, dimension( its:ite, kts:kte) :: thetavUPDmoist, wupd_dry real, dimension( its:ite, kts:kte) :: B, Bmoist real, dimension( its:ite, kts:kte) :: zm, zt, dzm, dzt real, dimension( its:ite, kts:kte) :: dthUPDdz, dqtup_temfxdz, dwUPDdz real, dimension( its:ite, kts:kte) :: dwUPDmoistdz real, dimension( its:ite, kts:kte) :: dUUPDdz, dVUPDdz, dTEUPDdz real, dimension( its:ite, kts:kte) :: TUPD, rstUPD, rUPD, rlUPD, qstUPD real, dimension( its:ite, kts:kte) :: N2, S, Ri, beta, ftau, fth, ratio real, dimension( its:ite, kts:kte) :: TKE, TE2 real, dimension( its:ite, kts:kte) :: ustrtilde, linv, leps real, dimension( its:ite, kts:kte) :: km, kh real, dimension( its:ite, kts:kte) :: Fz, QFK, uwk, vwk real, dimension( its:ite, kts:kte) :: km_conv, kh_conv, lconv real, dimension( its:ite, kts:kte) :: alpha2, beta2 ! For thetav flux calculation real, dimension( its:ite, kts:kte) :: THVF, buoy_src, srcs real, dimension( its:ite, kts:kte) :: u_new, v_new real, dimension( its:ite, kts:kte) :: thx_new, qvx_new, qcx_new real, dimension( its:ite, kts:kte) :: thup_new, qvup_new real, dimension( its:ite, kts:kte) :: beta1 ! For saturation humidity calculations real Cepsmf ! Prefactor for entrainment rate real red_fact ! for reducing MF components logical is_convective ! Vars for cloud fraction calculation real, dimension( its:ite, kts:kte) :: au, sigq, qst, satdef real sigq2, rst !---------------------------------------------------------------------- ! Grid staggering: Matlab version has mass and turbulence levels. ! WRF has full levels (with w) and half levels (u,v,theta,q*). Both ! sets of levels use the same indices (kts:kte). See pbl_driver or ! WRF Physics doc for (a few) details. ! So *mass levels correspond to half levels.* ! WRF full levels are ignored, we define our own turbulence levels ! in order to put the first one below the first half level. ! Another difference is that ! the Matlab version (and the Mauritsen et al. paper) consider the ! first mass level to be at z0 (effectively the surface). WRF considers ! the first half level to be above the effective surface. The first half ! level, at k=1, has nonzero values of u,v for example. Here we convert ! all incoming variables to internal ones with the correct indexing ! in order to make the code consistent with the Matlab version. We ! already had to do this for thetal and qt anyway, so the only additional ! overhead is for u and v. ! I use suffixes m for mass and t for turbulence as in Matlab for things ! like indices. ! Note that zsrf is the terrain height ASL, from Registry variable ht. ! Translations (Matlab to WRF): ! dzt -> calculated below ! dzm -> not supplied, calculated below ! k -> karman ! z0 -> znt ! z0t -> not in WRF, calculated below ! zt -> calculated below ! zm -> (z2d - zsrf) but NOTE zm(1) is now z0 (znt) and zm(2) is ! z2d(1) - zsrf ! ! WA I take the temperature at z0 to be ! TSK. This isn't exactly robust. ! WA 2/16/11 removed calculation of flhc, flqc which are not needed here. ! These should be removed from the calling sequence someday. ! ! Other notes: ! - I have often used 1 instead of kts below, because the scheme demands ! to know where the surface is. It won't work if kts .NE. 1. do i = its,ite ! Main loop ! Get incoming surface theta from TSK (WA for now) thetal(i,1) = tsk(i) / pi2d(i,1) ! WA really should use Exner func. at z0 ! WA TEST 7/27/17 use surface moisture, not first level ! qt(i,1) = qvx(i,1) qt(i,1) = qsfc(i) rv(i,1) = qt(i,1) / (1.-qt(i,1)) ! Water vapor rl(i,1) = 0. rt(i,1) = rv(i,1) + rl(i,1) ! Total water (without ice) chi_poisson(i,1) = rcp * (1.+rv(i,1)/ep2) / (1.+rv(i,1)*cpv/cp) gam(i,1) = rv(i,1) * r_v / (cp + rv(i,1)*cpv) thetav(i,1) = thetal(i,1) * (1. + 0.608*qt(i,1) - qcx(i,1)) ! WA 4/6/10 allow environment liquid z0t(i) = znt(i) ! Convert incoming theta to thetal and qv,qc to qt ! NOTE this is where the indexing gets changed from WRF to TEMF basis do k = kts+1,kte ! Convert specific humidities to mixing ratios rv(i,k) = qvx(i,k-1) / (1.-qvx(i,k-1)) ! Water vapor rl(i,k) = qcx(i,k-1) / (1.-qcx(i,k-1)) ! Liquid water rt(i,k) = rv(i,k) + rl(i,k) ! Total water (without ice) chi_poisson(i,k) = rcp * (1.+rv(i,k)/ep2) / (1.+rv(i,k)*cpv/cp) gam(i,k) = rt(i,k) * r_v / (cp + rt(i,k)*cpv) thetal(i,k) = thx(i,k-1) * & ((ep2+rv(i,k))/(ep2+rt(i,k)))**chi_poisson(i,k) * & (rv(i,k)/rt(i,k))**(-gam(i,k)) * exp( -xlv*rl(i,k) / & ((cp + rt(i,k)*cpv) * tx(i,k))) qt(i,k) = qvx(i,k-1) + qcx(i,k-1) thetav(i,k) = thetal(i,k) * (1. + 0.608*qt(i,k) - qcx(i,k-1)) ! WA 4/6/10 allow environment liquid end do ! Convert incoming u,v to internal u_temf, v_temf ! NOTE this is where the indexing gets changed from WRF to TEMF basis u_temf(i,1) = 0. ! zero winds at z0 v_temf(i,1) = 0. do k = kts+1,kte u_temf(i,k) = ux(i,k-1) v_temf(i,k) = vx(i,k-1) end do ! Get delta height at half (mass) levels zm(i,1) = znt(i) dzt(i,1) = z2d(i,1) - zsrf(i) - zm(i,1) ! Get height and delta at turbulence levels zt(i,1) = (z2d(i,1) - zsrf(i) - znt(i)) / 2. do kt = kts+1,kte zm(i,kt) = z2d(i,kt-1) - zsrf(i) ! Convert indexing from WRF to TEMF zt(i,kt) = (zm(i,kt) + z2d(i,kt) - zsrf(i)) / 2. dzm(i,kt) = zt(i,kt) - zt(i,kt-1) dzt(i,kt) = z2d(i,kt+1) - z2d(i,kt) end do dzm(i,1) = dzm(i,2) ! WA why? dzt(i,kte) = dzt(i,kte-1) ! WA 12/23/09 ! Gradients at first level dthdz(i,1) = (thetal(i,2)-thetal(i,1)) / (zt(i,1) * log10(zm(i,2)/z0t(i))) dqtdz(i,1) = (qt(i,2)-qt(i,1)) / (zt(i,1) * log10(zm(i,2)/z0t(i))) dudz(i,1) = (u_temf(i,2)-u_temf(i,1)) / (zt(i,1) * log10(zm(i,2)/znt(i))) dvdz(i,1) = (v_temf(i,2)-v_temf(i,1)) / (zt(i,1) * log10(zm(i,2)/znt(i))) ! Surface thetaV flux from Stull p.147 sfcTHVF(i) = hfx(i)/(rho(i,1)*cp) * (1.+0.608*(qvx(i,1)+qcx(i,1))) + 0.608*thetav(i,1)*qf_temfx(i,1) ! WA use hd_temf to calculate w* instead of finding h0 here???? ! Watch initialization! h0idx(i) = 1 h0(i) = zm(i,1) lepsmin(i,kts) = 0. ! WA 2/11/13 find index just above hmax for use below hmax_idx(i) = kte-1 do k = kts+1,kte-1 lepsmin(i,k) = 0. ! Mean gradients dthdz(i,k) = (thetal(i,k+1) - thetal(i,k)) / dzt(i,k) dqtdz(i,k) = (qt(i,k+1) - qt(i,k)) / dzt(i,k) dudz(i,k) = (u_temf(i,k+1) - u_temf(i,k)) / dzt(i,k) dvdz(i,k) = (v_temf(i,k+1) - v_temf(i,k)) / dzt(i,k) ! Find h0 (should eventually be interpolated for smoothness) ! WA TEST 10/30/17 reduce lower temperature to reduce h0 ! if (thetav(i,k) > thetav(i,1) .AND. h0idx(i) .EQ. 1) then if (thetav(i,k) > (thetav(i,1)+thetav(i,2))/2.0 .AND. h0idx(i) .EQ. 1) then ! WA 9/28/11 limit h0 as for hd and hct if (zm(i,k) < hmax) then h0idx(i) = k h0(i) = zm(i,k) else h0idx(i) = k h0(i) = hmax end if end if ! WA 2/11/13 find index just above hmax for use below if (zm(i,k) > hmax) then hmax_idx(i) = min(hmax_idx(i),k) end if end do ! Gradients at top level dthdz(i,kte) = dthdz(i,kte-1) dqtdz(i,kte) = dqtdz(i,kte-1) dudz(i,kte) = dudz(i,kte-1) dvdz(i,kte) = dvdz(i,kte-1) if ( hfx(i) > 0.) then ! wstr(i) = (g * h0(i) / thetav(i,2) * shf_temfx(i,1) ) ** (1./3.) wstr(i) = (g * h0(i) / thetav(i,2) * hfx(i)/(rho(i,1)*cp) ) ** (1./3.) else wstr(i) = 0. end if ! Set flag convective or not for use below is_convective = wstr(i) > 0. .AND. MFopt .AND. dthdz(i,1)<0. .AND. dthdz(i,2)<0. ! WA 12/16/09 require two levels of negative (unstable) gradient ! Find stability parameters and length scale (on turbulence levels) do kt = 1,kte-1 N2(i,kt) = 2. * g / (thetav(i,kt) + thetav(i,kt+1))*dthdz(i,kt) S(i,kt) = sqrt(dudz(i,kt)**2. + dvdz(i,kt)**2.) Ri(i,kt) = N2(i,kt) / S(i,kt)**2. if (S(i,kt) < 1e-15) then if (N2(i,kt) >= 0) then Ri(i,kt) = 10. else Ri(i,kt) = -1. end if end if beta(i,kt) = 2. * g / (thetav(i,kt)+thetav(i,kt+1)) if (Ri(i,kt) > 0) then ratio(i,kt) = Ri(i,kt)/(Cphi**2.*ftau0**2./(2.*Ceps**2.*fth0**2.)+3.*Ri(i,kt)) ftau(i,kt) = ftau0 * ((3./4.) / (1.+4.*Ri(i,kt)) + 1./4.) fth(i,kt) = fth0 / (1.+4.*Ri(i,kt)) TE2(i,kt) = 2. * te_temfx(i,kt) * ratio(i,kt) * N2(i,kt) / beta(i,kt)**2. else ratio(i,kt) = Ri(i,kt)/(Cphi**2.*ftau0**2./(-2.*Ceps**2.*fth0**2.)+2.*Ri(i,kt)) ftau(i,kt) = ftau0 fth(i,kt) = fth0 TE2(i,kt) = 0. end if TKE(i,kt) = te_temfx(i,kt) * (1. - ratio(i,kt)) ustrtilde(i,kt) = sqrt(ftau(i,kt) * TKE(i,kt)) if (N2(i,kt) > 0.) then linv(i,kt) = 1./karman / zt(i,kt) + abs(fCor(i)) / & (Cf*ustrtilde(i,kt)) + & sqrt(N2(i,kt))/(CN*ustrtilde(i,kt)) + 1./lasymp else linv(i,kt) = 1./karman / zt(i,kt) + abs(fCor(i)) / & (Cf*ustrtilde(i,kt)) + 1./lasymp end if leps(i,kt) = 1./linv(i,kt) leps(i,kt) = max(leps(i,kt),lepsmin(i,kt)) end do S(i,kte) = 0.0 N2(i,kte) = 0.0 TKE(i,kte) = 0.0 linv(i,kte) = linv(i,kte-1) leps(i,kte) = leps(i,kte-1) ! Find diffusion coefficients ! First use basic formulae for stable and neutral cases, ! then for convective conditions, and finally choose the larger ! WA 12/23/09 use convective form up to hd/2 always ! WA 12/28/09 after restructuring, this block is above MF block, ! so hd is not yet available for this timestep, must use h0, ! or use hd from previous timestep but be careful about initialization. do kt = 1,kte-1 ! WA 12/22/09 ! WA 4/8/10 remove beta term to avoid negative and huge values ! of km due to very small denominator. This is an interim fix ! until we find something better (more theoretically sound). ! km(i,kt) = TKE(i,kt)**1.5 * ftau(i,kt)**2. / (-beta(i,kt) * fth(i,kt) * sqrt(TE2(i,kt)) + Ceps * sqrt(TKE(i,kt)*te_temfx(i,kt)) / leps(i,kt)) km(i,kt) = TKE(i,kt)**1.5 * ftau(i,kt)**2. / (Ceps * sqrt(TKE(i,kt)*te_temfx(i,kt)) / leps(i,kt)) kh(i,kt) = 2. * leps(i,kt) * fth(i,kt)**2. * TKE(i,kt) / sqrt(te_temfx(i,kt)) / Cphi if ( is_convective) then ! WA 2/20/14 trap rare "equality" of h0 and zt (only when h0 = hmax) if (kt <= h0idx(i) .AND. h0(i)-zt(i,kt) > 1e-15) then lconv(i,kt) = 1. / (1. / (karman*zt(i,kt)) + Cc / (karman * (h0(i) - zt(i,kt)))) else lconv(i,kt) = 0. end if ! WA 12/15/09 use appropriate coeffs to match kh_conv and kh at neutral kh_conv(i,kt) = ftau0**2. / Ceps / PrT0 * sqrt(TKE(i,kt)) * lconv(i,kt) if (kh_conv(i,kt) < 0.) then kh_conv(i,kt) = 0. end if km_conv(i,kt) = PrT0 * kh_conv(i,kt) if (zt(i,kt) <= h0(i)/2.) then km(i,kt) = km_conv(i,kt) kh(i,kt) = kh_conv(i,kt) end if if (zt(i,kt) > h0(i)/2. .AND. kt <= h0idx(i)) then km(i,kt) = max(km(i,kt),km_conv(i,kt),visc_temf) kh(i,kt) = max(kh(i,kt),kh_conv(i,kt),conduc_temf) end if end if ! is_convective km(i,kt) = max(km(i,kt),visc_temf) kh(i,kt) = max(kh(i,kt),conduc_temf) Fz(i,kt) = -kh(i,kt) * dthdz(i,kt) ! Diffusive heat flux end do km(i,kte) = km(i,kte-1) kh(i,kte) = kh(i,kte-1) Fz(i,kte) = 0.0 !*** Mass flux block starts here *** if ( is_convective) then Cepsmf = 2.0 / max(200.,h0(i)) ! orig do k = kts,kte ! Calculate lateral entrainment fraction for subcloud layer ! epsilon and delta are defined on mass grid (half levels) epsmf(i,k) = Cepsmf ! WA TEST 8/3/17 make sure eps*dz<1 epsmf(i,k) = min(epsmf(i,k),0.9/dzt(i,k)) if (k<=hdidx(i) .AND. i==1 .AND. j==1 .AND. epsmf(i,k) 20.0) then te_temfx(i,1) = 20.0 ! WA 9/28/11 limit max TE end if sfcFTE(i) = -(km(i,1)+km(i,2)) / 2. * (te_temfx(i,2)-te_temfx(i,1)) / dzm(i,2) do kt = 1,kte if (THVF(i,kt) >= 0) then buoy_src(i,kt) = 2. * g / thetav(i,kt) * THVF(i,kt) else buoy_src(i,kt) = 0. ! Cancel buoyancy term when locally stable end if srcs(i,kt) = -uw_temfx(i,kt) * dudz(i,kt) - vw_temfx(i,kt) * dvdz(i,kt) - & Cgamma * te_temfx(i,kt)**1.5 * linv(i,kt) + buoy_src(i,kt) end do call solve_implicit_temf((km(i,kts:kte-1)+km(i,kts+1:kte))/2.0, & te_temfx(i,kts+1:kte),sfcFTE(i),dzt(i,kts+1:kte),dzt(i,kts:kte-1),kts,kte-1,dt,.false.) do kt = 2,kte-1 te_temfx(i,kt) = te_temfx(i,kt) + dt * srcs(i,kt) te_temfx(i,kt) = te_temfx(i,kt) + dt * (-(MFCTE(i,kt)-MFCTE(i,kt-1))) / dzt(i,kt) if (te_temfx(i,kt) < TEmin) te_temfx(i,kt) = TEmin end do te_temfx(i,kte) = 0.0 do kt = 2,kte-1 if (te_temfx(i,kt) > 20.0) then te_temfx(i,kt) = 20.0 ! WA 9/29/11 reduce limit max TE from 30 end if end do ! Done with updates, now convert internal variables back to WRF vars do k = kts,kte ! Populate kh_temfx, km_temfx from kh, km kh_temfx(i,k) = kh(i,k) km_temfx(i,k) = km(i,k) end do ! Convert thetal, qt back to theta, qv, qc ! See opposite conversion at top of subroutine ! WA this accounts for offset of indexing between ! WRF and TEMF, see notes at top of this routine. call thlqt2thqvqc(thetal(i,kts+1:kte),qt(i,kts+1:kte), & thx_new(i,kts:kte-1),qvx_new(i,kts:kte-1),qcx_new(i,kts:kte-1), & p2d(i,kts:kte-1),pi2d(i,kts:kte-1),kts,kte-1,ep2,xlv,cp) do k = kts,kte-1 ! Calculate tendency terms ! WA this accounts for offset of indexing between ! WRF and TEMF, see notes at top of this routine. rubltenx(i,k) = (u_new(i,k+1) - u_temf(i,k+1)) / dt rvbltenx(i,k) = (v_new(i,k+1) - v_temf(i,k+1)) / dt rthbltenx(i,k) = (thx_new(i,k) - thx(i,k)) / dt rqvbltenx(i,k) = (qvx_new(i,k) - qvx(i,k)) / dt rqcbltenx(i,k) = (qcx_new(i,k) - qcx(i,k)) / dt end do rubltenx(i,kte) = 0. rvbltenx(i,kte) = 0. rthbltenx(i,kte) = 0. rqvbltenx(i,kte) = 0. rqcbltenx(i,kte) = 0. ! Populate surface exchange coefficient variables to go back out ! for next time step of surface scheme ! WA 2/16/11 removed, not needed in BL ! Populate 10 m winds and 2 m temp ! WA Note this only works if first mass level is above 10 m u10(i) = u_new(i,2) * log(10.0/znt(i)) / log(zm(i,2)/znt(i)) v10(i) = v_new(i,2) * log(10.0/znt(i)) / log(zm(i,2)/znt(i)) t2(i) = (tsk(i)/pi2d(i,1) + (thx_new(i,1) - tsk(i)/pi2d(i,1)) * log(2.0/z0t(i)) / log(zm(i,2)/z0t(i))) * pi2d(i,1) ! WA this should also use pi at z0 ! Populate diagnostic variables hd_temfx(i) = hd(i) lcl_temfx(i) = lcl(i) hct_temfx(i) = hct(i) ! Send updraft liquid water back if ( is_convective) then do k = kts,kte-1 qlup_temfx(i,k) = qlUPD(i,k) end do else qlup_temfx(i,1) = qcx(i,1) do k = kts+1,kte-1 qlup_temfx(i,k) = qcx(i,k-1) end do end if qlup_temfx(i,kte) = qcx(i,kte) end do ! Main (i) loop end subroutine temf2d