From 5b52dfb23aeec100c385946e4b60541c379ba620 Mon Sep 17 00:00:00 2001 From: tackandr Date: Fri, 18 Jun 2021 14:51:00 +0300 Subject: [PATCH 1/3] add meps and ec support --- himan-scripts/bl-turbulence.lua | 147 ++++++++++++++++++++++---------- 1 file changed, 100 insertions(+), 47 deletions(-) diff --git a/himan-scripts/bl-turbulence.lua b/himan-scripts/bl-turbulence.lua index db284b363..bffefef5f 100644 --- a/himan-scripts/bl-turbulence.lua +++ b/himan-scripts/bl-turbulence.lua @@ -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) +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") +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 @@ -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") 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) + + -- 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 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 From b434fdf5e813ae3775682c96ba5b214f34d9f38d Mon Sep 17 00:00:00 2001 From: tackandr Date: Thu, 5 Aug 2021 14:56:33 +0300 Subject: [PATCH 2/3] incorporate feedback --- himan-scripts/bl-turbulence.lua | 31 +++++++++---------------------- 1 file changed, 9 insertions(+), 22 deletions(-) diff --git a/himan-scripts/bl-turbulence.lua b/himan-scripts/bl-turbulence.lua index bffefef5f..dbf6e6afc 100644 --- a/himan-scripts/bl-turbulence.lua +++ b/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,6 +12,7 @@ -- 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") @@ -37,26 +34,14 @@ end local wsms = luatool:FetchWithType(current_time, lev, ws, current_forecast_type) -- Wind gust -local gustlevel = {} +local gustlevel = level(HPLevelType.kHeight, 10) local gust = param("FFG-MS") -if currentProducerName == "ECGMTA" then - gustlevel = level(HPLevelType.kGround, 0) -else - gustlevel = level(HPLevelType.kHeight, 10) -end 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") -else - surface = level(HPLevelType.kHeight, 0) - landseamask = param("LC-0TO1") -end +local landseamask = param("LC-0TO1") +local surface = level(HPLevelType.kHeight, 0) local lcmask = luatool:FetchWithType(current_time, surface, landseamask, current_forecast_type) @@ -72,8 +57,10 @@ 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: + +-- 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_0 = luatool:FetchWithType(current_time, lev, U_HIR, current_forecast_type) @@ -171,7 +158,7 @@ for i = 1, #v_0 do if wskt>20 and (wgust>18 or maxshear>25) then turb = 4 end - if wskt>22 and (wgust>18 or maxshear>20) then + if wskt>22 and (wgust>16 or maxshear>20) then turb = 4 end From b36b03c0fc2b1faf04fa34d3f38a5c3f98e745ef Mon Sep 17 00:00:00 2001 From: tackandr Date: Fri, 13 Aug 2021 15:44:22 +0300 Subject: [PATCH 3/3] take Landcover from analysis time step --- himan-plugins/source/luatool.cpp | 1 + himan-scripts/bl-turbulence.lua | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/himan-plugins/source/luatool.cpp b/himan-plugins/source/luatool.cpp index 3b455c67e..d22048dee 100644 --- a/himan-plugins/source/luatool.cpp +++ b/himan-plugins/source/luatool.cpp @@ -1344,6 +1344,7 @@ void BindLib(lua_State* L) .def("Empty", &raw_time::Empty), class_("forecast_time") .def(constructor()) + .def(constructor()) .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)) diff --git a/himan-scripts/bl-turbulence.lua b/himan-scripts/bl-turbulence.lua index dbf6e6afc..bad360333 100644 --- a/himan-scripts/bl-turbulence.lua +++ b/himan-scripts/bl-turbulence.lua @@ -42,8 +42,9 @@ local gustms = luatool:FetchWithType(current_time, gustlevel, gust, current_fore -- 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)) -local lcmask = luatool:FetchWithType(current_time, surface, landseamask, current_forecast_type) +local lcmask = luatool:FetchWithType(mytime, surface, landseamask, current_forecast_type) if not (wsms and gustms and lcmask) then logger:Error("No data found")