diff --git a/scripts/simple_met_input/README.txt b/scripts/simple_met_input/README.txt index 12b1f383..80da4589 100644 --- a/scripts/simple_met_input/README.txt +++ b/scripts/simple_met_input/README.txt @@ -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. diff --git a/scripts/simple_met_input/main.py b/scripts/simple_met_input/main.py index 23033f8f..c08e52d9 100644 --- a/scripts/simple_met_input/main.py +++ b/scripts/simple_met_input/main.py @@ -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") diff --git a/scripts/simple_met_input/src/ISA.py b/scripts/simple_met_input/src/ISA.py index b263ae74..e4efd901 100644 --- a/scripts/simple_met_input/src/ISA.py +++ b/scripts/simple_met_input/src/ISA.py @@ -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) diff --git a/scripts/simple_met_input/src/met.py b/scripts/simple_met_input/src/met.py index 3e83cc85..e15c440b 100644 --- a/scripts/simple_met_input/src/met.py +++ b/scripts/simple_met_input/src/met.py @@ -3,13 +3,13 @@ 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. # @@ -17,21 +17,21 @@ def met_clean_slate(desired_altitudes, desired_RHis = None, RHi_background = 0., # 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,22 +140,22 @@ 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) @@ -163,9 +163,8 @@ def create_and_save_met(alts_m, RHis_PC, T_offset_K = 0, shear = 0.002, prefix = 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