Update of boundary layer turbulence algorithm

himan-plugins/source/luatool.cpp

@@ -1344,6 +1344,7 @@ void BindLib(lua_State* L)
 	              .def("Empty", &raw_time::Empty),
 	          class_<forecast_time>("forecast_time")
 	              .def(constructor<const raw_time&, const raw_time&>())
+                      .def(constructor<const raw_time&, const time_duration&>())
 	              .def("ClassName", &forecast_time::ClassName)
 	              .def("GetOriginDateTime", LUA_MEMFN(raw_time&, forecast_time, OriginDateTime, void))
 	              .def("GetValidDateTime", LUA_MEMFN(raw_time&, forecast_time, ValidDateTime, void))

himan-scripts/bl-turbulence.lua

@@ -1,8 +1,4 @@
--- Hirlam near-surface/boundary layer aircraft turbulence intensity.
--- Turbulent layer height may then be estimated by boundary layer height.
--- Requirements:
--- - TKE and wind on model levels (+model level heights) in the surface layer (0...500m)
--- - land/sea/lake differentiation
+-- Near-surface/boundary layer aircraft turbulence intensity.
 --
 -- The parameter gets (integer) values between 0-5:
 --   0 = nil turbulence/smooth
@@ -16,27 +12,41 @@
 -- https://wiki.fmi.fi/display/PROJEKTIT/Rajakerroksen+turbulenssi
 -- original smarttool script written by Simo Neiglick:
 -- junila 11/2015 new version 10/2016
+-- tack 06/2021 new version
 
 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
-end
+-- Wind gust
+local gustlevel = level(HPLevelType.kHeight, 10) 
+local gust = param("FFG-MS")
+
+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 = param("LC-0TO1")
+local surface = level(HPLevelType.kHeight, 0)
+local mytime = forecast_time(current_time:GetOriginDateTime(),time_duration(HPTimeResolution.kHourResolution,0))
 
--- 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 lcmask = luatool:FetchWithType(mytime, surface, landseamask, current_forecast_type)
 
-if not maxtke then
+if not (wsms and gustms and lcmask) then
   logger:Error("No data found")
   return
 end
@@ -48,116 +58,147 @@ logger:Info("Calculating wind shear")
 -- because only change in wind direction and speed is needed, not absolutely correct direction
 
 local dz = 304.8
+
+-- Use lowest model level wind for sfc shear in ***Hirlam/MEPS, 10m level at EC:
+-- Wind component differences between Hirlam lowest model level (L65 ~12m) and 1000ft (above it)
+
 local U_HIR = param("U-MS")
 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)
+
+-- Wind component at 1000ft
+local u1000 = hitool:VerticalValue(U_HIR,dz)
+local v1000 = hitool:VerticalValue(V_HIR,dz)
 
-if not u_dz or not u_0 or not v_dz or not v_0 then
+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
+-- Wind component at 2000ft
+local u2000 = hitool:VerticalValue(U_HIR,dz*2)
+local v2000 = hitool:VerticalValue(V_HIR,dz*2)
 
-  local u2 = u_dz[i]
-  local u1 = u_0[i]
-  local v2 = v_dz[i]
-  local v1 = v_0[i]
-
-  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)
 
+  -- 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>16 or maxshear>20) then
         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
     end
   end
   turbulence[i]= turb