From 5fd80015e9bca523ac52f22ca5b4fb662f76c2fb Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 6 Sep 2024 17:18:29 -0600 Subject: [PATCH 01/48] Update to build using scikit-build-core --- pyproject.toml | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 33bf3be94c..04ef962d24 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,9 +1,9 @@ [build-system] -requires = ["setuptools>=61", "wheel", "setuptools_scm[toml]>=3.4"] -build-backend = "setuptools.build_meta" +requires = ["setuptools>=61", "setuptools-scm[toml]>=3.4", "scikit-build-core", "pybind11"] +build-backend = "scikit_build_core.build" [project] -name = "MetPy" +name = "metpy" description = "Collection of tools for reading, visualizing and performing calculations with weather data." readme = "README.md" dynamic = ["version"] @@ -174,5 +174,9 @@ max-complexity = 61 [tool.ruff.lint.pydocstyle] convention = "numpy" +[tool.scikit-build] +metadata.version.provider = "scikit_build_core.metadata.setuptools_scm" +cmake.source-dir = "src" + [tool.setuptools_scm] version_scheme = "post-release" From 19399ff0d99b627c22aae3cbfc02a95a804150a9 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 6 Sep 2024 17:18:57 -0600 Subject: [PATCH 02/48] WIP: Add stub module that uses pybind11 --- src/CMakeLists.txt | 15 +++++++++++++++ src/calcmod.cpp | 10 ++++++++++ 2 files changed, 25 insertions(+) create mode 100644 src/CMakeLists.txt create mode 100644 src/calcmod.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000000..1080051ca0 --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 3.15...3.26) +project( + "${SKBUILD_PROJECT_NAME}" + LANGUAGES CXX + VERSION "${SKBUILD_PROJECT_VERSION}") + +find_package(Python COMPONENTS Interpreter Development.Module REQUIRED) +find_package(pybind11 CONFIG REQUIRED) + +pybind11_add_module(_calc_mod MODULE calcmod.cpp) + +target_compile_definitions(_calc_mod + PRIVATE VERSION_INFO=${PROJECT_VERSION}) + +install(TARGETS _calc_mod DESTINATION metpy) \ No newline at end of file diff --git a/src/calcmod.cpp b/src/calcmod.cpp new file mode 100644 index 0000000000..0d7837d01e --- /dev/null +++ b/src/calcmod.cpp @@ -0,0 +1,10 @@ +#include + +int add(int i, int j) { + return i + j; +} + +PYBIND11_MODULE(_calc_mod, m) { + m.doc() = "accelerator module docstring"; + m.def("add", &add, "Add two numbers"); +} \ No newline at end of file From e76004c4dd61246f8bccc2f7bd3fe03755448394 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 26 May 2025 14:44:51 -0600 Subject: [PATCH 03/48] add virtual_temperature.cpp and link it using pybind11 --- src/CMakeLists.txt | 11 ++++++++-- src/calcmod.cpp | 25 ++++++++++++++++++++-- src/virtual_temperature.cpp | 42 +++++++++++++++++++++++++++++++++++++ src/virtual_temperature.hpp | 24 +++++++++++++++++++++ 4 files changed, 98 insertions(+), 4 deletions(-) create mode 100644 src/virtual_temperature.cpp create mode 100644 src/virtual_temperature.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1080051ca0..59aeba7608 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,12 +4,19 @@ project( LANGUAGES CXX VERSION "${SKBUILD_PROJECT_VERSION}") +# Enforce C++14 or higher +set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + find_package(Python COMPONENTS Interpreter Development.Module REQUIRED) find_package(pybind11 CONFIG REQUIRED) -pybind11_add_module(_calc_mod MODULE calcmod.cpp) +pybind11_add_module(_calc_mod MODULE + calcmod.cpp + virtual_temperature.cpp +) target_compile_definitions(_calc_mod PRIVATE VERSION_INFO=${PROJECT_VERSION}) -install(TARGETS _calc_mod DESTINATION metpy) \ No newline at end of file +install(TARGETS _calc_mod DESTINATION metpy) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 0d7837d01e..3c48600549 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,10 +1,31 @@ #include +#include +#include "virtual_temperature.hpp" + +namespace py = pybind11; int add(int i, int j) { - return i + j; + return i - j; } + PYBIND11_MODULE(_calc_mod, m) { m.doc() = "accelerator module docstring"; + m.def("add", &add, "Add two numbers"); -} \ No newline at end of file + + m.def("virtual_temperature", + py::overload_cast&, const std::vector&, double>(&VirtualTemperature), + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), + "Compute virtual temperature from temperature and mixing ratio arrays"); + + m.def("virtual_temperature", + py::overload_cast&, double>(&VirtualTemperature), + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), + "Compute virtual temperature from scalar temperature and mixing ratio array"); + + m.def("virtual_temperature", + py::overload_cast&, double, double>(&VirtualTemperature), + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), + "Compute virtual temperature from temperature array and scalar mixing ratio"); +} diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp new file mode 100644 index 0000000000..b07236c872 --- /dev/null +++ b/src/virtual_temperature.cpp @@ -0,0 +1,42 @@ +#include "virtual_temperature.hpp" +#include + +std::vector VirtualTemperature( + const std::vector& temperature, + const std::vector& mixing_ratio, + double epsilon +) { + if (temperature.size() != mixing_ratio.size()) { + throw std::invalid_argument("Temperature and mixing ratio vectors must be the same size."); + } + + std::vector result; + result.reserve(temperature.size()); + + for (size_t i = 0; i < temperature.size(); ++i) { + double T = temperature[i]; + double w = mixing_ratio[i]; + result.push_back(T * (w + epsilon) / (epsilon * (1 + w))); + } + + return result; +} + +std::vector VirtualTemperature( + double temperature, + const std::vector& mixing_ratio, + double epsilon +) { + std::vector temperature_vec(mixing_ratio.size(), temperature); + return VirtualTemperature(temperature_vec, mixing_ratio, epsilon); +} + +std::vector VirtualTemperature( + const std::vector& temperature, + double mixing_ratio, + double epsilon +) { + std::vector mixing_ratio_vec(temperature.size(), mixing_ratio); + return VirtualTemperature(temperature, mixing_ratio_vec, epsilon); +} + diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp new file mode 100644 index 0000000000..aaceac9055 --- /dev/null +++ b/src/virtual_temperature.hpp @@ -0,0 +1,24 @@ +#ifndef VIRTUAL_TEMPERATURE_HPP // if not defined +#define VIRTUAL_TEMPERATURE_HPP // define the header file + +#include + +// Compute virtual temperature: vector + vector +std::vector VirtualTemperature( + const std::vector& temperature, + const std::vector& mixing_ratio, + double epsilon = 0.622); + +// Overload: scalar + vector +std::vector VirtualTemperature( + double temperature, + const std::vector& mixing_ratio, + double epsilon = 0.622); + +// Overload: vector + scalar +std::vector VirtualTemperature( + const std::vector& temperature, + double mixing_ratio, + double epsilon = 0.622); + +#endif // VIRTUAL_TEMPERATURE_HPP From 358161ada6f6483ca23bffa8719ebe8edc16793a Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 27 May 2025 11:50:53 -0600 Subject: [PATCH 04/48] debug 1.input data format compatibility 2.set default epsilon 0.622 --- src/calcmod.cpp | 60 +++++++++++++++++++++++++++---------- src/virtual_temperature.cpp | 35 ++++------------------ src/virtual_temperature.hpp | 14 +-------- 3 files changed, 51 insertions(+), 58 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 3c48600549..06d4d6098a 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,31 +1,59 @@ +#include "virtual_temperature.hpp" #include #include -#include "virtual_temperature.hpp" namespace py = pybind11; +// Flexible dispatcher that supports scalar/list combinations +std::vector virtual_temperature_py(py::object temperature, py::object mixing_ratio, double epsilon) { + std::vector temperature_vec; + std::vector mixing_ratio_vec; + + // Case 1: both are lists + if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { + temperature_vec = temperature.cast>(); + mixing_ratio_vec = mixing_ratio.cast>(); + if (temperature_vec.size() != mixing_ratio_vec.size()) { + throw std::invalid_argument("Temperature and mixing ratio lists must be the same length."); + } + } + // Case 2: temperature is float, mixing_ratio is list + else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { + mixing_ratio_vec = mixing_ratio.cast>(); + temperature_vec = std::vector(mixing_ratio_vec.size(), temperature.cast()); + } + // Case 3: temperature is list, mixing_ratio is float + else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { + temperature_vec = temperature.cast>(); + mixing_ratio_vec = std::vector(temperature_vec.size(), mixing_ratio.cast()); + } + // Case 4: both are floats + else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { + temperature_vec = {temperature.cast()}; + mixing_ratio_vec = {mixing_ratio.cast()}; + } + else { + throw std::invalid_argument("Inputs must be float or list."); + } + + return VirtualTemperature(temperature_vec, mixing_ratio_vec, epsilon); +} + int add(int i, int j) { return i - j; } - PYBIND11_MODULE(_calc_mod, m) { m.doc() = "accelerator module docstring"; m.def("add", &add, "Add two numbers"); - m.def("virtual_temperature", - py::overload_cast&, const std::vector&, double>(&VirtualTemperature), - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), - "Compute virtual temperature from temperature and mixing ratio arrays"); - - m.def("virtual_temperature", - py::overload_cast&, double>(&VirtualTemperature), - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), - "Compute virtual temperature from scalar temperature and mixing ratio array"); - - m.def("virtual_temperature", - py::overload_cast&, double, double>(&VirtualTemperature), - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), - "Compute virtual temperature from temperature array and scalar mixing ratio"); + // Unified binding with default epsilon + m.def("virtual_temperature", &virtual_temperature_py, + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622, + "Compute virtual temperature.\n" + "Accepts:\n" + " - two lists of equal length\n" + " - one scalar and one list\n" + "Defaults to epsilon = 0.622"); } diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index b07236c872..592365cce0 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,42 +1,19 @@ #include "virtual_temperature.hpp" -#include +//#include std::vector VirtualTemperature( const std::vector& temperature, const std::vector& mixing_ratio, double epsilon ) { - if (temperature.size() != mixing_ratio.size()) { - throw std::invalid_argument("Temperature and mixing ratio vectors must be the same size."); - } - - std::vector result; - result.reserve(temperature.size()); + std::vector result(temperature.size()); + double T, w; for (size_t i = 0; i < temperature.size(); ++i) { - double T = temperature[i]; - double w = mixing_ratio[i]; - result.push_back(T * (w + epsilon) / (epsilon * (1 + w))); + T = temperature[i]; + w = mixing_ratio[i]; + result[i] = T * (w + epsilon) / (epsilon * (1 + w)); } return result; } - -std::vector VirtualTemperature( - double temperature, - const std::vector& mixing_ratio, - double epsilon -) { - std::vector temperature_vec(mixing_ratio.size(), temperature); - return VirtualTemperature(temperature_vec, mixing_ratio, epsilon); -} - -std::vector VirtualTemperature( - const std::vector& temperature, - double mixing_ratio, - double epsilon -) { - std::vector mixing_ratio_vec(temperature.size(), mixing_ratio); - return VirtualTemperature(temperature, mixing_ratio_vec, epsilon); -} - diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index aaceac9055..2fc91564ef 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -7,18 +7,6 @@ std::vector VirtualTemperature( const std::vector& temperature, const std::vector& mixing_ratio, - double epsilon = 0.622); - -// Overload: scalar + vector -std::vector VirtualTemperature( - double temperature, - const std::vector& mixing_ratio, - double epsilon = 0.622); - -// Overload: vector + scalar -std::vector VirtualTemperature( - const std::vector& temperature, - double mixing_ratio, - double epsilon = 0.622); + double epsilon); #endif // VIRTUAL_TEMPERATURE_HPP From f18d93bd782766c8c11af6d089b8fa25f4903436 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 28 May 2025 15:55:49 -0600 Subject: [PATCH 05/48] remove setting c++14 in CMakeLists.txt --- src/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 59aeba7608..968e79ddc7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,8 +5,8 @@ project( VERSION "${SKBUILD_PROJECT_VERSION}") # Enforce C++14 or higher -set(CMAKE_CXX_STANDARD 14) -set(CMAKE_CXX_STANDARD_REQUIRED ON) +#set(CMAKE_CXX_STANDARD 14) +#set(CMAKE_CXX_STANDARD_REQUIRED ON) find_package(Python COMPONENTS Interpreter Development.Module REQUIRED) find_package(pybind11 CONFIG REQUIRED) From d12c77d9f3493b0cf1cd1627ef56bc713db27376 Mon Sep 17 00:00:00 2001 From: Linfeng Li <99068959+gitlinffff@users.noreply.github.com> Date: Tue, 3 Jun 2025 18:16:22 -0400 Subject: [PATCH 06/48] Pybind11_vectorize (#1) * vectorize the function using pybind11/numpy.h, achieving handling numpy array, float and dtype * add dewpoint, add constants.hpp, adding other functions WIP --- src/calcmod.cpp | 52 +++++++------------------------------ src/constants.hpp | 37 ++++++++++++++++++++++++++ src/virtual_temperature.cpp | 36 +++++++++++++++---------- src/virtual_temperature.hpp | 10 +++---- 4 files changed, 73 insertions(+), 62 deletions(-) create mode 100644 src/constants.hpp diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 06d4d6098a..b32aeb41be 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,46 +1,12 @@ #include "virtual_temperature.hpp" #include +#include #include namespace py = pybind11; -// Flexible dispatcher that supports scalar/list combinations -std::vector virtual_temperature_py(py::object temperature, py::object mixing_ratio, double epsilon) { - std::vector temperature_vec; - std::vector mixing_ratio_vec; - - // Case 1: both are lists - if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = temperature.cast>(); - mixing_ratio_vec = mixing_ratio.cast>(); - if (temperature_vec.size() != mixing_ratio_vec.size()) { - throw std::invalid_argument("Temperature and mixing ratio lists must be the same length."); - } - } - // Case 2: temperature is float, mixing_ratio is list - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - mixing_ratio_vec = mixing_ratio.cast>(); - temperature_vec = std::vector(mixing_ratio_vec.size(), temperature.cast()); - } - // Case 3: temperature is list, mixing_ratio is float - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = temperature.cast>(); - mixing_ratio_vec = std::vector(temperature_vec.size(), mixing_ratio.cast()); - } - // Case 4: both are floats - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = {temperature.cast()}; - mixing_ratio_vec = {mixing_ratio.cast()}; - } - else { - throw std::invalid_argument("Inputs must be float or list."); - } - - return VirtualTemperature(temperature_vec, mixing_ratio_vec, epsilon); -} - int add(int i, int j) { - return i - j; + return i + j; } PYBIND11_MODULE(_calc_mod, m) { @@ -49,11 +15,11 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("add", &add, "Add two numbers"); // Unified binding with default epsilon - m.def("virtual_temperature", &virtual_temperature_py, - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622, - "Compute virtual temperature.\n" - "Accepts:\n" - " - two lists of equal length\n" - " - one scalar and one list\n" - "Defaults to epsilon = 0.622"); + m.def("dewpoint", py::vectorize(DewPoint), + "Calculate dew point from water vapor partial pressure.", + py::arg("vapor_pressure")); + + m.def("virtual_temperature", py::vectorize(VirtualTemperature), + "Calculate virtual temperature from temperature and mixing ratio.", + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622); } diff --git a/src/constants.hpp b/src/constants.hpp new file mode 100644 index 0000000000..8e9d596a61 --- /dev/null +++ b/src/constants.hpp @@ -0,0 +1,37 @@ +#ifndef CONSTANTS_HPP +#define CONSTANTS_HPP + +namespace metpy_constants { +// Gas constants (J / kg / K) +constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) + +// Dry air +constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) +constexpr double Rd = R / Md; // Dry air (J / kg / K) +constexpr double dry_air_spec_heat_ratio = 1.4; +constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) + +// Water +constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) +constexpr double Rv = R / Mw; // Water vapor (J / kg / K) +constexpr double wv_spec_heat_ratio = 1.33; +constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) +constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) +constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) +constexpr double T0 = 273.16; // Triple point of water (K) + +// General Meteorological constants +constexpr double epsilon = Mw / Md; // ≈ 0.622 + + +// Standard gravity +constexpr double g = 9.80665; // m / s^2 + +// Reference pressure +constexpr double P0 = 100000.0; // Pa (hPa = 1000) + +} + +#endif // CONSTANTS_HPP diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 592365cce0..0204a07066 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,19 +1,29 @@ +#include +#include "constants.hpp" #include "virtual_temperature.hpp" //#include -std::vector VirtualTemperature( - const std::vector& temperature, - const std::vector& mixing_ratio, - double epsilon -) { +double water_latent_heat_vaporization(double temperature) { + // Calculate the latent heat of vaporization of water in J/kg at a given temperature. + using namespace metpy_constants; + return Lv - (Cp_l - Cp_v) * (temperature - T0); +} + +double _saturation_vapor_pressure(double temperature) { + // Calculate saturation (equilibrium) water vapor (partial) pressure over liquid water. + // Constants for the Magnus-Tetens approximation + //const double a = 17.67; + //const double b = 243.5; - std::vector result(temperature.size()); - double T, w; - for (size_t i = 0; i < temperature.size(); ++i) { - T = temperature[i]; - w = mixing_ratio[i]; - result[i] = T * (w + epsilon) / (epsilon * (1 + w)); - } + // Calculate saturation vapor pressure using the Magnus-Tetens formula + //return 6.112 * exp((a * temperature) / (b + temperature)); +} + +double DewPoint(double vapor_pressure) { + double val = log(vapor_pressure / 6.112); + return 243.5 * val / (17.67 - val); +} - return result; +double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { + return temperature * (mixing_ratio + epsilon) / (epsilon * (1. + mixing_ratio)); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 2fc91564ef..4b871bebd0 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -1,12 +1,10 @@ #ifndef VIRTUAL_TEMPERATURE_HPP // if not defined #define VIRTUAL_TEMPERATURE_HPP // define the header file -#include +//#include -// Compute virtual temperature: vector + vector -std::vector VirtualTemperature( - const std::vector& temperature, - const std::vector& mixing_ratio, - double epsilon); +double DewPoint(double vapor_pressure); + +double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); #endif // VIRTUAL_TEMPERATURE_HPP From 31dbf18284b5422ff0cc339764c023300492eb95 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 3 Jun 2025 18:16:47 -0600 Subject: [PATCH 07/48] wip --- src/calcmod.cpp | 3 +++ src/constants.cpp | 16 ++++++++++++++++ src/constants.hpp | 34 ++++------------------------------ src/constants_backup | 37 +++++++++++++++++++++++++++++++++++++ src/virtual_temperature.cpp | 14 +++++++++++--- 5 files changed, 71 insertions(+), 33 deletions(-) create mode 100644 src/constants.cpp create mode 100644 src/constants_backup diff --git a/src/calcmod.cpp b/src/calcmod.cpp index b32aeb41be..baf21ddc28 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,3 +1,4 @@ +#include "constants.hpp" #include "virtual_temperature.hpp" #include #include @@ -12,6 +13,8 @@ int add(int i, int j) { PYBIND11_MODULE(_calc_mod, m) { m.doc() = "accelerator module docstring"; + metpy_constants::load_constants_from_python(); + m.def("add", &add, "Add two numbers"); // Unified binding with default epsilon diff --git a/src/constants.cpp b/src/constants.cpp new file mode 100644 index 0000000000..7f54c29d35 --- /dev/null +++ b/src/constants.cpp @@ -0,0 +1,16 @@ +// constants.cpp +#include "constants.hpp" +#include + +namespace py = pybind11; +using namespace metpy_constants; + +// Definitions (required in C++11) +double sat_pressure_0c; + +void load_constants_from_python() { + py::object mod = py::module_::import("metpy.constants.default"); + + sat_pressure_0c = mod.attr("sat_pressure_0c").attr("magnitude").cast(); +} + diff --git a/src/constants.hpp b/src/constants.hpp index 8e9d596a61..fef19add59 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -2,36 +2,10 @@ #define CONSTANTS_HPP namespace metpy_constants { -// Gas constants (J / kg / K) -constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) - -// Dry air -constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) -constexpr double Rd = R / Md; // Dry air (J / kg / K) -constexpr double dry_air_spec_heat_ratio = 1.4; -constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) -constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) - -// Water -constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) -constexpr double Rv = R / Mw; // Water vapor (J / kg / K) -constexpr double wv_spec_heat_ratio = 1.33; -constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) -constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) -constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) -constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) -constexpr double T0 = 273.16; // Triple point of water (K) - -// General Meteorological constants -constexpr double epsilon = Mw / Md; // ≈ 0.622 - - -// Standard gravity -constexpr double g = 9.80665; // m / s^2 - -// Reference pressure -constexpr double P0 = 100000.0; // Pa (hPa = 1000) + extern double sat_pressure_0c; + void load_constants_from_python(); // call once in your PYBIND11_MODULE } -#endif // CONSTANTS_HPP +#endif + diff --git a/src/constants_backup b/src/constants_backup new file mode 100644 index 0000000000..8e9d596a61 --- /dev/null +++ b/src/constants_backup @@ -0,0 +1,37 @@ +#ifndef CONSTANTS_HPP +#define CONSTANTS_HPP + +namespace metpy_constants { +// Gas constants (J / kg / K) +constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) + +// Dry air +constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) +constexpr double Rd = R / Md; // Dry air (J / kg / K) +constexpr double dry_air_spec_heat_ratio = 1.4; +constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) + +// Water +constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) +constexpr double Rv = R / Mw; // Water vapor (J / kg / K) +constexpr double wv_spec_heat_ratio = 1.33; +constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) +constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) +constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) +constexpr double T0 = 273.16; // Triple point of water (K) + +// General Meteorological constants +constexpr double epsilon = Mw / Md; // ≈ 0.622 + + +// Standard gravity +constexpr double g = 9.80665; // m / s^2 + +// Reference pressure +constexpr double P0 = 100000.0; // Pa (hPa = 1000) + +} + +#endif // CONSTANTS_HPP diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 0204a07066..2cc5cccbed 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,12 +1,15 @@ #include +#include #include "constants.hpp" #include "virtual_temperature.hpp" //#include +namespace py = pybind11; + double water_latent_heat_vaporization(double temperature) { // Calculate the latent heat of vaporization of water in J/kg at a given temperature. - using namespace metpy_constants; - return Lv - (Cp_l - Cp_v) * (temperature - T0); +// using namespace metpy_constants; +// return Lv - (Cp_l - Cp_v) * (temperature - T0); } double _saturation_vapor_pressure(double temperature) { @@ -20,7 +23,12 @@ double _saturation_vapor_pressure(double temperature) { } double DewPoint(double vapor_pressure) { - double val = log(vapor_pressure / 6.112); + // fetch constants from python module + //py::object default_mod = py::module_::import("metpy.constants.default"); + // unit issue ignored for now + //double sat_pressure_0c = default_mod.attr("sat_pressure_0c").attr("magnitude").cast(); + + double val = log(vapor_pressure / sat_pressure_0c); return 243.5 * val / (17.67 - val); } From 8e1ecccb44c8f247882def973111a03856f10fd8 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 4 Jun 2025 17:55:31 -0600 Subject: [PATCH 08/48] create constants.cpp and .hpp to fetch all constants from python; converting saturation vapor pressure WIP --- src/CMakeLists.txt | 1 + src/calcmod.cpp | 8 ++++++++ src/constants.cpp | 25 ++++++++++++++++++------- src/constants.hpp | 31 +++++++++++++++++++++++++++++++ src/virtual_temperature.cpp | 25 +++++++++---------------- src/virtual_temperature.hpp | 2 ++ 6 files changed, 69 insertions(+), 23 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 968e79ddc7..c6d61218a2 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -13,6 +13,7 @@ find_package(pybind11 CONFIG REQUIRED) pybind11_add_module(_calc_mod MODULE calcmod.cpp + constants.cpp virtual_temperature.cpp ) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index baf21ddc28..573a1ad35d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -17,6 +17,14 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("add", &add, "Add two numbers"); + m.def("water_latent_heat_vaporization", py::vectorize(WaterLatentHeatVaporization), + "Calculate water latent heat vaporization from temperature.", + py::arg("temperature")); + + m.def("_saturation_vapor_pressure_liquid", py::vectorize(_SaturationVaporPressureLiquid), + "Calculate saturation vapor pressure from temperature.", + py::arg("temperature")); + // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), "Calculate dew point from water vapor partial pressure.", diff --git a/src/constants.cpp b/src/constants.cpp index 7f54c29d35..d99bfd3bba 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -3,14 +3,25 @@ #include namespace py = pybind11; -using namespace metpy_constants; -// Definitions (required in C++11) -double sat_pressure_0c; +namespace metpy_constants { + double Mw; + double Rv; + double Cp_v; + double Cp_l; + double Lv; + double sat_pressure_0c; + double T0; -void load_constants_from_python() { - py::object mod = py::module_::import("metpy.constants.default"); + void load_constants_from_python() { + py::object mod = py::module_::import("metpy.constants.default"); - sat_pressure_0c = mod.attr("sat_pressure_0c").attr("magnitude").cast(); + Mw = mod.attr("Mw").attr("to")("kg / mol").attr("magnitude").cast(); + Rv = mod.attr("Rv").attr("magnitude").cast(); + Cp_v = mod.attr("Cp_v").attr("magnitude").cast(); + Cp_l = mod.attr("Cp_l").attr("magnitude").cast(); + Lv = mod.attr("Lv").attr("magnitude").cast(); + sat_pressure_0c = mod.attr("sat_pressure_0c").attr("magnitude").cast(); + T0 = mod.attr("T0").attr("magnitude").cast(); + } } - diff --git a/src/constants.hpp b/src/constants.hpp index fef19add59..95ad4cb856 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -2,7 +2,38 @@ #define CONSTANTS_HPP namespace metpy_constants { + + // Gas constants (J / kg / K) +// extern double R; + + // Dry air +// extern double Md; +// extern double Rd; +// extern double dry_air_spec_heat_ratio; +// extern double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) +// extern double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) + + // Water + extern double Mw; + extern double Rv; +// extern double wv_spec_heat_ratio = 1.33; + extern double Cp_v; +// extern double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) + extern double Cp_l; + extern double Lv; extern double sat_pressure_0c; + extern double T0; + + // General Meteorological constants +// extern double epsilon = Mw / Md; // ≈ 0.622 + + + // Standard gravity +// extern double g = 9.80665; // m / s^2 + + // Reference pressure +// extern double P0 = 100000.0; // Pa (hPa = 1000) + void load_constants_from_python(); // call once in your PYBIND11_MODULE } diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 2cc5cccbed..b296043144 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -5,30 +5,23 @@ //#include namespace py = pybind11; +namespace mc = metpy_constants; -double water_latent_heat_vaporization(double temperature) { - // Calculate the latent heat of vaporization of water in J/kg at a given temperature. -// using namespace metpy_constants; -// return Lv - (Cp_l - Cp_v) * (temperature - T0); +double WaterLatentHeatVaporization(double temperature) { + return mc::Lv - (mc::Cp_l - mc::Cp_v) * (temperature - mc::T0); } -double _saturation_vapor_pressure(double temperature) { - // Calculate saturation (equilibrium) water vapor (partial) pressure over liquid water. - // Constants for the Magnus-Tetens approximation - //const double a = 17.67; - //const double b = 243.5; +double _SaturationVaporPressureLiquid(double temperature) { + double latent_heat = WaterLatentHeatVaporization(temperature); + double heat_power = (mc::Cp_l - mc::Cp_v) / mc::Rv; + double exp_term = (mc::Lv / mc::T0 - latent_heat / temperature) / mc::Rv; - // Calculate saturation vapor pressure using the Magnus-Tetens formula - //return 6.112 * exp((a * temperature) / (b + temperature)); + return mc::sat_pressure_0c * exp(exp_term) * pow(mc::T0 / temperature, heat_power); } double DewPoint(double vapor_pressure) { - // fetch constants from python module - //py::object default_mod = py::module_::import("metpy.constants.default"); - // unit issue ignored for now - //double sat_pressure_0c = default_mod.attr("sat_pressure_0c").attr("magnitude").cast(); - double val = log(vapor_pressure / sat_pressure_0c); + double val = log(vapor_pressure / mc::sat_pressure_0c); return 243.5 * val / (17.67 - val); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 4b871bebd0..230b3737f6 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -2,7 +2,9 @@ #define VIRTUAL_TEMPERATURE_HPP // define the header file //#include +double WaterLatentHeatVaporization(double temperature); +double _SaturationVaporPressureLiquid(double temperature); double DewPoint(double vapor_pressure); double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); From 6b81f0e89468e6dc6db079238898ff01e35b0fac Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 5 Jun 2025 16:17:53 -0600 Subject: [PATCH 09/48] c++ fetch nounit constants from nounit.py --- src/constants.cpp | 25 +++++++++++++++++-------- src/metpy/constants/nounit.py | 1 + src/virtual_temperature.cpp | 2 +- 3 files changed, 19 insertions(+), 9 deletions(-) diff --git a/src/constants.cpp b/src/constants.cpp index d99bfd3bba..55ebbbe656 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -14,14 +14,23 @@ namespace metpy_constants { double T0; void load_constants_from_python() { - py::object mod = py::module_::import("metpy.constants.default"); + py::object mod = py::module_::import("metpy.constants.nounit"); - Mw = mod.attr("Mw").attr("to")("kg / mol").attr("magnitude").cast(); - Rv = mod.attr("Rv").attr("magnitude").cast(); - Cp_v = mod.attr("Cp_v").attr("magnitude").cast(); - Cp_l = mod.attr("Cp_l").attr("magnitude").cast(); - Lv = mod.attr("Lv").attr("magnitude").cast(); - sat_pressure_0c = mod.attr("sat_pressure_0c").attr("magnitude").cast(); - T0 = mod.attr("T0").attr("magnitude").cast(); +// Mw = mod.attr("Mw").attr("magnitude").cast(); +// Rv = mod.attr("Rv").attr("magnitude").cast(); +// Cp_v = mod.attr("Cp_v").attr("magnitude").cast(); +// Cp_l = mod.attr("Cp_l").attr("magnitude").cast(); +// Lv = mod.attr("Lv").attr("magnitude").cast(); +// sat_pressure_0c = mod.attr("sat_pressure_0c").cast(); +// T0 = mod.attr("T0").attr("magnitude").cast(); + + + Mw = mod.attr("Mw").cast(); + Rv = mod.attr("Rv").cast(); + Cp_v = mod.attr("Cp_v").cast(); + Cp_l = mod.attr("Cp_l").cast(); + Lv = mod.attr("Lv").cast(); + sat_pressure_0c = mod.attr("sat_pressure_0c").cast(); + T0 = mod.attr("T0").cast(); } } diff --git a/src/metpy/constants/nounit.py b/src/metpy/constants/nounit.py index 2fd6e8a78f..d1b5e83730 100644 --- a/src/metpy/constants/nounit.py +++ b/src/metpy/constants/nounit.py @@ -6,6 +6,7 @@ from . import default from ..units import units +Mw = default.Mw.m_as('kg / mol') Rd = default.Rd.m_as('m**2 / K / s**2') Rv = default.Rv.m_as('m**2 / K / s**2') Lv = default.Lv.m_as('m**2 / s**2') diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index b296043144..273e1858cd 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -22,7 +22,7 @@ double _SaturationVaporPressureLiquid(double temperature) { double DewPoint(double vapor_pressure) { double val = log(vapor_pressure / mc::sat_pressure_0c); - return 243.5 * val / (17.67 - val); + return 243.5 * val / (17.67 - val); // use SI units instead } double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { From 0cb07ba3462b829aa97730407320475d0565c260 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 6 Jun 2025 11:15:17 -0600 Subject: [PATCH 10/48] finish converting saturation vapor pressure --- src/calcmod.cpp | 10 +++++++++- src/constants.cpp | 4 ++++ src/constants.hpp | 2 ++ src/virtual_temperature.cpp | 32 +++++++++++++++++++++++++++++++- src/virtual_temperature.hpp | 6 ++++++ 5 files changed, 52 insertions(+), 2 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 573a1ad35d..2732d1b1cc 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -20,11 +20,19 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("water_latent_heat_vaporization", py::vectorize(WaterLatentHeatVaporization), "Calculate water latent heat vaporization from temperature.", py::arg("temperature")); - + + m.def("water_latent_heat_sublimation", py::vectorize(WaterLatentHeatSublimation), + "Calculate water latent heat sublimation from temperature.", + py::arg("temperature")); + m.def("_saturation_vapor_pressure_liquid", py::vectorize(_SaturationVaporPressureLiquid), "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); + m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), + "Calculate saturation vapor pressure from temperature.", + py::arg("temperature"), py::arg("phase") = "liquid"); + // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), "Calculate dew point from water vapor partial pressure.", diff --git a/src/constants.cpp b/src/constants.cpp index 55ebbbe656..978eb93f9c 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -12,6 +12,8 @@ namespace metpy_constants { double Lv; double sat_pressure_0c; double T0; + double Ls; + double Cp_i; void load_constants_from_python() { py::object mod = py::module_::import("metpy.constants.nounit"); @@ -32,5 +34,7 @@ namespace metpy_constants { Lv = mod.attr("Lv").cast(); sat_pressure_0c = mod.attr("sat_pressure_0c").cast(); T0 = mod.attr("T0").cast(); + Ls = mod.attr("Ls").cast(); + Cp_i = mod.attr("Cp_i").cast(); } } diff --git a/src/constants.hpp b/src/constants.hpp index 95ad4cb856..08d76e07e8 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -23,6 +23,8 @@ namespace metpy_constants { extern double Lv; extern double sat_pressure_0c; extern double T0; + extern double Ls; + extern double Cp_i; // General Meteorological constants // extern double epsilon = Mw / Md; // ≈ 0.622 diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 273e1858cd..1c20b8670c 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,8 +1,9 @@ #include +#include #include #include "constants.hpp" #include "virtual_temperature.hpp" -//#include +#include namespace py = pybind11; namespace mc = metpy_constants; @@ -11,6 +12,10 @@ double WaterLatentHeatVaporization(double temperature) { return mc::Lv - (mc::Cp_l - mc::Cp_v) * (temperature - mc::T0); } +double WaterLatentHeatSublimation(double temperature) { + return mc::Ls - (mc::Cp_i - mc::Cp_v) * (temperature - mc::T0); +} + double _SaturationVaporPressureLiquid(double temperature) { double latent_heat = WaterLatentHeatVaporization(temperature); double heat_power = (mc::Cp_l - mc::Cp_v) / mc::Rv; @@ -19,6 +24,31 @@ double _SaturationVaporPressureLiquid(double temperature) { return mc::sat_pressure_0c * exp(exp_term) * pow(mc::T0 / temperature, heat_power); } +double _SaturationVaporPressureSolid(double temperature) { + double latent_heat = WaterLatentHeatSublimation(temperature); + double heat_power = (mc::Cp_i - mc::Cp_v) / mc::Rv; + double exp_term = (mc::Ls / mc::T0 - latent_heat / temperature) / mc::Rv; + + return mc::sat_pressure_0c * exp(exp_term) * pow(mc::T0 / temperature, heat_power); +} + +double SaturationVaporPressure(double temperature, std::string phase) { + if (phase == "liquid") { + return _SaturationVaporPressureLiquid(temperature); + } else if (phase == "solid") { + return _SaturationVaporPressureSolid(temperature); + } else if (phase == "auto") { + if (temperature > mc::T0) { + return _SaturationVaporPressureLiquid(temperature); + } else { + return _SaturationVaporPressureSolid(temperature); + } + } else { + throw std::invalid_argument("'" + phase + "' is not a valid option for phase. " + "Valid options are 'liquid', 'solid', or 'auto'."); + } +} + double DewPoint(double vapor_pressure) { double val = log(vapor_pressure / mc::sat_pressure_0c); diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 230b3737f6..d8182c2d0c 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -2,9 +2,15 @@ #define VIRTUAL_TEMPERATURE_HPP // define the header file //#include +#include + double WaterLatentHeatVaporization(double temperature); +double WaterLatentHeatSublimation(double temperature); double _SaturationVaporPressureLiquid(double temperature); +double _SaturationVaporPressureSolid(double temperature); +double SaturationVaporPressure(double temperature, std::string phase); + double DewPoint(double vapor_pressure); double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); From 5514401f542347feb2779abb1c9e9d5c4980bdd8 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 6 Jun 2025 15:08:06 -0600 Subject: [PATCH 11/48] calc.thermo calls c++ function WIP --- src/calcmod.cpp | 2 +- src/metpy/calc/thermo.py | 12 +++++++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 2732d1b1cc..c3a6ad866d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -40,5 +40,5 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("virtual_temperature", py::vectorize(VirtualTemperature), "Calculate virtual temperature from temperature and mixing ratio.", - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622); + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon")); } diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 40274192b2..d1ee282828 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -16,6 +16,7 @@ import scipy.optimize as so from scipy.special import lambertw import xarray as xr +import metpy._calc_mod as _calc_mod from .exceptions import InvalidSoundingError from .tools import (_greater_or_close, _less_or_close, _remove_nans, find_bounding_indices, @@ -2123,9 +2124,14 @@ def virtual_temperature( Renamed ``mixing`` parameter to ``mixing_ratio`` """ - return temperature * ((mixing_ratio + molecular_weight_ratio) - / (molecular_weight_ratio * (1 + mixing_ratio))) - + print(temperature) + print(type(temperature)) + exit() + T_mag = temperature.to('K').magnitude + w_mag = mixing_ratio.to('').magnitude + + return _calc_mod.virtual_temperature( + T_mag, w_mag, molecular_weight_ratio) * temperature.units @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', From 697afd6eccb8ea1f87d9e8f1f95f14e398d154a9 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 6 Jun 2025 16:23:55 -0600 Subject: [PATCH 12/48] apply c++ functions in thermo.py --- src/calcmod.cpp | 4 ++-- src/metpy/calc/thermo.py | 47 +++++++++------------------------------- 2 files changed, 12 insertions(+), 39 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index c3a6ad866d..d773b9fb85 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -29,9 +29,9 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); - m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), + m.def("_saturation_vapor_pressure_solid", py::vectorize(_SaturationVaporPressureSolid), "Calculate saturation vapor pressure from temperature.", - py::arg("temperature"), py::arg("phase") = "liquid"); + py::arg("temperature")); // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index d1ee282828..566b7d6629 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -186,10 +186,8 @@ def water_latent_heat_vaporization(temperature): Eq 15, [Ambaum2020]_, using MetPy-defined constants in place of cited values. """ - return (mpconsts.nounit.Lv - - (mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v) - * (temperature - mpconsts.nounit.T0)) - + # Calling c++ calculation module + return _calc_mod.water_latent_heat_vaporization(temperature) @exporter.export @preprocess_and_wrap(wrap_like='temperature') @@ -228,10 +226,8 @@ def water_latent_heat_sublimation(temperature): Eq 18, [Ambaum2020]_, using MetPy-defined constants in place of cited values. """ - return (mpconsts.nounit.Ls - - (mpconsts.nounit.Cp_i - mpconsts.nounit.Cp_v) - * (temperature - mpconsts.nounit.T0)) - + # Calling c++ calculation module + return _calc_mod.water_latent_heat_sublimation(temperature) @exporter.export @preprocess_and_wrap(wrap_like='temperature') @@ -1624,17 +1620,8 @@ def _saturation_vapor_pressure_liquid(temperature): .. math:: e = e_{s0} \frac{T_0}{T}^{(c_{pl} - c_{pv}) / R_v} \exp{ \frac{L_0}{R_v T_0} - \frac{L}{R_v T}} """ - latent_heat = water_latent_heat_vaporization._nounit(temperature) - heat_power = (mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v) / mpconsts.nounit.Rv - exp_term = ((mpconsts.nounit.Lv / mpconsts.nounit.T0 - latent_heat / temperature) - / mpconsts.nounit.Rv) - - return ( - mpconsts.nounit.sat_pressure_0c - * (mpconsts.nounit.T0 / temperature) ** heat_power - * np.exp(exp_term) - ) - + # Calling c++ calculation module + return _calc_mod._saturation_vapor_pressure_liquid(temperature) @preprocess_and_wrap(wrap_like='temperature') @process_units({'temperature': '[temperature]'}, '[pressure]') @@ -1661,17 +1648,8 @@ def _saturation_vapor_pressure_solid(temperature): .. math:: e_i = e_{i0} \frac{T_0}{T}^{(c_{pi} - c_{pv}) / R_v} \exp{ \frac{L_{s0}}{R_v T_0} - \frac{L_s}{R_v T}} """ - latent_heat = water_latent_heat_sublimation._nounit(temperature) - heat_power = (mpconsts.nounit.Cp_i - mpconsts.nounit.Cp_v) / mpconsts.nounit.Rv - exp_term = ((mpconsts.nounit.Ls / mpconsts.nounit.T0 - latent_heat / temperature) - / mpconsts.nounit.Rv) - - return ( - mpconsts.nounit.sat_pressure_0c - * (mpconsts.nounit.T0 / temperature) ** heat_power - * np.exp(exp_term) - ) - + # Calling c++ calculation module + return _calc_mod._saturation_vapor_pressure_solid(temperature) @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'relative_humidity')) @@ -2124,14 +2102,9 @@ def virtual_temperature( Renamed ``mixing`` parameter to ``mixing_ratio`` """ - print(temperature) - print(type(temperature)) - exit() - T_mag = temperature.to('K').magnitude - w_mag = mixing_ratio.to('').magnitude - + # Calling c++ calculation module return _calc_mod.virtual_temperature( - T_mag, w_mag, molecular_weight_ratio) * temperature.units + temperature, mixing_ratio, molecular_weight_ratio) @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', From 5072b95bb1e30a3472a29d61b43ec3480c56d1e3 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 6 Jun 2025 16:28:55 -0600 Subject: [PATCH 13/48] wip --- src/calcmod.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index d773b9fb85..7d4f51e94c 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -32,6 +32,10 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("_saturation_vapor_pressure_solid", py::vectorize(_SaturationVaporPressureSolid), "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); + + m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), + "Calculate saturation vapor pressure from temperature.", + py::arg("temperature"), py::arg("phase") = "liquid"); // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), From 40c94190be70681f336dcd3f1a78ce971677ae8c Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 9 Jun 2025 09:34:55 -0600 Subject: [PATCH 14/48] wip --- src/constants.cpp | 2 ++ src/constants.hpp | 1 + src/metpy/calc/thermo.py | 3 +-- src/virtual_temperature.cpp | 2 +- 4 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/constants.cpp b/src/constants.cpp index 978eb93f9c..1a64786a30 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -14,6 +14,7 @@ namespace metpy_constants { double T0; double Ls; double Cp_i; + double zero_degc; void load_constants_from_python() { py::object mod = py::module_::import("metpy.constants.nounit"); @@ -36,5 +37,6 @@ namespace metpy_constants { T0 = mod.attr("T0").cast(); Ls = mod.attr("Ls").cast(); Cp_i = mod.attr("Cp_i").cast(); + zero_degc = mod.attr("zero_degc").cast(); } } diff --git a/src/constants.hpp b/src/constants.hpp index 08d76e07e8..4d2f0bc7a4 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -25,6 +25,7 @@ namespace metpy_constants { extern double T0; extern double Ls; extern double Cp_i; + extern double zero_degc; // General Meteorological constants // extern double epsilon = Mw / Md; // ≈ 0.622 diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 566b7d6629..0b0dc786f8 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1729,9 +1729,8 @@ def dewpoint(vapor_pressure): Renamed ``e`` parameter to ``vapor_pressure`` """ - val = np.log(vapor_pressure / mpconsts.nounit.sat_pressure_0c) - return mpconsts.nounit.zero_degc + 243.5 * val / (17.67 - val) + return _calc_mod.dewpoint(vapor_pressure) @exporter.export @preprocess_and_wrap(wrap_like='partial_press', broadcast=('partial_press', 'total_press')) diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 1c20b8670c..29a94a8a3e 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -52,7 +52,7 @@ double SaturationVaporPressure(double temperature, std::string phase) { double DewPoint(double vapor_pressure) { double val = log(vapor_pressure / mc::sat_pressure_0c); - return 243.5 * val / (17.67 - val); // use SI units instead + return mc::zero_degc + 243.5 * val / (17.67 - val); } double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { From f913ea7ab587ae5830399817e763df688f934ef8 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 9 Jun 2025 15:47:40 -0600 Subject: [PATCH 15/48] wip --- src/calcmod.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 7d4f51e94c..7f3491da79 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -25,6 +25,10 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate water latent heat sublimation from temperature.", py::arg("temperature")); + m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), + "Calculate saturation vapor pressure from temperature.", + py::arg("temperature"), py::arg("phase") = "liquid"); + m.def("_saturation_vapor_pressure_liquid", py::vectorize(_SaturationVaporPressureLiquid), "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); @@ -33,11 +37,6 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); - m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), - "Calculate saturation vapor pressure from temperature.", - py::arg("temperature"), py::arg("phase") = "liquid"); - - // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), "Calculate dew point from water vapor partial pressure.", py::arg("vapor_pressure")); From e7d870ef9f09c9622a1aff2c9a4f7cdcf2f438fb Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 11 Jun 2025 16:27:44 -0600 Subject: [PATCH 16/48] add c++ functions: mixing_ratio, saturation_mixing_ratio, specific_humidity_from_mixing_ratio, specific_humidity_from_dewpoint --- src/calcmod.cpp | 17 +++++++++++++++++ src/constants.cpp | 2 ++ src/constants.hpp | 1 + src/metpy/calc/thermo.py | 5 +++-- src/virtual_temperature.cpp | 25 +++++++++++++++++++++++++ src/virtual_temperature.hpp | 7 +++++++ 6 files changed, 55 insertions(+), 2 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 7f3491da79..7fe49686b7 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -41,6 +41,23 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate dew point from water vapor partial pressure.", py::arg("vapor_pressure")); + m.def("mixing_ratio", py::vectorize(MixingRatio), + "Calculate the mixing ratio of a gas.", + py::arg("partial_press"), py::arg("total_press"), py::arg("epsilon")); + + m.def("saturation_mixing_ratio", py::vectorize(SaturationMixingRatio), + "Calculate the saturation mixing ratio of water vapor given total atmospheric pressure and temperature.", + py::arg("total_press"), py::arg("temperature"), py::arg("phase")); + + m.def("specific_humidity_from_mixing_ratio", py::vectorize(SpecificHumidityFromMixingRatio), + "Calculate the specific humidity from the mixing ratio.", + py::arg("mixing_ratio")); + + m.def("specific_humidity_from_dewpoint", py::vectorize(SpecificHumidityFromDewPoint), + "Calculate the specific humidity from the dewpoint temperature and pressure.", + py::arg("pressure"), py::arg("dewpoint"), py::arg("phase")); + + m.def("virtual_temperature", py::vectorize(VirtualTemperature), "Calculate virtual temperature from temperature and mixing ratio.", py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon")); diff --git a/src/constants.cpp b/src/constants.cpp index 1a64786a30..a83a7ed737 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -15,6 +15,7 @@ namespace metpy_constants { double Ls; double Cp_i; double zero_degc; + double epsilon; void load_constants_from_python() { py::object mod = py::module_::import("metpy.constants.nounit"); @@ -38,5 +39,6 @@ namespace metpy_constants { Ls = mod.attr("Ls").cast(); Cp_i = mod.attr("Cp_i").cast(); zero_degc = mod.attr("zero_degc").cast(); + epsilon = mod.attr("epsilon").cast(); } } diff --git a/src/constants.hpp b/src/constants.hpp index 4d2f0bc7a4..4707c7a046 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -26,6 +26,7 @@ namespace metpy_constants { extern double Ls; extern double Cp_i; extern double zero_degc; + extern double epsilon; // General Meteorological constants // extern double epsilon = Mw / Md; // ≈ 0.622 diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 0b0dc786f8..e36796cc50 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1729,7 +1729,7 @@ def dewpoint(vapor_pressure): Renamed ``e`` parameter to ``vapor_pressure`` """ - + # Calling c++ calculation module return _calc_mod.dewpoint(vapor_pressure) @exporter.export @@ -1790,7 +1790,8 @@ def mixing_ratio(partial_press, total_press, molecular_weight_ratio=mpconsts.nou Renamed ``part_press``, ``tot_press`` parameters to ``partial_press``, ``total_press`` """ - return molecular_weight_ratio * partial_press / (total_press - partial_press) + # Calling c++ calculation module + return _calc_mod.mixing_ratio(partial_press, total_press, molecular_weight_ratio) @exporter.export diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 29a94a8a3e..6255455a45 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -4,6 +4,8 @@ #include "constants.hpp" #include "virtual_temperature.hpp" #include +#include // for std::cerr +#include // for std::numeric_limits namespace py = pybind11; namespace mc = metpy_constants; @@ -55,6 +57,29 @@ double DewPoint(double vapor_pressure) { return mc::zero_degc + 243.5 * val / (17.67 - val); } +double MixingRatio(double partial_press, double total_press, double epsilon) { + return epsilon * partial_press / (total_press - partial_press); +} + +double SaturationMixingRatio(double total_press, double temperature, std::string phase) { + double e_s = SaturationVaporPressure(temperature, phase); + if (e_s >= total_press) { + std::cerr << "Total pressure must be greater than the saturation vapor pressure " + << "for liquid water to be in equilibrium.\n"; + return std::numeric_limits::quiet_NaN(); + } + return MixingRatio(e_s, total_press); +} + +double SpecificHumidityFromMixingRatio(double mixing_ratio) { + return mixing_ratio / (mixing_ratio + 1.0); +} + +double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::string phase) { + double mixing_ratio = SaturationMixingRatio(pressure, dew_point, phase); + return SpecificHumidityFromMixingRatio(mixing_ratio); +} + double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { return temperature * (mixing_ratio + epsilon) / (epsilon * (1. + mixing_ratio)); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index d8182c2d0c..32c6de0915 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -3,6 +3,9 @@ //#include #include +#include "constants.hpp" + +namespace mc = metpy_constants; double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); @@ -12,6 +15,10 @@ double _SaturationVaporPressureSolid(double temperature); double SaturationVaporPressure(double temperature, std::string phase); double DewPoint(double vapor_pressure); +double MixingRatio(double partial_press, double total_press, double epsilon=mc::epsilon); +double SaturationMixingRatio(double total_press, double temperature, std::string phase); +double SpecificHumidityFromMixingRatio(double mixing_ratio); +double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::string phase); double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); From ad18ce35fdc67a08b46c065a338a07b52f759bdf Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 11 Jun 2025 17:50:46 -0600 Subject: [PATCH 17/48] add c++ functions: moist_air_gas_constant, moist_air_specific_heat_pressure; WIP: LCL --- src/calcmod.cpp | 12 ++++++++++++ src/constants.cpp | 4 ++++ src/constants.hpp | 2 ++ src/virtual_temperature.cpp | 28 +++++++++++++++++++++++++--- src/virtual_temperature.hpp | 5 +++++ 5 files changed, 48 insertions(+), 3 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 7fe49686b7..16b4b13f7c 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -16,6 +16,14 @@ PYBIND11_MODULE(_calc_mod, m) { metpy_constants::load_constants_from_python(); m.def("add", &add, "Add two numbers"); + + m.def("moist_air_gas_constant", py::vectorize(MoistAirGasConstant), + "Calculate R_m, the gas constant for moist air.", + py::arg("specific_humidity")); + + m.def("moist_air_specific_heat_pressure", py::vectorize(MoistAirSpecificHeatPressure), + "Calculate C_pm, the specific heat of moist air at constant pressure.", + py::arg("specific_humidity")); m.def("water_latent_heat_vaporization", py::vectorize(WaterLatentHeatVaporization), "Calculate water latent heat vaporization from temperature.", @@ -24,6 +32,10 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("water_latent_heat_sublimation", py::vectorize(WaterLatentHeatSublimation), "Calculate water latent heat sublimation from temperature.", py::arg("temperature")); + + m.def("lcl", py::vectorize(LCL), + "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), "Calculate saturation vapor pressure from temperature.", diff --git a/src/constants.cpp b/src/constants.cpp index a83a7ed737..92c6327bea 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -6,7 +6,9 @@ namespace py = pybind11; namespace metpy_constants { double Mw; + double Rd; double Rv; + double Cp_d; double Cp_v; double Cp_l; double Lv; @@ -30,7 +32,9 @@ namespace metpy_constants { Mw = mod.attr("Mw").cast(); + Rd = mod.attr("Rd").cast(); Rv = mod.attr("Rv").cast(); + Cp_d = mod.attr("Cp_d").cast(); Cp_v = mod.attr("Cp_v").cast(); Cp_l = mod.attr("Cp_l").cast(); Lv = mod.attr("Lv").cast(); diff --git a/src/constants.hpp b/src/constants.hpp index 4707c7a046..a1c182623a 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -15,8 +15,10 @@ namespace metpy_constants { // Water extern double Mw; + extern double Rd; extern double Rv; // extern double wv_spec_heat_ratio = 1.33; + extern double Cp_d; extern double Cp_v; // extern double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) extern double Cp_l; diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 6255455a45..63e188e2af 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -10,6 +10,14 @@ namespace py = pybind11; namespace mc = metpy_constants; +double MoistAirGasConstant(double specific_humidity) { + return (1.0 - specific_humidity) * mc::Rd + specific_humidity * mc::Rv; +} + +double MoistAirSpecificHeatPressure(double specific_humidity) { + return (1.0 - specific_humidity) * mc::Cp_d + specific_humidity * mc::Cp_v; +} + double WaterLatentHeatVaporization(double temperature) { return mc::Lv - (mc::Cp_l - mc::Cp_v) * (temperature - mc::T0); } @@ -18,6 +26,21 @@ double WaterLatentHeatSublimation(double temperature) { return mc::Ls - (mc::Cp_i - mc::Cp_v) * (temperature - mc::T0); } +double LCL(double pressure, double temperature, double dewpoint) { + if (temperature <= dewpoint) { + std::cerr << "Temperature must be greater than dew point for LCL calculation.\n"; + return std::numeric_limits::quiet_NaN(); + } + + double q = SpecificHumidityFromDewPoint(pressure, dewpoint, "liquid"); + double moist_heat_ratio = MoistAirSpecificHeatPressure(q) / MoistAirGasConstant(q); + double spec_heat_diff = mc::Cp_l - mc::Cp_v; + + double a = moist_heat_ratio + spec_heat_diff / mc::Rv; + + return ; +} + double _SaturationVaporPressureLiquid(double temperature) { double latent_heat = WaterLatentHeatVaporization(temperature); double heat_power = (mc::Cp_l - mc::Cp_v) / mc::Rv; @@ -52,7 +75,6 @@ double SaturationVaporPressure(double temperature, std::string phase) { } double DewPoint(double vapor_pressure) { - double val = log(vapor_pressure / mc::sat_pressure_0c); return mc::zero_degc + 243.5 * val / (17.67 - val); } @@ -75,8 +97,8 @@ double SpecificHumidityFromMixingRatio(double mixing_ratio) { return mixing_ratio / (mixing_ratio + 1.0); } -double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::string phase) { - double mixing_ratio = SaturationMixingRatio(pressure, dew_point, phase); +double SpecificHumidityFromDewPoint(double pressure, double dewpoint, std::string phase) { + double mixing_ratio = SaturationMixingRatio(pressure, dewpoint, phase); return SpecificHumidityFromMixingRatio(mixing_ratio); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 32c6de0915..3352012e1c 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -7,9 +7,14 @@ namespace mc = metpy_constants; +double MoistAirGasConstant(double specific_humidity); +double MoistAirSpecificHeatPressure(double specific_humidity); + double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); +double LCL(double pressure, double temperature, double dewpoint); + double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); double SaturationVaporPressure(double temperature, std::string phase); From 89a17abd57124b9cee2e5d680840c12ca9c99604 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 11 Jun 2025 17:52:35 -0600 Subject: [PATCH 18/48] wip --- src/calcmod.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 16b4b13f7c..87c79c954e 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -50,7 +50,7 @@ PYBIND11_MODULE(_calc_mod, m) { py::arg("temperature")); m.def("dewpoint", py::vectorize(DewPoint), - "Calculate dew point from water vapor partial pressure.", + "Calculate dewpoint from water vapor partial pressure.", py::arg("vapor_pressure")); m.def("mixing_ratio", py::vectorize(MixingRatio), From 0af355b0657ca587008a335116522da8c8a3c028 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 12 Jun 2025 11:15:30 -0600 Subject: [PATCH 19/48] add function lambert_wm1: Lambert W function -1 section; LCL c++ WIP --- src/CMakeLists.txt | 1 + src/calcmod.cpp | 6 +++++- src/constants_backup | 37 ------------------------------------- src/math.cpp | 27 +++++++++++++++++++++++++++ src/math.hpp | 7 +++++++ src/virtual_temperature.cpp | 16 +++++++++++++++- src/virtual_temperature.hpp | 2 ++ 7 files changed, 57 insertions(+), 39 deletions(-) delete mode 100644 src/constants_backup create mode 100644 src/math.cpp create mode 100644 src/math.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c6d61218a2..aa59e45f65 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,6 +14,7 @@ find_package(pybind11 CONFIG REQUIRED) pybind11_add_module(_calc_mod MODULE calcmod.cpp constants.cpp + math.cpp virtual_temperature.cpp ) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 87c79c954e..3de2918ca9 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -32,7 +32,11 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("water_latent_heat_sublimation", py::vectorize(WaterLatentHeatSublimation), "Calculate water latent heat sublimation from temperature.", py::arg("temperature")); - + + m.def("relative_humidity_from_dewpoint", py::vectorize(RelativeHumidityFromDewPoint), + "Calculate relative humidity from temperature and dewpoint.", + py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); + m.def("lcl", py::vectorize(LCL), "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); diff --git a/src/constants_backup b/src/constants_backup deleted file mode 100644 index 8e9d596a61..0000000000 --- a/src/constants_backup +++ /dev/null @@ -1,37 +0,0 @@ -#ifndef CONSTANTS_HPP -#define CONSTANTS_HPP - -namespace metpy_constants { -// Gas constants (J / kg / K) -constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) - -// Dry air -constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) -constexpr double Rd = R / Md; // Dry air (J / kg / K) -constexpr double dry_air_spec_heat_ratio = 1.4; -constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) -constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) - -// Water -constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) -constexpr double Rv = R / Mw; // Water vapor (J / kg / K) -constexpr double wv_spec_heat_ratio = 1.33; -constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) -constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) -constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) -constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) -constexpr double T0 = 273.16; // Triple point of water (K) - -// General Meteorological constants -constexpr double epsilon = Mw / Md; // ≈ 0.622 - - -// Standard gravity -constexpr double g = 9.80665; // m / s^2 - -// Reference pressure -constexpr double P0 = 100000.0; // Pa (hPa = 1000) - -} - -#endif // CONSTANTS_HPP diff --git a/src/math.cpp b/src/math.cpp new file mode 100644 index 0000000000..816cca8f47 --- /dev/null +++ b/src/math.cpp @@ -0,0 +1,27 @@ +#include +#include +#include "math.hpp" + +double lambert_wm1(double x, double tol, int max_iter) { + if (x >= 0 || x < -1.0 / std::exp(1.0)) { + throw std::domain_error("W₋₁(x) is only defined for -1/e < x < 0"); + } + + // Initial guess for W₋₁, from Corless et al. (1996) // need to check if this is + // correct + double L1 = std::log(-x); + double L2 = std::log(-L1); + double w = L1 - L2 + (L2 / L1); + + for (int i = 0; i < max_iter; ++i) { + double ew = std::exp(w); + double wew = w * ew; + double diff = (wew - x) / (ew * (w + 1) - ((w + 2) * (wew - x)) / (2 * w + 2)); + w -= diff; + if (std::abs(diff) < tol) { + return w; + } + } + + throw std::runtime_error("lambert_wm1 did not converge"); +} diff --git a/src/math.hpp b/src/math.hpp new file mode 100644 index 0000000000..c4c99ea255 --- /dev/null +++ b/src/math.hpp @@ -0,0 +1,7 @@ +#ifndef MATH_HPP +#define MATH_HPP + +double lambert_wm1(double x, double tol = 1e-10, int max_iter = 100); + +#endif // MATH_HPP + diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 63e188e2af..5fbbab35c9 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,6 +1,7 @@ #include #include #include +#include "math.hpp" #include "constants.hpp" #include "virtual_temperature.hpp" #include @@ -26,6 +27,12 @@ double WaterLatentHeatSublimation(double temperature) { return mc::Ls - (mc::Cp_i - mc::Cp_v) * (temperature - mc::T0); } +double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase) { + double e_s = SaturationVaporPressure(temperature, phase); + double e = SaturationVaporPressure(dewpoint, phase); + return e / e_s; +} + double LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Temperature must be greater than dew point for LCL calculation.\n"; @@ -37,8 +44,15 @@ double LCL(double pressure, double temperature, double dewpoint) { double spec_heat_diff = mc::Cp_l - mc::Cp_v; double a = moist_heat_ratio + spec_heat_diff / mc::Rv; + double b = -(mc::Lv + spec_heat_diff * mc::T0) / (mc::Rv * temperature); + double c = b / a; + + double rh = RelativeHumidityFromDewPoint(temperature, dewpoint, "liquid"); + double w_minus1 = lambert_wm1(pow(rh, 1.0 / a) * c * exp(c)); + double t_lcl = c / w_minus1 * temperature; + double p_lcl = pressure * pow(t_lcl / temperature, moist_heat_ratio); - return ; + return t_lcl; // returning t_lcl and p_lcl together is needed } double _SaturationVaporPressureLiquid(double temperature) { diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 3352012e1c..c360006c5b 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -13,6 +13,8 @@ double MoistAirSpecificHeatPressure(double specific_humidity); double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); +double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase="liquid"); + double LCL(double pressure, double temperature, double dewpoint); double _SaturationVaporPressureLiquid(double temperature); From 94803bd748c754e9ade0af275afb7fde8b4c9284 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Sun, 15 Jun 2025 15:38:37 -0600 Subject: [PATCH 20/48] achieved manually vectorization of function LCL, support np.array of multiple dimensions --- src/calcmod.cpp | 64 +++++++++++++++++++++++++++++++++++-- src/virtual_temperature.cpp | 7 ++-- src/virtual_temperature.hpp | 4 +-- 3 files changed, 67 insertions(+), 8 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 3de2918ca9..a91f8b2c3d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -3,6 +3,8 @@ #include #include #include +#include // For std::pair +#include // For std::make_tuple namespace py = pybind11; @@ -37,9 +39,65 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate relative humidity from temperature and dewpoint.", py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); - m.def("lcl", py::vectorize(LCL), - "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", - py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); +// m.def("lcl", py::vectorize(LCL), +// "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", +// py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); + + m.def("lcl", [](py::array_t pressure, + py::array_t temperature, + py::array_t dewpoint) { + + // This helper ensures the arrays are in C-style contiguous memory. + // If an input array is already contiguous, it's a zero-cost operation. + // If it's a slice or has a different memory layout, it creates a copy. + // This makes the subsequent looping simple and safe. + auto p_contig = py::array::ensure(pressure, py::array::c_style); + auto t_contig = py::array::ensure(temperature, py::array::c_style); + auto d_contig = py::array::ensure(dewpoint, py::array::c_style); + + // --- Step 1: Check that all input arrays have the same shape --- + if (p_contig.ndim() != t_contig.ndim() || p_contig.ndim() != d_contig.ndim()) { + throw std::runtime_error("Input arrays must have the same number of dimensions."); + } + for (int i = 0; i < p_contig.ndim(); ++i) { + if (p_contig.shape(i) != t_contig.shape(i) || p_contig.shape(i) != d_contig.shape(i)) { + throw std::runtime_error("Input arrays must have the same shape."); + } + } + + // --- Step 2: Create output arrays with the exact same N-D shape as the inputs --- + auto p_lcl = py::array_t(p_contig.request().shape); + auto t_lcl = py::array_t(p_contig.request().shape); + + // --- Step 3: Get the total number of elements to loop over --- + ssize_t size = p_contig.size(); + + // --- Step 4: Get direct pointers to the (now contiguous) data buffers --- + const double* p_ptr = static_cast(p_contig.request().ptr); + const double* t_ptr = static_cast(t_contig.request().ptr); + const double* d_ptr = static_cast(d_contig.request().ptr); + double* p_lcl_ptr = p_lcl.mutable_data(); + double* t_lcl_ptr = t_lcl.mutable_data(); + + // --- Step 5: Loop through the data as if it were a single flat 1D array --- + for (ssize_t i = 0; i < size; i++) { + // Call the scalar c++ function for each element + std::pair result = LCL(p_ptr[i], t_ptr[i], d_ptr[i]); + + p_lcl_ptr[i] = result.first; + t_lcl_ptr[i] = result.second; + } + + // --- Step 6: Return a tuple of the two new, N-dimensional arrays --- + return std::make_tuple(p_lcl, t_lcl); + + }, "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); + + + + + m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), "Calculate saturation vapor pressure from temperature.", diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 5fbbab35c9..d523cbbb79 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,5 +1,6 @@ #include #include +#include // For std::pair #include #include "math.hpp" #include "constants.hpp" @@ -33,10 +34,10 @@ double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::st return e / e_s; } -double LCL(double pressure, double temperature, double dewpoint) { +std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Temperature must be greater than dew point for LCL calculation.\n"; - return std::numeric_limits::quiet_NaN(); + return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; } double q = SpecificHumidityFromDewPoint(pressure, dewpoint, "liquid"); @@ -52,7 +53,7 @@ double LCL(double pressure, double temperature, double dewpoint) { double t_lcl = c / w_minus1 * temperature; double p_lcl = pressure * pow(t_lcl / temperature, moist_heat_ratio); - return t_lcl; // returning t_lcl and p_lcl together is needed + return {p_lcl, t_lcl}; // returning t_lcl and p_lcl together is needed } double _SaturationVaporPressureLiquid(double temperature) { diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index c360006c5b..51488eeabb 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -1,8 +1,8 @@ #ifndef VIRTUAL_TEMPERATURE_HPP // if not defined #define VIRTUAL_TEMPERATURE_HPP // define the header file -//#include #include +#include // For std::pair #include "constants.hpp" namespace mc = metpy_constants; @@ -15,7 +15,7 @@ double WaterLatentHeatSublimation(double temperature); double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase="liquid"); -double LCL(double pressure, double temperature, double dewpoint); +std::pair LCL(double pressure, double temperature, double dewpoint); double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); From 8f50f44a5605aa1d38bbd06533b4a07e12568a97 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 16 Jun 2025 09:41:54 -0600 Subject: [PATCH 21/48] add reference to Lambert W function calculation --- src/math.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/math.cpp b/src/math.cpp index 816cca8f47..5b89a5f4cc 100644 --- a/src/math.cpp +++ b/src/math.cpp @@ -3,12 +3,13 @@ #include "math.hpp" double lambert_wm1(double x, double tol, int max_iter) { + // Corless, R.M., Gonnet, G.H., Hare, D.E.G. et al. On the LambertW function. + // Adv Comput Math 5, 329–359 (1996). https://doi.org/10.1007/BF02124750 + if (x >= 0 || x < -1.0 / std::exp(1.0)) { throw std::domain_error("W₋₁(x) is only defined for -1/e < x < 0"); } - // Initial guess for W₋₁, from Corless et al. (1996) // need to check if this is - // correct double L1 = std::log(-x); double L2 = std::log(-L1); double w = L1 - L2 + (L2 / L1); From 76ce4c2952c3aba6b2f971f53e6f14de61a5e011 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 16 Jun 2025 16:52:51 -0600 Subject: [PATCH 22/48] set warning messages and return NAN value --- src/math.cpp | 10 +++++++--- src/virtual_temperature.cpp | 8 +++++--- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/math.cpp b/src/math.cpp index 5b89a5f4cc..2c0dc4f39d 100644 --- a/src/math.cpp +++ b/src/math.cpp @@ -1,5 +1,5 @@ #include -#include +#include // for std::cerr #include "math.hpp" double lambert_wm1(double x, double tol, int max_iter) { @@ -7,7 +7,9 @@ double lambert_wm1(double x, double tol, int max_iter) { // Adv Comput Math 5, 329–359 (1996). https://doi.org/10.1007/BF02124750 if (x >= 0 || x < -1.0 / std::exp(1.0)) { - throw std::domain_error("W₋₁(x) is only defined for -1/e < x < 0"); + std::cerr << "Warning in function '" << __func__ + << "': lambert_wm1 is only defined for -1/e < x < 0\n"; + return std::numeric_limits::quiet_NaN(); } double L1 = std::log(-x); @@ -24,5 +26,7 @@ double lambert_wm1(double x, double tol, int max_iter) { } } - throw std::runtime_error("lambert_wm1 did not converge"); + std::cerr << "Warning in function '" << __func__ + << "': lambert_wm1 did not converge.\n"; + return std::numeric_limits::quiet_NaN(); } diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index d523cbbb79..5222b14bf0 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -36,7 +36,8 @@ double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::st std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { - std::cerr << "Temperature must be greater than dew point for LCL calculation.\n"; + std::cerr << "Warning in function '" << __func__ + << "': Temperature must be greater than dew point for LCL calculation.\n"; return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; } @@ -101,8 +102,9 @@ double MixingRatio(double partial_press, double total_press, double epsilon) { double SaturationMixingRatio(double total_press, double temperature, std::string phase) { double e_s = SaturationVaporPressure(temperature, phase); if (e_s >= total_press) { - std::cerr << "Total pressure must be greater than the saturation vapor pressure " - << "for liquid water to be in equilibrium.\n"; + std::cerr << "Warning in function '" << __func__ + << "': Total pressure must be greater than the saturation vapor pressure " + << "for liquid water to be in equilibrium.\n"; return std::numeric_limits::quiet_NaN(); } return MixingRatio(e_s, total_press); From 097844add4b2fd36f7f5b00bb3051589a8a8636f Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 18 Jun 2025 17:21:15 -0600 Subject: [PATCH 23/48] add c++ dry_lapse --- src/calcmod.cpp | 6 +++--- src/constants.cpp | 2 ++ src/constants.hpp | 1 + src/metpy/calc/thermo.py | 24 ++++++++++++++++++++++++ src/virtual_temperature.cpp | 35 ++++++++++++++++++++++++++++++++--- src/virtual_temperature.hpp | 4 ++++ 6 files changed, 66 insertions(+), 6 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index a91f8b2c3d..24e007f239 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -39,9 +39,9 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate relative humidity from temperature and dewpoint.", py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); -// m.def("lcl", py::vectorize(LCL), -// "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", -// py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); + m.def("dry_lapse", py::vectorize(DryLapse), + "Calculate the temperature at a level assuming only dry adiabatic process.", + py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); m.def("lcl", [](py::array_t pressure, py::array_t temperature, diff --git a/src/constants.cpp b/src/constants.cpp index 92c6327bea..6cc7bed626 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -18,6 +18,7 @@ namespace metpy_constants { double Cp_i; double zero_degc; double epsilon; + double kappa; void load_constants_from_python() { py::object mod = py::module_::import("metpy.constants.nounit"); @@ -44,5 +45,6 @@ namespace metpy_constants { Cp_i = mod.attr("Cp_i").cast(); zero_degc = mod.attr("zero_degc").cast(); epsilon = mod.attr("epsilon").cast(); + kappa = mod.attr("kappa").cast(); } } diff --git a/src/constants.hpp b/src/constants.hpp index a1c182623a..02271e28f6 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -29,6 +29,7 @@ namespace metpy_constants { extern double Cp_i; extern double zero_degc; extern double epsilon; + extern double kappa; // General Meteorological constants // extern double epsilon = Mw / Md; // ≈ 0.622 diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index e36796cc50..a1bc96913e 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -513,6 +513,30 @@ def dry_lapse(pressure, temperature, reference_pressure=None, vertical_dim=0): return temperature * (pressure / reference_pressure)**mpconsts.kappa + +@exporter.export +@preprocess_and_wrap( + wrap_like='temperature', + broadcast=('pressure', 'temperature', 'reference_pressure') +) +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'reference_pressure': '[pressure]' + }, + '[temperature]' +) +def dry_lapse_linfel(pressure, temperature, reference_pressure=None, vertical_dim=0): + """ + Linfeng's version of 'dry_lapse'. Added on Jun18 2025 + """ + if reference_pressure is None: + reference_pressure = pressure[0] + return _calc_mod.dry_lapse(pressure, temperature, reference_pressure) + + + @exporter.export @preprocess_and_wrap( wrap_like='temperature', diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 5222b14bf0..12d835019f 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,13 +1,14 @@ #include #include +#include #include // For std::pair #include -#include "math.hpp" -#include "constants.hpp" -#include "virtual_temperature.hpp" #include #include // for std::cerr #include // for std::numeric_limits +#include "math.hpp" +#include "constants.hpp" +#include "virtual_temperature.hpp" namespace py = pybind11; namespace mc = metpy_constants; @@ -34,6 +35,12 @@ double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::st return e / e_s; } +double DryLapse(double pressure, double ref_temperature, double ref_pressure) { + // calculate temperature at pressure given reference temperature and pressure + // assuming dry adiabatic process + return ref_temperature * pow(pressure / ref_pressure, mc::kappa); +} + std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Warning in function '" << __func__ @@ -57,6 +64,28 @@ std::pair LCL(double pressure, double temperature, double dewpoi return {p_lcl, t_lcl}; // returning t_lcl and p_lcl together is needed } +bool _CheckPressure(const std::vector& pressure) { + for (size_t i = 0; i + 1 < pressure.size(); ++i) { + if (pressure[i] < pressure[i + 1]) { + return false; // pressure increases (invalid) + } + } + return true; // strictly non-increasing +} + +void _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { + // Check that pressure does not increase. + if (!_CheckPressure(pressure)) { + throw std::runtime_error( + "Pressure increases between at least two points in your sounding. " + "Using a smoothing filter (e.g., scipy.signal.medfilt) may fix this."); + } + + // Find the LCL + //press_lcl, temp_lcl = LCL(pressure[0], temperature, dewpoint); + +} + double _SaturationVaporPressureLiquid(double temperature) { double latent_heat = WaterLatentHeatVaporization(temperature); double heat_power = (mc::Cp_l - mc::Cp_v) / mc::Rv; diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 51488eeabb..8620dee62a 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -2,6 +2,7 @@ #define VIRTUAL_TEMPERATURE_HPP // define the header file #include +#include #include // For std::pair #include "constants.hpp" @@ -14,9 +15,12 @@ double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase="liquid"); +double DryLapse(double pressure, double ref_temperature, double ref_pressure); std::pair LCL(double pressure, double temperature, double dewpoint); +bool _CheckPressure(const std::vector& pressure); + double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); double SaturationVaporPressure(double temperature, std::string phase); From 94fde62b07404f355bc487ba4b92293dae77996e Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 18 Jun 2025 18:10:11 -0600 Subject: [PATCH 24/48] add DryLapseProfile for c++ internal vectorization of DryLapse; _ParcelProfileHelper wip --- src/virtual_temperature.cpp | 25 +++++++++++++++++++++++-- src/virtual_temperature.hpp | 4 ++++ 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 12d835019f..0d846413bb 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -41,6 +41,22 @@ double DryLapse(double pressure, double ref_temperature, double ref_pressure) { return ref_temperature * pow(pressure / ref_pressure, mc::kappa); } +std::vector DryLapseProfile(const std::vector& pressure_profile, + double ref_temperature, + double ref_pressure) { + // Vectorized version of DryLapse for a pressure profile. C++ internally use. + + // calculate temperature profile of an air parcel lifting dry adiabatically + // through the given pressure profile. + std::vector temperature_profile; + temperature_profile.reserve(pressure_profile.size()); + + for (double p : pressure_profile) { + temperature_profile.push_back(DryLapse(p, ref_temperature, ref_pressure)); + } + return temperature_profile; +} + std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Warning in function '" << __func__ @@ -61,7 +77,7 @@ std::pair LCL(double pressure, double temperature, double dewpoi double t_lcl = c / w_minus1 * temperature; double p_lcl = pressure * pow(t_lcl / temperature, moist_heat_ratio); - return {p_lcl, t_lcl}; // returning t_lcl and p_lcl together is needed + return {p_lcl, t_lcl}; } bool _CheckPressure(const std::vector& pressure) { @@ -82,7 +98,12 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur } // Find the LCL - //press_lcl, temp_lcl = LCL(pressure[0], temperature, dewpoint); + std::pair result = LCL(pressure[0], temperature, dewpoint); + double press_lcl = result.first; + double temp_lcl = result.second; + + std::vector press_lower; + std::vector temp_lower = DryLapseProfile(press_lower, temperature, press_lower[0]); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 8620dee62a..2c28fcf283 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -15,7 +15,11 @@ double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase="liquid"); + double DryLapse(double pressure, double ref_temperature, double ref_pressure); +std::vector DryLapseProfile(const std::vector& pressure_profile, + double ref_temperature, + double ref_pressure); std::pair LCL(double pressure, double temperature, double dewpoint); From e9feef389308af479cd8f2fd4d05c0a360df9ccb Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 20 Jun 2025 22:28:59 -0600 Subject: [PATCH 25/48] add moist lapse, use rk4 to integrate dlnT_dlnP --- src/calcmod.cpp | 6 ++- src/virtual_temperature.cpp | 83 +++++++++++++++++++++++++++++++++++-- src/virtual_temperature.hpp | 7 ++++ 3 files changed, 92 insertions(+), 4 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 24e007f239..c430380f39 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -40,9 +40,13 @@ PYBIND11_MODULE(_calc_mod, m) { py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); m.def("dry_lapse", py::vectorize(DryLapse), - "Calculate the temperature at a level assuming only dry adiabatic process.", + "Calculate the temperature at pressure levels assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); + m.def("moist_lapse", py::vectorize(MoistLapse), + "Calculate the temperature at pressure levels assuming saturated adiabatic process.", + py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); + m.def("lcl", [](py::array_t pressure, py::array_t temperature, py::array_t dewpoint) { diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 0d846413bb..a967c66926 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -57,6 +57,61 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, return temperature_profile; } +double CaldlnTdlnP(double temperature, double pressure) { + // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. + double rs = SaturationMixingRatio(pressure, temperature, "liquid"); + double dlnT_dlnP_linfel = (mc::Rd + rs * mc::Rv) / (mc::Cp_d + rs * mc::Cp_v + + (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); + + double dlnT_dlnP_Bakhshaii2013 = (mc::Rd + mc::Lv * rs / temperature) / (mc::Cp_d + + (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); + + return dlnT_dlnP_Bakhshaii2013; +} + +double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep) { + // calculate temperature at pressure given reference temperature and pressure + // assuming moist adiabatic expansion (vapor condenses and removed from the air + // parcel) + + double dlnP = log(pressure / ref_pressure) / (double)nstep; + double T1 = ref_temperature; + double P1 = ref_pressure; + double k[4]; + + for (int i = 0; i < nstep; ++i) { + k[0] = CaldlnTdlnP(T1, P1); + k[1] = CaldlnTdlnP(T1 * exp(k[0] * dlnP/2.), P1 * exp(dlnP/2.)); + k[2] = CaldlnTdlnP(T1 * exp(k[1] * dlnP/2.), P1 * exp(dlnP/2.)); + k[3] = CaldlnTdlnP(T1 * exp(k[2] * dlnP), P1 * exp(dlnP)); + + T1 = T1 * exp((k[0] + 2.0 * k[1] + 2.0 * k[2] + k[3]) * dlnP / 6.0); + P1 = P1 * exp(dlnP); + } + + return T1; // check final T1 P1 +} + +std::vector MoistLapseProfile(const std::vector& pressure_profile, + double ref_temperature, + double ref_pressure) { + // Vectorized version of MoistLapse for a pressure profile. C++ internally use. + + // calculate temperature profile of an air parcel lifting saturated adiabatically + // through the given pressure profile. + std::vector temperature_profile; + temperature_profile.reserve(pressure_profile.size()); + + //double T1 = ref_temperature; + //double P1 = ref_pressure; + double T; + for (size_t i = 0; i < pressure_profile.size(); ++i) { + T = MoistLapse(pressure_profile[i], ref_temperature, ref_pressure, 50); + temperature_profile.push_back(T); + } + return temperature_profile; +} + std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Warning in function '" << __func__ @@ -98,13 +153,35 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur } // Find the LCL - std::pair result = LCL(pressure[0], temperature, dewpoint); - double press_lcl = result.first; - double temp_lcl = result.second; + auto [press_lcl, temp_lcl] = LCL(pressure[0], temperature, dewpoint); + // Establish profile below LCL std::vector press_lower; + for (double p : pressure) { + if (p >= press_lcl) { + press_lower.push_back(p); + } + } + press_lower.push_back(press_lcl); std::vector temp_lower = DryLapseProfile(press_lower, temperature, press_lower[0]); + // Early return if profile ends before reaching above LCL + if (pressure.back() >= press_lcl) { + press_lower.pop_back(); + temp_lower.pop_back(); +// return {press_lower, {}, press_lcl, temp_lower, {}, temp_lcl}; + } + + // Establish profile above LCL + std::vector press_upper; + press_upper.push_back(press_lcl); + for (double p : pressure) { + if (p < press_lcl) { + press_upper.push_back(p); + } + } + + } double _SaturationVaporPressureLiquid(double temperature) { diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 2c28fcf283..e362246d7a 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -21,6 +21,13 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, double ref_temperature, double ref_pressure); +double CaldlnTdlnP(double temperature, double pressure); +double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep); +std::vector MoistLapseProfile(const std::vector& pressure_profile, + double ref_temperature, + double ref_pressure); + + std::pair LCL(double pressure, double temperature, double dewpoint); bool _CheckPressure(const std::vector& pressure); From 569c92f2352b2a83eea3379bddf1709e21f0dfac Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 13:38:53 -0600 Subject: [PATCH 26/48] use Bakhshaii2013 for c++ moist_lapse --- src/metpy/calc/thermo.py | 24 ++++++++++++++++++++++++ src/virtual_temperature.cpp | 6 +++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a1bc96913e..3859fb1b05 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -675,6 +675,30 @@ def dt(p, t): return ret.squeeze() +@exporter.export +@preprocess_and_wrap( + wrap_like='temperature', + broadcast=('pressure', 'temperature', 'reference_pressure') +) +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'reference_pressure': '[pressure]' + }, + '[temperature]' +) +def moist_lapse_linfel(pressure, temperature, reference_pressure=None): + """ + Linfeng's version of 'moist_lapse'. Added on Jun23 2025 + """ + if reference_pressure is None: + reference_pressure = pressure[0] + # nstep for RK4 is set to 30 + return _calc_mod.moist_lapse(pressure, temperature, reference_pressure, 30) + + + @exporter.export @preprocess_and_wrap() @process_units( diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index a967c66926..cfcaea4178 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -60,9 +60,9 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, double CaldlnTdlnP(double temperature, double pressure) { // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. double rs = SaturationMixingRatio(pressure, temperature, "liquid"); - double dlnT_dlnP_linfel = (mc::Rd + rs * mc::Rv) / (mc::Cp_d + rs * mc::Cp_v + - (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); - + + //double dlnT_dlnP_linfel = (mc::Rd + rs * mc::Rv) / (mc::Cp_d + rs * mc::Cp_v + + // (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); double dlnT_dlnP_Bakhshaii2013 = (mc::Rd + mc::Lv * rs / temperature) / (mc::Cp_d + (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); From ca222aa26e0a02a4318a416d91e845f724c2dbd7 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 16:02:24 -0600 Subject: [PATCH 27/48] import c++ lcl into metpy thermo.py --- src/metpy/calc/thermo.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 3859fb1b05..f9fafbc431 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -784,6 +784,21 @@ def lcl(pressure, temperature, dewpoint, max_iters=None, eps=None): return p_lcl, t_lcl +@exporter.export +@preprocess_and_wrap() +@process_units( + {'pressure': '[pressure]', 'temperature': '[temperature]', 'dewpoint': '[temperature]'}, + ('[pressure]', '[temperature]') +) +def lcl_linfel(pressure, temperature, dewpoint, max_iters=None, eps=None): + """ + Linfeng's version of 'lcl'. Added on Jun23 2025 + """ + p_lcl, t_lcl = _calc_mod.lcl(pressure, temperature, dewpoint) + return p_lcl, t_lcl + + + @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]') From 7986a8acf931ca929cd0e5b13da642e3b8ce75d4 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 16:36:53 -0600 Subject: [PATCH 28/48] change file name 'virtual_temperature.cpp' -> 'thermo.cpp' --- src/CMakeLists.txt | 2 +- src/calcmod.cpp | 2 +- src/{virtual_temperature.cpp => thermo.cpp} | 2 +- src/{virtual_temperature.hpp => thermo.hpp} | 6 +++--- 4 files changed, 6 insertions(+), 6 deletions(-) rename src/{virtual_temperature.cpp => thermo.cpp} (99%) rename src/{virtual_temperature.hpp => thermo.hpp} (93%) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index aa59e45f65..d9be3d2855 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -15,7 +15,7 @@ pybind11_add_module(_calc_mod MODULE calcmod.cpp constants.cpp math.cpp - virtual_temperature.cpp + thermo.cpp ) target_compile_definitions(_calc_mod diff --git a/src/calcmod.cpp b/src/calcmod.cpp index c430380f39..5ffd42afaa 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,5 +1,5 @@ #include "constants.hpp" -#include "virtual_temperature.hpp" +#include "thermo.hpp" #include #include #include diff --git a/src/virtual_temperature.cpp b/src/thermo.cpp similarity index 99% rename from src/virtual_temperature.cpp rename to src/thermo.cpp index cfcaea4178..dbcb70be13 100644 --- a/src/virtual_temperature.cpp +++ b/src/thermo.cpp @@ -8,7 +8,7 @@ #include // for std::numeric_limits #include "math.hpp" #include "constants.hpp" -#include "virtual_temperature.hpp" +#include "thermo.hpp" namespace py = pybind11; namespace mc = metpy_constants; diff --git a/src/virtual_temperature.hpp b/src/thermo.hpp similarity index 93% rename from src/virtual_temperature.hpp rename to src/thermo.hpp index e362246d7a..415ce469a1 100644 --- a/src/virtual_temperature.hpp +++ b/src/thermo.hpp @@ -1,5 +1,5 @@ -#ifndef VIRTUAL_TEMPERATURE_HPP // if not defined -#define VIRTUAL_TEMPERATURE_HPP // define the header file +#ifndef THERMO_HPP // if not defined +#define THERMO_HPP // define the header file #include #include @@ -44,4 +44,4 @@ double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::stri double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); -#endif // VIRTUAL_TEMPERATURE_HPP +#endif // THERMO_HPP From 1f127493c54c1fc358e891cd6bb89d887d19cb9e Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 17:13:09 -0600 Subject: [PATCH 29/48] parcel profile wip --- src/thermo.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/thermo.cpp b/src/thermo.cpp index dbcb70be13..8ad53e52f2 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -180,6 +180,9 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur press_upper.push_back(p); } } + std::vector temp_lower = MoistLapseProfile(press_upper, temp_lcl, press_lcl); + + return //return type } From 8eb21f0f2de1f19e7d77d46700f2d5d508113d63 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 18:08:46 -0600 Subject: [PATCH 30/48] add a return struct for _ParcelProfileHelper wip --- src/thermo.cpp | 14 +++++++++----- src/thermo.hpp | 10 ++++++++++ 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/src/thermo.cpp b/src/thermo.cpp index 8ad53e52f2..4451132fd9 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -144,7 +144,8 @@ bool _CheckPressure(const std::vector& pressure) { return true; // strictly non-increasing } -void _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { + +ParcelProfile _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { // Check that pressure does not increase. if (!_CheckPressure(pressure)) { throw std::runtime_error( @@ -169,7 +170,7 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur if (pressure.back() >= press_lcl) { press_lower.pop_back(); temp_lower.pop_back(); -// return {press_lower, {}, press_lcl, temp_lower, {}, temp_lcl}; + return {press_lower, {}, press_lcl, temp_lower, {}, temp_lcl}; } // Establish profile above LCL @@ -180,11 +181,14 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur press_upper.push_back(p); } } - std::vector temp_lower = MoistLapseProfile(press_upper, temp_lcl, press_lcl); - - return //return type + std::vector temp_upper = MoistLapseProfile(press_upper, temp_lower.back(), press_lcl); + press_lower.pop_back(); + temp_lower.pop_back(); + press_upper.erase(press_upper.begin()); + temp_upper.erase(temp_upper.begin()); + return {press_lower, press_upper, press_lcl, temp_lower, temp_upper, temp_lcl}; } double _SaturationVaporPressureLiquid(double temperature) { diff --git a/src/thermo.hpp b/src/thermo.hpp index 415ce469a1..f100468e52 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -32,6 +32,16 @@ std::pair LCL(double pressure, double temperature, double dewpoi bool _CheckPressure(const std::vector& pressure); +// Return struct for _ParcelProfileHelper +struct ParcelProfile { + std::vector press_lower, press_upper; + double press_lcl; + std::vector temp_lower, temp_upper; + double temp_lcl; +}; + +ParcelProfile _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint); + double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); double SaturationVaporPressure(double temperature, std::string phase); From e33ff90c11753587f4879597d3f1300352a3f5a3 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 24 Jun 2025 16:24:31 -0600 Subject: [PATCH 31/48] c++ ParcelProfile finished, assembled dry and moist part of the profile --- src/calcmod.cpp | 11 +++++++++-- src/metpy/calc/thermo.py | 17 +++++++++++++++++ src/thermo.cpp | 19 ++++++++++++++++++- src/thermo.hpp | 8 ++++++-- 4 files changed, 50 insertions(+), 5 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 5ffd42afaa..95398e765d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -99,8 +99,15 @@ PYBIND11_MODULE(_calc_mod, m) { py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); - - + m.def("parcel_profile", + [](py::array_t pressure, double temperature, double dewpoint) { + // pressure.data() gives the beginning pointer of the data buffer + std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); + std::vector temp_prof = ParcelProfile(pressure_vec, temperature, dewpoint); + return py::array_t(temp_prof.size(), temp_prof.data()); + }, + "Compute the parcel temperature profile as it rises from a given pressure and temperature.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index f9fafbc431..a19d40f4fe 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1312,6 +1312,23 @@ def parcel_profile(pressure, temperature, dewpoint): return concatenate((t_l, t_u)) +@exporter.export +@preprocess_and_wrap(wrap_like='pressure') +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'dewpoint': '[temperature]' + }, + '[temperature]' +) # process units because no unit should be passed to c++ function +def parcel_profile_linfel(pressure, temperature, dewpoint): + """ + Linfeng's version of 'parcel_profile'. Added on Jun 24 2025 + """ + return _calc_mod.parcel_profile(pressure, temperature, dewpoint) + + @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]') diff --git a/src/thermo.cpp b/src/thermo.cpp index 4451132fd9..f79825f27d 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -145,7 +145,24 @@ bool _CheckPressure(const std::vector& pressure) { } -ParcelProfile _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { +std::vector ParcelProfile(const std::vector& pressure, + double temperature, + double dewpoint) { + // Returns a vector of temperatures corresponding to the input pressure levels. + + ParProStruct profile = _ParcelProfileHelper(pressure, temperature, dewpoint); + + // Combine lower and upper temperature profiles + std::vector combined_temp; + combined_temp.reserve(pressure.size()); + + combined_temp.insert(combined_temp.end(), profile.temp_lower.begin(), profile.temp_lower.end()); + combined_temp.insert(combined_temp.end(), profile.temp_upper.begin(), profile.temp_upper.end()); + + return combined_temp; +} + +ParProStruct _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { // Check that pressure does not increase. if (!_CheckPressure(pressure)) { throw std::runtime_error( diff --git a/src/thermo.hpp b/src/thermo.hpp index f100468e52..f4f3c0379e 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -33,14 +33,18 @@ std::pair LCL(double pressure, double temperature, double dewpoi bool _CheckPressure(const std::vector& pressure); // Return struct for _ParcelProfileHelper -struct ParcelProfile { +struct ParProStruct { std::vector press_lower, press_upper; double press_lcl; std::vector temp_lower, temp_upper; double temp_lcl; }; -ParcelProfile _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint); +std::vector ParcelProfile(const std::vector& pressure, + double temperature, + double dewpoint); + +ParProStruct _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint); double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); From 44f1dddf364a398559229b1c1f1b1fecf585b452 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 25 Jun 2025 18:07:10 -0600 Subject: [PATCH 32/48] c++ version moist_lapse can handle 2D multiple starting temperature now --- src/calcmod.cpp | 49 +++++++++++++++++++++++++++++++++++++--- src/metpy/calc/thermo.py | 5 +++- src/thermo.cpp | 23 ++++++++----------- src/thermo.hpp | 6 ++--- 4 files changed, 63 insertions(+), 20 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 95398e765d..70aa270230 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -5,6 +5,7 @@ #include #include // For std::pair #include // For std::make_tuple +#include namespace py = pybind11; @@ -43,9 +44,51 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate the temperature at pressure levels assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); - m.def("moist_lapse", py::vectorize(MoistLapse), - "Calculate the temperature at pressure levels assuming saturated adiabatic process.", - py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); + m.def("moist_lapse", [](py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure, + int rk_nstep) { + // This function calculates the moist adiabatic profile for multiple starting + // temperatures (2D surface) and a single communal starting pressure, along a + // 1D pressure profile. + + // --- Step 1: Prepare the C++ vector for pressure levels --- + if (pressure.ndim() != 1) { + throw std::runtime_error("Input 'pressure' array must be 1D."); + } + std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); + + // --- Step 2: Ensure the reference temperature array is contiguous --- + auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); + + // --- Step 3: Define the shape of the output array: (N+1) dimension--- + // Create a vector to hold the shape dimensions. + std::vector out_shape; + for(int i = 0; i < ref_temp_contig.ndim(); ++i) { + out_shape.push_back(ref_temp_contig.shape(i)); + } + ssize_t profile_len = pressure_vec.size(); + out_shape.push_back(profile_len); + + // Create the final output array with the correct N+1 dimensional shape. + auto out_array = py::array_t(out_shape); + + // --- Step 4: Get direct pointers to data buffers for fast access --- + const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); + double* out_array_ptr = out_array.mutable_data(); + ssize_t num_profiles = ref_temp_contig.size(); // Total number of profiles to calculate + + // --- Step 5: Loop through each reference temperature --- + for (ssize_t i = 0; i < num_profiles; ++i) { + for (ssize_t j = 0; j < profile_len; ++j) { + out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); + } + } + + return out_array; + }, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", + py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); + m.def("lcl", [](py::array_t pressure, py::array_t temperature, diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a19d40f4fe..774c15240a 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -690,7 +690,10 @@ def dt(p, t): ) def moist_lapse_linfel(pressure, temperature, reference_pressure=None): """ - Linfeng's version of 'moist_lapse'. Added on Jun23 2025 + Linfeng's version of 'moist_lapse'. Added on Jun 25 2025 + This function calculates the moist adiabatic profile for multiple starting + temperatures (2D surface) and a single communal starting pressure, along a + 1D pressure profile. """ if reference_pressure is None: reference_pressure = pressure[0] diff --git a/src/thermo.cpp b/src/thermo.cpp index f79825f27d..ba1c7bd0f7 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -92,24 +92,21 @@ double MoistLapse(double pressure, double ref_temperature, double ref_pressure, return T1; // check final T1 P1 } -std::vector MoistLapseProfile(const std::vector& pressure_profile, +std::vector MoistLapseProfile(const std::vector& press_profile, double ref_temperature, - double ref_pressure) { - // Vectorized version of MoistLapse for a pressure profile. C++ internally use. + double ref_pressure, + int nstep) { + // MoistLapse for one full pressure profile given one ref_temperature. C++ internally use. // calculate temperature profile of an air parcel lifting saturated adiabatically // through the given pressure profile. - std::vector temperature_profile; - temperature_profile.reserve(pressure_profile.size()); + std::vector temp_profile; + temp_profile.reserve(press_profile.size()); - //double T1 = ref_temperature; - //double P1 = ref_pressure; - double T; - for (size_t i = 0; i < pressure_profile.size(); ++i) { - T = MoistLapse(pressure_profile[i], ref_temperature, ref_pressure, 50); - temperature_profile.push_back(T); + for (double p : press_profile) { + temp_profile.push_back(MoistLapse(p, ref_temperature, ref_pressure, nstep)); } - return temperature_profile; + return temp_profile; } std::pair LCL(double pressure, double temperature, double dewpoint) { @@ -198,7 +195,7 @@ ParProStruct _ParcelProfileHelper(const std::vector& pressure, double te press_upper.push_back(p); } } - std::vector temp_upper = MoistLapseProfile(press_upper, temp_lower.back(), press_lcl); + std::vector temp_upper = MoistLapseProfile(press_upper, temp_lower.back(), press_lcl, 30); press_lower.pop_back(); temp_lower.pop_back(); diff --git a/src/thermo.hpp b/src/thermo.hpp index f4f3c0379e..2623bc366b 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -23,10 +23,10 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, double CaldlnTdlnP(double temperature, double pressure); double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep); -std::vector MoistLapseProfile(const std::vector& pressure_profile, +std::vector MoistLapseProfile(const std::vector& press_profile, double ref_temperature, - double ref_pressure); - + double ref_pressure, + int nstep); std::pair LCL(double pressure, double temperature, double dewpoint); From e0d1f1e8904bd5c294c2b9f975b3f2a25925d62f Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 26 Jun 2025 08:24:00 -0600 Subject: [PATCH 33/48] put MoistLapseVectorized in thermo.cpp, no need to vectorize in pybind11 module --- src/calcmod.cpp | 90 ++++++++++++++++++++++++++----------------------- src/thermo.cpp | 53 ++++++++++++++++++++++++++--- src/thermo.hpp | 11 ++++-- 3 files changed, 104 insertions(+), 50 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 70aa270230..5e8013be14 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -44,51 +44,55 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate the temperature at pressure levels assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); - m.def("moist_lapse", [](py::array_t pressure, - py::array_t ref_temperature, - double ref_pressure, - int rk_nstep) { - // This function calculates the moist adiabatic profile for multiple starting - // temperatures (2D surface) and a single communal starting pressure, along a - // 1D pressure profile. - - // --- Step 1: Prepare the C++ vector for pressure levels --- - if (pressure.ndim() != 1) { - throw std::runtime_error("Input 'pressure' array must be 1D."); - } - std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); - - // --- Step 2: Ensure the reference temperature array is contiguous --- - auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); - - // --- Step 3: Define the shape of the output array: (N+1) dimension--- - // Create a vector to hold the shape dimensions. - std::vector out_shape; - for(int i = 0; i < ref_temp_contig.ndim(); ++i) { - out_shape.push_back(ref_temp_contig.shape(i)); - } - ssize_t profile_len = pressure_vec.size(); - out_shape.push_back(profile_len); - - // Create the final output array with the correct N+1 dimensional shape. - auto out_array = py::array_t(out_shape); - - // --- Step 4: Get direct pointers to data buffers for fast access --- - const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); - double* out_array_ptr = out_array.mutable_data(); - ssize_t num_profiles = ref_temp_contig.size(); // Total number of profiles to calculate - - // --- Step 5: Loop through each reference temperature --- - for (ssize_t i = 0; i < num_profiles; ++i) { - for (ssize_t j = 0; j < profile_len; ++j) { - out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); - } - } - - return out_array; - }, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", + m.def("moist_lapse", &MoistLapseVectorized, + "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); +// m.def("moist_lapse", [](py::array_t pressure, +// py::array_t ref_temperature, +// double ref_pressure, +// int rk_nstep) { +// // This function calculates the moist adiabatic profile for multiple starting +// // temperatures (2D surface) and a single communal starting pressure, along a +// // 1D pressure profile. +// +// // --- Step 1: Prepare the C++ vector for pressure levels --- +// if (pressure.ndim() > 1) { +// throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); +// } +// std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); +// +// // --- Step 2: Ensure the reference temperature array is contiguous --- +// auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); +// +// // --- Step 3: Define the shape of the output array: (N+1) dimension--- +// // Create a vector to hold the shape dimensions. +// std::vector out_shape; +// for(int i = 0; i < ref_temp_contig.ndim(); ++i) { +// out_shape.push_back(ref_temp_contig.shape(i)); +// } +// ssize_t profile_len = pressure_vec.size(); +// out_shape.push_back(profile_len); +// +// // Create the final output array with the correct N+1 dimensional shape. +// auto out_array = py::array_t(out_shape); +// +// // --- Step 4: Get direct pointers to data buffers for fast access --- +// const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); +// double* out_array_ptr = out_array.mutable_data(); +// ssize_t num_profiles = ref_temp_contig.size(); // Total number of profiles to calculate +// +// // --- Step 5: Loop through each reference temperature --- +// for (ssize_t i = 0; i < num_profiles; ++i) { +// for (ssize_t j = 0; j < profile_len; ++j) { +// out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); +// } +// } +// +// return out_array; +// }, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", +// py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); + m.def("lcl", [](py::array_t pressure, py::array_t temperature, diff --git a/src/thermo.cpp b/src/thermo.cpp index ba1c7bd0f7..0ef694251f 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -3,6 +3,7 @@ #include #include // For std::pair #include +#include #include #include // for std::cerr #include // for std::numeric_limits @@ -69,17 +70,17 @@ double CaldlnTdlnP(double temperature, double pressure) { return dlnT_dlnP_Bakhshaii2013; } -double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep) { +double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep) { // calculate temperature at pressure given reference temperature and pressure // assuming moist adiabatic expansion (vapor condenses and removed from the air // parcel) - double dlnP = log(pressure / ref_pressure) / (double)nstep; + double dlnP = log(pressure / ref_pressure) / (double)rk_nstep; double T1 = ref_temperature; double P1 = ref_pressure; double k[4]; - for (int i = 0; i < nstep; ++i) { + for (int i = 0; i < rk_nstep; ++i) { k[0] = CaldlnTdlnP(T1, P1); k[1] = CaldlnTdlnP(T1 * exp(k[0] * dlnP/2.), P1 * exp(dlnP/2.)); k[2] = CaldlnTdlnP(T1 * exp(k[1] * dlnP/2.), P1 * exp(dlnP/2.)); @@ -95,7 +96,7 @@ double MoistLapse(double pressure, double ref_temperature, double ref_pressure, std::vector MoistLapseProfile(const std::vector& press_profile, double ref_temperature, double ref_pressure, - int nstep) { + int rk_nstep) { // MoistLapse for one full pressure profile given one ref_temperature. C++ internally use. // calculate temperature profile of an air parcel lifting saturated adiabatically @@ -104,11 +105,53 @@ std::vector MoistLapseProfile(const std::vector& press_profile, temp_profile.reserve(press_profile.size()); for (double p : press_profile) { - temp_profile.push_back(MoistLapse(p, ref_temperature, ref_pressure, nstep)); + temp_profile.push_back(MoistLapse(p, ref_temperature, ref_pressure, rk_nstep)); } return temp_profile; } +py::array_t MoistLapseVectorized(py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure, + int rk_nstep) { + // This function calculates the moist adiabatic profile for multiple starting + // temperatures (2D surface) and a single communal starting pressure, along a + // 1D pressure profile. + // --- Step 1: Prepare the C++ vector for pressure levels --- + if (pressure.ndim() > 1) { + throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); + } + std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); + + // --- Step 2: Ensure the reference temperature array is contiguous --- + auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); + + // --- Step 3: Define the shape of the output array: (N+1) dimension--- + std::vector out_shape; + for(int i = 0; i < ref_temp_contig.ndim(); ++i) { + out_shape.push_back(ref_temp_contig.shape(i)); + } + ssize_t profile_len = pressure_vec.size(); + out_shape.push_back(profile_len); + + auto out_array = py::array_t(out_shape); + + // --- Step 4: Get direct pointers to data buffers for fast access --- + const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); + double* out_array_ptr = out_array.mutable_data(); + ssize_t num_profiles = ref_temp_contig.size(); + + // --- Step 5: Loop through each reference temperature --- + for (ssize_t i = 0; i < num_profiles; ++i) { + for (ssize_t j = 0; j < profile_len; ++j) { + out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); + } + } + + return out_array; +} + + std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Warning in function '" << __func__ diff --git a/src/thermo.hpp b/src/thermo.hpp index 2623bc366b..b1632333de 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -4,8 +4,11 @@ #include #include #include // For std::pair +#include +#include #include "constants.hpp" +namespace py = pybind11; namespace mc = metpy_constants; double MoistAirGasConstant(double specific_humidity); @@ -22,11 +25,15 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, double ref_pressure); double CaldlnTdlnP(double temperature, double pressure); -double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep); +double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep); std::vector MoistLapseProfile(const std::vector& press_profile, double ref_temperature, double ref_pressure, - int nstep); + int rk_nstep); +py::array_t MoistLapseVectorized(py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure, + int rk_nstep); std::pair LCL(double pressure, double temperature, double dewpoint); From 853a213634afdd9b811b6b2984254890a0d150b8 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 26 Jun 2025 08:58:11 -0600 Subject: [PATCH 34/48] add DryLapseVectorized in thermo.cpp, no need to vectorize in pybind11 module --- src/calcmod.cpp | 49 ++----------------------------------------------- src/thermo.cpp | 41 +++++++++++++++++++++++++++++++++++++++++ src/thermo.hpp | 3 +++ 3 files changed, 46 insertions(+), 47 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 5e8013be14..f6f1c1aca0 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -40,59 +40,14 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate relative humidity from temperature and dewpoint.", py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); - m.def("dry_lapse", py::vectorize(DryLapse), - "Calculate the temperature at pressure levels assuming dry adiabatic process.", + m.def("dry_lapse", &DryLapseVectorized, + "Calculate the temperature along a pressure profile assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); m.def("moist_lapse", &MoistLapseVectorized, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); -// m.def("moist_lapse", [](py::array_t pressure, -// py::array_t ref_temperature, -// double ref_pressure, -// int rk_nstep) { -// // This function calculates the moist adiabatic profile for multiple starting -// // temperatures (2D surface) and a single communal starting pressure, along a -// // 1D pressure profile. -// -// // --- Step 1: Prepare the C++ vector for pressure levels --- -// if (pressure.ndim() > 1) { -// throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); -// } -// std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); -// -// // --- Step 2: Ensure the reference temperature array is contiguous --- -// auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); -// -// // --- Step 3: Define the shape of the output array: (N+1) dimension--- -// // Create a vector to hold the shape dimensions. -// std::vector out_shape; -// for(int i = 0; i < ref_temp_contig.ndim(); ++i) { -// out_shape.push_back(ref_temp_contig.shape(i)); -// } -// ssize_t profile_len = pressure_vec.size(); -// out_shape.push_back(profile_len); -// -// // Create the final output array with the correct N+1 dimensional shape. -// auto out_array = py::array_t(out_shape); -// -// // --- Step 4: Get direct pointers to data buffers for fast access --- -// const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); -// double* out_array_ptr = out_array.mutable_data(); -// ssize_t num_profiles = ref_temp_contig.size(); // Total number of profiles to calculate -// -// // --- Step 5: Loop through each reference temperature --- -// for (ssize_t i = 0; i < num_profiles; ++i) { -// for (ssize_t j = 0; j < profile_len; ++j) { -// out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); -// } -// } -// -// return out_array; -// }, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", -// py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); - m.def("lcl", [](py::array_t pressure, py::array_t temperature, diff --git a/src/thermo.cpp b/src/thermo.cpp index 0ef694251f..122c48a935 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -58,6 +58,47 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, return temperature_profile; } + +py::array_t DryLapseVectorized(py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure) { + // This function calculates the dry adiabatic profile for multiple starting + // temperatures (2D surface) and a single communal starting pressure, along a + // 1D pressure profile. + // --- Step 1: Prepare the C++ vector for pressure levels --- + if (pressure.ndim() > 1) { + throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); + } + std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); + + // --- Step 2: Ensure the reference temperature array is contiguous --- + auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); + + // --- Step 3: Define the shape of the output array: (N+1) dimension--- + std::vector out_shape; + for(int i = 0; i < ref_temp_contig.ndim(); ++i) { + out_shape.push_back(ref_temp_contig.shape(i)); + } + ssize_t profile_len = pressure_vec.size(); + out_shape.push_back(profile_len); + + auto out_array = py::array_t(out_shape); + + // --- Step 4: Get direct pointers to data buffers for fast access --- + const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); + double* out_array_ptr = out_array.mutable_data(); + ssize_t num_profiles = ref_temp_contig.size(); + + // --- Step 5: Loop through each reference temperature --- + for (ssize_t i = 0; i < num_profiles; ++i) { + for (ssize_t j = 0; j < profile_len; ++j) { + out_array_ptr[i * profile_len + j] = DryLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure); + } + } + + return out_array; +} + double CaldlnTdlnP(double temperature, double pressure) { // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. double rs = SaturationMixingRatio(pressure, temperature, "liquid"); diff --git a/src/thermo.hpp b/src/thermo.hpp index b1632333de..ea48f2ae0f 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -23,6 +23,9 @@ double DryLapse(double pressure, double ref_temperature, double ref_pressure); std::vector DryLapseProfile(const std::vector& pressure_profile, double ref_temperature, double ref_pressure); +py::array_t DryLapseVectorized(py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure); double CaldlnTdlnP(double temperature, double pressure); double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep); From e7679607451694f546f1cac38011212bccbaed49 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 27 Jun 2025 06:53:42 -0600 Subject: [PATCH 35/48] modify conditions for LCL to execute --- src/thermo.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/thermo.cpp b/src/thermo.cpp index 122c48a935..6aa15c91ac 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -194,7 +194,7 @@ py::array_t MoistLapseVectorized(py::array_t pressure, std::pair LCL(double pressure, double temperature, double dewpoint) { - if (temperature <= dewpoint) { + if (temperature < dewpoint) { std::cerr << "Warning in function '" << __func__ << "': Temperature must be greater than dew point for LCL calculation.\n"; return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; From 0afe0cc22533280243a9f0e98b8ed269aa66074e Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 30 Jun 2025 09:35:17 -0600 Subject: [PATCH 36/48] add c++ virtual_temperature_from_dewpoint --- src/calcmod.cpp | 6 ++++- src/metpy/calc/thermo.py | 53 +++++++++++++++++++++++++++++++++++++++- src/thermo.cpp | 7 ++++++ src/thermo.hpp | 4 ++- 4 files changed, 67 insertions(+), 3 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index f6f1c1aca0..fd99061c15 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -143,9 +143,13 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("specific_humidity_from_dewpoint", py::vectorize(SpecificHumidityFromDewPoint), "Calculate the specific humidity from the dewpoint temperature and pressure.", py::arg("pressure"), py::arg("dewpoint"), py::arg("phase")); - m.def("virtual_temperature", py::vectorize(VirtualTemperature), "Calculate virtual temperature from temperature and mixing ratio.", py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon")); + + m.def("virtual_temperature_from_dewpoint", py::vectorize(VirtualTemperatureFromDewpoint), + "Calculate virtual temperature from dewpoint.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint"), py::arg("epsilon"), py::arg("phase")); + } diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 774c15240a..ed832f9dbb 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2185,10 +2185,29 @@ def virtual_temperature( Renamed ``mixing`` parameter to ``mixing_ratio`` """ - # Calling c++ calculation module + return temperature * ((mixing_ratio + molecular_weight_ratio) + / (molecular_weight_ratio * (1 + mixing_ratio))) + +@exporter.export +@preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'mixing_ratio')) +@process_units( + { + 'temperature': '[temperature]', + 'mixing_ratio': '[dimensionless]', + 'molecular_weight_ratio': '[dimensionless]' + }, + '[temperature]', + ignore_inputs_for_output=('molecular_weight_ratio',) +) +def virtual_temperature_linfel( + temperature, mixing_ratio, molecular_weight_ratio=mpconsts.nounit.epsilon): + """ + Linfeng's version of 'virtual_temperature'. Added on Jun 30 2025 + """ return _calc_mod.virtual_temperature( temperature, mixing_ratio, molecular_weight_ratio) + @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', 'temperature', @@ -2257,6 +2276,38 @@ def virtual_temperature_from_dewpoint( return virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio) +@exporter.export +@preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', + 'temperature', + 'dewpoint')) +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'dewpoint': '[temperature]', + 'molecular_weight_ratio': '[dimensionless]' + }, + '[temperature]', + ignore_inputs_for_output=('molecular_weight_ratio',) +) +def virtual_temperature_from_dewpoint_linfel( + pressure, + temperature, + dewpoint, + molecular_weight_ratio=mpconsts.nounit.epsilon, + *, + phase='liquid' +): + """ + Linfeng's version of 'virtual_temperature_from_dewpoint'. Added on Jun 30 2025 + """ + return _calc_mod.virtual_temperature_from_dewpoint(pressure, + temperature, + dewpoint, + molecular_weight_ratio, + phase) + + @exporter.export @preprocess_and_wrap( wrap_like='temperature', diff --git a/src/thermo.cpp b/src/thermo.cpp index 6aa15c91ac..c96528b526 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -354,3 +354,10 @@ double SpecificHumidityFromDewPoint(double pressure, double dewpoint, std::strin double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { return temperature * (mixing_ratio + epsilon) / (epsilon * (1. + mixing_ratio)); } + +double VirtualTemperatureFromDewpoint(double pressure, double temperature, + double dewpoint, double epsilon, + std::string phase) { + double mixing_ratio = SaturationMixingRatio(pressure, dewpoint, phase); + return VirtualTemperature(temperature, mixing_ratio, epsilon); +} diff --git a/src/thermo.hpp b/src/thermo.hpp index ea48f2ae0f..1333c3a028 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -67,5 +67,7 @@ double SpecificHumidityFromMixingRatio(double mixing_ratio); double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::string phase); double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); - +double VirtualTemperatureFromDewpoint(double pressure, double temperature, + double dewpoint, double epsilon, + std::string phase); #endif // THERMO_HPP From 47845b53f317ae5b727a50e4cca2022aa195e661 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 7 Jul 2025 09:50:54 -0400 Subject: [PATCH 37/48] change ssize_t to size_t, ssize_t is not supported on windows --- src/calcmod.cpp | 4 ++-- src/thermo.cpp | 20 ++++++++++---------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index fd99061c15..f742dfa974 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -76,7 +76,7 @@ PYBIND11_MODULE(_calc_mod, m) { auto t_lcl = py::array_t(p_contig.request().shape); // --- Step 3: Get the total number of elements to loop over --- - ssize_t size = p_contig.size(); + size_t size = p_contig.size(); // --- Step 4: Get direct pointers to the (now contiguous) data buffers --- const double* p_ptr = static_cast(p_contig.request().ptr); @@ -86,7 +86,7 @@ PYBIND11_MODULE(_calc_mod, m) { double* t_lcl_ptr = t_lcl.mutable_data(); // --- Step 5: Loop through the data as if it were a single flat 1D array --- - for (ssize_t i = 0; i < size; i++) { + for (size_t i = 0; i < size; i++) { // Call the scalar c++ function for each element std::pair result = LCL(p_ptr[i], t_ptr[i], d_ptr[i]); diff --git a/src/thermo.cpp b/src/thermo.cpp index c96528b526..578f1b6e5b 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -75,11 +75,11 @@ py::array_t DryLapseVectorized(py::array_t pressure, auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); // --- Step 3: Define the shape of the output array: (N+1) dimension--- - std::vector out_shape; + std::vector out_shape; for(int i = 0; i < ref_temp_contig.ndim(); ++i) { out_shape.push_back(ref_temp_contig.shape(i)); } - ssize_t profile_len = pressure_vec.size(); + size_t profile_len = pressure_vec.size(); out_shape.push_back(profile_len); auto out_array = py::array_t(out_shape); @@ -87,11 +87,11 @@ py::array_t DryLapseVectorized(py::array_t pressure, // --- Step 4: Get direct pointers to data buffers for fast access --- const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); double* out_array_ptr = out_array.mutable_data(); - ssize_t num_profiles = ref_temp_contig.size(); + size_t num_profiles = ref_temp_contig.size(); // --- Step 5: Loop through each reference temperature --- - for (ssize_t i = 0; i < num_profiles; ++i) { - for (ssize_t j = 0; j < profile_len; ++j) { + for (size_t i = 0; i < num_profiles; ++i) { + for (size_t j = 0; j < profile_len; ++j) { out_array_ptr[i * profile_len + j] = DryLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure); } } @@ -168,11 +168,11 @@ py::array_t MoistLapseVectorized(py::array_t pressure, auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); // --- Step 3: Define the shape of the output array: (N+1) dimension--- - std::vector out_shape; + std::vector out_shape; for(int i = 0; i < ref_temp_contig.ndim(); ++i) { out_shape.push_back(ref_temp_contig.shape(i)); } - ssize_t profile_len = pressure_vec.size(); + size_t profile_len = pressure_vec.size(); out_shape.push_back(profile_len); auto out_array = py::array_t(out_shape); @@ -180,11 +180,11 @@ py::array_t MoistLapseVectorized(py::array_t pressure, // --- Step 4: Get direct pointers to data buffers for fast access --- const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); double* out_array_ptr = out_array.mutable_data(); - ssize_t num_profiles = ref_temp_contig.size(); + size_t num_profiles = ref_temp_contig.size(); // --- Step 5: Loop through each reference temperature --- - for (ssize_t i = 0; i < num_profiles; ++i) { - for (ssize_t j = 0; j < profile_len; ++j) { + for (size_t i = 0; i < num_profiles; ++i) { + for (size_t j = 0; j < profile_len; ++j) { out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); } } From 87b31eada8d4762b1b065dc840639d385564f1f6 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 7 Jul 2025 18:41:21 -0400 Subject: [PATCH 38/48] move lcl python binding from calcmod.cpp to thermo.cpp; debug syntax 'structured binding' --- src/calcmod.cpp | 54 +++-------------------------------------------- src/thermo.cpp | 56 +++++++++++++++++++++++++++++++++++++++++++++++-- src/thermo.hpp | 4 ++++ 3 files changed, 61 insertions(+), 53 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index f742dfa974..58ec3c7243 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -49,57 +49,9 @@ PYBIND11_MODULE(_calc_mod, m) { py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); - m.def("lcl", [](py::array_t pressure, - py::array_t temperature, - py::array_t dewpoint) { - - // This helper ensures the arrays are in C-style contiguous memory. - // If an input array is already contiguous, it's a zero-cost operation. - // If it's a slice or has a different memory layout, it creates a copy. - // This makes the subsequent looping simple and safe. - auto p_contig = py::array::ensure(pressure, py::array::c_style); - auto t_contig = py::array::ensure(temperature, py::array::c_style); - auto d_contig = py::array::ensure(dewpoint, py::array::c_style); - - // --- Step 1: Check that all input arrays have the same shape --- - if (p_contig.ndim() != t_contig.ndim() || p_contig.ndim() != d_contig.ndim()) { - throw std::runtime_error("Input arrays must have the same number of dimensions."); - } - for (int i = 0; i < p_contig.ndim(); ++i) { - if (p_contig.shape(i) != t_contig.shape(i) || p_contig.shape(i) != d_contig.shape(i)) { - throw std::runtime_error("Input arrays must have the same shape."); - } - } - - // --- Step 2: Create output arrays with the exact same N-D shape as the inputs --- - auto p_lcl = py::array_t(p_contig.request().shape); - auto t_lcl = py::array_t(p_contig.request().shape); - - // --- Step 3: Get the total number of elements to loop over --- - size_t size = p_contig.size(); - - // --- Step 4: Get direct pointers to the (now contiguous) data buffers --- - const double* p_ptr = static_cast(p_contig.request().ptr); - const double* t_ptr = static_cast(t_contig.request().ptr); - const double* d_ptr = static_cast(d_contig.request().ptr); - double* p_lcl_ptr = p_lcl.mutable_data(); - double* t_lcl_ptr = t_lcl.mutable_data(); - - // --- Step 5: Loop through the data as if it were a single flat 1D array --- - for (size_t i = 0; i < size; i++) { - // Call the scalar c++ function for each element - std::pair result = LCL(p_ptr[i], t_ptr[i], d_ptr[i]); - - p_lcl_ptr[i] = result.first; - t_lcl_ptr[i] = result.second; - } - - // --- Step 6: Return a tuple of the two new, N-dimensional arrays --- - return std::make_tuple(p_lcl, t_lcl); - - }, "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", - py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); - + m.def("lcl", &LCLVectorized, + "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); m.def("parcel_profile", [](py::array_t pressure, double temperature, double dewpoint) { diff --git a/src/thermo.cpp b/src/thermo.cpp index 578f1b6e5b..4bf70e6932 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -1,6 +1,7 @@ #include #include #include +#include // For std::make_tuple #include // For std::pair #include #include @@ -212,10 +213,59 @@ std::pair LCL(double pressure, double temperature, double dewpoi double w_minus1 = lambert_wm1(pow(rh, 1.0 / a) * c * exp(c)); double t_lcl = c / w_minus1 * temperature; double p_lcl = pressure * pow(t_lcl / temperature, moist_heat_ratio); - + return {p_lcl, t_lcl}; } + +std::tuple, py::array_t> LCLVectorized(py::array_t pressure, + py::array_t temperature, + py::array_t dewpoint) { + + // This helper ensures the arrays are in C-style contiguous memory. + // If an input array is already contiguous, it's a zero-cost operation. + // If it's a slice or has a different memory layout, it creates a copy. + // This makes the subsequent looping simple and safe. + auto p_contig = py::array::ensure(pressure, py::array::c_style); + auto t_contig = py::array::ensure(temperature, py::array::c_style); + auto d_contig = py::array::ensure(dewpoint, py::array::c_style); + + // --- Step 1: Check that all input arrays have the same shape --- + if (p_contig.ndim() != t_contig.ndim() || p_contig.ndim() != d_contig.ndim()) { + throw std::runtime_error("Input arrays must have the same number of dimensions."); + } + for (int i = 0; i < p_contig.ndim(); ++i) { + if (p_contig.shape(i) != t_contig.shape(i) || p_contig.shape(i) != d_contig.shape(i)) { + throw std::runtime_error("Input arrays must have the same shape."); + } + } + + // --- Step 2: Create output arrays with the exact same N-D shape as the inputs --- + auto p_lcl = py::array_t(p_contig.request().shape); + auto t_lcl = py::array_t(p_contig.request().shape); + + // --- Step 3: Get the total number of elements to loop over --- + size_t size = p_contig.size(); + + // --- Step 4: Get direct pointers to the (now contiguous) data buffers --- + const double* p_ptr = static_cast(p_contig.request().ptr); + const double* t_ptr = static_cast(t_contig.request().ptr); + const double* d_ptr = static_cast(d_contig.request().ptr); + double* p_lcl_ptr = p_lcl.mutable_data(); + double* t_lcl_ptr = t_lcl.mutable_data(); + + // --- Step 5: Loop through the data as if it were a single flat 1D array --- + for (size_t i = 0; i < size; i++) { + // Call the scalar c++ function for each element + std::pair result = LCL(p_ptr[i], t_ptr[i], d_ptr[i]); + p_lcl_ptr[i] = result.first; + t_lcl_ptr[i] = result.second; + } + + // --- Step 6: Return a tuple of the two new, N-dimensional arrays --- + return std::make_tuple(p_lcl, t_lcl); +} + bool _CheckPressure(const std::vector& pressure) { for (size_t i = 0; i + 1 < pressure.size(); ++i) { if (pressure[i] < pressure[i + 1]) { @@ -252,7 +302,9 @@ ParProStruct _ParcelProfileHelper(const std::vector& pressure, double te } // Find the LCL - auto [press_lcl, temp_lcl] = LCL(pressure[0], temperature, dewpoint); + std::pair lcl_result = LCL(pressure[0], temperature, dewpoint); + double press_lcl = lcl_result.first; + double temp_lcl = lcl_result.second; // Establish profile below LCL std::vector press_lower; diff --git a/src/thermo.hpp b/src/thermo.hpp index 1333c3a028..7af7f2c45d 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -4,6 +4,7 @@ #include #include #include // For std::pair +#include // For std::tuple #include #include #include "constants.hpp" @@ -39,6 +40,9 @@ py::array_t MoistLapseVectorized(py::array_t pressure, int rk_nstep); std::pair LCL(double pressure, double temperature, double dewpoint); +std::tuple, py::array_t> LCLVectorized(py::array_t pressure, + py::array_t temperature, + py::array_t dewpoint); bool _CheckPressure(const std::vector& pressure); From c93e6aa47014ac31da42b8517a7ef7e2f9b72748 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 11 Jul 2025 10:40:18 -0600 Subject: [PATCH 39/48] add DryLapseVectorized_3D --- src/calcmod.cpp | 4 ++++ src/thermo.cpp | 49 +++++++++++++++++++++++++++++++++++++++++++++++++ src/thermo.hpp | 3 +++ 3 files changed, 56 insertions(+) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 58ec3c7243..50929f8b59 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -44,6 +44,10 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate the temperature along a pressure profile assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); + m.def("dry_lapse_3d", &DryLapseVectorized_3D, + "Calculate the temperature along multiple pressure profiles assuming dry adiabatic process.", + py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); + m.def("moist_lapse", &MoistLapseVectorized, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); diff --git a/src/thermo.cpp b/src/thermo.cpp index 4bf70e6932..6b6e1ed43b 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -100,6 +100,55 @@ py::array_t DryLapseVectorized(py::array_t pressure, return out_array; } +py::array_t DryLapseVectorized_3D(py::array_t pressure, + py::array_t ref_temperature, + py::array_t ref_pressure) { + // --- Step 1: Ensure all input arrays are using a contiguous memory layout --- + auto p_contig = py::array::ensure(pressure, py::array::c_style); + auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); + auto ref_press_contig = py::array::ensure(ref_pressure, py::array::c_style); + + // --- Step 2: Perform comprehensive shape validation --- + if (ref_temp_contig.ndim() != ref_press_contig.ndim()) { + throw std::runtime_error("Input 'ref_temperature' and 'ref_pressure' must have the same number of dimensions."); + } + if (p_contig.ndim() != ref_temp_contig.ndim() + 1) { + throw std::runtime_error("Input 'pressure' must have one more dimension than 'ref_temperature'."); + } + for (int i = 0; i < ref_temp_contig.ndim(); ++i) { + if (ref_temp_contig.shape(i) != ref_press_contig.shape(i) || + p_contig.shape(i+1) != ref_temp_contig.shape(i)) { + throw std::runtime_error("The horizontal dimensions of all input arrays must match."); + } + } + + // --- Step 3: Define the shape of the output array --- + // The output shape will be identical to the input pressure array's shape. + auto out_array = py::array_t(p_contig.request().shape); + + // --- Step 4: Get direct pointers to data buffers for fast access --- + const double* pressure_ptr = static_cast(p_contig.request().ptr); + const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); + const double* ref_press_ptr = static_cast(ref_press_contig.request().ptr); + double* out_array_ptr = out_array.mutable_data(); + + // --- Step 5: Define loop boundaries --- + size_t num_profiles = ref_temp_contig.size(); // Total number of horizontal points + size_t profile_len = p_contig.shape(0); // Length of the vertical dimension + + // --- Step 6: Loop through each horizontal point and its vertical profile --- + for (int j = 0; j < profile_len; ++j) { + for (int i = 0; i < num_profiles; ++i) { + // Calculate the index for the current point in the flattened arrays. + out_array_ptr[i+j*num_profiles] = DryLapse(pressure_ptr[i+j*num_profiles], + ref_temp_ptr[i], + ref_press_ptr[i]); + } + } + + return out_array; +} + double CaldlnTdlnP(double temperature, double pressure) { // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. double rs = SaturationMixingRatio(pressure, temperature, "liquid"); diff --git a/src/thermo.hpp b/src/thermo.hpp index 7af7f2c45d..656a6b37c6 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -27,6 +27,9 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, py::array_t DryLapseVectorized(py::array_t pressure, py::array_t ref_temperature, double ref_pressure); +py::array_t DryLapseVectorized_3D(py::array_t pressure, + py::array_t ref_temperature, + py::array_t ref_pressure); double CaldlnTdlnP(double temperature, double pressure); double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep); From 57995ebabc25f7c2a11403d3947af2afd06f4323 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 16 Jul 2025 09:36:31 -0600 Subject: [PATCH 40/48] debug: number of dimension inconsistency when lcl has inputs of a scaler and an array with one element --- src/metpy/calc/thermo.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index ed832f9dbb..e387bf0a43 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -797,6 +797,7 @@ def lcl_linfel(pressure, temperature, dewpoint, max_iters=None, eps=None): """ Linfeng's version of 'lcl'. Added on Jun23 2025 """ + pressure, temperature, dewpoint = np.atleast_1d(pressure, temperature, dewpoint) p_lcl, t_lcl = _calc_mod.lcl(pressure, temperature, dewpoint) return p_lcl, t_lcl From 9dbc34871cd423e82073858cc048f835b260f95c Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 21 Jul 2025 09:53:33 -0600 Subject: [PATCH 41/48] code convention check --- src/calcmod.cpp | 4 ++-- src/math.cpp | 4 ++-- src/thermo.cpp | 4 ++-- src/thermo.hpp | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 50929f8b59..8c2cdbe64d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,8 +1,8 @@ -#include "constants.hpp" -#include "thermo.hpp" #include #include #include +#include "constants.hpp" +#include "thermo.hpp" #include // For std::pair #include // For std::make_tuple #include diff --git a/src/math.cpp b/src/math.cpp index 2c0dc4f39d..12c89469a0 100644 --- a/src/math.cpp +++ b/src/math.cpp @@ -18,8 +18,8 @@ double lambert_wm1(double x, double tol, int max_iter) { for (int i = 0; i < max_iter; ++i) { double ew = std::exp(w); - double wew = w * ew; - double diff = (wew - x) / (ew * (w + 1) - ((w + 2) * (wew - x)) / (2 * w + 2)); + double w_ew = w * ew; + double diff = (w_ew - x) / (ew * (w + 1) - ((w + 2) * (w_ew - x)) / (2 * w + 2)); w -= diff; if (std::abs(diff) < tol) { return w; diff --git a/src/thermo.cpp b/src/thermo.cpp index 6b6e1ed43b..9f04ab82c8 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -1,10 +1,10 @@ +#include +#include #include #include #include #include // For std::make_tuple #include // For std::pair -#include -#include #include #include // for std::cerr #include // for std::numeric_limits diff --git a/src/thermo.hpp b/src/thermo.hpp index 656a6b37c6..555bbc437a 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -1,12 +1,12 @@ #ifndef THERMO_HPP // if not defined #define THERMO_HPP // define the header file +#include +#include #include #include #include // For std::pair #include // For std::tuple -#include -#include #include "constants.hpp" namespace py = pybind11; From 2076d8edc63b5bec51a13a051dd1fa4b068ccc37 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 1 Jul 2025 14:06:32 -0600 Subject: [PATCH 42/48] update: let python version _parcel_profile_help call underlying c++ computing functions --- src/metpy/calc/thermo.py | 154 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 144 insertions(+), 10 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index e387bf0a43..d62d2eef53 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1318,19 +1318,14 @@ def parcel_profile(pressure, temperature, dewpoint): @exporter.export @preprocess_and_wrap(wrap_like='pressure') -@process_units( - { - 'pressure': '[pressure]', - 'temperature': '[temperature]', - 'dewpoint': '[temperature]' - }, - '[temperature]' -) # process units because no unit should be passed to c++ function +@check_units('[pressure]', '[temperature]', '[temperature]') def parcel_profile_linfel(pressure, temperature, dewpoint): """ - Linfeng's version of 'parcel_profile'. Added on Jun 24 2025 + Linfeng's version of 'parcel_profile'. Added on Jul 1 2025 """ - return _calc_mod.parcel_profile(pressure, temperature, dewpoint) + + _, _, _, t_l, _, t_u = _parcel_profile_helper_linfel(pressure, temperature, dewpoint) + return concatenate((t_l, t_u)) @exporter.export @@ -1424,6 +1419,100 @@ def parcel_profile_with_lcl(pressure, temperature, dewpoint): return new_press, new_temp, new_dewp, prof_temp +@exporter.export +@preprocess_and_wrap() +@check_units('[pressure]', '[temperature]', '[temperature]') +def parcel_profile_with_lcl_linfel(pressure, temperature, dewpoint): + """ + Linfeng's version of 'parcel_profile_with_lcl'. Added on Jul 1 2025 + """ + r"""Calculate the profile a parcel takes through the atmosphere. + + The parcel starts at `temperature`, and `dewpoint`, lifted up + dry adiabatically to the LCL, and then moist adiabatically from there. + `pressure` specifies the pressure levels for the profile. This function returns + a profile that includes the LCL. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure level(s) of interest. This array must be from + high to low pressure. + + temperature : `pint.Quantity` + Atmospheric temperature at the levels in `pressure`. The first entry should be at + the same level as the first `pressure` data point. + + dewpoint : `pint.Quantity` + Atmospheric dewpoint at the levels in `pressure`. The first entry should be at + the same level as the first `pressure` data point. + + Returns + ------- + pressure : `pint.Quantity` + The parcel profile pressures, which includes the specified levels and the LCL + + ambient_temperature : `pint.Quantity` + Atmospheric temperature values, including the value interpolated to the LCL level + + ambient_dew_point : `pint.Quantity` + Atmospheric dewpoint values, including the value interpolated to the LCL level + + profile_temperature : `pint.Quantity` + The parcel profile temperatures at all of the levels in the returned pressures array, + including the LCL + + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, parcel_profile_with_lcl + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compute parcel temperature + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> p_wLCL, T_wLCL, Td_wLCL, prof_wLCL = parcel_profile_with_lcl(p, T, Td) + + See Also + -------- + lcl, moist_lapse, dry_lapse, parcel_profile, parcel_profile_with_lcl_as_dataset + + Notes + ----- + Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). + Duplicate pressure levels return duplicate parcel temperatures. Consider preprocessing + low-precision, high frequency profiles with tools like `scipy.medfilt`, + `pandas.drop_duplicates`, or `numpy.unique`. + + Will only return Pint Quantities, even when given xarray DataArray profiles. To + obtain a xarray Dataset instead, use `parcel_profile_with_lcl_as_dataset` instead. + + .. versionchanged:: 1.0 + Renamed ``dewpt`` parameter to ``dewpoint`` + + """ + p_l, p_lcl, p_u, t_l, t_lcl, t_u = _parcel_profile_helper_linfel(pressure, temperature[0], + dewpoint[0]) + new_press = concatenate((p_l, p_lcl, p_u)) + prof_temp = concatenate((t_l, t_lcl, t_u)) + new_temp = _insert_lcl_level(pressure, temperature, p_lcl) + new_dewp = _insert_lcl_level(pressure, dewpoint, p_lcl) + return new_press, new_temp, new_dewp, prof_temp + + @exporter.export def parcel_profile_with_lcl_as_dataset(pressure, temperature, dewpoint): r"""Calculate the profile a parcel takes through the atmosphere, returning a Dataset. @@ -1556,6 +1645,51 @@ def _parcel_profile_helper(pressure, temperature, dewpoint): temp_lower[:-1], temp_lcl, temp_upper[1:]) +def _parcel_profile_helper_linfel(pressure, temperature, dewpoint): + """ Linfeng's version of _parcel_profile_helper. Added on Jul 1 2025 + """ + """Help calculate parcel profiles. + + Returns the temperature and pressure, above, below, and including the LCL. The + other calculation functions decide what to do with the pieces. + + """ + # Check that pressure does not increase. + _check_pressure_error(pressure) + + # Find the LCL + press_lcl, temp_lcl = lcl_linfel(pressure[0], temperature, dewpoint) + press_lcl = press_lcl.to(pressure.units) + + # Find the dry adiabatic profile, *including* the LCL. We need >= the LCL in case the + # LCL is included in the levels. It's slightly redundant in that case, but simplifies + # the logic for removing it later. + press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl)) + temp_lower = dry_lapse_linfel(press_lower, temperature) + + # If the pressure profile doesn't make it to the lcl, we can stop here + if _greater_or_close(np.nanmin(pressure), press_lcl): + return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units), + temp_lower[:-1], temp_lcl, units.Quantity(np.array([]), temp_lower.units)) + + # Establish profile above LCL + press_upper = concatenate((press_lcl, pressure[pressure < press_lcl])) + + # Check duplicate pressure values from remaining profile. + # unique will return remaining values sorted ascending. + unique, indices, counts = np.unique(press_upper.m, return_inverse=True, return_counts=True) + unique = units.Quantity(unique, press_upper.units) + if np.any(counts > 1): + _warnings.warn(f'Duplicate pressure(s) {unique[counts > 1]:~P} provided. ' + 'Output profile includes duplicate temperatures as a result.') + + # Find moist pseudo-adiabatic profile starting at the LCL + temp_upper = moist_lapse_linfel(press_upper, temp_lower[-1]).to(temp_lower.units) + + # Return profile pieces + return (press_lower[:-1], press_lcl, press_upper[1:], + temp_lower[:-1], temp_lcl, temp_upper[1:]) + def _insert_lcl_level(pressure, temperature, lcl_pressure): """Insert the LCL pressure into the profile.""" interp_temp = interpolate_1d(lcl_pressure, pressure, temperature) From 64d7247b1e026552ea0c73142f453caeeeb15b9b Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 1 Jul 2025 17:30:29 -0600 Subject: [PATCH 43/48] change lfc and el to use the underlying c++ function to compute parcel profile and lcl --- src/metpy/calc/thermo.py | 268 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 268 insertions(+) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index d62d2eef53..66853909a5 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1059,6 +1059,161 @@ def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoi dewpoint, intersect_type='LFC') +@exporter.export +@preprocess_and_wrap() +@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') +def lfc_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoint_start=None, + which='top'): + """ + Linfeng's version of 'lfc'. Added on Jul 1 2025 + """ + r"""Calculate the level of free convection (LFC). + + This works by finding the first intersection of the ideal parcel path and + the measured parcel temperature. If this intersection occurs below the LCL, + the LFC is determined to be the same as the LCL, based upon the conditions + set forth in [USAF1990]_, pg 4-14, where a parcel must be lifted dry adiabatically + to saturation before it can freely rise. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure profile. This array must be from high to low pressure. + + temperature : `pint.Quantity` + Temperature at the levels given by `pressure` + + dewpoint : `pint.Quantity` + Dewpoint at the levels given by `pressure` + + parcel_temperature_profile: `pint.Quantity`, optional + The parcel's temperature profile from which to calculate the LFC. Defaults to the + surface parcel profile. + + dewpoint_start: `pint.Quantity`, optional + Dewpoint of the parcel for which to calculate the LFC. Defaults to the surface + dewpoint. + + which: str, optional + Pick which LFC to return. Options are 'top', 'bottom', 'wide', 'most_cape', and 'all'; + 'top' returns the lowest-pressure LFC (default), + 'bottom' returns the highest-pressure LFC, + 'wide' returns the LFC whose corresponding EL is farthest away, + 'most_cape' returns the LFC that results in the most CAPE in the profile. + + Returns + ------- + `pint.Quantity` + LFC pressure, or array of same if which='all' + + `pint.Quantity` + LFC temperature, or array of same if which='all' + + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, lfc + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # calculate LFC + >>> lfc(p, T, Td) + (, ) + + See Also + -------- + parcel_profile + + Notes + ----- + Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). + Since this function returns scalar values when given a profile, this will return Pint + Quantities even when given xarray DataArray profiles. + + .. versionchanged:: 1.0 + Renamed ``dewpt``,``dewpoint_start`` parameters to ``dewpoint``, ``dewpoint_start`` + + """ + # Default to surface parcel if no profile or starting pressure level is given + if parcel_temperature_profile is None: + pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint) + new_profile = parcel_profile_with_lcl_linfel(pressure, temperature, dewpoint) + pressure, temperature, dewpoint, parcel_temperature_profile = new_profile + parcel_temperature_profile = parcel_temperature_profile.to(temperature.units) + else: + new_profile = _remove_nans(pressure, temperature, dewpoint, parcel_temperature_profile) + pressure, temperature, dewpoint, parcel_temperature_profile = new_profile + + if dewpoint_start is None: + dewpoint_start = dewpoint[0] + + # The parcel profile and data may have the same first data point. + # If that is the case, ignore that point to get the real first + # intersection for the LFC calculation. Use logarithmic interpolation. + if np.isclose(parcel_temperature_profile[0].to(temperature.units).m, temperature[0].m): + x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], + temperature[1:], direction='increasing', log_x=True) + else: + x, y = find_intersections(pressure, parcel_temperature_profile, + temperature, direction='increasing', log_x=True) + + # Compute LCL for this parcel for future comparisons + this_lcl = lcl_linfel(pressure[0], parcel_temperature_profile[0], dewpoint_start) + + # The LFC could: + # 1) Not exist + # 2) Exist but be equal to the LCL + # 3) Exist and be above the LCL + + # LFC does not exist or is LCL + if len(x) == 0: + # Is there any positive area above the LCL? + mask = pressure < this_lcl[0] + if np.all(_less_or_close(parcel_temperature_profile[mask], temperature[mask])): + # LFC doesn't exist + x = units.Quantity(np.nan, pressure.units) + y = units.Quantity(np.nan, temperature.units) + else: # LFC = LCL + x, y = this_lcl + return x, y + + # LFC exists. Make sure it is no lower than the LCL + else: + idx = x < this_lcl[0] + # LFC height < LCL height, so set LFC = LCL + if not any(idx): + el_pressure, _ = find_intersections(pressure[1:], parcel_temperature_profile[1:], + temperature[1:], direction='decreasing', + log_x=True) + if el_pressure.size and np.min(el_pressure) > this_lcl[0]: + # EL exists and it is below the LCL + x = units.Quantity(np.nan, pressure.units) + y = units.Quantity(np.nan, temperature.units) + else: + # EL exists and it is above the LCL or the EL does not exist + x, y = this_lcl + return x, y + # Otherwise, find all LFCs that exist above the LCL + # What is returned depends on which flag as described in the docstring + else: + return _multiple_el_lfc_options(x, y, idx, which, pressure, + parcel_temperature_profile, temperature, + dewpoint, intersect_type='LFC') + + def _multiple_el_lfc_options(intersect_pressures, intersect_temperatures, valid_x, which, pressure, parcel_temperature_profile, temperature, dewpoint, intersect_type): @@ -1238,6 +1393,119 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which=' units.Quantity(np.nan, temperature.units)) +@exporter.export +@preprocess_and_wrap() +@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') +def el_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='top'): + """ + Linfeng's version of 'el'. Added on Jul 1 2025 + """ + r"""Calculate the equilibrium level. + + This works by finding the last intersection of the ideal parcel path and + the measured environmental temperature. If there is one or fewer intersections, there is + no equilibrium level. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure profile. This array must be from high to low pressure. + + temperature : `pint.Quantity` + Temperature at the levels given by `pressure` + + dewpoint : `pint.Quantity` + Dewpoint at the levels given by `pressure` + + parcel_temperature_profile: `pint.Quantity`, optional + The parcel's temperature profile from which to calculate the EL. Defaults to the + surface parcel profile. + + which: str, optional + Pick which EL to return. Options are 'top', 'bottom', 'wide', 'most_cape', and 'all'. + 'top' returns the lowest-pressure EL, default. + 'bottom' returns the highest-pressure EL. + 'wide' returns the EL whose corresponding LFC is farthest away. + 'most_cape' returns the EL that results in the most CAPE in the profile. + + Returns + ------- + `pint.Quantity` + EL pressure, or array of same if which='all' + + `pint.Quantity` + EL temperature, or array of same if which='all' + + Examples + -------- + >>> from metpy.calc import el, dewpoint_from_relative_humidity, parcel_profile + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compute parcel profile temperature + >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') + >>> # calculate EL + >>> el(p, T, Td, prof) + (, ) + + See Also + -------- + parcel_profile + + Notes + ----- + Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). + Since this function returns scalar values when given a profile, this will return Pint + Quantities even when given xarray DataArray profiles. + + .. versionchanged:: 1.0 + Renamed ``dewpt`` parameter to ``dewpoint`` + + """ + # Default to surface parcel if no profile or starting pressure level is given + if parcel_temperature_profile is None: + pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint) + new_profile = parcel_profile_with_lcl_linfel(pressure, temperature, dewpoint) + pressure, temperature, dewpoint, parcel_temperature_profile = new_profile + parcel_temperature_profile = parcel_temperature_profile.to(temperature.units) + else: + new_profile = _remove_nans(pressure, temperature, dewpoint, parcel_temperature_profile) + pressure, temperature, dewpoint, parcel_temperature_profile = new_profile + + # If the top of the sounding parcel is warmer than the environment, there is no EL + if parcel_temperature_profile[-1] > temperature[-1]: + return (units.Quantity(np.nan, pressure.units), + units.Quantity(np.nan, temperature.units)) + + # Interpolate in log space to find the appropriate pressure - units have to be stripped + # and reassigned to allow np.log() to function properly. + x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], temperature[1:], + direction='decreasing', log_x=True) + lcl_p, _ = lcl_linfel(pressure[0], temperature[0], dewpoint[0]) + if len(x) > 0 and x[-1] < lcl_p: + idx = x < lcl_p + return _multiple_el_lfc_options(x, y, idx, which, pressure, + parcel_temperature_profile, temperature, dewpoint, + intersect_type='EL') + else: + return (units.Quantity(np.nan, pressure.units), + units.Quantity(np.nan, temperature.units)) + + @exporter.export @preprocess_and_wrap(wrap_like='pressure') @check_units('[pressure]', '[temperature]', '[temperature]') From 387bf1df6040b4fcbd314cf788169267435d60f6 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 2 Jul 2025 13:03:48 -0600 Subject: [PATCH 44/48] modified python cape_cin to use c++ functions to compute --- src/metpy/calc/thermo.py | 227 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 225 insertions(+), 2 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 66853909a5..3dba2e93f6 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2343,6 +2343,65 @@ def saturation_mixing_ratio(total_press, temperature, *, phase='liquid'): 'equilibrium.') return np.where(undefined, np.nan, mixing_ratio._nounit(e_s, total_press)) +@exporter.export +@preprocess_and_wrap(wrap_like='temperature', broadcast=('total_press', 'temperature')) +@process_units( + {'total_press': '[pressure]', 'temperature': '[temperature]'}, + '[dimensionless]' +) +def saturation_mixing_ratio_linfel(total_press, temperature, *, phase='liquid'): + """ + Linfeng's version of 'saturation_mixing_ratio'. Added on Jul 2 2025 + """ + r"""Calculate the saturation mixing ratio of water vapor. + + This calculation is given total atmospheric pressure and air temperature. + + Parameters + ---------- + total_press: `pint.Quantity` + Total atmospheric pressure + + temperature: `pint.Quantity` + Air temperature + + phase : {'liquid', 'solid', 'auto'} + Where applicable, adjust assumptions and constants to make calculation valid in + ``'liquid'`` water (default) or ``'solid'`` ice regimes. ``'auto'`` will change regime + based on determination of phase boundaries, eg `temperature` relative to freezing. + + Returns + ------- + `pint.Quantity` + Saturation mixing ratio, dimensionless + + Examples + -------- + >>> from metpy.calc import saturation_mixing_ratio + >>> from metpy.units import units + >>> saturation_mixing_ratio(983 * units.hPa, 25 * units.degC).to('g/kg') + + + Notes + ----- + This function is a straightforward implementation of the equation given in many places, + such as [Hobbs1977]_ pg.73: + + .. math:: r_s = \epsilon \frac{e_s}{p - e_s} + + By definition, this value is only defined for conditions where the saturation vapor + pressure (:math:`e_s`) for the given temperature is less than the given total pressure + (:math:`p`). Otherwise, liquid phase water cannot exist in equilibrium and there is only + water vapor present. For any value pairs that fall under this condition, the function will + warn and return NaN. + + .. versionchanged:: 1.0 + Renamed ``tot_press`` parameter to ``total_press`` + + """ + validate_choice({'liquid', 'solid', 'auto'}, phase=phase) + return _calc_mod.saturation_mixing_ratio(total_press, temperature, phase) + @exporter.export @preprocess_and_wrap( @@ -3234,7 +3293,6 @@ def relative_humidity_from_specific_humidity( mixing_ratio_from_specific_humidity(specific_humidity), phase=phase) - @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') @@ -3341,7 +3399,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature, dewpoint, parcel_profile) - pressure_lcl, _ = lcl(pressure[0], temperature[0], dewpoint[0]) + pressure_lcl, _ = lcl_linfel(pressure[0], temperature[0], dewpoint[0]) below_lcl = pressure > pressure_lcl # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated @@ -3400,6 +3458,171 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' cin = units.Quantity(0, 'J/kg') return cape, cin +@exporter.export +@preprocess_and_wrap() +@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') +def cape_cin_linfel(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom', + which_el='top'): + r"""Calculate CAPE and CIN. + + Calculate the convective available potential energy (CAPE) and convective inhibition (CIN) + of a given upper air profile and parcel path. CIN is integrated between the surface and + LFC, CAPE is integrated between the LFC and EL (or top of sounding). Intersection points + of the measured temperature profile and parcel profile are logarithmically interpolated. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure level(s) of interest, in order from highest to + lowest pressure + + temperature : `pint.Quantity` + Atmospheric temperature corresponding to pressure + + dewpoint : `pint.Quantity` + Atmospheric dewpoint corresponding to pressure + + parcel_profile : `pint.Quantity` + Temperature profile of the parcel + + which_lfc : str + Choose which LFC to integrate from. Valid options are 'top', 'bottom', 'wide', + and 'most_cape'. Default is 'bottom'. + + which_el : str + Choose which EL to integrate to. Valid options are 'top', 'bottom', 'wide', + and 'most_cape'. Default is 'top'. + + Returns + ------- + `pint.Quantity` + Convective Available Potential Energy (CAPE) + + `pint.Quantity` + Convective Inhibition (CIN) + + Examples + -------- + >>> from metpy.calc import cape_cin, dewpoint_from_relative_humidity, parcel_profile + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compture parcel temperature + >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') + >>> # calculate surface based CAPE/CIN + >>> cape_cin(p, T, Td, prof) + (, ) + + See Also + -------- + lfc, el + + Notes + ----- + Formula adopted from [Hobbs1977]_. + + .. math:: \text{CAPE} = -R_d \int_{LFC}^{EL} + (T_{{v}_{parcel}} - T_{{v}_{env}}) d\text{ln}(p) + + .. math:: \text{CIN} = -R_d \int_{SFC}^{LFC} + (T_{{v}_{parcel}} - T_{{v}_{env}}) d\text{ln}(p) + + + * :math:`CAPE` is convective available potential energy + * :math:`CIN` is convective inhibition + * :math:`LFC` is pressure of the level of free convection + * :math:`EL` is pressure of the equilibrium level + * :math:`SFC` is the level of the surface or beginning of parcel path + * :math:`R_d` is the gas constant + * :math:`g` is gravitational acceleration + * :math:`T_{{v}_{parcel}}` is the parcel virtual temperature + * :math:`T_{{v}_{env}}` is environment virtual temperature + * :math:`p` is atmospheric pressure + + Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). + Since this function returns scalar values when given a profile, this will return Pint + Quantities even when given xarray DataArray profiles. + + .. versionchanged:: 1.0 + Renamed ``dewpt`` parameter to ``dewpoint`` + + """ + pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature, + dewpoint, parcel_profile) + + pressure_lcl, _ = lcl_linfel(pressure[0], temperature[0], dewpoint[0]) + below_lcl = pressure > pressure_lcl + + # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated + # based on the temperature above the LCL + parcel_mixing_ratio = np.where( + below_lcl, + saturation_mixing_ratio_linfel(pressure[0], dewpoint[0]), + saturation_mixing_ratio_linfel(pressure, parcel_profile) + ) + + # Convert the temperature/parcel profile to virtual temperature + temperature = virtual_temperature_from_dewpoint_linfel(pressure, temperature, dewpoint) + parcel_profile = virtual_temperature_linfel(parcel_profile, parcel_mixing_ratio) + + # Calculate LFC limit of integration + lfc_pressure, _ = lfc_linfel(pressure, temperature, dewpoint, + parcel_temperature_profile=parcel_profile, which=which_lfc) + + # If there is no LFC, no need to proceed. + if np.isnan(lfc_pressure): + return units.Quantity(0, 'J/kg'), units.Quantity(0, 'J/kg') + else: + lfc_pressure = lfc_pressure.magnitude + + # Calculate the EL limit of integration + el_pressure, _ = el_linfel(pressure, temperature, dewpoint, + parcel_temperature_profile=parcel_profile, which=which_el) + + # No EL and we use the top reading of the sounding. + el_pressure = pressure[-1].magnitude if np.isnan(el_pressure) else el_pressure.magnitude + + # Difference between the parcel path and measured temperature profiles + y = (parcel_profile - temperature).to(units.degK) + + # Estimate zero crossings + x, y = _find_append_zero_crossings(np.copy(pressure), y) + + # CAPE + # Only use data between the LFC and EL for calculation + p_mask = _less_or_close(x.m, lfc_pressure) & _greater_or_close(x.m, el_pressure) + x_clipped = x[p_mask].magnitude + y_clipped = y[p_mask].magnitude + cape = (mpconsts.Rd + * units.Quantity(trapezoid(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg')) + + # CIN + # Only use data between the surface and LFC for calculation + p_mask = _greater_or_close(x.m, lfc_pressure) + x_clipped = x[p_mask].magnitude + y_clipped = y[p_mask].magnitude + cin = (mpconsts.Rd + * units.Quantity(trapezoid(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg')) + + # Set CIN to 0 if it's returned as a positive value (#1190) + if cin > units.Quantity(0, 'J/kg'): + cin = units.Quantity(0, 'J/kg') + return cape, cin + def _find_append_zero_crossings(x, y): r""" From 35ed80352922f60eac5429fbe048478f22bce651 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 2 Jul 2025 22:07:27 -0600 Subject: [PATCH 45/48] wip --- src/metpy/calc/thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 3dba2e93f6..8f2d21f9ec 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -3399,7 +3399,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature, dewpoint, parcel_profile) - pressure_lcl, _ = lcl_linfel(pressure[0], temperature[0], dewpoint[0]) + pressure_lcl, _ = lcl(pressure[0], temperature[0], dewpoint[0]) below_lcl = pressure > pressure_lcl # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated From b18fd24897b45552a48923a0c874880df9b8652c Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 15 Jul 2025 21:41:42 -0600 Subject: [PATCH 46/48] change to pybind11 vectorization for dry_lapse --- src/calcmod.cpp | 6 +--- src/thermo.cpp | 89 ------------------------------------------------- src/thermo.hpp | 6 ---- 3 files changed, 1 insertion(+), 100 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 8c2cdbe64d..903e19b164 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -40,14 +40,10 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate relative humidity from temperature and dewpoint.", py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); - m.def("dry_lapse", &DryLapseVectorized, + m.def("dry_lapse", py::vectorize(DryLapse), "Calculate the temperature along a pressure profile assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); - m.def("dry_lapse_3d", &DryLapseVectorized_3D, - "Calculate the temperature along multiple pressure profiles assuming dry adiabatic process.", - py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); - m.def("moist_lapse", &MoistLapseVectorized, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); diff --git a/src/thermo.cpp b/src/thermo.cpp index 9f04ab82c8..1ad00ee53b 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -60,95 +60,6 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, } -py::array_t DryLapseVectorized(py::array_t pressure, - py::array_t ref_temperature, - double ref_pressure) { - // This function calculates the dry adiabatic profile for multiple starting - // temperatures (2D surface) and a single communal starting pressure, along a - // 1D pressure profile. - // --- Step 1: Prepare the C++ vector for pressure levels --- - if (pressure.ndim() > 1) { - throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); - } - std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); - - // --- Step 2: Ensure the reference temperature array is contiguous --- - auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); - - // --- Step 3: Define the shape of the output array: (N+1) dimension--- - std::vector out_shape; - for(int i = 0; i < ref_temp_contig.ndim(); ++i) { - out_shape.push_back(ref_temp_contig.shape(i)); - } - size_t profile_len = pressure_vec.size(); - out_shape.push_back(profile_len); - - auto out_array = py::array_t(out_shape); - - // --- Step 4: Get direct pointers to data buffers for fast access --- - const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); - double* out_array_ptr = out_array.mutable_data(); - size_t num_profiles = ref_temp_contig.size(); - - // --- Step 5: Loop through each reference temperature --- - for (size_t i = 0; i < num_profiles; ++i) { - for (size_t j = 0; j < profile_len; ++j) { - out_array_ptr[i * profile_len + j] = DryLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure); - } - } - - return out_array; -} - -py::array_t DryLapseVectorized_3D(py::array_t pressure, - py::array_t ref_temperature, - py::array_t ref_pressure) { - // --- Step 1: Ensure all input arrays are using a contiguous memory layout --- - auto p_contig = py::array::ensure(pressure, py::array::c_style); - auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); - auto ref_press_contig = py::array::ensure(ref_pressure, py::array::c_style); - - // --- Step 2: Perform comprehensive shape validation --- - if (ref_temp_contig.ndim() != ref_press_contig.ndim()) { - throw std::runtime_error("Input 'ref_temperature' and 'ref_pressure' must have the same number of dimensions."); - } - if (p_contig.ndim() != ref_temp_contig.ndim() + 1) { - throw std::runtime_error("Input 'pressure' must have one more dimension than 'ref_temperature'."); - } - for (int i = 0; i < ref_temp_contig.ndim(); ++i) { - if (ref_temp_contig.shape(i) != ref_press_contig.shape(i) || - p_contig.shape(i+1) != ref_temp_contig.shape(i)) { - throw std::runtime_error("The horizontal dimensions of all input arrays must match."); - } - } - - // --- Step 3: Define the shape of the output array --- - // The output shape will be identical to the input pressure array's shape. - auto out_array = py::array_t(p_contig.request().shape); - - // --- Step 4: Get direct pointers to data buffers for fast access --- - const double* pressure_ptr = static_cast(p_contig.request().ptr); - const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); - const double* ref_press_ptr = static_cast(ref_press_contig.request().ptr); - double* out_array_ptr = out_array.mutable_data(); - - // --- Step 5: Define loop boundaries --- - size_t num_profiles = ref_temp_contig.size(); // Total number of horizontal points - size_t profile_len = p_contig.shape(0); // Length of the vertical dimension - - // --- Step 6: Loop through each horizontal point and its vertical profile --- - for (int j = 0; j < profile_len; ++j) { - for (int i = 0; i < num_profiles; ++i) { - // Calculate the index for the current point in the flattened arrays. - out_array_ptr[i+j*num_profiles] = DryLapse(pressure_ptr[i+j*num_profiles], - ref_temp_ptr[i], - ref_press_ptr[i]); - } - } - - return out_array; -} - double CaldlnTdlnP(double temperature, double pressure) { // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. double rs = SaturationMixingRatio(pressure, temperature, "liquid"); diff --git a/src/thermo.hpp b/src/thermo.hpp index 555bbc437a..0693ca4834 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -24,12 +24,6 @@ double DryLapse(double pressure, double ref_temperature, double ref_pressure); std::vector DryLapseProfile(const std::vector& pressure_profile, double ref_temperature, double ref_pressure); -py::array_t DryLapseVectorized(py::array_t pressure, - py::array_t ref_temperature, - double ref_pressure); -py::array_t DryLapseVectorized_3D(py::array_t pressure, - py::array_t ref_temperature, - py::array_t ref_pressure); double CaldlnTdlnP(double temperature, double pressure); double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep); From d24613d05a2105d7ee02d706691ddcc82f075b71 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 3 Jul 2025 11:50:49 -0600 Subject: [PATCH 47/48] replace original python version with functions using c++ version, clean all the function names in thermo.py --- src/metpy/calc/thermo.py | 1008 ++++---------------------------------- 1 file changed, 101 insertions(+), 907 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 8f2d21f9ec..67820ee4cd 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -455,14 +455,23 @@ def temperature_from_potential_temperature(pressure, potential_temperature): """ return potential_temperature * exner_function(pressure) - @exporter.export @preprocess_and_wrap( wrap_like='temperature', broadcast=('pressure', 'temperature', 'reference_pressure') ) -@check_units('[pressure]', '[temperature]', '[pressure]') +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'reference_pressure': '[pressure]' + }, + '[temperature]' +) def dry_lapse(pressure, temperature, reference_pressure=None, vertical_dim=0): + """ + Linfeng's version of 'dry_lapse'. Added on Jun18 2025 + """ r"""Calculate the temperature at a level assuming only dry processes. This function lifts a parcel starting at ``temperature``, conserving @@ -507,36 +516,11 @@ def dry_lapse(pressure, temperature, reference_pressure=None, vertical_dim=0): .. versionchanged:: 1.0 Renamed ``ref_pressure`` parameter to ``reference_pressure`` - """ - if reference_pressure is None: - reference_pressure = pressure[0] - return temperature * (pressure / reference_pressure)**mpconsts.kappa - - - -@exporter.export -@preprocess_and_wrap( - wrap_like='temperature', - broadcast=('pressure', 'temperature', 'reference_pressure') -) -@process_units( - { - 'pressure': '[pressure]', - 'temperature': '[temperature]', - 'reference_pressure': '[pressure]' - }, - '[temperature]' -) -def dry_lapse_linfel(pressure, temperature, reference_pressure=None, vertical_dim=0): - """ - Linfeng's version of 'dry_lapse'. Added on Jun18 2025 """ if reference_pressure is None: reference_pressure = pressure[0] return _calc_mod.dry_lapse(pressure, temperature, reference_pressure) - - @exporter.export @preprocess_and_wrap( wrap_like='temperature', @@ -551,6 +535,12 @@ def dry_lapse_linfel(pressure, temperature, reference_pressure=None, vertical_di '[temperature]' ) def moist_lapse(pressure, temperature, reference_pressure=None): + """ + Linfeng's version of 'moist_lapse'. Added on Jun 25 2025 + This function calculates the moist adiabatic profile for multiple starting + temperatures (2D surface) and a single communal starting pressure, along a + 1D pressure profile. + """ r"""Calculate the temperature at a level assuming liquid saturation processes. This function lifts a parcel starting at `temperature`. The starting pressure can @@ -604,104 +594,12 @@ def moist_lapse(pressure, temperature, reference_pressure=None): .. versionchanged:: 1.0 Renamed ``ref_pressure`` parameter to ``reference_pressure`` - """ - def dt(p, t): - rs = saturation_mixing_ratio._nounit(p, t) - frac = ( - (mpconsts.nounit.Rd * t + mpconsts.nounit.Lv * rs) - / (mpconsts.nounit.Cp_d + ( - mpconsts.nounit.Lv * mpconsts.nounit.Lv * rs * mpconsts.nounit.epsilon - / (mpconsts.nounit.Rd * t**2) - )) - ) - return frac / p - - temperature = np.atleast_1d(temperature) - pressure = np.atleast_1d(pressure) - if reference_pressure is None: - reference_pressure = pressure[0] - - if np.isnan(reference_pressure) or np.all(np.isnan(temperature)): - return np.full((temperature.size, pressure.size), np.nan) - - pres_decreasing = (pressure[0] > pressure[-1]) - if pres_decreasing: - # Everything is easier if pressures are in increasing order - pressure = pressure[::-1] - - # It would be preferable to use a regular solver like RK45, but as of scipy 1.8.0 - # anything other than LSODA goes into an infinite loop when given NaNs for y0. - solver_args = {'fun': dt, 'y0': temperature, - 'method': 'LSODA', 'atol': 1e-7, 'rtol': 1.5e-8} - - # Need to handle close points to avoid an error in the solver - close = np.isclose(pressure, reference_pressure) - if np.any(close): - ret = np.broadcast_to(temperature[:, np.newaxis], (temperature.size, np.sum(close))) - else: - ret = np.empty((temperature.size, 0), dtype=temperature.dtype) - - # Do we have any points above the reference pressure - points_above = (pressure < reference_pressure) & ~close - if np.any(points_above): - # Integrate upward--need to flip so values are properly ordered from ref to min - press_side = pressure[points_above][::-1] - - # Flip on exit so t values correspond to increasing pressure - result = si.solve_ivp(t_span=(reference_pressure, press_side[-1]), - t_eval=press_side, **solver_args) - if result.success: - ret = np.concatenate((result.y[..., ::-1], ret), axis=-1) - else: - raise ValueError('ODE Integration failed. This is likely due to trying to ' - 'calculate at too small values of pressure.') - - # Do we have any points below the reference pressure - points_below = ~points_above & ~close - if np.any(points_below): - # Integrate downward - press_side = pressure[points_below] - result = si.solve_ivp(t_span=(reference_pressure, press_side[-1]), - t_eval=press_side, **solver_args) - if result.success: - ret = np.concatenate((ret, result.y), axis=-1) - else: - raise ValueError('ODE Integration failed. This is likely due to trying to ' - 'calculate at too small values of pressure.') - - if pres_decreasing: - ret = ret[..., ::-1] - - return ret.squeeze() - - -@exporter.export -@preprocess_and_wrap( - wrap_like='temperature', - broadcast=('pressure', 'temperature', 'reference_pressure') -) -@process_units( - { - 'pressure': '[pressure]', - 'temperature': '[temperature]', - 'reference_pressure': '[pressure]' - }, - '[temperature]' -) -def moist_lapse_linfel(pressure, temperature, reference_pressure=None): - """ - Linfeng's version of 'moist_lapse'. Added on Jun 25 2025 - This function calculates the moist adiabatic profile for multiple starting - temperatures (2D surface) and a single communal starting pressure, along a - 1D pressure profile. """ if reference_pressure is None: reference_pressure = pressure[0] # nstep for RK4 is set to 30 return _calc_mod.moist_lapse(pressure, temperature, reference_pressure, 30) - - @exporter.export @preprocess_and_wrap() @process_units( @@ -709,6 +607,9 @@ def moist_lapse_linfel(pressure, temperature, reference_pressure=None): ('[pressure]', '[temperature]') ) def lcl(pressure, temperature, dewpoint, max_iters=None, eps=None): + """ + Linfeng's version of 'lcl'. Added on Jun23 2025 + """ r"""Calculate the lifted condensation level (LCL) from the starting point. The starting state for the parcel is defined by `temperature`, `dewpoint`, @@ -761,41 +662,6 @@ def lcl(pressure, temperature, dewpoint, max_iters=None, eps=None): .. versionchanged:: 1.0 Renamed ``dewpt`` parameter to ``dewpoint`` - """ - if max_iters or eps: - _warnings.warn( - 'max_iters, eps arguments unused and will be deprecated in a future version.', - PendingDeprecationWarning) - - q = specific_humidity_from_dewpoint._nounit(pressure, dewpoint, phase='liquid') - moist_heat_ratio = (moist_air_specific_heat_pressure._nounit(q) - / moist_air_gas_constant._nounit(q)) - spec_heat_diff = mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v - - a = moist_heat_ratio + spec_heat_diff / mpconsts.nounit.Rv - b = (-(mpconsts.nounit.Lv + spec_heat_diff * mpconsts.nounit.T0) - / (mpconsts.nounit.Rv * temperature)) - c = b / a - - w_minus1 = lambertw( - (relative_humidity_from_dewpoint._nounit(temperature, dewpoint, phase='liquid') - ** (1 / a) * c * np.exp(c)), k=-1).real - - t_lcl = c / w_minus1 * temperature - p_lcl = pressure * (t_lcl / temperature) ** moist_heat_ratio - - return p_lcl, t_lcl - - -@exporter.export -@preprocess_and_wrap() -@process_units( - {'pressure': '[pressure]', 'temperature': '[temperature]', 'dewpoint': '[temperature]'}, - ('[pressure]', '[temperature]') -) -def lcl_linfel(pressure, temperature, dewpoint, max_iters=None, eps=None): - """ - Linfeng's version of 'lcl'. Added on Jun23 2025 """ pressure, temperature, dewpoint = np.atleast_1d(pressure, temperature, dewpoint) p_lcl, t_lcl = _calc_mod.lcl(pressure, temperature, dewpoint) @@ -906,164 +772,11 @@ def ccl(pressure, temperature, dewpoint, height=None, mixed_layer_depth=None, wh x, y = x.to(pressure.units), y.to(temperature.units) return x, y, dry_lapse(pressure[0], y, x).to(temperature.units) - @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoint_start=None, which='top'): - r"""Calculate the level of free convection (LFC). - - This works by finding the first intersection of the ideal parcel path and - the measured parcel temperature. If this intersection occurs below the LCL, - the LFC is determined to be the same as the LCL, based upon the conditions - set forth in [USAF1990]_, pg 4-14, where a parcel must be lifted dry adiabatically - to saturation before it can freely rise. - - Parameters - ---------- - pressure : `pint.Quantity` - Atmospheric pressure profile. This array must be from high to low pressure. - - temperature : `pint.Quantity` - Temperature at the levels given by `pressure` - - dewpoint : `pint.Quantity` - Dewpoint at the levels given by `pressure` - - parcel_temperature_profile: `pint.Quantity`, optional - The parcel's temperature profile from which to calculate the LFC. Defaults to the - surface parcel profile. - - dewpoint_start: `pint.Quantity`, optional - Dewpoint of the parcel for which to calculate the LFC. Defaults to the surface - dewpoint. - - which: str, optional - Pick which LFC to return. Options are 'top', 'bottom', 'wide', 'most_cape', and 'all'; - 'top' returns the lowest-pressure LFC (default), - 'bottom' returns the highest-pressure LFC, - 'wide' returns the LFC whose corresponding EL is farthest away, - 'most_cape' returns the LFC that results in the most CAPE in the profile. - - Returns - ------- - `pint.Quantity` - LFC pressure, or array of same if which='all' - - `pint.Quantity` - LFC temperature, or array of same if which='all' - - Examples - -------- - >>> from metpy.calc import dewpoint_from_relative_humidity, lfc - >>> from metpy.units import units - >>> # pressure - >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., - ... 550., 500., 450., 400., 350., 300., 250., 200., - ... 175., 150., 125., 100., 80., 70., 60., 50., - ... 40., 30., 25., 20.] * units.hPa - >>> # temperature - >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, - ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, - ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, - ... -56.3, -51.7, -50.7, -47.5] * units.degC - >>> # relative humidity - >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, - ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, - ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless - >>> # calculate dewpoint - >>> Td = dewpoint_from_relative_humidity(T, rh) - >>> # calculate LFC - >>> lfc(p, T, Td) - (, ) - - See Also - -------- - parcel_profile - - Notes - ----- - Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). - Since this function returns scalar values when given a profile, this will return Pint - Quantities even when given xarray DataArray profiles. - - .. versionchanged:: 1.0 - Renamed ``dewpt``,``dewpoint_start`` parameters to ``dewpoint``, ``dewpoint_start`` - - """ - # Default to surface parcel if no profile or starting pressure level is given - if parcel_temperature_profile is None: - pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint) - new_profile = parcel_profile_with_lcl(pressure, temperature, dewpoint) - pressure, temperature, dewpoint, parcel_temperature_profile = new_profile - parcel_temperature_profile = parcel_temperature_profile.to(temperature.units) - else: - new_profile = _remove_nans(pressure, temperature, dewpoint, parcel_temperature_profile) - pressure, temperature, dewpoint, parcel_temperature_profile = new_profile - - if dewpoint_start is None: - dewpoint_start = dewpoint[0] - - # The parcel profile and data may have the same first data point. - # If that is the case, ignore that point to get the real first - # intersection for the LFC calculation. Use logarithmic interpolation. - if np.isclose(parcel_temperature_profile[0].to(temperature.units).m, temperature[0].m): - x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], - temperature[1:], direction='increasing', log_x=True) - else: - x, y = find_intersections(pressure, parcel_temperature_profile, - temperature, direction='increasing', log_x=True) - - # Compute LCL for this parcel for future comparisons - this_lcl = lcl(pressure[0], parcel_temperature_profile[0], dewpoint_start) - - # The LFC could: - # 1) Not exist - # 2) Exist but be equal to the LCL - # 3) Exist and be above the LCL - - # LFC does not exist or is LCL - if len(x) == 0: - # Is there any positive area above the LCL? - mask = pressure < this_lcl[0] - if np.all(_less_or_close(parcel_temperature_profile[mask], temperature[mask])): - # LFC doesn't exist - x = units.Quantity(np.nan, pressure.units) - y = units.Quantity(np.nan, temperature.units) - else: # LFC = LCL - x, y = this_lcl - return x, y - - # LFC exists. Make sure it is no lower than the LCL - else: - idx = x < this_lcl[0] - # LFC height < LCL height, so set LFC = LCL - if not any(idx): - el_pressure, _ = find_intersections(pressure[1:], parcel_temperature_profile[1:], - temperature[1:], direction='decreasing', - log_x=True) - if el_pressure.size and np.min(el_pressure) > this_lcl[0]: - # EL exists and it is below the LCL - x = units.Quantity(np.nan, pressure.units) - y = units.Quantity(np.nan, temperature.units) - else: - # EL exists and it is above the LCL or the EL does not exist - x, y = this_lcl - return x, y - # Otherwise, find all LFCs that exist above the LCL - # What is returned depends on which flag as described in the docstring - else: - return _multiple_el_lfc_options(x, y, idx, which, pressure, - parcel_temperature_profile, temperature, - dewpoint, intersect_type='LFC') - - -@exporter.export -@preprocess_and_wrap() -@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') -def lfc_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoint_start=None, - which='top'): """ Linfeng's version of 'lfc'. Added on Jul 1 2025 """ @@ -1150,7 +863,7 @@ def lfc_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, # Default to surface parcel if no profile or starting pressure level is given if parcel_temperature_profile is None: pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint) - new_profile = parcel_profile_with_lcl_linfel(pressure, temperature, dewpoint) + new_profile = parcel_profile_with_lcl(pressure, temperature, dewpoint) pressure, temperature, dewpoint, parcel_temperature_profile = new_profile parcel_temperature_profile = parcel_temperature_profile.to(temperature.units) else: @@ -1171,7 +884,7 @@ def lfc_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, temperature, direction='increasing', log_x=True) # Compute LCL for this parcel for future comparisons - this_lcl = lcl_linfel(pressure[0], parcel_temperature_profile[0], dewpoint_start) + this_lcl = lcl(pressure[0], parcel_temperature_profile[0], dewpoint_start) # The LFC could: # 1) Not exist @@ -1249,154 +962,46 @@ def _wide_option(intersect_type, p_list, t_list, pressure, parcel_temperature_pr log_x=True) else: # intersect_type == 'EL' el_p_list = p_list - # Find LFC intersection pressure values - lfc_p_list, _ = find_intersections(pressure, parcel_temperature_profile, - temperature, direction='increasing', - log_x=True) - diff = [lfc_p.m - el_p.m for lfc_p, el_p in zip(lfc_p_list, el_p_list, strict=False)] - return (p_list[np.where(diff == np.max(diff))][0], - t_list[np.where(diff == np.max(diff))][0]) - - -def _most_cape_option(intersect_type, p_list, t_list, pressure, temperature, dewpoint, - parcel_temperature_profile): - """Calculate the LFC or EL that produces the most CAPE in the profile.""" - # Need to loop through all possible combinations of cape, find greatest cape profile - cape_list, pair_list = [], [] - for which_lfc in ['top', 'bottom']: - for which_el in ['top', 'bottom']: - cape, _ = cape_cin(pressure, temperature, dewpoint, parcel_temperature_profile, - which_lfc=which_lfc, which_el=which_el) - cape_list.append(cape.m) - pair_list.append([which_lfc, which_el]) - (lfc_chosen, el_chosen) = pair_list[np.where(cape_list == np.max(cape_list))[0][0]] - if intersect_type == 'LFC': - if lfc_chosen == 'top': - x, y = p_list[-1], t_list[-1] - else: # 'bottom' is returned - x, y = p_list[0], t_list[0] - else: # EL is returned - if el_chosen == 'top': - x, y = p_list[-1], t_list[-1] - else: - x, y = p_list[0], t_list[0] - return x, y - - -@exporter.export -@preprocess_and_wrap() -@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') -def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='top'): - r"""Calculate the equilibrium level. - - This works by finding the last intersection of the ideal parcel path and - the measured environmental temperature. If there is one or fewer intersections, there is - no equilibrium level. - - Parameters - ---------- - pressure : `pint.Quantity` - Atmospheric pressure profile. This array must be from high to low pressure. - - temperature : `pint.Quantity` - Temperature at the levels given by `pressure` - - dewpoint : `pint.Quantity` - Dewpoint at the levels given by `pressure` - - parcel_temperature_profile: `pint.Quantity`, optional - The parcel's temperature profile from which to calculate the EL. Defaults to the - surface parcel profile. - - which: str, optional - Pick which EL to return. Options are 'top', 'bottom', 'wide', 'most_cape', and 'all'. - 'top' returns the lowest-pressure EL, default. - 'bottom' returns the highest-pressure EL. - 'wide' returns the EL whose corresponding LFC is farthest away. - 'most_cape' returns the EL that results in the most CAPE in the profile. - - Returns - ------- - `pint.Quantity` - EL pressure, or array of same if which='all' - - `pint.Quantity` - EL temperature, or array of same if which='all' - - Examples - -------- - >>> from metpy.calc import el, dewpoint_from_relative_humidity, parcel_profile - >>> from metpy.units import units - >>> # pressure - >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., - ... 550., 500., 450., 400., 350., 300., 250., 200., - ... 175., 150., 125., 100., 80., 70., 60., 50., - ... 40., 30., 25., 20.] * units.hPa - >>> # temperature - >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, - ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, - ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, - ... -56.3, -51.7, -50.7, -47.5] * units.degC - >>> # relative humidity - >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, - ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, - ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless - >>> # calculate dewpoint - >>> Td = dewpoint_from_relative_humidity(T, rh) - >>> # compute parcel profile temperature - >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') - >>> # calculate EL - >>> el(p, T, Td, prof) - (, ) - - See Also - -------- - parcel_profile - - Notes - ----- - Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). - Since this function returns scalar values when given a profile, this will return Pint - Quantities even when given xarray DataArray profiles. - - .. versionchanged:: 1.0 - Renamed ``dewpt`` parameter to ``dewpoint`` - - """ - # Default to surface parcel if no profile or starting pressure level is given - if parcel_temperature_profile is None: - pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint) - new_profile = parcel_profile_with_lcl(pressure, temperature, dewpoint) - pressure, temperature, dewpoint, parcel_temperature_profile = new_profile - parcel_temperature_profile = parcel_temperature_profile.to(temperature.units) - else: - new_profile = _remove_nans(pressure, temperature, dewpoint, parcel_temperature_profile) - pressure, temperature, dewpoint, parcel_temperature_profile = new_profile + # Find LFC intersection pressure values + lfc_p_list, _ = find_intersections(pressure, parcel_temperature_profile, + temperature, direction='increasing', + log_x=True) + diff = [lfc_p.m - el_p.m for lfc_p, el_p in zip(lfc_p_list, el_p_list, strict=False)] + return (p_list[np.where(diff == np.max(diff))][0], + t_list[np.where(diff == np.max(diff))][0]) + + +def _most_cape_option(intersect_type, p_list, t_list, pressure, temperature, dewpoint, + parcel_temperature_profile): + """Calculate the LFC or EL that produces the most CAPE in the profile.""" + # Need to loop through all possible combinations of cape, find greatest cape profile + cape_list, pair_list = [], [] + for which_lfc in ['top', 'bottom']: + for which_el in ['top', 'bottom']: + cape, _ = cape_cin(pressure, temperature, dewpoint, parcel_temperature_profile, + which_lfc=which_lfc, which_el=which_el) + cape_list.append(cape.m) + pair_list.append([which_lfc, which_el]) + (lfc_chosen, el_chosen) = pair_list[np.where(cape_list == np.max(cape_list))[0][0]] + if intersect_type == 'LFC': + if lfc_chosen == 'top': + x, y = p_list[-1], t_list[-1] + else: # 'bottom' is returned + x, y = p_list[0], t_list[0] + else: # EL is returned + if el_chosen == 'top': + x, y = p_list[-1], t_list[-1] + else: + x, y = p_list[0], t_list[0] + return x, y - # If the top of the sounding parcel is warmer than the environment, there is no EL - if parcel_temperature_profile[-1] > temperature[-1]: - return (units.Quantity(np.nan, pressure.units), - units.Quantity(np.nan, temperature.units)) - # Interpolate in log space to find the appropriate pressure - units have to be stripped - # and reassigned to allow np.log() to function properly. - x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], temperature[1:], - direction='decreasing', log_x=True) - lcl_p, _ = lcl(pressure[0], temperature[0], dewpoint[0]) - if len(x) > 0 and x[-1] < lcl_p: - idx = x < lcl_p - return _multiple_el_lfc_options(x, y, idx, which, pressure, - parcel_temperature_profile, temperature, dewpoint, - intersect_type='EL') - else: - return (units.Quantity(np.nan, pressure.units), - units.Quantity(np.nan, temperature.units)) @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') -def el_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='top'): +def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='top'): """ Linfeng's version of 'el'. Added on Jul 1 2025 """ @@ -1479,7 +1084,7 @@ def el_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, # Default to surface parcel if no profile or starting pressure level is given if parcel_temperature_profile is None: pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint) - new_profile = parcel_profile_with_lcl_linfel(pressure, temperature, dewpoint) + new_profile = parcel_profile_with_lcl(pressure, temperature, dewpoint) pressure, temperature, dewpoint, parcel_temperature_profile = new_profile parcel_temperature_profile = parcel_temperature_profile.to(temperature.units) else: @@ -1495,7 +1100,7 @@ def el_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, # and reassigned to allow np.log() to function properly. x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], temperature[1:], direction='decreasing', log_x=True) - lcl_p, _ = lcl_linfel(pressure[0], temperature[0], dewpoint[0]) + lcl_p, _ = lcl(pressure[0], temperature[0], dewpoint[0]) if len(x) > 0 and x[-1] < lcl_p: idx = x < lcl_p return _multiple_el_lfc_options(x, y, idx, which, pressure, @@ -1506,10 +1111,15 @@ def el_linfel(pressure, temperature, dewpoint, parcel_temperature_profile=None, units.Quantity(np.nan, temperature.units)) + + @exporter.export @preprocess_and_wrap(wrap_like='pressure') @check_units('[pressure]', '[temperature]', '[temperature]') def parcel_profile(pressure, temperature, dewpoint): + """ + Linfeng's version of 'parcel_profile'. Added on Jul 1 2025 + """ r"""Calculate the profile a parcel takes through the atmosphere. The parcel starts at `temperature`, and `dewpoint`, lifted up @@ -1580,117 +1190,14 @@ def parcel_profile(pressure, temperature, dewpoint): Renamed ``dewpt`` parameter to ``dewpoint`` """ - _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint) - return concatenate((t_l, t_u)) - - -@exporter.export -@preprocess_and_wrap(wrap_like='pressure') -@check_units('[pressure]', '[temperature]', '[temperature]') -def parcel_profile_linfel(pressure, temperature, dewpoint): - """ - Linfeng's version of 'parcel_profile'. Added on Jul 1 2025 - """ - _, _, _, t_l, _, t_u = _parcel_profile_helper_linfel(pressure, temperature, dewpoint) + _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint) return concatenate((t_l, t_u)) - @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]') def parcel_profile_with_lcl(pressure, temperature, dewpoint): - r"""Calculate the profile a parcel takes through the atmosphere. - - The parcel starts at `temperature`, and `dewpoint`, lifted up - dry adiabatically to the LCL, and then moist adiabatically from there. - `pressure` specifies the pressure levels for the profile. This function returns - a profile that includes the LCL. - - Parameters - ---------- - pressure : `pint.Quantity` - Atmospheric pressure level(s) of interest. This array must be from - high to low pressure. - - temperature : `pint.Quantity` - Atmospheric temperature at the levels in `pressure`. The first entry should be at - the same level as the first `pressure` data point. - - dewpoint : `pint.Quantity` - Atmospheric dewpoint at the levels in `pressure`. The first entry should be at - the same level as the first `pressure` data point. - - Returns - ------- - pressure : `pint.Quantity` - The parcel profile pressures, which includes the specified levels and the LCL - - ambient_temperature : `pint.Quantity` - Atmospheric temperature values, including the value interpolated to the LCL level - - ambient_dew_point : `pint.Quantity` - Atmospheric dewpoint values, including the value interpolated to the LCL level - - profile_temperature : `pint.Quantity` - The parcel profile temperatures at all of the levels in the returned pressures array, - including the LCL - - Examples - -------- - >>> from metpy.calc import dewpoint_from_relative_humidity, parcel_profile_with_lcl - >>> from metpy.units import units - >>> # pressure - >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., - ... 550., 500., 450., 400., 350., 300., 250., 200., - ... 175., 150., 125., 100., 80., 70., 60., 50., - ... 40., 30., 25., 20.] * units.hPa - >>> # temperature - >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, - ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, - ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, - ... -56.3, -51.7, -50.7, -47.5] * units.degC - >>> # relative humidity - >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, - ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, - ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless - >>> # calculate dewpoint - >>> Td = dewpoint_from_relative_humidity(T, rh) - >>> # compute parcel temperature - >>> Td = dewpoint_from_relative_humidity(T, rh) - >>> p_wLCL, T_wLCL, Td_wLCL, prof_wLCL = parcel_profile_with_lcl(p, T, Td) - - See Also - -------- - lcl, moist_lapse, dry_lapse, parcel_profile, parcel_profile_with_lcl_as_dataset - - Notes - ----- - Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). - Duplicate pressure levels return duplicate parcel temperatures. Consider preprocessing - low-precision, high frequency profiles with tools like `scipy.medfilt`, - `pandas.drop_duplicates`, or `numpy.unique`. - - Will only return Pint Quantities, even when given xarray DataArray profiles. To - obtain a xarray Dataset instead, use `parcel_profile_with_lcl_as_dataset` instead. - - .. versionchanged:: 1.0 - Renamed ``dewpt`` parameter to ``dewpoint`` - - """ - p_l, p_lcl, p_u, t_l, t_lcl, t_u = _parcel_profile_helper(pressure, temperature[0], - dewpoint[0]) - new_press = concatenate((p_l, p_lcl, p_u)) - prof_temp = concatenate((t_l, t_lcl, t_u)) - new_temp = _insert_lcl_level(pressure, temperature, p_lcl) - new_dewp = _insert_lcl_level(pressure, dewpoint, p_lcl) - return new_press, new_temp, new_dewp, prof_temp - - -@exporter.export -@preprocess_and_wrap() -@check_units('[pressure]', '[temperature]', '[temperature]') -def parcel_profile_with_lcl_linfel(pressure, temperature, dewpoint): """ Linfeng's version of 'parcel_profile_with_lcl'. Added on Jul 1 2025 """ @@ -1772,7 +1279,7 @@ def parcel_profile_with_lcl_linfel(pressure, temperature, dewpoint): Renamed ``dewpt`` parameter to ``dewpoint`` """ - p_l, p_lcl, p_u, t_l, t_lcl, t_u = _parcel_profile_helper_linfel(pressure, temperature[0], + p_l, p_lcl, p_u, t_l, t_lcl, t_u = _parcel_profile_helper(pressure, temperature[0], dewpoint[0]) new_press = concatenate((p_l, p_lcl, p_u)) prof_temp = concatenate((t_l, t_lcl, t_u)) @@ -1868,52 +1375,9 @@ def _check_pressure_error(pressure): 'your sounding. Using scipy.signal.medfilt may fix this.') -def _parcel_profile_helper(pressure, temperature, dewpoint): - """Help calculate parcel profiles. - - Returns the temperature and pressure, above, below, and including the LCL. The - other calculation functions decide what to do with the pieces. - - """ - # Check that pressure does not increase. - _check_pressure_error(pressure) - - # Find the LCL - press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpoint) - press_lcl = press_lcl.to(pressure.units) - - # Find the dry adiabatic profile, *including* the LCL. We need >= the LCL in case the - # LCL is included in the levels. It's slightly redundant in that case, but simplifies - # the logic for removing it later. - press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl)) - temp_lower = dry_lapse(press_lower, temperature) - - # If the pressure profile doesn't make it to the lcl, we can stop here - if _greater_or_close(np.nanmin(pressure), press_lcl): - return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units), - temp_lower[:-1], temp_lcl, units.Quantity(np.array([]), temp_lower.units)) - - # Establish profile above LCL - press_upper = concatenate((press_lcl, pressure[pressure < press_lcl])) - - # Remove duplicate pressure values from remaining profile. Needed for solve_ivp in - # moist_lapse. unique will return remaining values sorted ascending. - unique, indices, counts = np.unique(press_upper.m, return_inverse=True, return_counts=True) - unique = units.Quantity(unique, press_upper.units) - if np.any(counts > 1): - _warnings.warn(f'Duplicate pressure(s) {unique[counts > 1]:~P} provided. ' - 'Output profile includes duplicate temperatures as a result.') - - # Find moist pseudo-adiabatic profile starting at the LCL, reversing above sorting - temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units) - temp_upper = temp_upper[::-1][indices] - - # Return profile pieces - return (press_lower[:-1], press_lcl, press_upper[1:], - temp_lower[:-1], temp_lcl, temp_upper[1:]) -def _parcel_profile_helper_linfel(pressure, temperature, dewpoint): +def _parcel_profile_helper(pressure, temperature, dewpoint): """ Linfeng's version of _parcel_profile_helper. Added on Jul 1 2025 """ """Help calculate parcel profiles. @@ -1926,14 +1390,14 @@ def _parcel_profile_helper_linfel(pressure, temperature, dewpoint): _check_pressure_error(pressure) # Find the LCL - press_lcl, temp_lcl = lcl_linfel(pressure[0], temperature, dewpoint) + press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpoint) press_lcl = press_lcl.to(pressure.units) # Find the dry adiabatic profile, *including* the LCL. We need >= the LCL in case the # LCL is included in the levels. It's slightly redundant in that case, but simplifies # the logic for removing it later. press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl)) - temp_lower = dry_lapse_linfel(press_lower, temperature) + temp_lower = dry_lapse(press_lower, temperature) # If the pressure profile doesn't make it to the lcl, we can stop here if _greater_or_close(np.nanmin(pressure), press_lcl): @@ -1952,7 +1416,7 @@ def _parcel_profile_helper_linfel(pressure, temperature, dewpoint): 'Output profile includes duplicate temperatures as a result.') # Find moist pseudo-adiabatic profile starting at the LCL - temp_upper = moist_lapse_linfel(press_upper, temp_lower[-1]).to(temp_lower.units) + temp_upper = moist_lapse(press_upper, temp_lower[-1]).to(temp_lower.units) # Return profile pieces return (press_lower[:-1], press_lcl, press_upper[1:], @@ -2241,107 +1705,43 @@ def mixing_ratio(partial_press, total_press, molecular_weight_ratio=mpconsts.nou partial_press : `pint.Quantity` Partial pressure of the constituent gas - total_press : `pint.Quantity` - Total air pressure - - molecular_weight_ratio : `pint.Quantity` or float, optional - The ratio of the molecular weight of the constituent gas to that assumed - for air. Defaults to the ratio for water vapor to dry air - (:math:`\epsilon\approx0.622`). - - Returns - ------- - `pint.Quantity` - The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g) - - Examples - -------- - >>> from metpy.calc import mixing_ratio - >>> from metpy.units import units - >>> mixing_ratio(25 * units.hPa, 1000 * units.hPa).to('g/kg') - - - See Also - -------- - saturation_mixing_ratio, vapor_pressure - - Notes - ----- - This function is a straightforward implementation of the equation given in many places, - such as [Hobbs1977]_ pg.73: - - .. math:: r = \epsilon \frac{e}{p - e} - - .. versionchanged:: 1.0 - Renamed ``part_press``, ``tot_press`` parameters to ``partial_press``, ``total_press`` - - """ - # Calling c++ calculation module - return _calc_mod.mixing_ratio(partial_press, total_press, molecular_weight_ratio) - - -@exporter.export -@preprocess_and_wrap(wrap_like='temperature', broadcast=('total_press', 'temperature')) -@process_units( - {'total_press': '[pressure]', 'temperature': '[temperature]'}, - '[dimensionless]' -) -def saturation_mixing_ratio(total_press, temperature, *, phase='liquid'): - r"""Calculate the saturation mixing ratio of water vapor. - - This calculation is given total atmospheric pressure and air temperature. - - Parameters - ---------- - total_press: `pint.Quantity` - Total atmospheric pressure - - temperature: `pint.Quantity` - Air temperature + total_press : `pint.Quantity` + Total air pressure - phase : {'liquid', 'solid', 'auto'} - Where applicable, adjust assumptions and constants to make calculation valid in - ``'liquid'`` water (default) or ``'solid'`` ice regimes. ``'auto'`` will change regime - based on determination of phase boundaries, eg `temperature` relative to freezing. + molecular_weight_ratio : `pint.Quantity` or float, optional + The ratio of the molecular weight of the constituent gas to that assumed + for air. Defaults to the ratio for water vapor to dry air + (:math:`\epsilon\approx0.622`). Returns ------- `pint.Quantity` - Saturation mixing ratio, dimensionless + The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g) Examples -------- - >>> from metpy.calc import saturation_mixing_ratio + >>> from metpy.calc import mixing_ratio >>> from metpy.units import units - >>> saturation_mixing_ratio(983 * units.hPa, 25 * units.degC).to('g/kg') - + >>> mixing_ratio(25 * units.hPa, 1000 * units.hPa).to('g/kg') + + + See Also + -------- + saturation_mixing_ratio, vapor_pressure Notes ----- This function is a straightforward implementation of the equation given in many places, such as [Hobbs1977]_ pg.73: - .. math:: r_s = \epsilon \frac{e_s}{p - e_s} - - By definition, this value is only defined for conditions where the saturation vapor - pressure (:math:`e_s`) for the given temperature is less than the given total pressure - (:math:`p`). Otherwise, liquid phase water cannot exist in equilibrium and there is only - water vapor present. For any value pairs that fall under this condition, the function will - warn and return NaN. + .. math:: r = \epsilon \frac{e}{p - e} .. versionchanged:: 1.0 - Renamed ``tot_press`` parameter to ``total_press`` + Renamed ``part_press``, ``tot_press`` parameters to ``partial_press``, ``total_press`` """ - validate_choice({'liquid', 'solid', 'auto'}, phase=phase) - e_s = saturation_vapor_pressure._nounit(temperature, phase=phase) - undefined = e_s >= total_press - if np.any(undefined): - _warnings.warn('Saturation mixing ratio is undefined for some requested pressure/' - 'temperature combinations. Total pressure must be greater than the ' - 'water vapor saturation pressure for liquid water to be in ' - 'equilibrium.') - return np.where(undefined, np.nan, mixing_ratio._nounit(e_s, total_press)) + # Calling c++ calculation module + return _calc_mod.mixing_ratio(partial_press, total_press, molecular_weight_ratio) @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('total_press', 'temperature')) @@ -2349,7 +1749,7 @@ def saturation_mixing_ratio(total_press, temperature, *, phase='liquid'): {'total_press': '[pressure]', 'temperature': '[temperature]'}, '[dimensionless]' ) -def saturation_mixing_ratio_linfel(total_press, temperature, *, phase='liquid'): +def saturation_mixing_ratio(total_press, temperature, *, phase='liquid'): """ Linfeng's version of 'saturation_mixing_ratio'. Added on Jul 2 2025 """ @@ -2594,7 +1994,6 @@ def wet_bulb_potential_temperature(pressure, temperature, dewpoint): theta_w = units.Quantity(theta_e.m_as('kelvin') - np.exp(a / b), 'kelvin') return np.where(theta_e <= units.Quantity(173.15, 'kelvin'), theta_e, theta_w) - @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'mixing_ratio')) @process_units( @@ -2607,8 +2006,10 @@ def wet_bulb_potential_temperature(pressure, temperature, dewpoint): ignore_inputs_for_output=('molecular_weight_ratio',) ) def virtual_temperature( - temperature, mixing_ratio, molecular_weight_ratio=mpconsts.nounit.epsilon -): + temperature, mixing_ratio, molecular_weight_ratio=mpconsts.nounit.epsilon): + """ + Linfeng's version of 'virtual_temperature'. Added on Jun 30 2025 + """ r"""Calculate virtual temperature. This calculation must be given an air parcel's temperature and mixing ratio. @@ -2647,34 +2048,23 @@ def virtual_temperature( Renamed ``mixing`` parameter to ``mixing_ratio`` """ - return temperature * ((mixing_ratio + molecular_weight_ratio) - / (molecular_weight_ratio * (1 + mixing_ratio))) + return _calc_mod.virtual_temperature( + temperature, mixing_ratio, molecular_weight_ratio) @exporter.export -@preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'mixing_ratio')) +@preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', + 'temperature', + 'dewpoint')) @process_units( { + 'pressure': '[pressure]', 'temperature': '[temperature]', - 'mixing_ratio': '[dimensionless]', + 'dewpoint': '[temperature]', 'molecular_weight_ratio': '[dimensionless]' }, '[temperature]', ignore_inputs_for_output=('molecular_weight_ratio',) ) -def virtual_temperature_linfel( - temperature, mixing_ratio, molecular_weight_ratio=mpconsts.nounit.epsilon): - """ - Linfeng's version of 'virtual_temperature'. Added on Jun 30 2025 - """ - return _calc_mod.virtual_temperature( - temperature, mixing_ratio, molecular_weight_ratio) - - -@exporter.export -@preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', - 'temperature', - 'dewpoint')) -@check_units('[pressure]', '[temperature]', '[temperature]') def virtual_temperature_from_dewpoint( pressure, temperature, @@ -2683,6 +2073,9 @@ def virtual_temperature_from_dewpoint( *, phase='liquid' ): + """ + Linfeng's version of 'virtual_temperature_from_dewpoint'. Added on Jun 30 2025 + """ r"""Calculate virtual temperature. This calculation must be given an air parcel's temperature and mixing ratio. @@ -2730,46 +2123,12 @@ def virtual_temperature_from_dewpoint( """ validate_choice({'liquid', 'solid', 'auto'}, phase=phase) - - # Convert dewpoint to mixing ratio - mixing_ratio = saturation_mixing_ratio(pressure, dewpoint, phase=phase) - - # Calculate virtual temperature with given parameters - return virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio) - - -@exporter.export -@preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', - 'temperature', - 'dewpoint')) -@process_units( - { - 'pressure': '[pressure]', - 'temperature': '[temperature]', - 'dewpoint': '[temperature]', - 'molecular_weight_ratio': '[dimensionless]' - }, - '[temperature]', - ignore_inputs_for_output=('molecular_weight_ratio',) -) -def virtual_temperature_from_dewpoint_linfel( - pressure, - temperature, - dewpoint, - molecular_weight_ratio=mpconsts.nounit.epsilon, - *, - phase='liquid' -): - """ - Linfeng's version of 'virtual_temperature_from_dewpoint'. Added on Jun 30 2025 - """ return _calc_mod.virtual_temperature_from_dewpoint(pressure, temperature, dewpoint, molecular_weight_ratio, phase) - @exporter.export @preprocess_and_wrap( wrap_like='temperature', @@ -3458,171 +2817,6 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' cin = units.Quantity(0, 'J/kg') return cape, cin -@exporter.export -@preprocess_and_wrap() -@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') -def cape_cin_linfel(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom', - which_el='top'): - r"""Calculate CAPE and CIN. - - Calculate the convective available potential energy (CAPE) and convective inhibition (CIN) - of a given upper air profile and parcel path. CIN is integrated between the surface and - LFC, CAPE is integrated between the LFC and EL (or top of sounding). Intersection points - of the measured temperature profile and parcel profile are logarithmically interpolated. - - Parameters - ---------- - pressure : `pint.Quantity` - Atmospheric pressure level(s) of interest, in order from highest to - lowest pressure - - temperature : `pint.Quantity` - Atmospheric temperature corresponding to pressure - - dewpoint : `pint.Quantity` - Atmospheric dewpoint corresponding to pressure - - parcel_profile : `pint.Quantity` - Temperature profile of the parcel - - which_lfc : str - Choose which LFC to integrate from. Valid options are 'top', 'bottom', 'wide', - and 'most_cape'. Default is 'bottom'. - - which_el : str - Choose which EL to integrate to. Valid options are 'top', 'bottom', 'wide', - and 'most_cape'. Default is 'top'. - - Returns - ------- - `pint.Quantity` - Convective Available Potential Energy (CAPE) - - `pint.Quantity` - Convective Inhibition (CIN) - - Examples - -------- - >>> from metpy.calc import cape_cin, dewpoint_from_relative_humidity, parcel_profile - >>> from metpy.units import units - >>> # pressure - >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., - ... 550., 500., 450., 400., 350., 300., 250., 200., - ... 175., 150., 125., 100., 80., 70., 60., 50., - ... 40., 30., 25., 20.] * units.hPa - >>> # temperature - >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, - ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, - ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, - ... -56.3, -51.7, -50.7, -47.5] * units.degC - >>> # relative humidity - >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, - ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, - ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless - >>> # calculate dewpoint - >>> Td = dewpoint_from_relative_humidity(T, rh) - >>> # compture parcel temperature - >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') - >>> # calculate surface based CAPE/CIN - >>> cape_cin(p, T, Td, prof) - (, ) - - See Also - -------- - lfc, el - - Notes - ----- - Formula adopted from [Hobbs1977]_. - - .. math:: \text{CAPE} = -R_d \int_{LFC}^{EL} - (T_{{v}_{parcel}} - T_{{v}_{env}}) d\text{ln}(p) - - .. math:: \text{CIN} = -R_d \int_{SFC}^{LFC} - (T_{{v}_{parcel}} - T_{{v}_{env}}) d\text{ln}(p) - - - * :math:`CAPE` is convective available potential energy - * :math:`CIN` is convective inhibition - * :math:`LFC` is pressure of the level of free convection - * :math:`EL` is pressure of the equilibrium level - * :math:`SFC` is the level of the surface or beginning of parcel path - * :math:`R_d` is the gas constant - * :math:`g` is gravitational acceleration - * :math:`T_{{v}_{parcel}}` is the parcel virtual temperature - * :math:`T_{{v}_{env}}` is environment virtual temperature - * :math:`p` is atmospheric pressure - - Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). - Since this function returns scalar values when given a profile, this will return Pint - Quantities even when given xarray DataArray profiles. - - .. versionchanged:: 1.0 - Renamed ``dewpt`` parameter to ``dewpoint`` - - """ - pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature, - dewpoint, parcel_profile) - - pressure_lcl, _ = lcl_linfel(pressure[0], temperature[0], dewpoint[0]) - below_lcl = pressure > pressure_lcl - - # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated - # based on the temperature above the LCL - parcel_mixing_ratio = np.where( - below_lcl, - saturation_mixing_ratio_linfel(pressure[0], dewpoint[0]), - saturation_mixing_ratio_linfel(pressure, parcel_profile) - ) - - # Convert the temperature/parcel profile to virtual temperature - temperature = virtual_temperature_from_dewpoint_linfel(pressure, temperature, dewpoint) - parcel_profile = virtual_temperature_linfel(parcel_profile, parcel_mixing_ratio) - - # Calculate LFC limit of integration - lfc_pressure, _ = lfc_linfel(pressure, temperature, dewpoint, - parcel_temperature_profile=parcel_profile, which=which_lfc) - - # If there is no LFC, no need to proceed. - if np.isnan(lfc_pressure): - return units.Quantity(0, 'J/kg'), units.Quantity(0, 'J/kg') - else: - lfc_pressure = lfc_pressure.magnitude - - # Calculate the EL limit of integration - el_pressure, _ = el_linfel(pressure, temperature, dewpoint, - parcel_temperature_profile=parcel_profile, which=which_el) - - # No EL and we use the top reading of the sounding. - el_pressure = pressure[-1].magnitude if np.isnan(el_pressure) else el_pressure.magnitude - - # Difference between the parcel path and measured temperature profiles - y = (parcel_profile - temperature).to(units.degK) - - # Estimate zero crossings - x, y = _find_append_zero_crossings(np.copy(pressure), y) - - # CAPE - # Only use data between the LFC and EL for calculation - p_mask = _less_or_close(x.m, lfc_pressure) & _greater_or_close(x.m, el_pressure) - x_clipped = x[p_mask].magnitude - y_clipped = y[p_mask].magnitude - cape = (mpconsts.Rd - * units.Quantity(trapezoid(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg')) - - # CIN - # Only use data between the surface and LFC for calculation - p_mask = _greater_or_close(x.m, lfc_pressure) - x_clipped = x[p_mask].magnitude - y_clipped = y[p_mask].magnitude - cin = (mpconsts.Rd - * units.Quantity(trapezoid(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg')) - - # Set CIN to 0 if it's returned as a positive value (#1190) - if cin > units.Quantity(0, 'J/kg'): - cin = units.Quantity(0, 'J/kg') - return cape, cin - def _find_append_zero_crossings(x, y): r""" From 90e1c85fe8afc6b6717c0fecf3d7b05321bc7ddd Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 3 Jul 2025 12:44:19 -0600 Subject: [PATCH 48/48] wip --- src/metpy/calc/thermo.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 67820ee4cd..3be12dd192 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1483,6 +1483,9 @@ def vapor_pressure(pressure, mixing_ratio): @preprocess_and_wrap(wrap_like='temperature') @process_units({'temperature': '[temperature]'}, '[pressure]') def saturation_vapor_pressure(temperature, *, phase='liquid'): + """ + Linfeng's version of 'saturation_vapor_pressure'. Added on Jul 3 2025 + """ r"""Calculate the saturation (equilibrium) water vapor (partial) pressure. Parameters @@ -1548,6 +1551,9 @@ def saturation_vapor_pressure(temperature, *, phase='liquid'): @preprocess_and_wrap(wrap_like='temperature') @process_units({'temperature': '[temperature]'}, '[pressure]') def _saturation_vapor_pressure_liquid(temperature): + """ + Linfeng's version of '_saturation_vapor_pressure_liquid'. Added on Jul 3 2025 + """ r"""Calculate saturation (equilibrium) water vapor (partial) pressure over liquid water. Parameters @@ -1576,6 +1582,9 @@ def _saturation_vapor_pressure_liquid(temperature): @preprocess_and_wrap(wrap_like='temperature') @process_units({'temperature': '[temperature]'}, '[pressure]') def _saturation_vapor_pressure_solid(temperature): + """ + Linfeng's version of '_saturation_vapor_pressure_solid'. Added on Jul 3 2025 + """ r"""Calculate the saturation water vapor (partial) pressure over solid water (ice). Parameters @@ -1694,6 +1703,9 @@ def dewpoint(vapor_pressure): ignore_inputs_for_output=('molecular_weight_ratio',) ) def mixing_ratio(partial_press, total_press, molecular_weight_ratio=mpconsts.nounit.epsilon): + """ + Linfeng's version of 'mixing_ratio'. Added on Jul 3 2025 + """ r"""Calculate the mixing ratio of a gas. This calculates mixing ratio given its partial pressure and the total pressure of