Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update of boundary layer turbulence algorithm #10

Closed
wants to merge 5 commits into from
Closed

Update of boundary layer turbulence algorithm #10

Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5b52dfb
add meps and ec support
tackandr Jun 18, 2021
b434fdf
incorporate feedback
tackandr Aug 5, 2021
1b4672c
Merge branch 'development' into STU-16161
tackandr Aug 5, 2021
bb122fe
Merge branch 'development' into STU-16161
tackandr Aug 12, 2021
b36b03c
take Landcover from analysis time step
tackandr Aug 13, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
147 changes: 100 additions & 47 deletions himan-scripts/bl-turbulence.lua
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,24 +19,48 @@

logger:Info("Calculating boundary layer turbulence")

-- fetch producer information
local currentProducer = configuration:GetTargetProducer()
local currentProducerName = currentProducer.GetName(currentProducer)

local turbparam = param("BLTURB-N")

-- Surface (10m) wind speed in m/s = lowest model level = 65
local lev = {}
local ws = param("FF-MS")
local lev = level(HPLevelType.kHybrid, 65)
if currentProducerName == "ECGMTA" then
lev = level(HPLevelType.kHeight, 10)
else
lev = level(HPLevelType.kHybrid, 65)
end

local wsms = luatool:FetchWithType(current_time, lev, ws, current_forecast_type)

if not wsms then
logger:Error("No data found")
return
-- Wind gust
local gustlevel = {}
local gust = param("FFG-MS")
if currentProducerName == "ECGMTA" then
gustlevel = level(HPLevelType.kGround, 0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just did this:

radon=> insert into param_level_transform values (default, 131, (select id from param where name = 'FFG-MS'), 1, 0, 6, 10);
INSERT 0 1

So Himan automatically knows to fetch FFG-MS from ground/0 if it's requested for height/10.

else
gustlevel = level(HPLevelType.kHeight, 10)
end

-- Max tke in the layer 0-500m (note: max tke is generally found at the 1st model level above the surface)
local tke = param("TKEN-JKG")
local maxtke = hitool:VerticalMaximum(tke,0,500)
local gustms = luatool:FetchWithType(current_time, gustlevel, gust, current_forecast_type)

-- land-sea-mask, proportion from 0 to 1 where 1=land, 0=sea
local landseamask = {}
local surface = {}
if currentProducerName == "ECGMTA" then
surface = level(HPLevelType.kGround, 0)
landseamask = param("N-0TO1")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why cloud cover for ECMWF?

else
surface = level(HPLevelType.kHeight, 0)
landseamask = param("LC-0TO1")
end

local lcmask = luatool:FetchWithType(current_time, surface, landseamask, current_forecast_type)

if not maxtke then
if not (wsms and gustms and lcmask) then
logger:Error("No data found")
return
end
Expand All @@ -48,116 +72,145 @@ logger:Info("Calculating wind shear")
-- because only change in wind direction and speed is needed, not absolutely correct direction

local dz = 304.8
-- 2) Use lowest model level wind for sfc shear in ***Hirlam:
-- Wind component differences between Hirlam lowest model level (L65 ~12m) and 1000ft (above it)
local U_HIR = param("U-MS")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are these two variables called "HIR" ?

local V_HIR = param("V-MS")
local u_dz = hitool:VerticalValue(U_HIR,dz)
local u_0 = luatool:FetchWithType(current_time, lev, U_HIR, current_forecast_type)
local v_dz = hitool:VerticalValue(V_HIR,dz)
local v_0 = luatool:FetchWithType(current_time, lev, V_HIR, current_forecast_type)
local u = hitool:VerticalValue(U_HIR,dz+12)
local v = hitool:VerticalValue(V_HIR,dz+12)

if not u_dz or not u_0 or not v_dz or not v_0 then
-- Wind component at 1000ft
local u1000 = hitool:VerticalValue(U_HIR,dz)
local v1000 = hitool:VerticalValue(V_HIR,dz)

if not (u and v and u1000 and v1000 and u_0 and v_0) then
return
end

local windshear = {}

for i = 1, #v_0 do

local u2 = u_dz[i]
local u1 = u_0[i]
local v2 = v_dz[i]
local v1 = v_0[i]
-- Wind component at 2000ft
local u2000 = hitool:VerticalValue(U_HIR,dz*2)
local v2000 = hitool:VerticalValue(V_HIR,dz*2)

