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
a = addfile("./wrfout_d01_2015-05-01_00:00:00.nc","r")
type = "x11"
wks = gsn_open_wks(type,"plt_Surface1")
res = True
res@MainTitle = "Real-Time WRF"
pltres = True
mpres = True
Times = wrf_user_getvar(a,"times",-1)
ntimes = dimsizes(Times)
local dims, itr_rs, ws, conv2ftmin, cbouy
; 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*"
;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
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"
pblh = wrf_user_getvar(a,"PBLH",it)
zi = pblh+ter
zi@description="Boundary layer top height"
wstar = calcws(qsurf,pblh)
opts = res
opts@cnFillOn = True
;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)