I am trying to calculate the Deardorff velocity w*

more fully described http://glossary.ametsoc.org/wiki/Convective_velocity_scale

A few questions have come up. (I am a pilot/engineer trying to create local convective layer forecasts in my area, similar to http://www.drjack.info/BLIP/NAM/INFO/parameters.html#wfpm and am learning the meteorological nitty gritty as I go )

I am planning to use an arbitrary constant Tv of approximately 300 kelvin, but should I be calculating/using a potential temperature instead?

I am using the PBLH parameter from my wrf output as Zi -- as I understand it, PBLH is 'over' terrain, and if I wanted to convert it to an altitude over sea level, I would compute ter+PBLH?

I am struggling to calculate the kinematic vertical turbulent flux of virtual potential temperature.

At present, I have found

HFX -- wrf output field defined as the surface sensible heat flux in W/m^2

LH -- latent heat flux defined as the surface latent heat flux in W/m^2

I found equations that describe surface virtual heat flux (as part of calculating bouyancy production of turbulent kinetic energy), which come down to something nice and simple

HFX + 0.07 * LH

I believe that I should be able to put this into the formula for w*, but when I do, I get numbers that are off by almost an order of magnitude. w* should be on the order of 1 m/s certainly not more than 10, so when I calculate something that gives me w* > 25 m/s I am sure that I am doing something wrong.

My code for calculating and plotting w* is included below. Any thoughts would be appreciated.

- Code: Select all

load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

;load "plotdefaults.ncl"

a = addfile("./wrfout_d01_2015-05-01_00:00:00.nc","r")

type = "x11"

wks = gsn_open_wks(type,"plt_Surface1")

gsn_define_colormap(wks,"WhViBlGrYeOrReWh")

res = True

res@MainTitle = "Real-Time WRF"

pltres = True

mpres = True

Times = wrf_user_getvar(a,"times",-1)

ntimes = dimsizes(Times)

function calcws(q,z)

local dims, itr_rs, ws, conv2ftmin, cbouy

begin

; wstar = ((g/Tavg) * Qsurf * D) ^ (1/3)

; g := gravitational field = ~ 9.8062 m/s^2

; Tavg := average temperature ~300 k

; Qsurf := virtual temperature heat flux at surface into air W/m^2

; D := boundary layer thickness

dims = dimsizes(z)

ws = new (dims,float)

ws@description="Thermal updraft velocity W*"

ws@units="m/s"

;conv2ftmin = 196.85 ; const conversion factor from m/s to ft/min

cbouy = 9.8062 / 300.0 ; hard code average temp as 300 K

do i = 0, dims(0) -1

do j = 0, dims(1)-1

itr_rs = cbouy*q(i,j)*z(i,j)

if (itr_rs .lt. 0.05)

; don't waste time calculating small values of wstar

ws(i,j) = 0.0

else

ws(i,j)=(itr_rs^0.333)

end if

end do

end do

return(ws)

end

ter = wrf_user_getvar(a,"HGT_M",0)

do it = 6,ntimes-1,1

print("Working on time: " + Times(it) )

ts = wrf_user_getvar(a,"TSK",it)

hfx = wrf_user_getvar(a,"HFX",it)

tc = wrf_user_getvar(a,"tc",it)

vhf = wrf_user_getvar(a,"LH",it)

qsurf = hfx + 0.07*vhf

qsurf@description="Surface virtual heat flux"

qsurf@units="W/m^2"

pblh = wrf_user_getvar(a,"PBLH",it)

zi = pblh+ter

zi@description="Boundary layer top height"

zi@units="meters"

wstar = calcws(qsurf,pblh)

opts = res

opts@cnFillOn = True

;con_ts=wrf_contour(a,wks,ts,opts)

;con_hfx=wrf_contour(a,wks,hfx,opts)

;con_vhf=wrf_contour(a,wks,qsurf,opts)

con_pblh=wrf_contour(a,wks,pblh,opts)

con_qsurf=wrf_contour(a,wks,qsurf,opts)

con_wstar=wrf_contour(a,wks,wstar,opts)

;plot = wrf_map_overlays(a,wks,(/con_ts/),pltres,mpres)

;plot = wrf_map_overlays(a,wks,(/con_hfx/),pltres,mpres)

plot = wrf_map_overlays(a,wks,(/con_qsurf/),pltres,mpres)

plot = wrf_map_overlays(a,wks,(/con_pblh/),pltres,mpres)

plot = wrf_map_overlays(a,wks,(/con_wstar/),pltres,mpres)

end do