local shear = math.sqrt(math.pow(((u2-u1)/(dz-10)),2)+ math.pow(((v2-v1)/(dz-10)),2))*(dz-10)*1.946
if not (u2000 and v2000) then
return
end

windshear[i]= shear
-- Wind component at 3000ft
local u3000 = hitool:VerticalValue(U_HIR,dz*3)
local v3000 = hitool:VerticalValue(V_HIR,dz*3)

if not (u3000 and v3000) then
return
end

-- land-sea-mask, proportion from 0 to 1 where 1=land, 0=sea
local landseamask = param("LC-0TO1")
local surface = level(HPLevelType.kHeight, 0)
local lcmask = luatool:FetchWithType(current_time, surface, landseamask, current_forecast_type)
-- Wind component at 4000ft
local u4000 = hitool:VerticalValue(U_HIR,dz*4)
local v4000 = hitool:VerticalValue(V_HIR,dz*4)

if not lcmask then
if not (u4000 and v4000) then
return
end

-- Turbulence intensity

local turbulence = {}
local Missing = missing

for i = 1, #maxtke do
for i = 1, #v_0 do

local du = u[i] - u_0[i]
local dv = v[i] - v_0[i]
local du12 = u1000[i] - u2000[i]
local dv12 = v1000[i] - v2000[i]
local du23 = u2000[i] - u3000[i]
local dv23 = v2000[i] - v3000[i]
local du34 = u3000[i] - u4000[i]
local dv34 = v3000[i] - v4000[i]

local shear = math.sqrt(math.pow(du/dz,2) + math.pow(dv/dz,2))*dz*1.946
local shear12 = math.sqrt(math.pow(du12/dz,2) + math.pow(dv12/dz,2))*dz*1.946
local shear23 = math.sqrt(math.pow(du23/dz,2) + math.pow(dv23/dz,2))*dz*1.946
local shear34 = math.sqrt(math.pow(du34/dz,2) + math.pow(dv34/dz,2))*dz*1.946

local maxshear = math.max(shear, shear12, shear23, shear34)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice way to reduce the verbose if-statements in the original macro 👍


-- Turbulence intensity
local wskt = wsms[i]*1.946
local lc = lcmask[i]
local maxTKE = maxtke[i]
local shear = windshear[i]
local turb = Missing
local wgust = (gustms[i] - wsms[i])*1.946

if wskt == wskt and maxTKE == maxTKE and shear == shear then
if IsValid(wskt) and IsValid(shear) then

turb = 0

-- Turbulence over land areas

if lc > 0.5 then
if lc > 0.4 then

-- FBL-MOD over land

if wskt>11 and maxTKE>4 then
if wskt>12 and shear >20 then
turb = 2
end
if wskt>14 and maxTKE>3.5 then
if wskt>14 and (wgust>14 or maxshear>15) then
turb = 2
end

-- MOD over land

if wskt>12 and shear*maxTKE >80 then
turb = 3
end
if wskt>14 and (maxTKE>4 or shear>25) then
if wskt>16 and (wgust>18 or maxshear>20) then
turb = 3
end
if wskt>17 and maxTKE>2.5 then
if wskt>18 and (wgust>14 or maxshear>15) then
turb = 3
end

-- MOD-SEV over land

if wskt>20 and maxTKE>4.5 then
if wskt>20 and (wgust>18 or maxshear>25) then
turb = 4
end
if wskt>22 and maxTKE>3.5 then
if wskt>22 and (wgust>18 or maxshear>20) then
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

wgust-conditional should be:

wgust>16

turb = 4
end

-- SEV over land

if wskt>=26 and maxTKE>4 then
if wskt>=24 and wgust>20 then
turb = 5
end
if wskt>26 and (wgust>18 or maxshear>25) then
turb = 5
end

end

-- Turbulence at sea (or over lakes)

if lc <= 0.5 then
if lc <= 0.4 then

-- FBL-MOD over sea

if wskt>25 and wgust>14 then
turb = 2
end

-- MOD over sea

if wskt>=35 and maxTKE>3 then -- MOD over sea
if wskt>=30 and wgust>15 then
turb = 3
end

-- MOD-SEV over sea

if wskt>=50 and maxTKE>4 then -- SEV over sea
if wskt>=35 and wgust>20 then
turb = 4
end

-- SEV over sea
if wskt>40 and wgust>23 then
turb = 5
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Overall looks good to me and I can't find any problems, good work.

end
end
turbulence[i]= turb
Expand Down