Skip to content

Commit

Permalink
Merge pull request #34 from ca525/main
Browse files Browse the repository at this point in the history
Improve code and comment readibility in the met file scripts
sdeastham authored Aug 28, 2024
2 parents 4b68ac8 + 33b3b6c commit 67f9541
Showing 4 changed files with 97 additions and 96 deletions.
2 changes: 2 additions & 0 deletions scripts/simple_met_input/README.txt
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
This folder contains python scripts to produce idealised met files containing trapezoidal moist layers (in an altitude
vs RHi plot).

TO RUN THE "main.py" SCRIPT, PLEASE OPEN THE FOLDER "simple_met_input" IN YOUR EDITOR/SET IT AS YOUR PYTHON RUNDIR.

The trapezoid is symmetrical about its centre, and is parametrised as follows:
- centre_altitude_m: The centre altitude in meters.
- moist_depth_m: The moist layer depth in meters; the vertical extent of the moist layer.
6 changes: 3 additions & 3 deletions scripts/simple_met_input/main.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from src.met import *
from src.met import make_idealised_met
import matplotlib.pyplot as plt

######################### MAIN FUNCTION (EXAMPLE) #########################
@@ -12,13 +12,13 @@
grad_RHi_PC_per_m = None, # % RHi / m; None indicates a rectangular moist region
alt_resolution_m = 100, # m
upper_alt_m = 11001, # m
shear = 2e-3, # 1/s
shear_over_s = 2e-3, # 1/s
T_offset_K = 0, # K
filename_prefix = "example"
)

RHi_PC = met.relative_humidity_ice.to_numpy()
alt_m = met.altitude.to_numpy() * 1e3
alt_m = met.altitude.to_numpy() * 1e3 # km to m

plt.plot(RHi_PC, alt_m, color = 'b', lw = 1.5)
plt.ylabel("Altitude, m")
116 changes: 58 additions & 58 deletions scripts/simple_met_input/src/ISA.py
Original file line number Diff line number Diff line change
@@ -1,21 +1,19 @@
import numpy as np

######
## ISA Equations and constants from https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere
## ISA equations and constants from https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere
######

# TODO: Include humidity in these calculations...

class ISA:
# ISA Class needed to store the constants without making them global

# ISA Temperature Lapse Rate Magnitude (Troposphere)
# ISA temperature lapse rate magnitude at the troposphere
alpha = 6.5e-3 # K / m

# Perfect Gas Constant for Air
# Perfect gas constant for air
R = 287.053 # J / kg K

# Acceleration Due to Gravity
# Acceleration due to gravity
g = 9.80665 # m / s^2

# Sea level quantities
@@ -26,108 +24,110 @@ class ISA:
# Sea level temperature offfset
T_offset = 0 # K

# Quantities at the boundary between the troposphere and the stratosphere (11000m)
h_11 = 11000 # m
p_11 = np.nan # Calculated during the object init
T_11 = np.nan # Calculated during the object init
# Quantities at the tropopause (11000m AMSL)
h_11 = 11000 # m AMSL
p_11 = np.nan # Calculated during the object initialisation
T_11 = np.nan # Calculated during the object initialisation

def __init__(self, T_offset_K = 0):
# Evaluate the quantities at the boundary between the troposphere and the stratosphere (11000m) upon initialisation
# Evaluate the quantities at the tropopause (11000m AMSL) upon initialisation
# according to the International Standard Atmosphere equations from the following link:
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere

self.T_offset = T_offset_K

self.p_11 = self.p_0 * np.power(1 - self.alpha / self.T_0 * self.h_11,
self.g / self.R / self.alpha)
T_11 = self.T_0 - self.alpha * self.h_11 + self.T_offset

def compute_p_from_h(self, h):
# Returns the pressure (in Pa) at a given altitude h (in metres)
# according to the International Standard Atmosphere from the following link:
def compute_p_from_h(self, h_m):
# Returns the pressure (in Pa) at a given altitude h (in metres AMSL)
# according to the International Standard Atmosphere equations from the following link:
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere

