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

CBTCU v1.9 #15

Merged
merged 2 commits into from
May 3, 2024
Merged

CBTCU v1.9 #15

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3088ac0
CBTCU v1.9
tackandr Apr 30, 2024
5b1ac55
Apply code review feedback
tackandr May 3, 2024
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
100 changes: 56 additions & 44 deletions himan-scripts/CB-TCU-cloud.lua
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,31 +10,39 @@ local MU = level(HPLevelType.kMaximumThetaE,0)
local HL = level(HPLevelType.kHeightLayer,500,0)
local HG = level(HPLevelType.kHeight,0)

Copy link
Member

Choose a reason for hiding this comment

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

"local" could be used for these variables

Copy link
Contributor Author

Choose a reason for hiding this comment

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

What does "local" actually do?

Copy link
Member

Choose a reason for hiding this comment

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

Well it declares the variable as a local variable instead of global, as you may have guessed :D

Basically it is the way lua recommends to be done: always use "local" unless you have a reason not to. It also has some performance implications as fetching a variable from global scope is more expensive. I have not benchmarked this though and I guess it's rather small, unless running on a very tight loop.

EL500 = luatool:Fetch(current_time, HL, param("EL-LAST-HPA"), current_forecast_type)
LCL500 = luatool:Fetch(current_time, HL, param("LCL-HPA"), current_forecast_type)
LFC500 = luatool:Fetch(current_time, HL, param("LFC-HPA"), current_forecast_type)
CAPE500 = luatool:Fetch(current_time, HL, param("CAPE-JKG"), current_forecast_type)
CIN500 = luatool:Fetch(current_time, HL, param("CIN-JKG"), current_forecast_type)

LCLmu = luatool:Fetch(current_time, MU, param("LCL-HPA"), current_forecast_type)
LFCmu = luatool:Fetch(current_time, MU, param("LFC-HPA"), current_forecast_type)
ELmu = luatool:Fetch(current_time, MU, param("EL-LAST-HPA"), current_forecast_type)
CINmu = luatool:Fetch(current_time, MU, param("CIN-JKG"), current_forecast_type)
CAPEmu = luatool:Fetch(current_time, MU, param("CAPE-JKG"), current_forecast_type)

Ttop = luatool:Fetch(current_time, HL, param("EL-K"), current_forecast_type)
Tbase = luatool:Fetch(current_time, HL, param("LCL-K"), current_forecast_type)
TtopMU = luatool:Fetch(current_time, MU, param("EL-K"), current_forecast_type)
--LFC probably better than LCL for elev conv. base
TbaseMU = luatool:Fetch(current_time, MU, param("LFC-K"), current_forecast_type)

if not EL500 or not LCLmu or not Ttop then
local EL500 = luatool:Fetch(current_time, HL, param("EL-LAST-M"), current_forecast_type)
local pEL500 = luatool:Fetch(current_time, HL, param("EL-LAST-HPA"), current_forecast_type)
local LCL500 = luatool:Fetch(current_time, HL, param("LCL-M"), current_forecast_type)
local CAPE500 = luatool:Fetch(current_time, HL, param("CAPE-JKG"), current_forecast_type)
local CIN500 = luatool:Fetch(current_time, HL, param("CIN-JKG"), current_forecast_type)
local LFCmu = luatool:Fetch(current_time, MU, param("LFC-M"), current_forecast_type)
local pLFCmu = luatool:Fetch(current_time, MU, param("LFC-HPA"), current_forecast_type)
local ELmu = luatool:Fetch(current_time, MU, param("EL-LAST-M"), current_forecast_type)
local pELmu = luatool:Fetch(current_time, MU, param("EL-LAST-HPA"), current_forecast_type)
local CINmu = luatool:Fetch(current_time, MU, param("CIN-JKG"), current_forecast_type)
local CAPEmu = luatool:Fetch(current_time, MU, param("CAPE-JKG"), current_forecast_type)
local Ttop = luatool:Fetch(current_time, HL, param("EL-K"), current_forecast_type)
local TtopMU = luatool:Fetch(current_time, MU, param("EL-K"), current_forecast_type)