p_troposphere = self.p_0 * np.power(1 - self.alpha / self.T_0 * h,
p_troposphere_Pa = self.p_0 * np.power(1 - self.alpha / self.T_0 * h_m,
self.g / self.R / self.alpha)
p_stratosphere = self.p_11 * np.exp(-self.g / self.R / self.T_11 * (h - self.h_11))
p_stratosphere_Pa = self.p_11 * np.exp(-self.g / self.R / self.T_11 * (h_m - self.h_11))

p_out = np.where(h <= self.h_11, p_troposphere, p_stratosphere)
p_out_Pa = np.where(h_m <= self.h_11, p_troposphere_Pa, p_stratosphere_Pa)

# Enforce altitudes beneath SL to be equal to SL values.
return np.where(p_out > self.p_0, self.p_0, p_out)
return np.where(p_out_Pa > self.p_0, self.p_0, p_out_Pa)

def compute_T_from_h(self, h):
# Returns the temperature (in degrees K) at a given altitude h (in metres)
# according to the International Standard Atmosphere from the following link:
def compute_T_from_h(self, h_m):
# Returns the temperature (in degrees K) at a given altitude h (in metres AMSL)
# according to the International Standard Atmosphere equations from the following link:
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere

T_troposphere = self.T_0 - self.alpha * h
T_stratosphere = self.T_11
T_troposphere_K = self.T_0 - self.alpha * h_m
T_stratosphere_K = self.T_11

T_out = np.where(h <= self.h_11, T_troposphere, T_stratosphere)
T_out_K = np.where(h_m <= self.h_11, T_troposphere_K, T_stratosphere_K)

# Enforce altitudes beneath SL to be equal to SL values. Also apply the T offset.
return self.T_offset + np.where(T_out > self.T_0, self.T_0, T_out)
return self.T_offset + np.where(T_out_K > self.T_0, self.T_0, T_out_K)

def compute_rho_from_h(self, h):
# Returns the air density (in kg / m^3) at a given altitude h (in metres)
# according to the International Standard Atmosphere from the following link:
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere
def compute_rho_from_h(self, h_m):
# Returns the air density (in kg / m^3) at a given altitude h (in metres AMSL)
# according to the ideal gas law.
#
# NOTE: NO DISTINCTION IS MADE BETWEEN DRY AND HUMID AIR DENSITY HERE!

p = self.compute_p_from_h(h)
T = self.compute_T_from_h(h)
p_Pa = self.compute_p_from_h(h_m)
T_K = self.compute_T_from_h(h_m)

return p / (self.R * T)
return p_Pa / (self.R * T_K)

def compute_h_from_p(self, p):
# Returns the pressure (in kPa) at a given altitude h (in metres)
# according to the International Standard Atmosphere from the following link:
def compute_h_from_p(self, p_Pa):
# Returns the altitude (in metres AMSL) at a given altitude ambient pressure (in Pa)
# according to the International Standard Atmosphere equations from the following link:
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere

h_troposphere = self.T_0 / self.alpha * ( 1 - np.power(p / self.p_0,
h_troposphere_m = self.T_0 / self.alpha * ( 1 - np.power(p_Pa / self.p_0,
self.R * self.alpha / self.g) )
h_stratosphere = self.h_11 - self.R * self.T_11 / self.g * np.log(p/self.p_11)
h_stratosphere_m = self.h_11 - self.R * self.T_11 / self.g * np.log(p_Pa / self.p_11)

h_out = np.where(p >= self.p_11, h_troposphere, h_stratosphere)
h_out_m = np.where(p_Pa >= self.p_11, h_troposphere_m, h_stratosphere_m)

# Enforce altitudes beneath SL to be equal to the SL.
return np.where(h_out < 0, 0, h_out)
return np.where(h_out_m < 0, 0, h_out_m)

def compute_p_ISA(h):
# Returns the pressure (in Pa) at a given altitude h (in metres)
# according to the International Standard Atmosphere from the following link:
def compute_p_ISA(h_m):
# Returns the pressure (in Pa) at a given altitude h (in metres AMSL)
# according to the International Standard Atmosphere equations from the following link:
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere

ISA_computer = ISA()
return ISA_computer.compute_p_from_h(h)
return ISA_computer.compute_p_from_h(h_m)

def compute_T_ISA(h, T_offset_K = 0):
# Returns the temperature (in degrees K) at a given altitude h (in metres)
# according to the International Standard Atmosphere from the following link:
def compute_T_ISA(h_m, T_offset_K = 0):
# Returns the temperature (in degrees K) at a given altitude h (in metres AMSL)
# according to the International Standard Atmosphere equations from the following link:
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere

ISA_computer = ISA(T_offset_K)
return ISA_computer.compute_T_from_h(h)
return ISA_computer.compute_T_from_h(h_m)

def compute_rho_ISA(h, T_offset_K = 0):
# Returns the air density (in kg / m^3) at a given altitude h (in metres)
# according to the International Standard Atmosphere from the following link:
def compute_rho_ISA(h_m, T_offset_K = 0):
# Returns the air density (in kg / m^3) at a given altitude h (in metres AMSL)
# according to the ideal gas law.
#
# NOTE: NO DISTINCTION IS MADE BETWEEN DRY AND HUMID AIR DENSITY HERE!
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere

ISA_computer = ISA(T_offset_K)
return ISA_computer.compute_rho_from_h(h)
return ISA_computer.compute_rho_from_h(h_m)

def compute_h_from_p_ISA(p):
# Returns the pressure (in kPa) at a given altitude h (in metres)
# according to the International Standard Atmosphere from the following link:
def compute_h_from_p_ISA(p_Pa):
# Returns the altitude (in metres AMSL) at a given altitude ambient pressure (in Pa)
# according to the International Standard Atmosphere equations from the following link:
#
# https://eng.libretexts.org/Bookshelves/Aerospace_Engineering/Fundamentals_of_Aerospace_Engineering_(Arnedo)/02%3A_Generalities/2.03%3A_Standard_atmosphere

ISA_computer = ISA()
return ISA_computer.compute_h_from_p(p)
return ISA_computer.compute_h_from_p(p_Pa)
69 changes: 34 additions & 35 deletions scripts/simple_met_input/src/met.py
Original file line number Diff line number Diff line change
@@ -3,35 +3,35 @@
import matplotlib.pyplot as plt

from pathlib import Path
from src.ISA import *
from src.moist import *
from src.thermo import *
from src.ISA import compute_T_ISA, compute_p_ISA
from src.moist import trapezium_moist_layer
from src.thermo import convert_RHi_to_RH

def met_clean_slate(desired_altitudes, desired_RHis = None, RHi_background = 0., shear = 0.002, T_offset_K = 0):
def met_clean_slate(desired_altitudes_m, desired_RHis_PC = None, RHi_background_PC = 0., shear_over_s = 0.002, T_offset_K = 0):
# Creates a clean slate xr.Dataset object with the weather parameters required by APCEMM at
# each altitude. The vertical resolution is determined by the "desired_altitudes" input.
# each altitude. The vertical resolution is determined by the "desired_altitudes_m" input.
#
# The met created by this function is time-invariant and is only defined for 24 hours.
#
# This function uses the ISA. However, the temperature offset ("T_offset_K") does not affect the
# ISA calculations (which assume zero temperature offset), except for the temperature and density.
#
# The following weather parameters can be changed:
# - RHi: Must be in %. Set by the use of a vertical profile (set in "desired_RHis") or by the
# use of a constant RHi (set in "RHi_background"; only used when "desired_RHis" is "None").
# - Shear: In units of 1/s. A uniform shear field is assumed. (Set by the "shear" input.)
# - RHi: Must be in %. Set by the use of a vertical profile (set in "desired_RHis_PC") or by the
# use of a constant RHi (set in "RHi_background_PC"; only used when "desired_RHis_PC" is "None").
# - Shear: In units of 1/s. A uniform shear field is assumed. (Set by the "shear_over_s" input.)
# - Temperature offset: In units of K. A linear shift in the ISA temperature, can be either
# positive or negative.
#
# Returns a xr.Dataset object.

desired_altitudes = desired_altitudes / 1e3 # m to km
desired_altitudes_km = desired_altitudes_m / 1e3 # m to km

# Time vector, just in case it is needed
time_ip = np.concatenate((np.linspace(15, 23, 9), np.linspace(0, 14, 15))) # hours

# Initialise the met variables
len_alt = len(desired_altitudes)
len_alt = len(desired_altitudes_km)
pressure_ip = np.zeros(len_alt) # hPa
temperature_ip = np.zeros(len_alt) # Degrees Kelvin
relative_humidity_ip = np.zeros(len_alt) # %
@@ -48,37 +48,37 @@ def met_clean_slate(desired_altitudes, desired_RHis = None, RHi_background = 0.,
shear=(["altitude"], shear_ip, {'units':'s^-1'}),
),
coords=dict(
altitude=(["altitude"], desired_altitudes, {'units':'km'}),
altitude=(["altitude"], desired_altitudes_km, {'units':'km'}),
time=(["time"], time_ip, {'units':'h'}),
),
)

if desired_RHis is None:
if desired_RHis_PC is None:
met_data_ISA = met_from_ISA(
met = met_data_temp,
RHi_background = RHi_background,
RHi_background_PC = RHi_background_PC,
T_offset_K = T_offset_K,
shear = shear
shear_over_s = shear_over_s
)
else:
RHi_background = np.min(desired_RHis)
RHi_background_PC = np.min(desired_RHis_PC)
met_data_ISA = met_from_ISA(
met = met_data_temp,
RHi_background = RHi_background,
RHi_background_PC = RHi_background_PC,
T_offset_K = T_offset_K,
shear = shear
shear_over_s = shear_over_s
)
RHi = met_data_ISA.relative_humidity_ice
RH = met_data_ISA.relative_humidity
Ts_K = met_data_ISA.temperature.to_numpy()
desired_RHs = convert_RHi_to_RH(Ts_K, desired_RHis)
desired_RHs = convert_RHi_to_RH(Ts_K, desired_RHis_PC)

# Enforce the desired RHis to each altitude
for i in range(len(desired_altitudes)):
current_alt = desired_altitudes[i]
for i in range(len(desired_altitudes_km)):
current_alt = desired_altitudes_km[i]

mask = RHi.altitude != current_alt
RHi = RHi.where(mask, desired_RHis[i])
RHi = RHi.where(mask, desired_RHis_PC[i])
RH = RH.where(mask, desired_RHs[i])

met_data_ISA["relative_humidity_ice"] = RHi
@@ -89,7 +89,7 @@ def met_clean_slate(desired_altitudes, desired_RHis = None, RHi_background = 0.,
def truncate_met(met, truncation_alt_km = 20.):
# Truncates the APCEMM met data to remain beneath the tropopause.

altitudes = met.altitude
altitudes = met.altitude # km
truncation_idx = len(altitudes)

# Find the index at which the altitude surpasses the truncation threshold
@@ -100,7 +100,7 @@ def truncate_met(met, truncation_alt_km = 20.):

return met.isel(altitude = slice(0, truncation_idx))

def met_from_ISA(met, RHi_background = 0., shear = 0.002, T_offset_K = 0):
def met_from_ISA(met, RHi_background_PC = 0., shear_over_s = 0.002, T_offset_K = 0):
# Clears the existing temperature values of an input xr.Dataset ("met"), and
# enforces ISA conditions at each point.
#
@@ -110,8 +110,8 @@ def met_from_ISA(met, RHi_background = 0., shear = 0.002, T_offset_K = 0):
# assume zero temperature offset), except for the temperature itself.
#
# The following weather parameters can be changed:
# - Background RHi: Must be in %. Set by the use of "RHi_background".
# - Shear: In units of 1/s. A uniform shear field is assumed. (Set by the "shear" input.)
# - Background RHi: Must be in %. Set by the use of "RHi_background_PC".
# - Shear: In units of 1/s. A uniform shear field is assumed. (Set by the "shear_over_s" input.)
# - Temperature offset: In units of K. A shift in the ISA temperature, can be either
# positive or negative. Does not affect the other ISA parameters.
#
@@ -131,7 +131,7 @@ def met_from_ISA(met, RHi_background = 0., shear = 0.002, T_offset_K = 0):

mask = met.altitude != altitudes[i]
T = T.where(mask, T_current)
RH = RH.where(mask, convert_RHi_to_RH(T_current, RHi_background))
RH = RH.where(mask, convert_RHi_to_RH(T_current, RHi_background_PC))
p = p.where(mask, p_current)

met["temperature"] = T
@@ -140,32 +140,31 @@ def met_from_ISA(met, RHi_background = 0., shear = 0.002, T_offset_K = 0):

# Set the RHi of the met data to the RHi default value
RHi_met = met["relative_humidity_ice"]
RHi_met = RHi_met.where(False, RHi_background)
RHi_met = RHi_met.where(False, RHi_background_PC)
met["relative_humidity_ice"] = RHi_met

# Set the shear of the met data to the default value
shear_met = met["shear"]
shear_met = shear_met.where(False, shear)
shear_met = shear_met.where(False, shear_over_s)
met["shear"] = shear_met

return met

def create_and_save_met(alts_m, RHis_PC, T_offset_K = 0, shear = 0.002, prefix = "00"):
def create_and_save_met(alts_m, RHis_PC, T_offset_K = 0, shear_over_s = 0.002, prefix = "00"):
met = met_clean_slate(
desired_altitudes = alts_m,
desired_RHis = RHis_PC,
desired_altitudes_m = alts_m,
desired_RHis_PC = RHis_PC,
T_offset_K = T_offset_K,
shear = shear
shear_over_s = shear_over_s
)

Path("outputs/").mkdir(parents=True, exist_ok=True)
met.to_netcdf('outputs/' + prefix + '-met.nc', engine="netcdf4")
return met

def make_idealised_met(centre_altitude_m, moist_depth_m, RHi_background_PC = 0., RHi_peak_PC = 150,
grad_RHi_PC_per_m = None, alt_resolution_m = 100, upper_alt_m = 11001, shear = 2e-3,
grad_RHi_PC_per_m = None, alt_resolution_m = 100, upper_alt_m = 11001, shear_over_s = 2e-3,
T_offset_K = 0, filename_prefix = "default"):

# Creates an idealised met file with a trapezium moist layer.
#
# The altitude is sampled at intervals with spacing of alt_resolution_m, up until
@@ -183,6 +182,6 @@ def make_idealised_met(centre_altitude_m, moist_depth_m, RHi_background_PC = 0.,
RHi_peak_PC, grad_RHi_PC_per_m, alt_resolution_m,
upper_alt_m)

met = create_and_save_met(alts_m=altitudes,RHis_PC=RHis,T_offset_K=T_offset_K,shear=shear,prefix=filename_prefix)
met = create_and_save_met(alts_m=altitudes,RHis_PC=RHis,T_offset_K=T_offset_K,shear_over_s=shear_over_s,prefix=filename_prefix)

return met

0 comments on commit 67f9541

Please sign in to comment.