if not EL500 or
not pEL500 or
not LCL500 or
not CAPE500 or
not CIN500 or
not LFCmu or
not pLFCmu or
not ELmu or
not pELmu or
not CINmu or
not CAPEmu or
not Ttop or
not TtopMU then
logger:Error("Some data not found")
return
end

NL = luatool:Fetch(current_time, HG, param("NL-PRCNT"), current_forecast_type)
NM = luatool:Fetch(current_time, HG, param("NM-PRCNT"), current_forecast_type)
local NL = luatool:Fetch(current_time, HG, param("NL-PRCNT"), current_forecast_type)
local NM = luatool:Fetch(current_time, HG, param("NM-PRCNT"), current_forecast_type)

if not NL then
NL = luatool:Fetch(current_time, HG, param("NL-0TO1"), current_forecast_type)
Expand All @@ -44,16 +52,19 @@ if not NM then
NM = luatool:Fetch(current_time, HG, param("NM-0TO1"), current_forecast_type)
end

RR = luatool:Fetch(current_time, HG, param("RRR-KGM2"), current_forecast_type)
local RR = luatool:Fetch(current_time, HG, param("RRR-KGM2"), current_forecast_type)

CBlimit = 12 --required vertical thickness [degrees C] to consider a CB (tweak this..!)
TCUlimit = 9 --required vertical thickness [degrees C] to consider a TCU (tweak this..!)
CBtopLim = 263.15 --required top T [K] to consider a CB (tweakable!)
CINlimTCU = -1 --CIN limit for TCu
RRlimit = 0.1 -- precipitation limit [mm/h] to consider a Cb
if not NL or not NM or not RR then
logger:Error("Some data not found")
return
end

--Max height [FL] to check for top
TopLim = 650
local CBlimit = 2000 --required vertical thickness [m] to consider a CB (tweak this..!)
local TCUlimit = 1500 --required vertical thickness [m] to consider a TCU (tweak this..!)
local CBtopLim = -10 --required top T [K] to consider a CB (tweakable!)
local CINlimTCU = -1 --CIN limit for TCu
local RRlimit = 0.1 -- precipitation limit [mm/h] to consider a Cb
local CAPElimit = 2.71828 --euler constant

local i = 0
local res = {}
Expand All @@ -64,10 +75,11 @@ for i=1, #EL500 do
res[i] = Missing

--TCU
if ((Tbase[i]-Ttop[i]>TCUlimit) and (NL[i]>0) and (CIN500[i]>CINlimTCU)) then
res[i] = FlightLevel_(EL500[i]*100)
if ((EL500[i] - LCL500[i] > TCUlimit) and (NL[i] > 0) and (CIN500[i] > CINlimTCU) ) then
--we don't use vertical search for flight level of EL500 but calculate directly from EL500 pressure
res[i] = FlightLevel_(pEL500[i] * 100)
Copy link
Member

@mpartio mpartio May 2, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

This is clever 👍

In thing that caught my eye is that the smarttool function call is:

vertfl_findh(par3_meps,0,TopLim,EL500,1)

The search is capped to TopLim which is set to FL 650. Should this also be checked here, ie.

if res[i] > 650 then 
  res[i] = Missing
end

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I think it is just there because smarttool function requires some value and FL650 is so high that EL500 will always be below. At least in our latitudes. And if there is 20km tall convection I think it's better to put it as it is than to turn it into missing.

--Limit top value
if (CAPE500[i]>math.exp(1)) then
if (CAPE500[i] > CAPElimit) then
--Add for overshooting top based on CAPE
res[i] = -(res[i] + CAPE500[i] / (math.log(CAPE500[i]) * 10))
Copy link
Member

@mpartio mpartio May 2, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

In smarttool the last term (10) has a division operator, not multiplication. I tried to think if you are doing some kind of mathematical trick here but I couldn't get it :D

I don't understand why Simo is using exp(1) here, it's just a constant number.

Copy link
Contributor Author

@tackandr tackandr May 2, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

I put () around it which is the mathematical trick you missed. On a high nerd level I think multiplication is a cheaper operation than division.

exp(1) is used to generate the constant e that LUA doesn't provide. This was needed because we had some problems where CAPE was small and then the logarithm function goes to 0 for ln(1) and division by 0 is bad as well as for values below 1 it gives negative numbers that potentially turns the sign of the result around. Which by the logic of this macro turns CB into TCU. Probably it's not wise to have two things in one parameter and we should output CB and TCU as two separate fields but that's what it is now.

Copy link
Member

Choose a reason for hiding this comment

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

I was just wondering why e(1), does it have some special significance. Why doesn't he just declare "const CAPElimit = 2.78".

I agree it's a bit weird that this parameter has different interpretation depending on the sign of the data value.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah, I think I'll change it to that CAPElimit version then.

else
Expand All @@ -76,33 +88,33 @@ for i=1, #EL500 do
end

--CB
if ((Ttop[i]<CBtopLim) and (Tbase[i]-Ttop[i]>CBlimit) and (RR[i]>RRlimit)) then
res[i] = FlightLevel_(EL500[i]*100)
if ((Ttop[i] < CBtopLim) and (EL500[i] - LCL500[i] > Cblimit) and (RR[i] > RRlimit)) then
res[i] = FlightLevel_(pEL500[i] * 100)
--Limit top value
if (CAPE500[i]>math.exp(1)) then
if (CAPE500[i] > math.exp(1)) then
--Add for overshooting top based on CAPE
res[i] = res[i] + CAPE500[i] / (math.log(CAPE500[i]) * 10)
end
end

--If no TOP from above, check also with MU values, for elev. conv. only from blw 3,5km
if ( IsMissing(res[i]) and (TbaseMU[i]<Tbase[i]) and (LFCmu[i]>650)) then
-- TCU
if ((TbaseMU[i]-TtopMU[i]>TCUlimit) and ((NL[i]>0) or (NM[i]>0)) and (CINmu[i]>CINlimTCU)) then
res[i] = FlightLevel_(ELmu[i]*100)
if ( IsMissing(res[i]) and (LFCmu[i] > LCL500[i]) and (pLFCmu[i] > 650)) then
-- TCU elevated
if ((ELmu[i] - LFCmu[i] > TCUlimit) and ((NL[i] > 0) or (NM[i] > 0)) and (CINmu[i] > CINlimTCU)) then
res[i] = FlightLevel_(pELmu[i] * 100)
--Limit top value
if (CAPEmu[i]>math.exp(1)) then
if (CAPEmu[i] > CAPElimit) then
--Add for overshooting top based on CAPE, +1000ft/350J/kg (tweak this!)
res[i] = -(res[i] + CAPEmu[i] / (math.log(CAPEmu[i]) * 10))
else
res[i] = -res[i]
end
end
--CB
if ((TtopMU[i]<CBtopLim) and (TbaseMU[i]-TtopMU[i]>CBlimit) and (RR[i]>RRlimit)) then
res[i] = FlightLevel_(ELmu[i]*100)
--CB elevated
if ((TtopMU[i] < CBtopLim) and (ELmu[i] - LFCmu[i] > CBlimit) and (RR[i] > RRlimit)) then
res[i] = FlightLevel_(pELmu[i] * 100)
--Limit top value
if (CAPEmu[i]>math.exp(1)) then
if (CAPEmu[i] > CAPElimit) then
--Add for overshooting top based on CAPE, +1000ft/350J/kg (tweak this!)
res[i] = res[i] + CAPEmu[i] / (math.log(CAPEmu[i]) * 10)
end
Expand Down