diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 651ea33..7ee4b69 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,11 +2,8 @@ ecbuild_generate_config_headers( DESTINATION ${INSTALL_INCLUDE_DIR}/vader ) list( APPEND vader_src_files vader/DefaultCookbook.h -vader/vader.h -vader/RecipeBase.h vader/RecipeBase.cc -vader/vader.cc -vader/VaderParameters.h +vader/RecipeBase.h vader/recipes/AirPotentialTemperature_A.cc vader/recipes/AirPotentialTemperature_B.cc vader/recipes/AirPotentialTemperature.h @@ -46,14 +43,26 @@ vader/recipes/DryAirDensity.h vader/recipes/DryAirDensity_A.cc vader/recipes/DryAirDensityLevelsMinusOne.h vader/recipes/DryAirDensityLevelsMinusOne_A.cc +vader/recipes/GeopotentialAtInterface_A.cc +vader/recipes/GeopotentialAtInterface.h +vader/recipes/GeopotentialHeight_A.cc +vader/recipes/GeopotentialHeight.h +vader/recipes/GeopotentialHeightAtInterface_A.cc +vader/recipes/GeopotentialHeightAtInterface.h +vader/recipes/GeopotentialHeightAtSurface_A.cc +vader/recipes/GeopotentialHeightAtSurface.h +vader/recipes/HeightAboveMeanSeaLevelAtSurface_A.cc +vader/recipes/HeightAboveMeanSeaLevelAtSurface.h vader/recipes/HydrostaticExnerLevels.h vader/recipes/HydrostaticExnerLevels_A.cc vader/recipes/HydrostaticPressureLevels.h vader/recipes/HydrostaticPressureLevels_A.cc +vader/recipes/LnAirPressure_A.cc +vader/recipes/LnAirPressure.h vader/recipes/LnAirPressureAtInterface_A.cc vader/recipes/LnAirPressureAtInterface.h -vader/recipes/LogDerivativeSaturationVaporPressure.h vader/recipes/LogDerivativeSaturationVaporPressure_A.cc +vader/recipes/LogDerivativeSaturationVaporPressure.h vader/recipes/ParticulateMatter2p5.h vader/recipes/ParticulateMatter2p5_A.cc vader/recipes/ParticulateMatter2p5_B.cc @@ -61,10 +70,13 @@ vader/recipes/RainMixingRatio.h vader/recipes/RainMixingRatio_A.cc vader/recipes/RelativeHumidity.h vader/recipes/RelativeHumidity_A.cc +vader/recipes/RelativeHumidity_B.cc vader/recipes/SaturationSpecificHumidity.h vader/recipes/SaturationSpecificHumidity_A.cc +vader/recipes/SaturationSpecificHumidity_B.cc vader/recipes/SaturationVaporPressure.h vader/recipes/SaturationVaporPressure_A.cc +vader/recipes/SaturationVaporPressure_B.cc vader/recipes/SulfateMassFraction.h vader/recipes/SulfateMassFraction_A.cc vader/recipes/SurfaceAirPressure.h @@ -81,6 +93,30 @@ vader/recipes/TotalWaterMixingRatioWrtWetAir.h vader/recipes/TotalWaterMixingRatioWrtWetAir_A.cc vader/recipes/EastwardWindAt10m_A.cc vader/recipes/EastwardWindAt10m.h +vader/recipes/EffectiveRadiusOfCloudIceParticle.h +vader/recipes/EffectiveRadiusOfCloudIceParticle_A.cc +vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle.h +vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle_A.cc +vader/recipes/EffectiveRadiusOfGraupelParticle.h +vader/recipes/EffectiveRadiusOfGraupelParticle_A.cc +vader/recipes/EffectiveRadiusOfHailParticle.h +vader/recipes/EffectiveRadiusOfHailParticle_A.cc +vader/recipes/EffectiveRadiusOfRainParticle.h +vader/recipes/EffectiveRadiusOfRainParticle_A.cc +vader/recipes/EffectiveRadiusOfSnowParticle.h +vader/recipes/EffectiveRadiusOfSnowParticle_A.cc +vader/recipes/MassContentOfCloudIceInAtmosphereLayer.h +vader/recipes/MassContentOfCloudIceInAtmosphereLayer_A.cc +vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer.h +vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer_A.cc +vader/recipes/MassContentOfGraupelInAtmosphereLayer.h +vader/recipes/MassContentOfGraupelInAtmosphereLayer_A.cc +vader/recipes/MassContentOfHailInAtmosphereLayer.h +vader/recipes/MassContentOfHailInAtmosphereLayer_A.cc +vader/recipes/MassContentOfRainInAtmosphereLayer.h +vader/recipes/MassContentOfRainInAtmosphereLayer_A.cc +vader/recipes/MassContentOfSnowInAtmosphereLayer.h +vader/recipes/MassContentOfSnowInAtmosphereLayer_A.cc vader/recipes/VirtualPotentialTemperature_A.cc vader/recipes/VirtualPotentialTemperature_B.cc vader/recipes/VirtualPotentialTemperature.h @@ -96,8 +132,11 @@ vader/recipes/WaterVaporMixingRatioWrtWetAir.h vader/recipes/WaterVaporMixingRatioWrtWetAir_A.cc vader/recipes/WaterVaporMixingRatioWrtWetAir2m.h vader/recipes/WaterVaporMixingRatioWrtWetAir2m_A.cc -mo/common_varchange.h +vader/vader.cc +vader/vader.h +vader/VaderParameters.h mo/common_varchange.cc +mo/common_varchange.h mo/constants.h mo/functions.h mo/functions.cc @@ -159,12 +198,12 @@ mo/svpW_lookup.h if( gsw_FOUND ) list( APPEND vader_src_files - vader/recipes/SeaWaterPotentialTemperature.h - vader/recipes/SeaWaterPotentialTemperature_A.cc - vader/recipes/SeaWaterTemperature.h - vader/recipes/SeaWaterTemperature_A.cc OceanConversions/OceanConversions.interface.F90 OceanConversions/OceanConversions.interface.h + vader/recipes/SeaWaterPotentialTemperature_A.cc + vader/recipes/SeaWaterPotentialTemperature.h + vader/recipes/SeaWaterTemperature_A.cc + vader/recipes/SeaWaterTemperature.h ) endif( gsw_FOUND ) diff --git a/src/vader/DefaultCookbook.h b/src/vader/DefaultCookbook.h index 18765af..8647f35 100644 --- a/src/vader/DefaultCookbook.h +++ b/src/vader/DefaultCookbook.h @@ -21,10 +21,26 @@ #include "recipes/DryAirDensity.h" #include "recipes/DryAirDensityLevelsMinusOne.h" #include "recipes/EastwardWindAt10m.h" +#include "recipes/EffectiveRadiusOfCloudIceParticle.h" +#include "recipes/EffectiveRadiusOfCloudLiquidWaterParticle.h" +#include "recipes/EffectiveRadiusOfGraupelParticle.h" +#include "recipes/EffectiveRadiusOfHailParticle.h" +#include "recipes/EffectiveRadiusOfRainParticle.h" +#include "recipes/EffectiveRadiusOfSnowParticle.h" #include "recipes/HydrostaticExnerLevels.h" +#include "recipes/MassContentOfCloudIceInAtmosphereLayer.h" +#include "recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer.h" +#include "recipes/MassContentOfGraupelInAtmosphereLayer.h" +#include "recipes/MassContentOfHailInAtmosphereLayer.h" +#include "recipes/MassContentOfRainInAtmosphereLayer.h" +#include "recipes/MassContentOfSnowInAtmosphereLayer.h" #include "recipes/NorthwardWindAt10m.h" #include "recipes/ParticulateMatter2p5.h" #include "recipes/RainMixingRatio.h" +#include "recipes/RelativeHumidity.h" +#include "recipes/SaturationSpecificHumidity.h" +#include "recipes/SaturationVaporPressure.h" +#include "recipes/TotalRelativeHumidity.h" #include "recipes/TotalWater.h" #include "recipes/TotalWaterMixingRatioWrtDryAir.h" #include "recipes/VirtualPotentialTemperature.h" @@ -80,7 +96,37 @@ const cookbookConfigType Vader::defaultCookbookDefinition = { {vwind_at_10m_A::Name}}, {oops::Variable{"mass_density_of_particulate_matter_2p5_in_air"}, {ParticulateMatter2p5_A::Name, - ParticulateMatter2p5_B::Name}} + ParticulateMatter2p5_B::Name}}, + {oops::Variable{"mass_content_of_cloud_liquid_water_in_atmosphere_layer"}, + {MassContentOfCloudLiquidWaterInAtmosphereLayer_A::Name}}, + {oops::Variable{"mass_content_of_cloud_ice_in_atmosphere_layer"}, + {MassContentOfCloudIceInAtmosphereLayer_A::Name}}, + {oops::Variable{"mass_content_of_rain_in_atmosphere_layer"}, + {MassContentOfRainInAtmosphereLayer_A::Name}}, + {oops::Variable{"mass_content_of_snow_in_atmosphere_layer"}, + {MassContentOfSnowInAtmosphereLayer_A::Name}}, + {oops::Variable{"mass_content_of_graupel_in_atmosphere_layer"}, + {MassContentOfGraupelInAtmosphereLayer_A::Name}}, + {oops::Variable{"mass_content_of_hail_in_atmosphere_layer"}, + {MassContentOfHailInAtmosphereLayer_A::Name}}, + {oops::Variable{"effective_radius_of_cloud_liquid_water_particle"}, + {EffectiveRadiusOfCloudLiquidWaterParticle_A::Name}}, + {oops::Variable{"effective_radius_of_cloud_ice_particle"}, + {EffectiveRadiusOfCloudIceParticle_A::Name}}, + {oops::Variable{"effective_radius_of_rain_particle"}, + {EffectiveRadiusOfRainParticle_A::Name}}, + {oops::Variable{"effective_radius_of_snow_particle"}, + {EffectiveRadiusOfSnowParticle_A::Name}}, + {oops::Variable{"effective_radius_of_graupel_particle"}, + {EffectiveRadiusOfGraupelParticle_A::Name}}, + {oops::Variable{"effective_radius_of_hail_particle"}, + {EffectiveRadiusOfHailParticle_A::Name}}, + {oops::Variable{"svp"}, + {SaturationVaporPressure_B::Name}}, + {oops::Variable{"water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation"}, + {SaturationSpecificHumidity_B::Name}}, + {oops::Variable{"relative_humidity"}, + {RelativeHumidity_A::Name, RelativeHumidity_B::Name}} }; } // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfCloudIceParticle.h b/src/vader/recipes/EffectiveRadiusOfCloudIceParticle.h new file mode 100644 index 0000000..e356ba3 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfCloudIceParticle.h @@ -0,0 +1,65 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class EffectiveRadiusOfCloudIceParticle_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(EffectiveRadiusOfCloudIceParticle_AParameters, + RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'EffectiveRadiusOfCloudIceParticle_A' defines a recipe for + * effective radius of ice particle + * + * \details This instantiation of RecipeBase produces effective radius + * using a simplified diagnostic approximation with temperature-dependent + * relationship, based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature + * Output units: m + */ +class EffectiveRadiusOfCloudIceParticle_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef EffectiveRadiusOfCloudIceParticle_AParameters Parameters_; + + EffectiveRadiusOfCloudIceParticle_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfCloudIceParticle_A.cc b/src/vader/recipes/EffectiveRadiusOfCloudIceParticle_A.cc new file mode 100644 index 0000000..c0a0dfe --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfCloudIceParticle_A.cc @@ -0,0 +1,104 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/EffectiveRadiusOfCloudIceParticle.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char EffectiveRadiusOfCloudIceParticle_A::Name[] = + "EffectiveRadiusOfCloudIceParticle_A"; +const oops::Variables EffectiveRadiusOfCloudIceParticle_A::Ingredients{ + std::vector{"air_temperature"}}; + +// Register the maker +static RecipeMaker + makerEffectiveRadiusOfCloudIceParticle_A_( + EffectiveRadiusOfCloudIceParticle_A::Name); + +// ------------------------------------------------------------------------------------------------- + +EffectiveRadiusOfCloudIceParticle_A::EffectiveRadiusOfCloudIceParticle_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "EffectiveRadiusOfCloudIceParticle_A::" + << "EffectiveRadiusOfCloudIceParticle_A" << std::endl; +} + +std::string EffectiveRadiusOfCloudIceParticle_A::name() const { + return EffectiveRadiusOfCloudIceParticle_A::Name; +} + +oops::Variable EffectiveRadiusOfCloudIceParticle_A::product() const { + return oops::Variable{"effective_radius_of_cloud_ice_particle"}; +} + +oops::Variables EffectiveRadiusOfCloudIceParticle_A::ingredients() const { + return EffectiveRadiusOfCloudIceParticle_A::Ingredients; +} + +size_t EffectiveRadiusOfCloudIceParticle_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace EffectiveRadiusOfCloudIceParticle_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void EffectiveRadiusOfCloudIceParticle_A::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "EffectiveRadiusOfCloudIceParticle_A::executeNL starting" + << std::endl; + + // Get input field + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + + // Create output field + auto reff_field = afieldset.field("effective_radius_of_cloud_ice_particle"); + auto reff = atlas::array::make_view(reff_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Temperature-dependent effective radius for ice particles + // Typical range: 20-100 micrometers + // Colder temperatures -> smaller crystals + const double T_freeze = 273.15; // K + const double T_cold = 233.15; // -40C + const double r_min = 20.0e-6; // minimum radius (m) + const double r_max = 100.0e-6; // maximum radius (m) + + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + + if (T < T_freeze) { + // Warmer ice clouds have larger particles due to aggregation + const double temp_factor = (T - T_cold) / (T_freeze - T_cold); + reff(jn, jl) = r_min + (r_max - r_min) * std::max(0.0, std::min(1.0, temp_factor)); + } else { + // Above freezing, use minimum value + reff(jn, jl) = r_min; + } + } + } + + oops::Log::trace() << "EffectiveRadiusOfCloudIceParticle_A::executeNL done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle.h b/src/vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle.h new file mode 100644 index 0000000..0f8e1f2 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle.h @@ -0,0 +1,65 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class EffectiveRadiusOfCloudLiquidWaterParticle_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(EffectiveRadiusOfCloudLiquidWaterParticle_AParameters, + RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'EffectiveRadiusOfCloudLiquidWaterParticle_A' defines a recipe for + * effective radius of cloud liquid water particle + * + * \details This instantiation of RecipeBase produces effective radius + * using a simplified diagnostic approximation with temperature-dependent + * relationship, based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature + * Output units: m + */ +class EffectiveRadiusOfCloudLiquidWaterParticle_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef EffectiveRadiusOfCloudLiquidWaterParticle_AParameters Parameters_; + + EffectiveRadiusOfCloudLiquidWaterParticle_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle_A.cc b/src/vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle_A.cc new file mode 100644 index 0000000..462e7b6 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle_A.cc @@ -0,0 +1,108 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/EffectiveRadiusOfCloudLiquidWaterParticle.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char EffectiveRadiusOfCloudLiquidWaterParticle_A::Name[] = + "EffectiveRadiusOfCloudLiquidWaterParticle_A"; +const oops::Variables EffectiveRadiusOfCloudLiquidWaterParticle_A::Ingredients{ + std::vector{"air_temperature"}}; + +// Register the maker +static RecipeMaker + makerEffectiveRadiusOfCloudLiquidWaterParticle_A_( + EffectiveRadiusOfCloudLiquidWaterParticle_A::Name); + +// ------------------------------------------------------------------------------------------------- + +EffectiveRadiusOfCloudLiquidWaterParticle_A::EffectiveRadiusOfCloudLiquidWaterParticle_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "EffectiveRadiusOfCloudLiquidWaterParticle_A::" + << "EffectiveRadiusOfCloudLiquidWaterParticle_A" << std::endl; +} + +std::string EffectiveRadiusOfCloudLiquidWaterParticle_A::name() const { + return EffectiveRadiusOfCloudLiquidWaterParticle_A::Name; +} + +oops::Variable EffectiveRadiusOfCloudLiquidWaterParticle_A::product() const { + return oops::Variable{"effective_radius_of_cloud_liquid_water_particle"}; +} + +oops::Variables EffectiveRadiusOfCloudLiquidWaterParticle_A::ingredients() const { + return EffectiveRadiusOfCloudLiquidWaterParticle_A::Ingredients; +} + +size_t EffectiveRadiusOfCloudLiquidWaterParticle_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace EffectiveRadiusOfCloudLiquidWaterParticle_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void EffectiveRadiusOfCloudLiquidWaterParticle_A::executeNL( + atlas::FieldSet & afieldset) { + oops::Log::trace() << "EffectiveRadiusOfCloudLiquidWaterParticle_A::" + << "executeNL starting" << std::endl; + + // Get input field + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + + // Create output field + auto reff_field = afieldset.field( + "effective_radius_of_cloud_liquid_water_particle"); + auto reff = atlas::array::make_view(reff_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Temperature-dependent effective radius + // Following Martin et al. (1994) + // Typical range: 4-15 micrometers for liquid water clouds + const double T_freeze = 273.15; // K + const double r_min = 4.0e-6; // minimum radius (m) + const double r_max = 15.0e-6; // maximum radius (m) + + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + + // Linear temperature dependence + // Warmer clouds have larger droplets + if (T > T_freeze) { + const double T_range = 30.0; // Temperature range for scaling + const double temp_factor = std::min(1.0, (T - T_freeze) / T_range); + reff(jn, jl) = r_min + (r_max - r_min) * temp_factor; + } else { + // Below freezing, set to minimum + reff(jn, jl) = r_min; + } + } + } + + oops::Log::trace() << "EffectiveRadiusOfCloudLiquidWaterParticle_A::" + << "executeNL done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfGraupelParticle.h b/src/vader/recipes/EffectiveRadiusOfGraupelParticle.h new file mode 100644 index 0000000..439a439 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfGraupelParticle.h @@ -0,0 +1,63 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class EffectiveRadiusOfGraupelParticle_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(EffectiveRadiusOfGraupelParticle_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'EffectiveRadiusOfGraupelParticle_A' defines a recipe for + * effective radius of graupel particle + * + * \details This instantiation of RecipeBase produces effective radius + * using a simplified diagnostic approximation with fixed value, + * based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Output units: m + */ +class EffectiveRadiusOfGraupelParticle_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef EffectiveRadiusOfGraupelParticle_AParameters Parameters_; + + EffectiveRadiusOfGraupelParticle_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfGraupelParticle_A.cc b/src/vader/recipes/EffectiveRadiusOfGraupelParticle_A.cc new file mode 100644 index 0000000..0e0f6d3 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfGraupelParticle_A.cc @@ -0,0 +1,93 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/EffectiveRadiusOfGraupelParticle.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char EffectiveRadiusOfGraupelParticle_A::Name[] = + "EffectiveRadiusOfGraupelParticle_A"; +const oops::Variables EffectiveRadiusOfGraupelParticle_A::Ingredients{ + std::vector{"air_temperature"}}; + +// Register the maker +static RecipeMaker + makerEffectiveRadiusOfGraupelParticle_A_( + EffectiveRadiusOfGraupelParticle_A::Name); + +// ------------------------------------------------------------------------------------------------- + +EffectiveRadiusOfGraupelParticle_A::EffectiveRadiusOfGraupelParticle_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "EffectiveRadiusOfGraupelParticle_A::" + << "EffectiveRadiusOfGraupelParticle_A" << std::endl; +} + +std::string EffectiveRadiusOfGraupelParticle_A::name() const { + return EffectiveRadiusOfGraupelParticle_A::Name; +} + +oops::Variable EffectiveRadiusOfGraupelParticle_A::product() const { + return oops::Variable{"effective_radius_of_graupel_particle"}; +} + +oops::Variables EffectiveRadiusOfGraupelParticle_A::ingredients() const { + return EffectiveRadiusOfGraupelParticle_A::Ingredients; +} + +size_t EffectiveRadiusOfGraupelParticle_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace EffectiveRadiusOfGraupelParticle_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void EffectiveRadiusOfGraupelParticle_A::executeNL( + atlas::FieldSet & afieldset) { + oops::Log::trace() << "EffectiveRadiusOfGraupelParticle_A::" + << "executeNL starting" << std::endl; + + // Get input field for shape only + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + + // Create output field + auto reff_field = afieldset.field( + "effective_radius_of_graupel_particle"); + auto reff = atlas::array::make_view(reff_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Fixed effective radius for graupel + // Typical value: ~500-1500 micrometers for graupel particles + const double r_graupel = 750.0e-6; // 750 micrometers (m) + + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + reff(jn, jl) = r_graupel; + } + } + + oops::Log::trace() << "EffectiveRadiusOfGraupelParticle_A::" + << "executeNL done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfHailParticle.h b/src/vader/recipes/EffectiveRadiusOfHailParticle.h new file mode 100644 index 0000000..96e4d0b --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfHailParticle.h @@ -0,0 +1,63 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class EffectiveRadiusOfHailParticle_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(EffectiveRadiusOfHailParticle_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'EffectiveRadiusOfHailParticle_A' defines a recipe for + * effective radius of hail particle + * + * \details This instantiation of RecipeBase produces effective radius + * using a simplified diagnostic approximation with fixed value, + * based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Output units: m + */ +class EffectiveRadiusOfHailParticle_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef EffectiveRadiusOfHailParticle_AParameters Parameters_; + + EffectiveRadiusOfHailParticle_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfHailParticle_A.cc b/src/vader/recipes/EffectiveRadiusOfHailParticle_A.cc new file mode 100644 index 0000000..3a45b94 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfHailParticle_A.cc @@ -0,0 +1,93 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/EffectiveRadiusOfHailParticle.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char EffectiveRadiusOfHailParticle_A::Name[] = + "EffectiveRadiusOfHailParticle_A"; +const oops::Variables EffectiveRadiusOfHailParticle_A::Ingredients{ + std::vector{"air_temperature"}}; + +// Register the maker +static RecipeMaker + makerEffectiveRadiusOfHailParticle_A_( + EffectiveRadiusOfHailParticle_A::Name); + +// ------------------------------------------------------------------------------------------------- + +EffectiveRadiusOfHailParticle_A::EffectiveRadiusOfHailParticle_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "EffectiveRadiusOfHailParticle_A::" + << "EffectiveRadiusOfHailParticle_A" << std::endl; +} + +std::string EffectiveRadiusOfHailParticle_A::name() const { + return EffectiveRadiusOfHailParticle_A::Name; +} + +oops::Variable EffectiveRadiusOfHailParticle_A::product() const { + return oops::Variable{"effective_radius_of_hail_particle"}; +} + +oops::Variables EffectiveRadiusOfHailParticle_A::ingredients() const { + return EffectiveRadiusOfHailParticle_A::Ingredients; +} + +size_t EffectiveRadiusOfHailParticle_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace EffectiveRadiusOfHailParticle_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void EffectiveRadiusOfHailParticle_A::executeNL( + atlas::FieldSet & afieldset) { + oops::Log::trace() << "EffectiveRadiusOfHailParticle_A::" + << "executeNL starting" << std::endl; + + // Get input field for shape only + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + + // Create output field + auto reff_field = afieldset.field( + "effective_radius_of_hail_particle"); + auto reff = atlas::array::make_view(reff_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Fixed effective radius for hail + // Typical value: ~2-10 mm for hail stones + const double r_hail = 5.0e-3; // 5 mm (m) + + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + reff(jn, jl) = r_hail; + } + } + + oops::Log::trace() << "EffectiveRadiusOfHailParticle_A::" + << "executeNL done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfRainParticle.h b/src/vader/recipes/EffectiveRadiusOfRainParticle.h new file mode 100644 index 0000000..4045490 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfRainParticle.h @@ -0,0 +1,63 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class EffectiveRadiusOfRainParticle_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(EffectiveRadiusOfRainParticle_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'EffectiveRadiusOfRainParticle_A' defines a recipe for + * effective radius of rain particle + * + * \details This instantiation of RecipeBase produces effective radius + * using a simplified diagnostic approximation with fixed value, + * based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Output units: m + */ +class EffectiveRadiusOfRainParticle_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef EffectiveRadiusOfRainParticle_AParameters Parameters_; + + EffectiveRadiusOfRainParticle_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfRainParticle_A.cc b/src/vader/recipes/EffectiveRadiusOfRainParticle_A.cc new file mode 100644 index 0000000..82b8b92 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfRainParticle_A.cc @@ -0,0 +1,93 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/EffectiveRadiusOfRainParticle.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char EffectiveRadiusOfRainParticle_A::Name[] = + "EffectiveRadiusOfRainParticle_A"; +const oops::Variables EffectiveRadiusOfRainParticle_A::Ingredients{ + std::vector{"air_temperature"}}; + +// Register the maker +static RecipeMaker + makerEffectiveRadiusOfRainParticle_A_( + EffectiveRadiusOfRainParticle_A::Name); + +// ------------------------------------------------------------------------------------------------- + +EffectiveRadiusOfRainParticle_A::EffectiveRadiusOfRainParticle_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "EffectiveRadiusOfRainParticle_A::" + << "EffectiveRadiusOfRainParticle_A" << std::endl; +} + +std::string EffectiveRadiusOfRainParticle_A::name() const { + return EffectiveRadiusOfRainParticle_A::Name; +} + +oops::Variable EffectiveRadiusOfRainParticle_A::product() const { + return oops::Variable{"effective_radius_of_rain_particle"}; +} + +oops::Variables EffectiveRadiusOfRainParticle_A::ingredients() const { + return EffectiveRadiusOfRainParticle_A::Ingredients; +} + +size_t EffectiveRadiusOfRainParticle_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace EffectiveRadiusOfRainParticle_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void EffectiveRadiusOfRainParticle_A::executeNL( + atlas::FieldSet & afieldset) { + oops::Log::trace() << "EffectiveRadiusOfRainParticle_A::" + << "executeNL starting" << std::endl; + + // Get input field for shape only + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + + // Create output field + auto reff_field = afieldset.field( + "effective_radius_of_rain_particle"); + auto reff = atlas::array::make_view(reff_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Fixed effective radius for rain droplets + // Typical value: ~200-500 micrometers for rain + const double r_rain = 250.0e-6; // 250 micrometers (m) + + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + reff(jn, jl) = r_rain; + } + } + + oops::Log::trace() << "EffectiveRadiusOfRainParticle_A::" + << "executeNL done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfSnowParticle.h b/src/vader/recipes/EffectiveRadiusOfSnowParticle.h new file mode 100644 index 0000000..c794154 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfSnowParticle.h @@ -0,0 +1,64 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class EffectiveRadiusOfSnowParticle_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(EffectiveRadiusOfSnowParticle_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'EffectiveRadiusOfSnowParticle_A' defines a recipe for + * effective radius of snow particle + * + * \details This instantiation of RecipeBase produces effective radius + * using a simplified diagnostic approximation with temperature-dependent + * relationship, based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature + * Output units: m + */ +class EffectiveRadiusOfSnowParticle_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef EffectiveRadiusOfSnowParticle_AParameters Parameters_; + + EffectiveRadiusOfSnowParticle_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/EffectiveRadiusOfSnowParticle_A.cc b/src/vader/recipes/EffectiveRadiusOfSnowParticle_A.cc new file mode 100644 index 0000000..7d98a69 --- /dev/null +++ b/src/vader/recipes/EffectiveRadiusOfSnowParticle_A.cc @@ -0,0 +1,107 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/EffectiveRadiusOfSnowParticle.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char EffectiveRadiusOfSnowParticle_A::Name[] = + "EffectiveRadiusOfSnowParticle_A"; +const oops::Variables EffectiveRadiusOfSnowParticle_A::Ingredients{ + std::vector{"air_temperature"}}; + +// Register the maker +static RecipeMaker + makerEffectiveRadiusOfSnowParticle_A_( + EffectiveRadiusOfSnowParticle_A::Name); + +// ------------------------------------------------------------------------------------------------- + +EffectiveRadiusOfSnowParticle_A::EffectiveRadiusOfSnowParticle_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "EffectiveRadiusOfSnowParticle_A::" + << "EffectiveRadiusOfSnowParticle_A" << std::endl; +} + +std::string EffectiveRadiusOfSnowParticle_A::name() const { + return EffectiveRadiusOfSnowParticle_A::Name; +} + +oops::Variable EffectiveRadiusOfSnowParticle_A::product() const { + return oops::Variable{"effective_radius_of_snow_particle"}; +} + +oops::Variables EffectiveRadiusOfSnowParticle_A::ingredients() const { + return EffectiveRadiusOfSnowParticle_A::Ingredients; +} + +size_t EffectiveRadiusOfSnowParticle_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace EffectiveRadiusOfSnowParticle_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void EffectiveRadiusOfSnowParticle_A::executeNL( + atlas::FieldSet & afieldset) { + oops::Log::trace() << "EffectiveRadiusOfSnowParticle_A::" + << "executeNL starting" << std::endl; + + // Get input field + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + + // Create output field + auto reff_field = afieldset.field( + "effective_radius_of_snow_particle"); + auto reff = atlas::array::make_view(reff_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Temperature-dependent effective radius for snow aggregates + // Typical range: 100-1000 micrometers + // Warmer snow (near 0C) has larger aggregates + const double T_freeze = 273.15; // K + const double T_cold = 233.15; // -40C + const double r_min = 100.0e-6; // minimum radius (m) + const double r_max = 1000.0e-6; // maximum radius (m) + + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + + if (T < T_freeze) { + // Warmer temperatures allow more aggregation -> larger particles + const double temp_factor = (T - T_cold) / (T_freeze - T_cold); + reff(jn, jl) = r_min + (r_max - r_min) * std::max(0.0, std::min(1.0, temp_factor)); + } else { + // Above freezing, use minimum + reff(jn, jl) = r_min; + } + } + } + + oops::Log::trace() << "EffectiveRadiusOfSnowParticle_A::" + << "executeNL done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/GeopotentialAtInterface.h b/src/vader/recipes/GeopotentialAtInterface.h new file mode 100644 index 0000000..0b63a2d --- /dev/null +++ b/src/vader/recipes/GeopotentialAtInterface.h @@ -0,0 +1,64 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader { + +class GeopotentialAtInterface_A_Parameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(GeopotentialAtInterface_A_Parameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{ + "recipe name", + this}; +}; + +// ------------------------------------------------------------------------------------------------- + +/*! \brief GeopotentialAtInterface_A computes geopotential at model interfaces + * + * Uses hydrostatic integration with layer-averaged virtual temperature. + * Inputs: geopotential at full levels, virtual_temperature, ln_air_pressure at full levels + * and interfaces. Output: geopotential_at_interface (m^2 s^-2). + * Full TL/AD support. + */ +class GeopotentialAtInterface_A : public RecipeBase { + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef GeopotentialAtInterface_A_Parameters Parameters_; + + GeopotentialAtInterface_A(const Parameters_ &, const VaderConfigVars &); + + // Recipe base class overrides + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +} // namespace vader diff --git a/src/vader/recipes/GeopotentialAtInterface_A.cc b/src/vader/recipes/GeopotentialAtInterface_A.cc new file mode 100644 index 0000000..90d69ce --- /dev/null +++ b/src/vader/recipes/GeopotentialAtInterface_A.cc @@ -0,0 +1,247 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "atlas/util/Metadata.h" +#include "oops/util/abor1_cpp.h" +#include "oops/util/Logger.h" +#include "vader/recipes/GeopotentialAtInterface.h" + +namespace vader +{ +// ------------------------------------------------------------------------------------------------ + +// Static attribute initialization +const char GeopotentialAtInterface_A::Name[] = "GeopotentialAtInterface_A"; +const oops::Variables GeopotentialAtInterface_A::Ingredients{ + std::vector{ + "geopotential", + "virtual_temperature", + "ln_air_pressure", + "ln_air_pressure_at_interface"}}; + +// Register the maker +static RecipeMaker makerGeopotentialAtInterface_A_( + GeopotentialAtInterface_A::Name); + +GeopotentialAtInterface_A::GeopotentialAtInterface_A(const Parameters_ & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} +{ + oops::Log::trace() << "GeopotentialAtInterface_A::GeopotentialAtInterface_A(params)" + << std::endl; +} + +std::string GeopotentialAtInterface_A::name() const +{ + return GeopotentialAtInterface_A::Name; +} + +oops::Variable GeopotentialAtInterface_A::product() const +{ + return oops::Variable{"geopotential_levels"}; +} + +oops::Variables GeopotentialAtInterface_A::ingredients() const +{ + return GeopotentialAtInterface_A::Ingredients; +} + +size_t GeopotentialAtInterface_A::productLevels(const atlas::FieldSet & afieldset) const +{ + return afieldset.field("ln_air_pressure_at_interface").shape(1); +} +atlas::FunctionSpace GeopotentialAtInterface_A::productFunctionSpace(const atlas::FieldSet & + afieldset) const +{ + return afieldset.field("geopotential").functionspace(); +} + +// ------------------------------------------------------------------------------------------------- + +void GeopotentialAtInterface_A::executeNL(atlas::FieldSet & afieldset) +{ + oops::Log::trace() << "entering GeopotentialAtInterface_A::executeNL function" << std::endl; + + // Extract values from client config + const double Rd = configVariables_.getDouble("gas_constant_of_dry_air"); + + atlas::Field phi = afieldset.field("geopotential"); + atlas::Field tv = afieldset.field("virtual_temperature"); + atlas::Field ln_p = afieldset.field("ln_air_pressure"); + atlas::Field ln_p_int = afieldset.field("ln_air_pressure_at_interface"); + atlas::Field phi_int = afieldset.field("geopotential_levels"); + + auto phi_view = atlas::array::make_view(phi); + auto tv_view = atlas::array::make_view(tv); + auto ln_p_view = atlas::array::make_view(ln_p); + auto ln_p_int_view = atlas::array::make_view(ln_p_int); + auto phi_int_view = atlas::array::make_view(phi_int); + + const size_t npoint = phi.shape(0); + const int nlev = phi.shape(1); + const int nint = ln_p_int.shape(1); + + if (nlev < 2 || nint < 2) { + oops::Log::error() << "GeopotentialAtInterface_A::executeNL: need at least 2 full levels " + "and 2 interfaces" << std::endl; + ABORT("GeopotentialAtInterface_A::executeNL: need at least 2 full levels and 2 interfaces"); + } + + std::vector tv_layer(nlev - 1); + + for (size_t ip = 0; ip < npoint; ++ip) { + // Layer-mean Tv between full levels (top->bottom) + for (int k = 0; k < nlev - 1; ++k) { + tv_layer[k] = 0.5 * (tv_view(ip, k) + tv_view(ip, k + 1)); + } + + // Top interface: extrapolate half-layer using precomputed ln pressures + double dln_top = ln_p_view(ip, 0) - ln_p_int_view(ip, 0); + phi_int_view(ip, 0) = phi_view(ip, 0) + Rd * tv_layer[0] * dln_top; + + // Integrate downward using precomputed ln(p) at interfaces + for (int k = 0; k < nint - 1; ++k) { + double dln = ln_p_int_view(ip, k + 1) - ln_p_int_view(ip, k); + phi_int_view(ip, k + 1) = phi_int_view(ip, k) - Rd * tv_layer[k] * dln; + } + } + + oops::Log::trace() << "leaving GeopotentialAtInterface_A::executeNL function" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +void GeopotentialAtInterface_A::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & /*afieldsetTraj*/) { + oops::Log::trace() << "entering GeopotentialAtInterface_A::executeTL function" << std::endl; + + // Extract values from client config + const double Rd = configVariables_.getDouble("gas_constant_of_dry_air"); + + atlas::Field phi_tl = afieldsetTL.field("geopotential"); + atlas::Field tv_tl = afieldsetTL.field("virtual_temperature"); + atlas::Field ln_p = afieldsetTL.field("ln_air_pressure"); + atlas::Field ln_p_int = afieldsetTL.field("ln_air_pressure_at_interface"); + atlas::Field phi_int_tl = afieldsetTL.field("geopotential_levels"); + + auto phi_tl_view = atlas::array::make_view(phi_tl); + auto tv_tl_view = atlas::array::make_view(tv_tl); + auto ln_p_view = atlas::array::make_view(ln_p); + auto ln_p_int_view = atlas::array::make_view(ln_p_int); + auto phi_int_tl_view = atlas::array::make_view(phi_int_tl); + + const size_t npoint = phi_tl.shape(0); + const int nlev = phi_tl.shape(1); + const int nint = ln_p_int_view.shape(1); + + if (nlev < 2 || nint < 2) { + oops::Log::error() << "GeopotentialAtInterface_A::executeNL: need at least 2 full levels " + "and 2 interfaces" << std::endl; + ABORT("GeopotentialAtInterface_A::executeNL: need at least 2 full levels and 2 interfaces"); + } + + std::vector tv_layer_tl(nlev - 1); + + for (size_t ip = 0; ip < npoint; ++ip) { + // Layer-mean Tv' for TL + for (int k = 0; k < nlev - 1; ++k) { + tv_layer_tl[k] = 0.5 * (tv_tl_view(ip, k) + tv_tl_view(ip, k + 1)); + } + + // Top interface TL + double dln_top = ln_p_view(ip, 0) - ln_p_int_view(ip, 0); + phi_int_tl_view(ip, 0) = phi_tl_view(ip, 0) + Rd * tv_layer_tl[0] * dln_top; + + // Propagate TL downward + for (int k = 0; k < nint - 1; ++k) { + double dln = ln_p_int_view(ip, k + 1) - ln_p_int_view(ip, k); + phi_int_tl_view(ip, k + 1) = + phi_int_tl_view(ip, k) - Rd * tv_layer_tl[k] * dln; + } + } + + oops::Log::trace() << "leaving GeopotentialAtInterface_A::executeTL function" << std::endl; +} + +// ------------------------------------------------------------------------------------------------ + +void GeopotentialAtInterface_A::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & /*afieldsetTraj*/) { + oops::Log::trace() << "entering GeopotentialAtInterface_A::executeAD function" << std::endl; + + // Extract values from client config + const double Rd = configVariables_.getDouble("gas_constant_of_dry_air"); + + atlas::Field phi_ad = afieldsetAD.field("geopotential"); + atlas::Field tv_ad = afieldsetAD.field("virtual_temperature"); + atlas::Field ln_p = afieldsetAD.field("ln_air_pressure"); + atlas::Field ln_p_int = afieldsetAD.field("ln_air_pressure_at_interface"); + atlas::Field phi_int_ad = afieldsetAD.field("geopotential_levels"); + + auto phi_ad_view = atlas::array::make_view(phi_ad); + auto tv_ad_view = atlas::array::make_view(tv_ad); + auto ln_p_view = atlas::array::make_view(ln_p); + auto ln_p_int_view = atlas::array::make_view(ln_p_int); + auto phi_int_ad_view = atlas::array::make_view(phi_int_ad); + + const size_t npoint = phi_ad.shape(0); + const int nlev = phi_ad.shape(1); + const int nint = ln_p_int_view.shape(1); + + if (nlev < 2 || nint < 2) { + oops::Log::error() << "GeopotentialAtInterface_A::executeNL: need at least 2 full levels " + "and 2 interfaces" << std::endl; + ABORT("GeopotentialAtInterface_A::executeNL: need at least 2 full levels and 2 interfaces"); + } + + for (size_t ip = 0; ip < npoint; ++ip) { + // Propagate adjoint from bottom->top (reverse of TL forward order) + for (int k = nint - 1; k >= 1; --k) { + double lam = phi_int_ad_view(ip, k); + if (lam != 0.0) { + // Pass adjoint to previous interface + phi_int_ad_view(ip, k - 1) += lam; + + // Contribution to tv_ad at full levels k-1 and k (layer k-1) + double dln = ln_p_int_view(ip, k) - ln_p_int_view(ip, k - 1); + double coeff = -rdry * 0.5 * dln; + tv_ad_view(ip, k - 1) += coeff * lam; + tv_ad_view(ip, k) += coeff * lam; + + // Zero out phi_int_ad after processing + phi_int_ad_view(ip, k) = 0.0; + } + } + + // Handle top interface (index 0) + double lam0 = phi_int_ad_view(ip, 0); + if (lam0 != 0.0) { + // Contribution to phi_ad at top full level + phi_ad_view(ip, 0) += lam0; + + // Contribution to tv_ad from extrapolation term + double dln_top = ln_p_view(ip, 0) - ln_p_int_view(ip, 0); + double coeff0 = Rd * 0.5 * dln_top; + tv_ad_view(ip, 0) += coeff0 * lam0; + tv_ad_view(ip, 1) += coeff0 * lam0; + + // Zero out phi_int_ad after processing + phi_int_ad_view(ip, 0) = 0.0; + } + } + + oops::Log::trace() << "leaving GeopotentialAtInterface_A::executeAD function" << std::endl; +} + +// ------------------------------------------------------------------------------------------------ + +} // namespace vader diff --git a/src/vader/recipes/GeopotentialHeight.h b/src/vader/recipes/GeopotentialHeight.h new file mode 100644 index 0000000..8dd6f76 --- /dev/null +++ b/src/vader/recipes/GeopotentialHeight.h @@ -0,0 +1,64 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader { + +class GeopotentialHeight_A_Parameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(GeopotentialHeight_A_Parameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{ + "recipe name", + this}; +}; + +// ------------------------------------------------------------------------------------------------- + +/*! \brief GeopotentialHeight_A computes geopotential height from geopotential + * + * Formula: z = phi / g + * where phi is geopotential (m^2 s^-2), g is standard_gravitational_acceleration (m s^-2), + * and z is geopotential_height (m). + * Full TL/AD support. + */ +class GeopotentialHeight_A : public RecipeBase { + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef GeopotentialHeight_A_Parameters Parameters_; + + GeopotentialHeight_A(const Parameters_ &, const VaderConfigVars &); + + // Recipe base class overrides + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +} // namespace vader diff --git a/src/vader/recipes/GeopotentialHeightAtInterface.h b/src/vader/recipes/GeopotentialHeightAtInterface.h new file mode 100644 index 0000000..75b19c6 --- /dev/null +++ b/src/vader/recipes/GeopotentialHeightAtInterface.h @@ -0,0 +1,64 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader { + +class GeopotentialHeightAtInterface_A_Parameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(GeopotentialHeightAtInterface_A_Parameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{ + "recipe name", + this}; +}; + +// ------------------------------------------------------------------------------------------------- + +/*! \brief GeopotentialHeightAtInterface_A computes geopotential height at model interfaces + * + * Formula: z_int = phi_int / g + * where phi_int is geopotential_at_interface (m^2 s^-2), g is standard_gravitational_acceleration (m s^-2), + * and z_int is geopotential_height_levels (m). + * Full TL/AD support. + */ +class GeopotentialHeightAtInterface_A : public RecipeBase { + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef GeopotentialHeightAtInterface_A_Parameters Parameters_; + + GeopotentialHeightAtInterface_A(const Parameters_ &, const VaderConfigVars &); + + // Recipe base class overrides + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +} // namespace vader diff --git a/src/vader/recipes/GeopotentialHeightAtInterface_A.cc b/src/vader/recipes/GeopotentialHeightAtInterface_A.cc new file mode 100644 index 0000000..a7737c4 --- /dev/null +++ b/src/vader/recipes/GeopotentialHeightAtInterface_A.cc @@ -0,0 +1,150 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "atlas/util/Metadata.h" +#include "oops/util/Logger.h" +#include "vader/recipes/GeopotentialHeightAtInterface.h" + +namespace vader +{ +// ------------------------------------------------------------------------------------------------ + +// Static attribute initialization +const char GeopotentialHeightAtInterface_A::Name[] = "GeopotentialHeightAtInterface_A"; +const oops::Variables GeopotentialHeightAtInterface_A::Ingredients{std::vector{ + "geopotential_levels"}}; + +// Register the maker +static RecipeMaker makerGeopotentialHeightAtInterface_A_( + GeopotentialHeightAtInterface_A::Name); + +GeopotentialHeightAtInterface_A::GeopotentialHeightAtInterface_A(const Parameters_ & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} +{ + oops::Log::trace() << "GeopotentialHeightAtInterface_A::GeopotentialHeightAtInterface_A(params)" + << std::endl; +} + +std::string GeopotentialHeightAtInterface_A::name() const +{ + return GeopotentialHeightAtInterface_A::Name; +} + +oops::Variable GeopotentialHeightAtInterface_A::product() const +{ + return oops::Variable{"geopotential_height_levels"}; +} + +oops::Variables GeopotentialHeightAtInterface_A::ingredients() const +{ + return GeopotentialHeightAtInterface_A::Ingredients; +} + +size_t GeopotentialHeightAtInterface_A::productLevels(const atlas::FieldSet & afieldset) const +{ + return afieldset.field("geopotential_levels").shape(1); +} + +atlas::FunctionSpace GeopotentialHeightAtInterface_A::productFunctionSpace(const atlas::FieldSet & + afieldset) const +{ + return afieldset.field("geopotential_levels").functionspace(); +} + +void GeopotentialHeightAtInterface_A::executeNL(atlas::FieldSet & afieldset) +{ + oops::Log::trace() << "entering GeopotentialHeightAtInterface_A::executeNL function" + << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi_int = afieldset.field("geopotential_levels"); + atlas::Field z_int = afieldset.field("geopotential_height_levels"); + + auto phi_int_view = atlas::array::make_view(phi_int); + auto z_int_view = atlas::array::make_view(z_int); + + const size_t grid_size = phi_int.shape(0); + const int nlevels = phi_int.shape(1); + const double inv_g = 1.0 / grav; + + for (int level = 0; level < nlevels; ++level) { + for ( size_t jnode = 0; jnode < grid_size ; ++jnode ) { + z_int_view(jnode, level) = phi_int_view(jnode, level) * inv_g; + } + } + oops::Log::trace() << "leaving GeopotentialHeightAtInterface_A::executeNL function" << std::endl; +} + +void GeopotentialHeightAtInterface_A::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & /*afieldsetTraj*/) +{ + oops::Log::trace() << "entering GeopotentialHeightAtInterface_A::executeTL function" + << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi_int_tl = afieldsetTL.field("geopotential_levels"); + atlas::Field z_int_tl = afieldsetTL.field("geopotential_height_levels"); + + auto phi_int_tl_view = atlas::array::make_view(phi_int_tl); + auto z_int_tl_view = atlas::array::make_view(z_int_tl); + + const size_t grid_size = phi_int_tl.shape(0); + const int nlevels = phi_int_tl.shape(1); + const double inv_g = 1.0 / grav; + + for (int level = 0; level < nlevels; ++level) { + for (size_t jnode = 0; jnode < grid_size; ++jnode) { + z_int_tl_view(jnode, level) = phi_int_tl_view(jnode, level) * inv_g; + } + } + + oops::Log::trace() << "leaving GeopotentialHeightAtInterface_A::executeTL function" + << std::endl; +} + +void GeopotentialHeightAtInterface_A::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & /*afieldsetTraj*/) +{ + oops::Log::trace() << "entering GeopotentialHeightAtInterface_A::executeAD function" + << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi_int_ad = afieldsetAD.field("geopotential_levels"); + atlas::Field z_int_ad = afieldsetAD.field("geopotential_height_levels"); + + auto phi_int_ad_view = atlas::array::make_view(phi_int_ad); + auto z_int_ad_view = atlas::array::make_view(z_int_ad); + + const size_t grid_size = phi_int_ad.shape(0); + const int nlevels = phi_int_ad.shape(1); + const double inv_g = 1.0 / grav; + + for (int level = 0; level < nlevels; ++level) { + for (size_t jnode = 0; jnode < grid_size; ++jnode) { + phi_int_ad_view(jnode, level) += z_int_ad_view(jnode, level) * inv_g; + z_int_ad_view(jnode, level) = 0.0; + } + } + + oops::Log::trace() << "leaving GeopotentialHeightAtInterface_A::executeAD function" + << std::endl; +} + +} // namespace vader diff --git a/src/vader/recipes/GeopotentialHeightAtSurface.h b/src/vader/recipes/GeopotentialHeightAtSurface.h new file mode 100644 index 0000000..d4897f1 --- /dev/null +++ b/src/vader/recipes/GeopotentialHeightAtSurface.h @@ -0,0 +1,64 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader { + +class GeopotentialHeightAtSurface_A_Parameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(GeopotentialHeightAtSurface_A_Parameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{ + "recipe name", + this}; +}; + +// ------------------------------------------------------------------------------------------------- + +/*! \brief GeopotentialHeightAtSurface_A computes geopotential height at surface + * + * Formula: z_surf = phi_surf / g + * where phi_surf is geopotential_at_surface (m^2 s^-2), g is standard_gravitational_acceleration (m s^-2), + * and z_surf is geopotential_height_at_surface (m). + * Full TL/AD support. + */ +class GeopotentialHeightAtSurface_A : public RecipeBase { + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef GeopotentialHeightAtSurface_A_Parameters Parameters_; + + GeopotentialHeightAtSurface_A(const Parameters_ &, const VaderConfigVars &); + + // Recipe base class overrides + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +} // namespace vader diff --git a/src/vader/recipes/GeopotentialHeightAtSurface_A.cc b/src/vader/recipes/GeopotentialHeightAtSurface_A.cc new file mode 100644 index 0000000..8a830e3 --- /dev/null +++ b/src/vader/recipes/GeopotentialHeightAtSurface_A.cc @@ -0,0 +1,136 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "atlas/util/Metadata.h" +#include "oops/util/Logger.h" +#include "vader/recipes/GeopotentialHeightAtSurface.h" + +namespace vader +{ +// ------------------------------------------------------------------------------------------------ + +// Static attribute initialization +const char GeopotentialHeightAtSurface_A::Name[] = "GeopotentialHeightAtSurface_A"; +const oops::Variables GeopotentialHeightAtSurface_A::Ingredients{std::vector{ + "geopotential_at_surface"}}; + +// Register the maker +static RecipeMaker makerGeopotentialHeightAtSurface_A_( + GeopotentialHeightAtSurface_A::Name); + +GeopotentialHeightAtSurface_A::GeopotentialHeightAtSurface_A(const Parameters_ & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} +{ + oops::Log::trace() << "GeopotentialHeightAtSurface_A::GeopotentialHeightAtSurface_A(params)" + << std::endl; +} + +std::string GeopotentialHeightAtSurface_A::name() const +{ + return GeopotentialHeightAtSurface_A::Name; +} + +oops::Variable GeopotentialHeightAtSurface_A::product() const +{ + return oops::Variable{"geopotential_height_at_surface"}; +} + +oops::Variables GeopotentialHeightAtSurface_A::ingredients() const +{ + return GeopotentialHeightAtSurface_A::Ingredients; +} + +size_t GeopotentialHeightAtSurface_A::productLevels(const atlas::FieldSet & afieldset) const +{ + return 1; +} + +atlas::FunctionSpace GeopotentialHeightAtSurface_A::productFunctionSpace(const atlas::FieldSet & + afieldset) const +{ + return afieldset.field("geopotential_at_surface").functionspace(); +} + +void GeopotentialHeightAtSurface_A::executeNL(atlas::FieldSet & afieldset) +{ + oops::Log::trace() << "entering GeopotentialHeightAtSurface_A::executeNL" << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi_surf = afieldset.field("geopotential_at_surface"); + atlas::Field z_surf = afieldset.field("geopotential_height_at_surface"); + + auto phi_surf_view = atlas::array::make_view(phi_surf); + auto z_surf_view = atlas::array::make_view(z_surf); + + const size_t grid_size = phi_surf.shape(0); + const double inv_g = 1.0 / grav; + + for ( size_t jnode = 0; jnode < grid_size ; ++jnode ) { + z_surf_view(jnode, 0) = phi_surf_view(jnode, 0) * inv_g; + } + oops::Log::trace() << "leaving GeopotentialHeightAtSurface_A::executeNL" << std::endl; +} + +void GeopotentialHeightAtSurface_A::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & /*afieldsetTraj*/) +{ + oops::Log::trace() << "entering GeopotentialHeightAtSurface_A::executeTL function" << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi_surf_tl = afieldsetTL.field("geopotential_at_surface"); + atlas::Field z_surf_tl = afieldsetTL.field("geopotential_height_at_surface"); + + auto phi_surf_tl_view = atlas::array::make_view(phi_surf_tl); + auto z_surf_tl_view = atlas::array::make_view(z_surf_tl); + + const size_t grid_size = phi_surf_tl.shape(0); + const double inv_g = 1.0 / grav; + + for (size_t jnode = 0; jnode < grid_size; ++jnode) { + z_surf_tl_view(jnode, 0) = phi_surf_tl_view(jnode, 0) * inv_g; + } + + oops::Log::trace() << "leaving GeopotentialHeightAtSurface_A::executeTL function" << std::endl; +} + +void GeopotentialHeightAtSurface_A::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & /*afieldsetTraj*/) +{ + oops::Log::trace() << "entering GeopotentialHeightAtSurface_A::executeAD function" << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi_surf_ad = afieldsetAD.field("geopotential_at_surface"); + atlas::Field z_surf_ad = afieldsetAD.field("geopotential_height_at_surface"); + + auto phi_surf_ad_view = atlas::array::make_view(phi_surf_ad); + auto z_surf_ad_view = atlas::array::make_view(z_surf_ad); + + const size_t grid_size = phi_surf_ad.shape(0); + const double inv_g = 1.0 / grav; + + for (size_t jnode = 0; jnode < grid_size; ++jnode) { + phi_surf_ad_view(jnode, 0) += z_surf_ad_view(jnode, 0) * inv_g; + z_surf_ad_view(jnode, 0) = 0.0; + } + + oops::Log::trace() << "leaving GeopotentialHeightAtSurface_A::executeAD function" << std::endl; +} + +} // namespace vader diff --git a/src/vader/recipes/GeopotentialHeight_A.cc b/src/vader/recipes/GeopotentialHeight_A.cc new file mode 100644 index 0000000..b405923 --- /dev/null +++ b/src/vader/recipes/GeopotentialHeight_A.cc @@ -0,0 +1,144 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "atlas/util/Metadata.h" +#include "oops/util/Logger.h" +#include "vader/recipes/GeopotentialHeight.h" + +namespace vader +{ +// ------------------------------------------------------------------------------------------------ + +// Static attribute initialization +const char GeopotentialHeight_A::Name[] = "GeopotentialHeight_A"; +const oops::Variables GeopotentialHeight_A::Ingredients{std::vector{ + "geopotential"}}; + +// Register the maker +static RecipeMaker makerGeopotentialHeight_A_(GeopotentialHeight_A::Name); + +GeopotentialHeight_A::GeopotentialHeight_A(const Parameters_ & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} +{ + oops::Log::trace() << "GeopotentialHeight_A::GeopotentialHeight_A(params)" << std::endl; +} + +std::string GeopotentialHeight_A::name() const +{ + return GeopotentialHeight_A::Name; +} + +oops::Variable GeopotentialHeight_A::product() const +{ + return oops::Variable{"geopotential_height"}; +} + +oops::Variables GeopotentialHeight_A::ingredients() const +{ + return GeopotentialHeight_A::Ingredients; +} + +size_t GeopotentialHeight_A::productLevels(const atlas::FieldSet & afieldset) const +{ + return afieldset.field("geopotential").shape(1); +} + +atlas::FunctionSpace GeopotentialHeight_A::productFunctionSpace(const atlas::FieldSet & + afieldset) const +{ + return afieldset.field("geopotential").functionspace(); +} + +void GeopotentialHeight_A::executeNL(atlas::FieldSet & afieldset) +{ + oops::Log::trace() << "entering GeopotentialHeight_A::executeNL function" << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi = afieldset.field("geopotential"); + atlas::Field z = afieldset.field("geopotential_height"); + + auto phi_view = atlas::array::make_view(phi); + auto z_view = atlas::array::make_view(z); + + const size_t grid_size = phi.shape(0); + const int nlevels = phi.shape(1); + const double inv_g = 1.0 / grav; + + for (int level = 0; level < nlevels; ++level) { + for ( size_t jnode = 0; jnode < grid_size ; ++jnode ) { + z_view(jnode, level) = phi_view(jnode, level) * inv_g; + } + } + oops::Log::trace() << "leaving GeopotentialHeight_A::executeNL function" << std::endl; +} + +void GeopotentialHeight_A::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & /*afieldsetTraj*/) +{ + oops::Log::trace() << "entering GeopotentialHeight_A::executeTL function" << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi_tl = afieldsetTL.field("geopotential"); + atlas::Field z_tl = afieldsetTL.field("geopotential_height"); + + auto phi_tl_view = atlas::array::make_view(phi_tl); + auto z_tl_view = atlas::array::make_view(z_tl); + + const size_t grid_size = phi_tl.shape(0); + const int nlevels = phi_tl.shape(1); + const double inv_g = 1.0 / grav; + + for (int level = 0; level < nlevels; ++level) { + for (size_t jnode = 0; jnode < grid_size; ++jnode) { + z_tl_view(jnode, level) = + phi_tl_view(jnode, level) * inv_g; + } + } + + oops::Log::trace() << "leaving GeopotentialHeight_A::executeTL function" << std::endl; +} + +void GeopotentialHeight_A::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & /*afieldsetTraj*/) +{ + oops::Log::trace() << "entering GeopotentialHeight_A::executeAD function" << std::endl; + + // Extract values from client config + const double grav = configVariables_.getDouble("standard_gravitational_acceleration"); + + atlas::Field phi_ad = afieldsetAD.field("geopotential"); + atlas::Field z_ad = afieldsetAD.field("geopotential_height"); + + auto phi_ad_view = atlas::array::make_view(phi_ad); + auto z_ad_view = atlas::array::make_view(z_ad); + + const size_t grid_size = phi_ad.shape(0); + const int nlevels = phi_ad.shape(1); + const double inv_g = 1.0 / grav; + + for (int level = 0; level < nlevels; ++level) { + for (size_t jnode = 0; jnode < grid_size; ++jnode) { + phi_ad_view(jnode, level) += z_ad_view(jnode, level) * inv_g; + z_ad_view(jnode, level) = 0.0; + } + } + + oops::Log::trace() << "leaving GeopotentialHeight_A::executeAD function" << std::endl; +} + +} // namespace vader diff --git a/src/vader/recipes/HeightAboveMeanSeaLevelAtSurface.h b/src/vader/recipes/HeightAboveMeanSeaLevelAtSurface.h new file mode 100644 index 0000000..2cce856 --- /dev/null +++ b/src/vader/recipes/HeightAboveMeanSeaLevelAtSurface.h @@ -0,0 +1,63 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader { + +class HeightAboveMeanSeaLevelAtSurface_A_Parameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(HeightAboveMeanSeaLevelAtSurface_A_Parameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{ + "recipe name", + this}; +}; + +// ------------------------------------------------------------------------------------------------- + +/*! \brief HeightAboveMeanSeaLevelAtSurface_A computes height above mean sea level at surface + * + * Formula: z_msl = z_surf (direct copy) + * where z_surf is geopotential_height_at_surface (m) and z_msl is height_above_mean_sea_level_at_surface (m). + * Full TL/AD support. + */ +class HeightAboveMeanSeaLevelAtSurface_A : public RecipeBase { + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef HeightAboveMeanSeaLevelAtSurface_A_Parameters Parameters_; + + HeightAboveMeanSeaLevelAtSurface_A(const Parameters_ &, const VaderConfigVars &); + + // Recipe base class overrides + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +} // namespace vader diff --git a/src/vader/recipes/HeightAboveMeanSeaLevelAtSurface_A.cc b/src/vader/recipes/HeightAboveMeanSeaLevelAtSurface_A.cc new file mode 100644 index 0000000..2c0f745 --- /dev/null +++ b/src/vader/recipes/HeightAboveMeanSeaLevelAtSurface_A.cc @@ -0,0 +1,129 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "atlas/util/Metadata.h" +#include "oops/util/Logger.h" +#include "vader/recipes/HeightAboveMeanSeaLevelAtSurface.h" + +namespace vader +{ +// ------------------------------------------------------------------------------------------------ + +// Static attribute initialization +const char HeightAboveMeanSeaLevelAtSurface_A::Name[] = "HeightAboveMeanSeaLevelAtSurface_A"; +const oops::Variables HeightAboveMeanSeaLevelAtSurface_A::Ingredients{std::vector{ + "geopotential_height_at_surface"}}; + +// Register the maker +static RecipeMaker makerHeightAboveMeanSeaLevelAtSurface_A_( + HeightAboveMeanSeaLevelAtSurface_A::Name); + +HeightAboveMeanSeaLevelAtSurface_A::HeightAboveMeanSeaLevelAtSurface_A(const Parameters_ & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} +{ + oops::Log::trace() << "HeightAboveMeanSeaLevelAtSurface_A::" + << "HeightAboveMeanSeaLevelAtSurface_A(params)" + << std::endl; +} + +std::string HeightAboveMeanSeaLevelAtSurface_A::name() const +{ + return HeightAboveMeanSeaLevelAtSurface_A::Name; +} + +oops::Variable HeightAboveMeanSeaLevelAtSurface_A::product() const +{ + return oops::Variable{"height_above_mean_sea_level_at_surface"}; +} + +oops::Variables HeightAboveMeanSeaLevelAtSurface_A::ingredients() const +{ + return HeightAboveMeanSeaLevelAtSurface_A::Ingredients; +} + +size_t HeightAboveMeanSeaLevelAtSurface_A::productLevels(const atlas::FieldSet & afieldset) const +{ + return 1; +} + +atlas::FunctionSpace HeightAboveMeanSeaLevelAtSurface_A::productFunctionSpace(const + atlas::FieldSet & afieldset) const +{ + return afieldset.field("geopotential_height_at_surface").functionspace(); +} + +void HeightAboveMeanSeaLevelAtSurface_A::executeNL(atlas::FieldSet & afieldset) +{ + oops::Log::trace() << "entering HeightAboveMeanSeaLevelAtSurface_A::executeNL" << std::endl; + + atlas::Field phi_surf = afieldset.field("geopotential_height_at_surface"); + atlas::Field z_surf = afieldset.field("height_above_mean_sea_level_at_surface"); + + auto phi_surf_view = atlas::array::make_view(phi_surf); + auto z_surf_view = atlas::array::make_view(z_surf); + + const size_t grid_size = phi_surf.shape(0); + + for ( size_t jnode = 0; jnode < grid_size ; ++jnode ) { + z_surf_view(jnode, 0) = phi_surf_view(jnode, 0); + } + oops::Log::trace() << "leaving HeightAboveMeanSeaLevelAtSurface_A::executeNL" << std::endl; +} + +void HeightAboveMeanSeaLevelAtSurface_A::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & /*afieldsetTraj*/) +{ + oops::Log::trace() << "entering HeightAboveMeanSeaLevelAtSurface_A::executeTL function" + << std::endl; + + atlas::Field phi_surf_tl = afieldsetTL.field("geopotential_height_at_surface"); + atlas::Field z_surf_tl = afieldsetTL.field("height_above_mean_sea_level_at_surface"); + + auto phi_surf_tl_view = atlas::array::make_view(phi_surf_tl); + auto z_surf_tl_view = atlas::array::make_view(z_surf_tl); + + const size_t grid_size = phi_surf_tl.shape(0); + + for (size_t jnode = 0; jnode < grid_size; ++jnode) { + z_surf_tl_view(jnode, 0) = phi_surf_tl_view(jnode, 0); + } + + oops::Log::trace() << "leaving HeightAboveMeanSeaLevelAtSurface_A::executeTL function" + << std::endl; +} + +void HeightAboveMeanSeaLevelAtSurface_A::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & /*afieldsetTraj*/) +{ + oops::Log::trace() << "entering HeightAboveMeanSeaLevelAtSurface_A::executeAD function" + << std::endl; + + atlas::Field phi_surf_ad = afieldsetAD.field("geopotential_height_at_surface"); + atlas::Field z_surf_ad = afieldsetAD.field("height_above_mean_sea_level_at_surface"); + + auto phi_surf_ad_view = atlas::array::make_view(phi_surf_ad); + auto z_surf_ad_view = atlas::array::make_view(z_surf_ad); + + const size_t grid_size = phi_surf_ad.shape(0); + + for (size_t jnode = 0; jnode < grid_size; ++jnode) { + phi_surf_ad_view(jnode, 0) += z_surf_ad_view(jnode, 0); + z_surf_ad_view(jnode, 0) = 0.0; + } + + oops::Log::trace() << "leaving HeightAboveMeanSeaLevelAtSurface_A::executeAD function" + << std::endl; +} + +} // namespace vader diff --git a/src/vader/recipes/LnAirPressure.h b/src/vader/recipes/LnAirPressure.h new file mode 100644 index 0000000..20ecaa6 --- /dev/null +++ b/src/vader/recipes/LnAirPressure.h @@ -0,0 +1,62 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class LnAirPressure_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(LnAirPressure_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief LnAirPressure_A computes natural logarithm of air pressure + * + * Formula: ln_p = ln(p) + * where p is air_pressure (Pa) and ln_p is ln(air_pressure) (dimensionless). + * Full TL/AD support. + */ +class LnAirPressure_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef LnAirPressure_AParameters Parameters_; + + LnAirPressure_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + oops::Variables trajectoryVars() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/LnAirPressureAtInterface.h b/src/vader/recipes/LnAirPressureAtInterface.h index bd5adf1..d13b94a 100644 --- a/src/vader/recipes/LnAirPressureAtInterface.h +++ b/src/vader/recipes/LnAirPressureAtInterface.h @@ -1,5 +1,5 @@ /* - * (C) Copyright 203 UCAR + * (C) Copyright 2023 UCAR * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. @@ -26,31 +26,35 @@ class LnAirPressureAtInterface_AParameters : public RecipeParametersBase { OOPS_CONCRETE_PARAMETERS(LnAirPressureAtInterface_AParameters, RecipeParametersBase) public: - oops::RequiredParameter name{"recipe name", this}; + oops::RequiredParameter name{"recipe name", this}; }; -/*! \brief LnAirPressureAtInterface_A class defines a recipe for pressure levels from pressure - thickness. +/*! \brief LnAirPressureAtInterface_A computes natural logarithm of air pressure at interfaces * - * \details This recipe uses pressure at the interfaces, along with the Phillips method to - * compute pressure at the mid points. It does not provide TL/AD algorithms. + * Formula: ln_p_int = ln(p_int) + * where p_int is air_pressure_levels (Pa) and ln_p_int is ln_air_pressure_at_interface (dimensionless). + * Full TL/AD support. */ class LnAirPressureAtInterface_A : public RecipeBase { public: - static const char Name[]; - static const oops::Variables Ingredients; - - typedef LnAirPressureAtInterface_AParameters Parameters_; - - LnAirPressureAtInterface_A(const Parameters_ &, const VaderConfigVars &); - - std::string name() const override; - oops::Variable product() const override; - oops::Variables ingredients() const override; - size_t productLevels(const atlas::FieldSet &) const override; - atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; - void executeNL(atlas::FieldSet &) override; + static const char Name[]; + static const oops::Variables Ingredients; + + typedef LnAirPressureAtInterface_AParameters Parameters_; + + LnAirPressureAtInterface_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + oops::Variables trajectoryVars() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; }; // ------------------------------------------------------------------------------------------------- diff --git a/src/vader/recipes/LnAirPressureAtInterface_A.cc b/src/vader/recipes/LnAirPressureAtInterface_A.cc index 4767729..adfc704 100644 --- a/src/vader/recipes/LnAirPressureAtInterface_A.cc +++ b/src/vader/recipes/LnAirPressureAtInterface_A.cc @@ -1,11 +1,11 @@ /* - * (C) Copyright 2023 UCAR. + * (C) Copyright 2025 UCAR * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. */ -#include +#include #include #include @@ -17,88 +17,146 @@ namespace vader { - // ------------------------------------------------------------------------------------------------- // Static attribute initialization const char LnAirPressureAtInterface_A::Name[] = "LnAirPressureAtInterface_A"; const oops::Variables LnAirPressureAtInterface_A::Ingredients{ - std::vector{"air_pressure_levels"}}; + std::vector{"air_pressure_levels"}}; // ------------------------------------------------------------------------------------------------- // Register the maker static RecipeMaker - makerLnAirPressureAtInterface_A_(LnAirPressureAtInterface_A::Name); + makerLnAirPressureAtInterface_A_(LnAirPressureAtInterface_A::Name); // ------------------------------------------------------------------------------------------------- LnAirPressureAtInterface_A::LnAirPressureAtInterface_A(const - LnAirPressureAtInterface_AParameters & params, - const VaderConfigVars & configVariables) -{ - oops::Log::trace() << "LnAirPressureAtInterface_A::LnAirPressureAtInterface_A Starting" - << std::endl; - oops::Log::trace() << "LnAirPressureAtInterface_A::LnAirPressureAtInterface_A Done" - << std::endl; + LnAirPressureAtInterface_AParameters & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} { + oops::Log::trace() << "LnAirPressureAtInterface_A::LnAirPressureAtInterface_A Starting" + << std::endl; + oops::Log::trace() << "LnAirPressureAtInterface_A::LnAirPressureAtInterface_A Done" + << std::endl; } // ------------------------------------------------------------------------------------------------- std::string LnAirPressureAtInterface_A::name() const { - return LnAirPressureAtInterface_A::Name; + return LnAirPressureAtInterface_A::Name; } // ------------------------------------------------------------------------------------------------- oops::Variable LnAirPressureAtInterface_A::product() const { - return oops::Variable{"ln_air_pressure_at_interface"}; + return oops::Variable{"ln_air_pressure_at_interface"}; } // ------------------------------------------------------------------------------------------------- oops::Variables LnAirPressureAtInterface_A::ingredients() const { - return LnAirPressureAtInterface_A::Ingredients; + return LnAirPressureAtInterface_A::Ingredients; +} + +// ------------------------------------------------------------------------------------------------- + +oops::Variables LnAirPressureAtInterface_A::trajectoryVars() const { + return oops::Variables{std::vector{"air_pressure"}}; } // ------------------------------------------------------------------------------------------------- size_t LnAirPressureAtInterface_A::productLevels(const atlas::FieldSet & afieldset) const { - return afieldset.field("air_pressure_levels").shape(1); + return afieldset.field("air_pressure_levels").shape(1); } // ------------------------------------------------------------------------------------------------- atlas::FunctionSpace LnAirPressureAtInterface_A::productFunctionSpace(const - atlas::FieldSet & afieldset) const { - return afieldset.field("air_pressure_levels").functionspace(); + atlas::FieldSet & afieldset) const { + return afieldset.field("air_pressure_levels").functionspace(); } // ------------------------------------------------------------------------------------------------- void LnAirPressureAtInterface_A::executeNL(atlas::FieldSet & afieldset) { - // - oops::Log::trace() << "LnAirPressureAtInterface_A::executeNL Starting" << std::endl; - - // Get fields - atlas::Field airPressureLevelsF = afieldset.field("air_pressure_levels"); - atlas::Field lnAirPressureAtInterfaceF = afieldset.field("ln_air_pressure_at_interface"); - - // Get field views - auto airPressureLevels = atlas::array::make_view(airPressureLevelsF); - auto lnAirPressureAtInterface = atlas::array::make_view(lnAirPressureAtInterfaceF); - - // Grid dimensions - size_t h_size = airPressureLevelsF.shape(0); - int v_size = airPressureLevelsF.shape(1); - - // Calculate the output variable - for (int vv = 0; vv < v_size; ++vv) { - for ( size_t hh = 0; hh < h_size ; ++hh ) { - lnAirPressureAtInterface(hh, vv) = log(airPressureLevels(hh, vv)); - } + oops::Log::trace() << "LnAirPressureAtInterface_A::executeNL Starting" << std::endl; + + // Get fields + atlas::Field p_int = afieldset.field("air_pressure_levels"); + atlas::Field ln_p_int = afieldset.field("ln_air_pressure_at_interface"); + + // Get field views + auto p_int_view = atlas::array::make_view(p_int); + auto ln_p_int_view = atlas::array::make_view(ln_p_int); + + // Grid dimensions + const size_t npoints = p_int.shape(0); + const size_t nLevel = p_int.shape(1); + + // Calculate the output variable + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < gridSize; ++jn) { + ln_p_int_view(jn, jl) = std::log(p_int_view(jn, jl)); + } + } + oops::Log::trace() << "LnAirPressureAtInterface_A::executeNL Done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +void LnAirPressureAtInterface_A::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "LnAirPressureAtInterface_A::executeTL Starting" << std::endl; + + atlas::Field p_int_tl = afieldsetTL.field("air_pressure_levels"); + atlas::Field ln_p_int_tl = afieldsetTL.field("ln_air_pressure_at_interface"); + atlas::Field p_int = afieldsetTraj.field("air_pressure_levels"); + + auto p_int_tl_view = atlas::array::make_view(p_int_tl); + auto ln_p_int_tl_view = atlas::array::make_view(ln_p_int_tl); + auto p_int_view = atlas::array::make_view(p_int); + + const size_t npoints = p_int_tl.shape(0); + const size_t nLevel = p_int_tl.shape(1); + + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < gridSize; ++jn) { + ln_p_int_tl_view(jn, jl) = p_int_tl_view(jn, jl) / p_int_view(jn, jl); + } + } + oops::Log::trace() << "LnAirPressureAtInterface_A::executeTL Done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +void LnAirPressureAtInterface_A::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "LnAirPressureAtInterface_A::executeAD Starting" << std::endl; + + atlas::Field p_int_ad = afieldsetAD.field("air_pressure_levels"); + atlas::Field ln_p_int_ad = afieldsetAD.field("ln_air_pressure_at_interface"); + atlas::Field p_int = afieldsetTraj.field("air_pressure_levels"); + + auto p_int_ad_view = atlas::array::make_view(p_int_ad); + auto ln_p_int_ad_view = atlas::array::make_view(ln_p_int_ad); + auto p_int_view = atlas::array::make_view(p_int); + + const size_t npoints = p_int_ad.shape(0); + const size_t nLevel = p_int_ad.shape(1); + + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < gridSize; ++jn) { + const double lam = ln_p_int_ad_view(jn, jl); + if (lam != 0.0) { + p_int_ad_view(jn, jl) += lam / p_int_view(jn, jl); + ln_p_int_ad_view(jn, jl) = 0.0; } - oops::Log::trace() << "LnAirPressureAtInterface_A::executeNL Done" << std::endl; + } + } + oops::Log::trace() << "LnAirPressureAtInterface_A::executeAD Done" << std::endl; } // ------------------------------------------------------------------------------------------------- diff --git a/src/vader/recipes/LnAirPressure_A.cc b/src/vader/recipes/LnAirPressure_A.cc new file mode 100644 index 0000000..7d3ad9a --- /dev/null +++ b/src/vader/recipes/LnAirPressure_A.cc @@ -0,0 +1,160 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "atlas/util/Metadata.h" +#include "oops/util/Logger.h" +#include "vader/recipes/LnAirPressure.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- + +// Static attribute initialization +const char LnAirPressure_A::Name[] = "LnAirPressure_A"; +const oops::Variables LnAirPressure_A::Ingredients{ + std::vector{"air_pressure"}}; + +// ------------------------------------------------------------------------------------------------- + +// Register the maker +static RecipeMaker makerLnAirPressure_A_(LnAirPressure_A::Name); + +// ------------------------------------------------------------------------------------------------- + +LnAirPressure_A::LnAirPressure_A(const LnAirPressure_AParameters & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} { + oops::Log::trace() << "LnAirPressure_A::LnAirPressure_A Starting" << std::endl; + oops::Log::trace() << "LnAirPressure_A::LnAirPressure_A Done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +std::string LnAirPressure_A::name() const { + return LnAirPressure_A::Name; +} + +// ------------------------------------------------------------------------------------------------- + +oops::Variable LnAirPressure_A::product() const { + return oops::Variable{"ln_air_pressure"}; +} + +// ------------------------------------------------------------------------------------------------- + +oops::Variables LnAirPressure_A::ingredients() const { + return LnAirPressure_A::Ingredients; +} + +// ------------------------------------------------------------------------------------------------- + +oops::Variables LnAirPressure_A::trajectoryVars() const { + return oops::Variables{std::vector{"air_pressure"}}; +} + +// ------------------------------------------------------------------------------------------------- + +size_t LnAirPressure_A::productLevels(const atlas::FieldSet & afieldset) const { + return afieldset.field("air_pressure").shape(1); +} + +// ------------------------------------------------------------------------------------------------- + +atlas::FunctionSpace LnAirPressure_A::productFunctionSpace(const + atlas::FieldSet & afieldset) const { + return afieldset.field("air_pressure").functionspace(); +} + +// ------------------------------------------------------------------------------------------------- + +void LnAirPressure_A::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "LnAirPressure_A::executeNL Starting" << std::endl; + + // Get fields + atlas::Field p = afieldset.field("air_pressure"); + atlas::Field ln_p = afieldset.field("ln_air_pressure"); + + // Get field views + auto p_view = atlas::array::make_view(p); + auto ln_p_view = atlas::array::make_view(ln_p); + + // Grid dimensions + const size_t npoints = p.shape(0); + const size_t nlevels = p.shape(1); + + // Calculate the output variable + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < npoints; ++jn) { + ln_p_view(jn, jl) = std::log(p_view(jn, jl)); + } + } + oops::Log::trace() << "LnAirPressure_A::executeNL Done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +void LnAirPressure_A::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "LnAirPressure_A::executeTL Starting" << std::endl; + + atlas::Field p_tl = afieldsetTL.field("air_pressure"); + atlas::Field ln_p_tl = afieldsetTL.field("ln_air_pressure"); + atlas::Field p = afieldsetTraj.field("air_pressure"); + + auto p_tl_view = atlas::array::make_view(p_tl); + auto ln_p_tl_view = atlas::array::make_view(ln_p_tl); + auto p_view = atlas::array::make_view(p); + + const size_t npoints = p_tl.shape(0); + const size_t nlevels = p_tl.shape(1); + + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < npoints; ++jn) { + ln_p_tl_view(jn, jl) = p_tl_view(jn, jl) / p_view(jn, jl); + } + } + oops::Log::trace() << "LnAirPressure_A::executeTL Done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +void LnAirPressure_A::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "LnAirPressure_A::executeAD Starting" << std::endl; + + atlas::Field p_ad = afieldsetAD.field("air_pressure"); + atlas::Field ln_p_ad = afieldsetAD.field("ln_air_pressure"); + atlas::Field p = afieldsetTraj.field("air_pressure"); + + auto p_ad_view = atlas::array::make_view(p_ad); + auto ln_p_ad_view = atlas::array::make_view(ln_p_ad); + auto p_view = atlas::array::make_view(p); + + const size_t npoints = p_ad.shape(0); + const size_t nlevels = p_ad.shape(1); + + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < npoints; ++jn) { + const double lam = ln_p_ad_view(jn, jl); + if (lam != 0.0) { + p_ad_view(jn, jl) += lam / p_view(jn, jl); + ln_p_ad_view(jn, jl) = 0.0; + } + } + } + oops::Log::trace() << "LnAirPressure_A::executeAD Done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfCloudIceInAtmosphereLayer.h b/src/vader/recipes/MassContentOfCloudIceInAtmosphereLayer.h new file mode 100644 index 0000000..6275e27 --- /dev/null +++ b/src/vader/recipes/MassContentOfCloudIceInAtmosphereLayer.h @@ -0,0 +1,65 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class MassContentOfCloudIceInAtmosphereLayer_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(MassContentOfCloudIceInAtmosphereLayer_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'MassContentOfCloudIceInAtmosphereLayer_A' defines a recipe for + * mass content of cloud ice in atmosphere layer + * + * \details This instantiation of RecipeBase produces mass content of cloud ice + * using a simplified diagnostic approximation based on pressure-dependent + * RH thresholds and temperature, based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature, air_pressure, relative_humidity, + * air_pressure_thickness, air_pressure_at_surface + * Output units: kg m-2 + */ +class MassContentOfCloudIceInAtmosphereLayer_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef MassContentOfCloudIceInAtmosphereLayer_AParameters Parameters_; + + MassContentOfCloudIceInAtmosphereLayer_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfCloudIceInAtmosphereLayer_A.cc b/src/vader/recipes/MassContentOfCloudIceInAtmosphereLayer_A.cc new file mode 100644 index 0000000..c87214b --- /dev/null +++ b/src/vader/recipes/MassContentOfCloudIceInAtmosphereLayer_A.cc @@ -0,0 +1,123 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/MassContentOfCloudIceInAtmosphereLayer.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char MassContentOfCloudIceInAtmosphereLayer_A::Name[] = + "MassContentOfCloudIceInAtmosphereLayer_A"; +const oops::Variables MassContentOfCloudIceInAtmosphereLayer_A::Ingredients{ + std::vector{ + "air_temperature", "air_pressure", "relative_humidity", + "air_pressure_thickness", "air_pressure_at_surface"}}; + +// Register the maker +static RecipeMaker + makerMassContentOfCloudIceInAtmosphereLayer_A_( + MassContentOfCloudIceInAtmosphereLayer_A::Name); + +// ------------------------------------------------------------------------------------------------- + +MassContentOfCloudIceInAtmosphereLayer_A::MassContentOfCloudIceInAtmosphereLayer_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "MassContentOfCloudIceInAtmosphereLayer_A::" + << "MassContentOfCloudIceInAtmosphereLayer_A" << std::endl; +} + +std::string MassContentOfCloudIceInAtmosphereLayer_A::name() const { + return MassContentOfCloudIceInAtmosphereLayer_A::Name; +} + +oops::Variable MassContentOfCloudIceInAtmosphereLayer_A::product() const { + return oops::Variable{"mass_content_of_cloud_ice_in_atmosphere_layer"}; +} + +oops::Variables MassContentOfCloudIceInAtmosphereLayer_A::ingredients() const { + return MassContentOfCloudIceInAtmosphereLayer_A::Ingredients; +} + +size_t MassContentOfCloudIceInAtmosphereLayer_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace MassContentOfCloudIceInAtmosphereLayer_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void MassContentOfCloudIceInAtmosphereLayer_A::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "MassContentOfCloudIceInAtmosphereLayer_A::executeNL starting" + << std::endl; + + // Get input fields + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + auto pressure = atlas::array::make_view(afieldset.field("air_pressure")); + auto rh = atlas::array::make_view(afieldset.field("relative_humidity")); + auto delta_p = atlas::array::make_view(afieldset.field("air_pressure_thickness")); + auto ps = atlas::array::make_view(afieldset.field("air_pressure_at_surface")); + + // Create output field + auto ice_field = afieldset.field("mass_content_of_cloud_ice_in_atmosphere_layer"); + auto ice = atlas::array::make_view(ice_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Physical constants + const double g = configVariables_.getDouble("standard_gravitational_acceleration"); + // Diagnostic cloud scheme parameters + const double RH_crit_base = 0.75; // Base critical RH (pressure-dependent) + const double T_freeze = 273.15; // Freezing point (K) + const double T_cold = 233.15; // Cold cloud threshold (-40C) + + // Calculate mass content for each grid point and level + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + const double RH = rh(jn, jl); + + // Pressure-dependent critical RH + const double P = pressure(jn, jl); + const double p_ratio = P / ps(jn, 0); + const double RH_crit = RH_crit_base + 0.10 * p_ratio; + + // Only form ice clouds below freezing + if (T < T_freeze && RH > RH_crit) { + // Excess relative humidity + const double RH_excess = std::max(0.0, RH - RH_crit); + + // Ice water content (kg/kg) - simple diagnostic based on RH and temperature + // More ice at colder temperatures + const double temp_factor = std::max(0.0, (T_freeze - T) / (T_freeze - T_cold)); + const double q_ice = 1e-4 * RH_excess / (1.0 - RH_crit) * (1.0 + temp_factor); + // Mass content = mixing ratio * pressure thickness / g (kg/m2) + ice(jn, jl) = q_ice * delta_p(jn, jl) / g; + } else { + ice(jn, jl) = 0.0; + } + } + } + + oops::Log::trace() << "MassContentOfCloudIceInAtmosphereLayer_A::executeNL done" + << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer.h b/src/vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer.h new file mode 100644 index 0000000..e2cee66 --- /dev/null +++ b/src/vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer.h @@ -0,0 +1,66 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class MassContentOfCloudLiquidWaterInAtmosphereLayer_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(MassContentOfCloudLiquidWaterInAtmosphereLayer_AParameters, + RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'MassContentOfCloudLiquidWaterInAtmosphereLayer_A' defines a recipe for + * mass content of cloud liquid water in atmosphere layer + * + * \details This instantiation of RecipeBase produces mass content of cloud liquid water + * using a simplified diagnostic approximation based on pressure-dependent + * RH thresholds, based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature, air_pressure, relative_humidity, + * air_pressure_thickness, air_pressure_at_surface + * Output units: kg m-2 + */ +class MassContentOfCloudLiquidWaterInAtmosphereLayer_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef MassContentOfCloudLiquidWaterInAtmosphereLayer_AParameters Parameters_; + + MassContentOfCloudLiquidWaterInAtmosphereLayer_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer_A.cc b/src/vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer_A.cc new file mode 100644 index 0000000..ea07aa6 --- /dev/null +++ b/src/vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer_A.cc @@ -0,0 +1,122 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/MassContentOfCloudLiquidWaterInAtmosphereLayer.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char MassContentOfCloudLiquidWaterInAtmosphereLayer_A::Name[] = + "MassContentOfCloudLiquidWaterInAtmosphereLayer_A"; +const oops::Variables MassContentOfCloudLiquidWaterInAtmosphereLayer_A::Ingredients{ + std::vector{ + "air_temperature", "air_pressure", "relative_humidity", + "air_pressure_thickness", "air_pressure_at_surface"}}; + +// Register the maker +static RecipeMaker + makerMassContentOfCloudLiquidWaterInAtmosphereLayer_A_( + MassContentOfCloudLiquidWaterInAtmosphereLayer_A::Name); + +// ------------------------------------------------------------------------------------------------- + +MassContentOfCloudLiquidWaterInAtmosphereLayer_A::MassContentOfCloudLiquidWaterInAtmosphereLayer_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "MassContentOfCloudLiquidWaterInAtmosphereLayer_A::" + << "MassContentOfCloudLiquidWaterInAtmosphereLayer_A" << std::endl; +} + +std::string MassContentOfCloudLiquidWaterInAtmosphereLayer_A::name() const { + return MassContentOfCloudLiquidWaterInAtmosphereLayer_A::Name; +} + +oops::Variable MassContentOfCloudLiquidWaterInAtmosphereLayer_A::product() const { + return oops::Variable{"mass_content_of_cloud_liquid_water_in_atmosphere_layer"}; +} + +oops::Variables MassContentOfCloudLiquidWaterInAtmosphereLayer_A::ingredients() const { + return MassContentOfCloudLiquidWaterInAtmosphereLayer_A::Ingredients; +} + +size_t MassContentOfCloudLiquidWaterInAtmosphereLayer_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace MassContentOfCloudLiquidWaterInAtmosphereLayer_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void MassContentOfCloudLiquidWaterInAtmosphereLayer_A::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "MassContentOfCloudLiquidWaterInAtmosphereLayer_A::executeNL starting" + << std::endl; + + // Get input fields + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + auto pressure = atlas::array::make_view(afieldset.field("air_pressure")); + auto rh = atlas::array::make_view(afieldset.field("relative_humidity")); + auto delta_p = atlas::array::make_view(afieldset.field("air_pressure_thickness")); + auto ps = atlas::array::make_view(afieldset.field("air_pressure_at_surface")); + + // Create output field + auto clw_field = afieldset.field("mass_content_of_cloud_liquid_water_in_atmosphere_layer"); + auto clw = atlas::array::make_view(clw_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Physical constants + const double g = configVariables_.getDouble("standard_gravitational_acceleration"); + // Diagnostic cloud scheme parameters + const double RH_crit_base = 0.75; // Base critical RH (pressure-dependent) + const double T_freeze = 273.15; // Freezing point (K) + const double T_warm = 288.15; // Temperature for warm clouds (K) + + // Calculate mass content for each grid point and level + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + const double RH = rh(jn, jl); + + // Pressure-dependent critical RH + const double P = pressure(jn, jl); + const double p_ratio = P / ps(jn, 0); + const double RH_crit = RH_crit_base + 0.10 * p_ratio; + + // Only form liquid water clouds above freezing + if (T > T_freeze && RH > RH_crit) { + // Excess relative humidity + const double RH_excess = std::max(0.0, RH - RH_crit); + + // Liquid water content (kg/kg) - simple diagnostic based on RH + // Assumes cloud liquid water increases with RH above critical threshold + const double q_liq = 1e-4 * RH_excess / (1.0 - RH_crit); + // Mass content = mixing ratio * pressure thickness / g (kg/m2) + clw(jn, jl) = q_liq * delta_p(jn, jl) / g; + } else { + clw(jn, jl) = 0.0; + } + } + } + + oops::Log::trace() << "MassContentOfCloudLiquidWaterInAtmosphereLayer_A::executeNL done" + << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfGraupelInAtmosphereLayer.h b/src/vader/recipes/MassContentOfGraupelInAtmosphereLayer.h new file mode 100644 index 0000000..5893ecb --- /dev/null +++ b/src/vader/recipes/MassContentOfGraupelInAtmosphereLayer.h @@ -0,0 +1,65 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class MassContentOfGraupelInAtmosphereLayer_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(MassContentOfGraupelInAtmosphereLayer_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'MassContentOfGraupelInAtmosphereLayer_A' defines a recipe for + * mass content of graupel in atmosphere layer + * + * \details This instantiation of RecipeBase produces mass content of graupel + * using a simplified diagnostic approximation based on mixed-phase conditions, + * based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature, air_pressure, relative_humidity, + * air_pressure_thickness, air_pressure_at_surface + * Output units: kg m-2 + */ +class MassContentOfGraupelInAtmosphereLayer_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef MassContentOfGraupelInAtmosphereLayer_AParameters Parameters_; + + MassContentOfGraupelInAtmosphereLayer_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfGraupelInAtmosphereLayer_A.cc b/src/vader/recipes/MassContentOfGraupelInAtmosphereLayer_A.cc new file mode 100644 index 0000000..bcf8686 --- /dev/null +++ b/src/vader/recipes/MassContentOfGraupelInAtmosphereLayer_A.cc @@ -0,0 +1,127 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/MassContentOfGraupelInAtmosphereLayer.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char MassContentOfGraupelInAtmosphereLayer_A::Name[] = + "MassContentOfGraupelInAtmosphereLayer_A"; +const oops::Variables MassContentOfGraupelInAtmosphereLayer_A::Ingredients{ + std::vector{ + "air_temperature", "air_pressure", "relative_humidity", + "air_pressure_thickness", "air_pressure_at_surface"}}; + +// Register the maker +static RecipeMaker + makerMassContentOfGraupelInAtmosphereLayer_A_( + MassContentOfGraupelInAtmosphereLayer_A::Name); + +// ------------------------------------------------------------------------------------------------- + +MassContentOfGraupelInAtmosphereLayer_A::MassContentOfGraupelInAtmosphereLayer_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "MassContentOfGraupelInAtmosphereLayer_A::" + << "MassContentOfGraupelInAtmosphereLayer_A" << std::endl; +} + +std::string MassContentOfGraupelInAtmosphereLayer_A::name() const { + return MassContentOfGraupelInAtmosphereLayer_A::Name; +} + +oops::Variable MassContentOfGraupelInAtmosphereLayer_A::product() const { + return oops::Variable{"mass_content_of_graupel_in_atmosphere_layer"}; +} + +oops::Variables MassContentOfGraupelInAtmosphereLayer_A::ingredients() const { + return MassContentOfGraupelInAtmosphereLayer_A::Ingredients; +} + +size_t MassContentOfGraupelInAtmosphereLayer_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace MassContentOfGraupelInAtmosphereLayer_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void MassContentOfGraupelInAtmosphereLayer_A::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "MassContentOfGraupelInAtmosphereLayer_A::executeNL starting" + << std::endl; + + // Get input fields + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + auto pressure = atlas::array::make_view(afieldset.field("air_pressure")); + auto rh = atlas::array::make_view(afieldset.field("relative_humidity")); + auto delta_p = atlas::array::make_view(afieldset.field("air_pressure_thickness")); + auto ps = atlas::array::make_view(afieldset.field("air_pressure_at_surface")); + + // Create output field + auto graupel_field = afieldset.field("mass_content_of_graupel_in_atmosphere_layer"); + auto graupel = atlas::array::make_view(graupel_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Physical constants + const double g = configVariables_.getDouble("standard_gravitational_acceleration"); + // Diagnostic scheme parameters for graupel (rimed ice) + const double RH_crit_base = 0.88; // High RH for graupel formation + const double T_freeze = 273.15; // Freezing point (K) + const double T_graupel_min = 258.15; // Min temp for graupel (-15C) + const double T_graupel_max = 268.15; // Max temp for graupel (-5C) + + // Calculate mass content for each grid point and level + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + const double RH = rh(jn, jl); + + // Pressure-dependent critical RH + const double P = pressure(jn, jl); + const double p_ratio = P / ps(jn, 0); + const double RH_crit = RH_crit_base + 0.10 * p_ratio; + + // Graupel forms in mixed-phase clouds with high RH (riming region) + if (T < T_graupel_max && T > T_graupel_min && RH > RH_crit) { + // Excess relative humidity + const double RH_excess = std::max(0.0, RH - RH_crit); + + // Graupel content from riming (kg/kg) + // Optimal in mixed-phase region around -10C + const double T_optimal = 263.15; // -10C + const double temp_factor = 1.0 - std::abs(T - T_optimal) / 10.0; + const double q_graupel = 6e-5 * RH_excess / (1.0 - RH_crit) * + std::max(0.0, temp_factor); + + // Mass content (kg/m2) + graupel(jn, jl) = q_graupel * delta_p(jn, jl) / g; + } else { + graupel(jn, jl) = 0.0; + } + } + } + + oops::Log::trace() << "MassContentOfGraupelInAtmosphereLayer_A::executeNL done" + << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfHailInAtmosphereLayer.h b/src/vader/recipes/MassContentOfHailInAtmosphereLayer.h new file mode 100644 index 0000000..d61b292 --- /dev/null +++ b/src/vader/recipes/MassContentOfHailInAtmosphereLayer.h @@ -0,0 +1,65 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class MassContentOfHailInAtmosphereLayer_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(MassContentOfHailInAtmosphereLayer_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'MassContentOfHailInAtmosphereLayer_A' defines a recipe for + * mass content of hail in atmosphere layer + * + * \details This instantiation of RecipeBase produces mass content of hail + * using a simplified diagnostic approximation based on strong updraft conditions, + * based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature, air_pressure, relative_humidity, + * air_pressure_thickness, air_pressure_at_surface + * Output units: kg m-2 + */ +class MassContentOfHailInAtmosphereLayer_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef MassContentOfHailInAtmosphereLayer_AParameters Parameters_; + + MassContentOfHailInAtmosphereLayer_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfHailInAtmosphereLayer_A.cc b/src/vader/recipes/MassContentOfHailInAtmosphereLayer_A.cc new file mode 100644 index 0000000..bfebdb9 --- /dev/null +++ b/src/vader/recipes/MassContentOfHailInAtmosphereLayer_A.cc @@ -0,0 +1,126 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/MassContentOfHailInAtmosphereLayer.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char MassContentOfHailInAtmosphereLayer_A::Name[] = + "MassContentOfHailInAtmosphereLayer_A"; +const oops::Variables MassContentOfHailInAtmosphereLayer_A::Ingredients{ + std::vector{ + "air_temperature", "air_pressure", "relative_humidity", + "air_pressure_thickness", "air_pressure_at_surface"}}; + +// Register the maker +static RecipeMaker + makerMassContentOfHailInAtmosphereLayer_A_( + MassContentOfHailInAtmosphereLayer_A::Name); + +// ------------------------------------------------------------------------------------------------- + +MassContentOfHailInAtmosphereLayer_A::MassContentOfHailInAtmosphereLayer_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "MassContentOfHailInAtmosphereLayer_A::" + << "MassContentOfHailInAtmosphereLayer_A" << std::endl; +} + +std::string MassContentOfHailInAtmosphereLayer_A::name() const { + return MassContentOfHailInAtmosphereLayer_A::Name; +} + +oops::Variable MassContentOfHailInAtmosphereLayer_A::product() const { + return oops::Variable{"mass_content_of_hail_in_atmosphere_layer"}; +} + +oops::Variables MassContentOfHailInAtmosphereLayer_A::ingredients() const { + return MassContentOfHailInAtmosphereLayer_A::Ingredients; +} + +size_t MassContentOfHailInAtmosphereLayer_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace MassContentOfHailInAtmosphereLayer_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void MassContentOfHailInAtmosphereLayer_A::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "MassContentOfHailInAtmosphereLayer_A::executeNL starting" + << std::endl; + + // Get input fields + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + auto pressure = atlas::array::make_view(afieldset.field("air_pressure")); + auto rh = atlas::array::make_view(afieldset.field("relative_humidity")); + auto delta_p = atlas::array::make_view(afieldset.field("air_pressure_thickness")); + auto ps = atlas::array::make_view(afieldset.field("air_pressure_at_surface")); + + // Create output field + auto hail_field = afieldset.field("mass_content_of_hail_in_atmosphere_layer"); + auto hail = atlas::array::make_view(hail_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Physical constants + const double g = configVariables_.getDouble("standard_gravitational_acceleration"); + + // Diagnostic scheme parameters for hail + // Hail requires very strong convection - use very high RH threshold + const double RH_crit_base = 0.92; // Base critical RH for hail + const double T_freeze = 273.15; // Freezing point (K) + const double T_hail_min = 253.15; // Min temp for hail (-20C) + const double T_hail_max = 268.15; // Max temp for hail (-5C) + + // Calculate mass content for each grid point and level + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + const double RH = rh(jn, jl); + + // Pressure-dependent critical RH + const double P = pressure(jn, jl); + const double p_ratio = P / ps(jn, 0); + const double RH_crit = RH_crit_base + 0.10 * p_ratio; + + // Hail forms in strong convective mixed-phase clouds + // Require very high RH as proxy for strong updrafts + if (T < T_hail_max && T > T_hail_min && RH > RH_crit) { + // Excess relative humidity + const double RH_excess = std::max(0.0, RH - RH_crit); + + // Hail content (kg/kg) - less common than other hydrometeors + const double q_hail = 3e-5 * RH_excess / (1.0 - RH_crit); + + // Mass content (kg/m2) + hail(jn, jl) = q_hail * delta_p(jn, jl) / g; + } else { + hail(jn, jl) = 0.0; + } + } + } + + oops::Log::trace() << "MassContentOfHailInAtmosphereLayer_A::executeNL done" + << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfRainInAtmosphereLayer.h b/src/vader/recipes/MassContentOfRainInAtmosphereLayer.h new file mode 100644 index 0000000..6df7ace --- /dev/null +++ b/src/vader/recipes/MassContentOfRainInAtmosphereLayer.h @@ -0,0 +1,65 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class MassContentOfRainInAtmosphereLayer_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(MassContentOfRainInAtmosphereLayer_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'MassContentOfRainInAtmosphereLayer_A' defines a recipe for + * mass content of rain in atmosphere layer + * + * \details This instantiation of RecipeBase produces mass content of rain + * using a simplified diagnostic approximation based on cloud liquid water + * autoconversion, based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature, air_pressure, relative_humidity, + * air_pressure_thickness, air_pressure_at_surface + * Output units: kg m-2 + */ +class MassContentOfRainInAtmosphereLayer_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef MassContentOfRainInAtmosphereLayer_AParameters Parameters_; + + MassContentOfRainInAtmosphereLayer_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfRainInAtmosphereLayer_A.cc b/src/vader/recipes/MassContentOfRainInAtmosphereLayer_A.cc new file mode 100644 index 0000000..8ca6e0f --- /dev/null +++ b/src/vader/recipes/MassContentOfRainInAtmosphereLayer_A.cc @@ -0,0 +1,125 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/MassContentOfRainInAtmosphereLayer.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char MassContentOfRainInAtmosphereLayer_A::Name[] = + "MassContentOfRainInAtmosphereLayer_A"; +const oops::Variables MassContentOfRainInAtmosphereLayer_A::Ingredients{ + std::vector{ + "air_temperature", "air_pressure", "relative_humidity", + "air_pressure_thickness", "air_pressure_at_surface"}}; + +// Register the maker +static RecipeMaker + makerMassContentOfRainInAtmosphereLayer_A_( + MassContentOfRainInAtmosphereLayer_A::Name); + +// ------------------------------------------------------------------------------------------------- + +MassContentOfRainInAtmosphereLayer_A::MassContentOfRainInAtmosphereLayer_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "MassContentOfRainInAtmosphereLayer_A::" + << "MassContentOfRainInAtmosphereLayer_A" << std::endl; +} + +std::string MassContentOfRainInAtmosphereLayer_A::name() const { + return MassContentOfRainInAtmosphereLayer_A::Name; +} + +oops::Variable MassContentOfRainInAtmosphereLayer_A::product() const { + return oops::Variable{"mass_content_of_rain_in_atmosphere_layer"}; +} + +oops::Variables MassContentOfRainInAtmosphereLayer_A::ingredients() const { + return MassContentOfRainInAtmosphereLayer_A::Ingredients; +} + +size_t MassContentOfRainInAtmosphereLayer_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace MassContentOfRainInAtmosphereLayer_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void MassContentOfRainInAtmosphereLayer_A::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "MassContentOfRainInAtmosphereLayer_A::executeNL starting" + << std::endl; + + // Get input fields + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + auto pressure = atlas::array::make_view(afieldset.field("air_pressure")); + auto rh = atlas::array::make_view(afieldset.field("relative_humidity")); + auto delta_p = atlas::array::make_view(afieldset.field("air_pressure_thickness")); + auto ps = atlas::array::make_view(afieldset.field("air_pressure_at_surface")); + + // Create output field + auto rain_field = afieldset.field("mass_content_of_rain_in_atmosphere_layer"); + auto rain = atlas::array::make_view(rain_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Physical constants + const double g = configVariables_.getDouble("standard_gravitational_acceleration"); + + // Diagnostic scheme parameters + const double RH_crit_base = 0.85; // Base critical RH for rain formation + const double T_freeze = 273.15; // Freezing point (K) + const double T_warm = 283.15; // Warm rain threshold (10C) + + // Calculate mass content for each grid point and level + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + const double RH = rh(jn, jl); + + // Pressure-dependent critical RH + const double P = pressure(jn, jl); + const double p_ratio = P / ps(jn, 0); + const double RH_crit = RH_crit_base + 0.10 * p_ratio; + + // Rain forms from liquid clouds above freezing with high RH + if (T > T_freeze && RH > RH_crit) { + // Excess relative humidity + const double RH_excess = std::max(0.0, RH - RH_crit); + + // Rain water content from autoconversion (kg/kg) + // More rain in warmer, more saturated conditions + const double temp_factor = std::min(1.0, (T - T_freeze) / (T_warm - T_freeze)); + const double q_rain = 5e-5 * RH_excess / (1.0 - RH_crit) * temp_factor; + + // Mass content (kg/m2) + rain(jn, jl) = q_rain * delta_p(jn, jl) / g; + } else { + rain(jn, jl) = 0.0; + } + } + } + + oops::Log::trace() << "MassContentOfRainInAtmosphereLayer_A::executeNL done" + << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfSnowInAtmosphereLayer.h b/src/vader/recipes/MassContentOfSnowInAtmosphereLayer.h new file mode 100644 index 0000000..a1e8341 --- /dev/null +++ b/src/vader/recipes/MassContentOfSnowInAtmosphereLayer.h @@ -0,0 +1,65 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include + +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/FunctionSpace.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "vader/RecipeBase.h" + +namespace vader +{ + +// ------------------------------------------------------------------------------------------------- + +class MassContentOfSnowInAtmosphereLayer_AParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(MassContentOfSnowInAtmosphereLayer_AParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{"recipe name", this}; +}; + +/*! \brief The class 'MassContentOfSnowInAtmosphereLayer_A' defines a recipe for + * mass content of snow in atmosphere layer + * + * \details This instantiation of RecipeBase produces mass content of snow + * using a simplified diagnostic approximation based on ice cloud conditions, + * based on Martin et al. (1994, J. Atmos. Sci., 51, 1823-1842) + * Inputs: air_temperature, air_pressure, relative_humidity, + * air_pressure_thickness, air_pressure_at_surface + * Output units: kg m-2 + */ +class MassContentOfSnowInAtmosphereLayer_A : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef MassContentOfSnowInAtmosphereLayer_AParameters Parameters_; + + MassContentOfSnowInAtmosphereLayer_A(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + void executeNL(atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/MassContentOfSnowInAtmosphereLayer_A.cc b/src/vader/recipes/MassContentOfSnowInAtmosphereLayer_A.cc new file mode 100644 index 0000000..b205cd6 --- /dev/null +++ b/src/vader/recipes/MassContentOfSnowInAtmosphereLayer_A.cc @@ -0,0 +1,124 @@ +/* + * (C) Copyright 2025 UCAR. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/MassContentOfSnowInAtmosphereLayer.h" + +namespace vader { + +// ------------------------------------------------------------------------------------------------- +// Static attribute initialization +const char MassContentOfSnowInAtmosphereLayer_A::Name[] = + "MassContentOfSnowInAtmosphereLayer_A"; +const oops::Variables MassContentOfSnowInAtmosphereLayer_A::Ingredients{ + std::vector{ + "air_temperature", "air_pressure", "relative_humidity", + "air_pressure_thickness", "air_pressure_at_surface"}}; + +// Register the maker +static RecipeMaker + makerMassContentOfSnowInAtmosphereLayer_A_( + MassContentOfSnowInAtmosphereLayer_A::Name); + +// ------------------------------------------------------------------------------------------------- + +MassContentOfSnowInAtmosphereLayer_A::MassContentOfSnowInAtmosphereLayer_A( + const Parameters_ & params, const VaderConfigVars & configVariables): + configVariables_{configVariables} { + oops::Log::trace() << "MassContentOfSnowInAtmosphereLayer_A::" + << "MassContentOfSnowInAtmosphereLayer_A" << std::endl; +} + +std::string MassContentOfSnowInAtmosphereLayer_A::name() const { + return MassContentOfSnowInAtmosphereLayer_A::Name; +} + +oops::Variable MassContentOfSnowInAtmosphereLayer_A::product() const { + return oops::Variable{"mass_content_of_snow_in_atmosphere_layer"}; +} + +oops::Variables MassContentOfSnowInAtmosphereLayer_A::ingredients() const { + return MassContentOfSnowInAtmosphereLayer_A::Ingredients; +} + +size_t MassContentOfSnowInAtmosphereLayer_A::productLevels( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace MassContentOfSnowInAtmosphereLayer_A::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void MassContentOfSnowInAtmosphereLayer_A::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "MassContentOfSnowInAtmosphereLayer_A::executeNL starting" + << std::endl; + + // Get input fields + auto temp = atlas::array::make_view(afieldset.field("air_temperature")); + auto pressure = atlas::array::make_view(afieldset.field("air_pressure")); + auto rh = atlas::array::make_view(afieldset.field("relative_humidity")); + auto delta_p = atlas::array::make_view(afieldset.field("air_pressure_thickness")); + auto ps = atlas::array::make_view(afieldset.field("air_pressure_at_surface")); + + // Create output field + auto snow_field = afieldset.field("mass_content_of_snow_in_atmosphere_layer"); + auto snow = atlas::array::make_view(snow_field); + + const size_t npoints = temp.shape(0); + const size_t nlevels = temp.shape(1); + + // Physical constants + const double g = configVariables_.getDouble("standard_gravitational_acceleration"); + + // Diagnostic scheme parameters + const double RH_crit_base = 0.82; // Base critical RH for snow + const double T_freeze = 273.15; // Freezing point (K) + const double T_snow_max = 268.15; // Max temp for snow (-5C) + + // Calculate mass content for each grid point and level + for (size_t jn = 0; jn < npoints; ++jn) { + for (size_t jl = 0; jl < nlevels; ++jl) { + const double T = temp(jn, jl); + const double RH = rh(jn, jl); + + // Pressure-dependent critical RH + const double P = pressure(jn, jl); + const double p_ratio = P / ps(jn, 0); + const double RH_crit = RH_crit_base + 0.10 * p_ratio; + + // Snow forms from ice clouds below freezing + if (T < T_freeze && T > T_snow_max - 40.0 && RH > RH_crit) { + // Excess relative humidity + const double RH_excess = std::max(0.0, RH - RH_crit); + + // Snow content from ice aggregation (kg/kg) + const double temp_factor = std::max(0.0, (T_freeze - T) / 40.0); + const double q_snow = 8e-5 * RH_excess / (1.0 - RH_crit) * temp_factor; + + // Mass content (kg/m2) + snow(jn, jl) = q_snow * delta_p(jn, jl) / g; + } else { + snow(jn, jl) = 0.0; + } + } + } + + oops::Log::trace() << "MassContentOfSnowInAtmosphereLayer_A::executeNL done" + << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace vader diff --git a/src/vader/recipes/RelativeHumidity.h b/src/vader/recipes/RelativeHumidity.h index 6b44aab..f9451de 100644 --- a/src/vader/recipes/RelativeHumidity.h +++ b/src/vader/recipes/RelativeHumidity.h @@ -70,4 +70,51 @@ class RelativeHumidity_A : public RecipeBase { const VaderConfigVars & configVariables_; }; +// ------------------------------------------------------------------------------------------------ + +class RelativeHumidity_BParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(RelativeHumidity_BParameters, RecipeParametersBase) + + public: + oops::RequiredParameter name{ + "recipe name", + this}; +}; + +// ------------------------------------------------------------------------------------------------ +/*! \brief RelativeHumidity_B class defines a recipe for relative humidity + * + * \details This instantiation of RecipeBase produces relative humidity + * using modular ingredients. Works with any qsat source. + * Formula: rh = 100 * q / qsat + * Inputs: water_vapor_mixing_ratio_wrt_moist_air, + * water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation + * Output units: % + * Full TL/AD support + */ +class RelativeHumidity_B : public RecipeBase { + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef RelativeHumidity_BParameters Parameters_; + + RelativeHumidity_B(const Parameters_ &, const VaderConfigVars &); + + // Recipe base class overrides + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + oops::Variables trajectoryVars() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + } // namespace vader diff --git a/src/vader/recipes/RelativeHumidity_B.cc b/src/vader/recipes/RelativeHumidity_B.cc new file mode 100644 index 0000000..0b9c5f9 --- /dev/null +++ b/src/vader/recipes/RelativeHumidity_B.cc @@ -0,0 +1,192 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/RelativeHumidity.h" + +namespace vader { + +// Static attribute initialization +const char RelativeHumidity_B::Name[] = "RelativeHumidity_B"; +const oops::Variables RelativeHumidity_B::Ingredients{std::vector{ + "water_vapor_mixing_ratio_wrt_moist_air", + "water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation"}}; + +// Register the maker +static RecipeMaker makerRelativeHumidity_B_(RelativeHumidity_B::Name); + +RelativeHumidity_B::RelativeHumidity_B(const Parameters_ & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} { + oops::Log::trace() << "RelativeHumidity_B::RelativeHumidity_B" << std::endl; +} + +std::string RelativeHumidity_B::name() const { + return RelativeHumidity_B::Name; +} + +oops::Variable RelativeHumidity_B::product() const { + return oops::Variable{"relative_humidity"}; +} + +oops::Variables RelativeHumidity_B::ingredients() const { + return RelativeHumidity_B::Ingredients; +} + +size_t RelativeHumidity_B::productLevels(const atlas::FieldSet & afieldset) const { + return afieldset.field("water_vapor_mixing_ratio_wrt_moist_air").shape(1); +} + +atlas::FunctionSpace RelativeHumidity_B::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("water_vapor_mixing_ratio_wrt_moist_air").functionspace(); +} + +void RelativeHumidity_B::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "RelativeHumidity_B::executeNL Starting" << std::endl; + + // Get fields + auto qView = atlas::array::make_view( + afieldset.field("water_vapor_mixing_ratio_wrt_moist_air")); + auto qsatView = atlas::array::make_view( + afieldset.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + auto rhView = atlas::array::make_view( + afieldset.field("relative_humidity")); + + const size_t nx = qView.shape(0); + const size_t nz = qView.shape(1); + + // Simple formula: rh = 100 * q / qsat + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double q = qView(jn, jl); + const double qsat = qsatView(jn, jl); + + if (qsat > 1.0e-12) { + rhView(jn, jl) = 100.0 * q / qsat; + } else { + rhView(jn, jl) = 0.0; + } + + // Ensure non-negative (allows supersaturation > 100%) + rhView(jn, jl) = std::max(0.0, rhView(jn, jl)); + } + } + + oops::Log::trace() << "RelativeHumidity_B::executeNL Done" << std::endl; +} + +oops::Variables RelativeHumidity_B::trajectoryVars() const { + return oops::Variables{std::vector{ + "water_vapor_mixing_ratio_wrt_moist_air", + "water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation", + "relative_humidity"}}; +} + +void RelativeHumidity_B::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "RelativeHumidity_B::executeTL Starting" << std::endl; + + // Get trajectory fields + auto qView = atlas::array::make_view( + afieldsetTraj.field("water_vapor_mixing_ratio_wrt_moist_air")); + auto qsatView = atlas::array::make_view( + afieldsetTraj.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + auto rhView = atlas::array::make_view( + afieldsetTraj.field("relative_humidity")); + + // Get TL fields + auto qTLView = atlas::array::make_view( + afieldsetTL.field("water_vapor_mixing_ratio_wrt_moist_air")); + auto qsatTLView = atlas::array::make_view( + afieldsetTL.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + auto rhTLView = atlas::array::make_view( + afieldsetTL.field("relative_humidity")); + + const size_t nx = qView.shape(0); + const size_t nz = qView.shape(1); + + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double q = qView(jn, jl); + const double qsat = qsatView(jn, jl); + const double rh = rhView(jn, jl); + + const double dq = qTLView(jn, jl); + const double dqsat = qsatTLView(jn, jl); + + if (qsat > 1.0e-12 && rh >= 0.0) { + // Linearization of: rh = 100 * q / qsat + // drh/dq = 100 / qsat + // drh/dqsat = -100 * q / qsat² + const double drh_dq = 100.0 / qsat; + const double drh_dqsat = -100.0 * q / (qsat * qsat); + + rhTLView(jn, jl) = drh_dq * dq + drh_dqsat * dqsat; + } else { + rhTLView(jn, jl) = 0.0; + } + } + } + + oops::Log::trace() << "RelativeHumidity_B::executeTL Done" << std::endl; +} + +void RelativeHumidity_B::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "RelativeHumidity_B::executeAD Starting" << std::endl; + + // Get trajectory fields + auto qView = atlas::array::make_view( + afieldsetTraj.field("water_vapor_mixing_ratio_wrt_moist_air")); + auto qsatView = atlas::array::make_view( + afieldsetTraj.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + auto rhView = atlas::array::make_view( + afieldsetTraj.field("relative_humidity")); + + // Get AD fields + auto qADView = atlas::array::make_view( + afieldsetAD.field("water_vapor_mixing_ratio_wrt_moist_air")); + auto qsatADView = atlas::array::make_view( + afieldsetAD.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + auto rhADView = atlas::array::make_view( + afieldsetAD.field("relative_humidity")); + + const size_t nx = qView.shape(0); + const size_t nz = qView.shape(1); + + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double q = qView(jn, jl); + const double qsat = qsatView(jn, jl); + const double rh = rhView(jn, jl); + const double rh_ad = rhADView(jn, jl); + + if (qsat > 1.0e-12 && rh >= 0.0) { + // Adjoint of TL + const double drh_dq = 100.0 / qsat; + const double drh_dqsat = -100.0 * q / (qsat * qsat); + + qADView(jn, jl) += drh_dq * rh_ad; + qsatADView(jn, jl) += drh_dqsat * rh_ad; + } + + // Zero out rh_ad after processing + rhADView(jn, jl) = 0.0; + } + } + + oops::Log::trace() << "RelativeHumidity_B::executeAD Done" << std::endl; +} + +} // namespace vader diff --git a/src/vader/recipes/SaturationSpecificHumidity.h b/src/vader/recipes/SaturationSpecificHumidity.h index c5fdef4..21d2c3b 100644 --- a/src/vader/recipes/SaturationSpecificHumidity.h +++ b/src/vader/recipes/SaturationSpecificHumidity.h @@ -47,4 +47,48 @@ class SaturationSpecificHumidity_A : public RecipeBase // ------------------------------------------------------------------------------------------------- +class SaturationSpecificHumidity_BParameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(SaturationSpecificHumidity_BParameters, RecipeParametersBase) + public: +}; + +// ------------------------------------------------------------------------------------------------- + +/*! \brief The class 'SaturationSpecificHumidity_B' defines a recipe for qsat + * + * \details This instantiation of RecipeBase produces saturation specific humidity (qsat) + * using exact formula with TL/AD support. + * Inputs: svp, air_pressure + * Output: water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation (kg/kg) + * + * Formula: qsat = 0.622 * es / (p - 0.378 * es) + */ +class SaturationSpecificHumidity_B : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef SaturationSpecificHumidity_BParameters Parameters_; + + SaturationSpecificHumidity_B(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + oops::Variables trajectoryVars() const override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + } // namespace vader diff --git a/src/vader/recipes/SaturationSpecificHumidity_B.cc b/src/vader/recipes/SaturationSpecificHumidity_B.cc new file mode 100644 index 0000000..509b170 --- /dev/null +++ b/src/vader/recipes/SaturationSpecificHumidity_B.cc @@ -0,0 +1,211 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/SaturationSpecificHumidity.h" + +namespace vader { + +// Static attribute initialization +const char SaturationSpecificHumidity_B::Name[] = "SaturationSpecificHumidity_B"; +const oops::Variables SaturationSpecificHumidity_B::Ingredients{ + std::vector{"svp", "air_pressure"}}; + +// Register the maker +static RecipeMaker + makerSaturationSpecificHumidity_B_(SaturationSpecificHumidity_B::Name); + +SaturationSpecificHumidity_B::SaturationSpecificHumidity_B(const Parameters_ & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} { + oops::Log::trace() << "SaturationSpecificHumidity_B::SaturationSpecificHumidity_B" + << std::endl; +} + +std::string SaturationSpecificHumidity_B::name() const { + return SaturationSpecificHumidity_B::Name; +} + +oops::Variable SaturationSpecificHumidity_B::product() const { + return oops::Variable{"water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation"}; +} + +oops::Variables SaturationSpecificHumidity_B::ingredients() const { + return Ingredients; +} + +size_t SaturationSpecificHumidity_B::productLevels(const atlas::FieldSet & afieldset) const { + return afieldset.field("air_pressure").shape(1); +} + +atlas::FunctionSpace SaturationSpecificHumidity_B::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_pressure").functionspace(); +} + +void SaturationSpecificHumidity_B::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "SaturationSpecificHumidity_B::executeNL Starting" << std::endl; + + // Get fields + auto svpView = atlas::array::make_view( + afieldset.field("svp")); + auto pressureView = atlas::array::make_view( + afieldset.field("air_pressure")); + auto qsatView = atlas::array::make_view( + afieldset.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + + const size_t nx = pressureView.shape(0); + const size_t nz = pressureView.shape(1); + + // Exact formula: qsat = eps * es / (p - c * es) + const double eps = configVariables_.getDouble("epsilon"); + const double c = 1.0 - eps; + + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double es = svpView(jn, jl); // Pa + const double p = pressureView(jn, jl); // Pa + + // Compute qsat with protection against division by zero + const double denom = p - c * es; + if (denom > 1.0) { // Ensure denominator is positive and reasonable + qsatView(jn, jl) = eps * es / denom; + } else { + qsatView(jn, jl) = 0.0; + } + + // Ensure non-negative and bounded + qsatView(jn, jl) = std::max(0.0, std::min(1.0, qsatView(jn, jl))); + } + } + + oops::Log::trace() << "SaturationSpecificHumidity_B::executeNL Done" << std::endl; +} + +oops::Variables SaturationSpecificHumidity_B::trajectoryVars() const { + return oops::Variables{std::vector{"svp", "air_pressure", + "water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation"}}; +} + +void SaturationSpecificHumidity_B::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "SaturationSpecificHumidity_B::executeTL Starting" << std::endl; + + // Get trajectory fields + auto svpView = atlas::array::make_view( + afieldsetTraj.field("svp")); + auto pressureView = atlas::array::make_view( + afieldsetTraj.field("air_pressure")); + auto qsatView = atlas::array::make_view( + afieldsetTraj.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + + // Get TL fields + auto svpTLView = atlas::array::make_view( + afieldsetTL.field("svp")); + auto pressureTLView = atlas::array::make_view( + afieldsetTL.field("air_pressure")); + auto qsatTLView = atlas::array::make_view( + afieldsetTL.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + + const size_t nx = pressureView.shape(0); + const size_t nz = pressureView.shape(1); + + const double eps = configVariables_.getDouble("epsilon"); + const double c = 1.0 - eps; + + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double es = svpView(jn, jl); + const double p = pressureView(jn, jl); + const double qsat = qsatView(jn, jl); + + const double des = svpTLView(jn, jl); + const double dp = pressureTLView(jn, jl); + + const double denom = p - c * es; + + if (denom > 1.0 && qsat > 0.0 && qsat < 1.0) { + // Linearization of: qsat = eps * es / (p - c * es) + // dqsat/des = eps / (p - c*es) + eps*es*c / (p - c*es)² + // = eps * [1 + c*es/(p - c*es)] / (p - c*es) + // = eps * p / (p - c*es)² + // dqsat/dp = -eps * es / (p - c*es)² + + const double denom2 = denom * denom; + const double dqsat_des = eps * p / denom2; + const double dqsat_dp = -eps * es / denom2; + + qsatTLView(jn, jl) = dqsat_des * des + dqsat_dp * dp; + } else { + qsatTLView(jn, jl) = 0.0; + } + } + } + + oops::Log::trace() << "SaturationSpecificHumidity_B::executeTL Done" << std::endl; +} + +void SaturationSpecificHumidity_B::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "SaturationSpecificHumidity_B::executeAD Starting" << std::endl; + + // Get trajectory fields + auto svpView = atlas::array::make_view( + afieldsetTraj.field("svp")); + auto pressureView = atlas::array::make_view( + afieldsetTraj.field("air_pressure")); + auto qsatView = atlas::array::make_view( + afieldsetTraj.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + + // Get AD fields + auto svpADView = atlas::array::make_view( + afieldsetAD.field("svp")); + auto pressureADView = atlas::array::make_view( + afieldsetAD.field("air_pressure")); + auto qsatADView = atlas::array::make_view( + afieldsetAD.field("water_vapor_mixing_ratio_wrt_moist_air_assuming_saturation")); + + const size_t nx = pressureView.shape(0); + const size_t nz = pressureView.shape(1); + + const double eps = configVariables_.getDouble("epsilon"); + const double c = 1.0 - eps; + + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double es = svpView(jn, jl); + const double p = pressureView(jn, jl); + const double qsat = qsatView(jn, jl); + const double qsat_ad = qsatADView(jn, jl); + + const double denom = p - c * es; + + if (denom > 1.0 && qsat > 0.0 && qsat < 1.0) { + // Adjoint of TL + const double denom2 = denom * denom; + const double dqsat_des = eps * p / denom2; + const double dqsat_dp = -eps * es / denom2; + + svpADView(jn, jl) += dqsat_des * qsat_ad; + pressureADView(jn, jl) += dqsat_dp * qsat_ad; + } + + // Zero out qsat_ad after processing + qsatADView(jn, jl) = 0.0; + } + } + + oops::Log::trace() << "SaturationSpecificHumidity_B::executeAD Done" << std::endl; +} + +} // namespace vader diff --git a/src/vader/recipes/SaturationVaporPressure.h b/src/vader/recipes/SaturationVaporPressure.h index 3f91780..a036f56 100644 --- a/src/vader/recipes/SaturationVaporPressure.h +++ b/src/vader/recipes/SaturationVaporPressure.h @@ -47,4 +47,49 @@ class SaturationVaporPressure_A : public RecipeBase // ------------------------------------------------------------------------------------------------- +class SaturationVaporPressure_B_Parameters : public RecipeParametersBase { + OOPS_CONCRETE_PARAMETERS(SaturationVaporPressure_B_Parameters, RecipeParametersBase) + public: +}; + +// ------------------------------------------------------------------------------------------------- + +/*! \brief The class 'SaturationVaporPressure_B' defines a recipe for 'saturation vapor pressure' + * + * \details This instantiation of RecipeBase produces saturation vapor pressure (svp) + * using Murphy & Koop (2005) formula with TL/AD support. + * Handles both liquid water (T >= 273.15K) and ice (T < 273.15K). + * Input: air_temperature + * Output: svp (Pa) + * + * Reference: Murphy & Koop (2005), valid 110-332K with ~0.01% accuracy + */ +class SaturationVaporPressure_B : public RecipeBase +{ + public: + static const char Name[]; + static const oops::Variables Ingredients; + + typedef SaturationVaporPressure_B_Parameters Parameters_; + + SaturationVaporPressure_B(const Parameters_ &, const VaderConfigVars &); + + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + oops::Variables trajectoryVars() const override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; + + private: + const VaderConfigVars & configVariables_; +}; + +// ------------------------------------------------------------------------------------------------- + } // namespace vader diff --git a/src/vader/recipes/SaturationVaporPressure_B.cc b/src/vader/recipes/SaturationVaporPressure_B.cc new file mode 100644 index 0000000..6bb739d --- /dev/null +++ b/src/vader/recipes/SaturationVaporPressure_B.cc @@ -0,0 +1,249 @@ +/* + * (C) Copyright 2025 UCAR + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "oops/util/Logger.h" +#include "vader/recipes/SaturationVaporPressure.h" + +namespace vader { + +// Static attribute initialization +const char SaturationVaporPressure_B::Name[] = "SaturationVaporPressure_B"; +const oops::Variables SaturationVaporPressure_B::Ingredients{ + std::vector{"air_temperature"}}; + +// Register the maker +static RecipeMaker + makerSaturationVaporPressure_B_(SaturationVaporPressure_B::Name); + +SaturationVaporPressure_B::SaturationVaporPressure_B(const Parameters_ & params, + const VaderConfigVars & configVariables) : + configVariables_{configVariables} { + oops::Log::trace() << "SaturationVaporPressure_B::SaturationVaporPressure_B" << std::endl; +} + +std::string SaturationVaporPressure_B::name() const { + return SaturationVaporPressure_B::Name; +} + +oops::Variable SaturationVaporPressure_B::product() const { + return oops::Variable{"svp"}; +} + +oops::Variables SaturationVaporPressure_B::ingredients() const { + return Ingredients; +} + +size_t SaturationVaporPressure_B::productLevels(const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").shape(1); +} + +atlas::FunctionSpace SaturationVaporPressure_B::productFunctionSpace( + const atlas::FieldSet & afieldset) const { + return afieldset.field("air_temperature").functionspace(); +} + +void SaturationVaporPressure_B::executeNL(atlas::FieldSet & afieldset) { + oops::Log::trace() << "SaturationVaporPressure_B::executeNL Starting" << std::endl; + + // Get fields + auto temperatureView = atlas::array::make_view( + afieldset.field("air_temperature")); + auto svpView = atlas::array::make_view( + afieldset.field("svp")); + + // Murphy & Koop (2005) coefficients for saturation vapor pressure + // Reference: Murphy & Koop (2005), Q. J. R. Meteorol. Soc., 131, 1539-1565 + // Formula: ln(es) = A/T + B + C*T + D*T² + E*T³ + F*ln(T) + + // Liquid water coefficients (valid 123K to 332K, accuracy ~0.01%) + const double A_liq = -6096.9385; + const double B_liq = 21.2409642; + const double C_liq = -2.711193e-2; + const double D_liq = 1.673952e-5; + const double E_liq = 2.433502e-8; + const double F_liq = 0.0; + + // Ice coefficients (valid 110K to 273.16K, accuracy ~0.01%) + const double A_ice = -6024.5282; + const double B_ice = 29.32707; + const double C_ice = 1.0613868e-2; + const double D_ice = -1.3198825e-5; + const double E_ice = 0.0; + const double F_ice = -0.49382577; + + // Transition temperature (273.15K = 0°C) + const double T_transition = 273.15; + + const size_t nx = temperatureView.shape(0); + const size_t nz = temperatureView.shape(1); + + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double T = temperatureView(jn, jl); // Kelvin + + // Select coefficients based on temperature + double A, B, C, D, E, F; + if (T >= T_transition) { + // Use liquid water formulation above freezing + A = A_liq; B = B_liq; C = C_liq; D = D_liq; E = E_liq; F = F_liq; + } else { + // Use ice formulation below freezing + A = A_ice; B = B_ice; C = C_ice; D = D_ice; E = E_ice; F = F_ice; + } + + // Murphy & Koop (2005) formula + const double lnT = (F != 0.0) ? std::log(T) : 0.0; + const double lnes = A/T + B + C*T + D*T*T + E*T*T*T + F*lnT; + const double es = std::exp(lnes); // Pa + + // Ensure non-negative + svpView(jn, jl) = std::max(0.0, es); + } + } + + oops::Log::trace() << "SaturationVaporPressure_B::executeNL Done" << std::endl; +} + +oops::Variables SaturationVaporPressure_B::trajectoryVars() const { + return oops::Variables{std::vector{"air_temperature", "svp"}}; +} + +void SaturationVaporPressure_B::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "SaturationVaporPressure_B::executeTL Starting" << std::endl; + + // Get trajectory fields + auto temperatureView = atlas::array::make_view( + afieldsetTraj.field("air_temperature")); + auto svpView = atlas::array::make_view( + afieldsetTraj.field("svp")); + + // Get TL fields + auto temperatureTLView = atlas::array::make_view( + afieldsetTL.field("air_temperature")); + auto svpTLView = atlas::array::make_view( + afieldsetTL.field("svp")); + + // Murphy & Koop (2005) coefficients + // Liquid water coefficients + const double A_liq = -6096.9385; + const double C_liq = -2.711193e-2; + const double D_liq = 1.673952e-5; + const double E_liq = 2.433502e-8; + const double F_liq = 0.0; + + // Ice coefficients + const double A_ice = -6024.5282; + const double C_ice = 1.0613868e-2; + const double D_ice = -1.3198825e-5; + const double E_ice = 0.0; + const double F_ice = -0.49382577; + + // Transition temperature + const double T_transition = 273.15; + + const size_t nx = temperatureView.shape(0); + const size_t nz = temperatureView.shape(1); + + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double T = temperatureView(jn, jl); + const double es = svpView(jn, jl); + const double dT = temperatureTLView(jn, jl); + + // Select coefficients based on temperature + double A, C, D, E, F; + if (T >= T_transition) { + A = A_liq; C = C_liq; D = D_liq; E = E_liq; F = F_liq; + } else { + A = A_ice; C = C_ice; D = D_ice; E = E_ice; F = F_ice; + } + + // Derivative: d(es)/dT = es * d(ln(es))/dT + // d(ln(es))/dT = -A/T² + C + 2*D*T + 3*E*T² + F/T + const double dlnes_dT = -A/(T*T) + C + 2.0*D*T + 3.0*E*T*T + F/T; + const double des_dT = es * dlnes_dT; + + svpTLView(jn, jl) = des_dT * dT; + } + } + + oops::Log::trace() << "SaturationVaporPressure_B::executeTL Done" << std::endl; +} + +void SaturationVaporPressure_B::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & afieldsetTraj) { + oops::Log::trace() << "SaturationVaporPressure_B::executeAD Starting" << std::endl; + + // Get trajectory fields + auto temperatureView = atlas::array::make_view( + afieldsetTraj.field("air_temperature")); + auto svpView = atlas::array::make_view( + afieldsetTraj.field("svp")); + + // Get AD fields + auto temperatureADView = atlas::array::make_view( + afieldsetAD.field("air_temperature")); + auto svpADView = atlas::array::make_view( + afieldsetAD.field("svp")); + + // Murphy & Koop (2005) coefficients + // Liquid water coefficients + const double A_liq = -6096.9385; + const double C_liq = -2.711193e-2; + const double D_liq = 1.673952e-5; + const double E_liq = 2.433502e-8; + const double F_liq = 0.0; + + // Ice coefficients + const double A_ice = -6024.5282; + const double C_ice = 1.0613868e-2; + const double D_ice = -1.3198825e-5; + const double E_ice = 0.0; + const double F_ice = -0.49382577; + + // Transition temperature + const double T_transition = 273.15; + + const size_t nx = temperatureView.shape(0); + const size_t nz = temperatureView.shape(1); + + for (size_t jn = 0; jn < nx; ++jn) { + for (size_t jl = 0; jl < nz; ++jl) { + const double T = temperatureView(jn, jl); + const double es = svpView(jn, jl); + const double es_ad = svpADView(jn, jl); + + // Select coefficients based on temperature + double A, C, D, E, F; + if (T >= T_transition) { + A = A_liq; C = C_liq; D = D_liq; E = E_liq; F = F_liq; + } else { + A = A_ice; C = C_ice; D = D_ice; E = E_ice; F = F_ice; + } + + // Adjoint of: des = des_dT * dT + const double dlnes_dT = -A/(T*T) + C + 2.0*D*T + 3.0*E*T*T + F/T; + const double des_dT = es * dlnes_dT; + + temperatureADView(jn, jl) += des_dT * es_ad; + + // Zero out svp_ad after processing + svpADView(jn, jl) = 0.0; + } + } + + oops::Log::trace() << "SaturationVaporPressure_B::executeAD Done" << std::endl; +} + +} // namespace vader diff --git a/src/vader/recipes/SurfaceAirPressure.h b/src/vader/recipes/SurfaceAirPressure.h index efc2423..7a8254b 100644 --- a/src/vader/recipes/SurfaceAirPressure.h +++ b/src/vader/recipes/SurfaceAirPressure.h @@ -24,37 +24,37 @@ class SurfaceAirPressure_AParameters : public RecipeParametersBase { OOPS_CONCRETE_PARAMETERS(SurfaceAirPressure_AParameters, RecipeParametersBase) public: - oops::RequiredParameter name{"recipe name", this}; + oops::RequiredParameter name{"recipe name", this}; }; -// ------------------------------------------------------------------------------------------------- -/// Recipe base class - -/*! \brief SurfaceAirPressure_A class defines a recipe for surface pressure +/*! \brief SurfaceAirPressure_A computes surface pressure from pressure thickness layers * - * \details This recipe produces surface pressure from air pressure thickness (delp) by summing - * the pressure at the model top with all the delp values. It does not provide - * TL/AD algorithms. + * Formula: ps = ptop + sum(delp) + * where delp is air_pressure_thickness (Pa), ps is air_pressure_at_surface (Pa), + * and ptop is air_pressure_at_top_of_atmosphere_model (Pa). + * Full TL/AD support. */ class SurfaceAirPressure_A : public RecipeBase { public: - static const char Name[]; - static const oops::Variables Ingredients; + static const char Name[]; + static const oops::Variables Ingredients; - typedef SurfaceAirPressure_AParameters Parameters_; + typedef SurfaceAirPressure_AParameters Parameters_; - SurfaceAirPressure_A(const Parameters_ &, const VaderConfigVars &); + SurfaceAirPressure_A(const Parameters_ &, const VaderConfigVars &); - std::string name() const override; - oops::Variable product() const override; - oops::Variables ingredients() const override; - size_t productLevels(const atlas::FieldSet &) const override; - atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; - void executeNL(atlas::FieldSet &) override; + std::string name() const override; + oops::Variable product() const override; + oops::Variables ingredients() const override; + size_t productLevels(const atlas::FieldSet &) const override; + atlas::FunctionSpace productFunctionSpace(const atlas::FieldSet &) const override; + bool hasTLAD() const override { return true; } + void executeNL(atlas::FieldSet &) override; + void executeTL(atlas::FieldSet &, const atlas::FieldSet &) override; + void executeAD(atlas::FieldSet &, const atlas::FieldSet &) override; private: - std::map p0Defaults_; - const VaderConfigVars & configVariables_; + const VaderConfigVars & configVariables_; }; } // namespace vader diff --git a/src/vader/recipes/SurfaceAirPressure_A.cc b/src/vader/recipes/SurfaceAirPressure_A.cc index b2f7beb..e0c80fb 100644 --- a/src/vader/recipes/SurfaceAirPressure_A.cc +++ b/src/vader/recipes/SurfaceAirPressure_A.cc @@ -1,11 +1,11 @@ /* - * (C) Copyright 2021-2022 UCAR. + * (C) Copyright 2025 UCAR * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. */ -#include +#include #include #include @@ -21,7 +21,7 @@ namespace vader // Static attribute initialization const char SurfaceAirPressure_A::Name[] = "SurfaceAirPressure_A"; const oops::Variables SurfaceAirPressure_A::Ingredients{ - std::vector{"air_pressure_thickness"}}; + std::vector{"air_pressure_thickness"}}; // Register the maker static RecipeMaker makerSurfaceAirPressure_A_(SurfaceAirPressure_A::Name); @@ -29,86 +29,147 @@ static RecipeMaker makerSurfaceAirPressure_A_(SurfaceAirPr // ------------------------------------------------------------------------------------------------- SurfaceAirPressure_A::SurfaceAirPressure_A(const SurfaceAirPressure_AParameters & params, - const VaderConfigVars & configVariables) : - configVariables_{configVariables} { - oops::Log::trace() << "SurfaceAirPressure_A::SurfaceAirPressure_A Starting" << std::endl; - oops::Log::trace() << "SurfaceAirPressure_A::SurfaceAirPressure_A Done" << std::endl; + const VaderConfigVars & configVariables) : + configVariables_{configVariables} { + oops::Log::trace() << "SurfaceAirPressure_A::SurfaceAirPressure_A Starting" << std::endl; + oops::Log::trace() << "SurfaceAirPressure_A::SurfaceAirPressure_A Done" << std::endl; } // ------------------------------------------------------------------------------------------------- std::string SurfaceAirPressure_A::name() const { - return SurfaceAirPressure_A::Name; + return SurfaceAirPressure_A::Name; } // ------------------------------------------------------------------------------------------------- oops::Variable SurfaceAirPressure_A::product() const { - return oops::Variable{"air_pressure_at_surface"}; + return oops::Variable{"air_pressure_at_surface"}; } // ------------------------------------------------------------------------------------------------- oops::Variables SurfaceAirPressure_A::ingredients() const { - return SurfaceAirPressure_A::Ingredients; + return SurfaceAirPressure_A::Ingredients; } // ------------------------------------------------------------------------------------------------- size_t SurfaceAirPressure_A::productLevels(const atlas::FieldSet & afieldset) const { - return 1; + return 1; } // ------------------------------------------------------------------------------------------------- atlas::FunctionSpace SurfaceAirPressure_A::productFunctionSpace(const atlas::FieldSet & afieldset) const { - return afieldset.field("air_pressure_thickness").functionspace(); + return afieldset.field("air_pressure_thickness").functionspace(); } // ------------------------------------------------------------------------------------------------- void SurfaceAirPressure_A::executeNL(atlas::FieldSet & afieldset) { - // - oops::Log::trace() << "SurfaceAirPressure_A::executeNL Starting" << std::endl; - - // Get the fields - atlas::Field delp = afieldset.field("air_pressure_thickness"); - atlas::Field ps = afieldset.field("air_pressure_at_surface"); - - // Get the units - std::string delp_units, ps_units, prsi_units; - delp.metadata().get("units", delp_units); - ps.metadata().get("units", ps_units); - - // Assert that the units match - ASSERT_MSG(ps_units == delp_units, "In Vader::SurfaceAirPressure_A::executeNL the units for " - "surface pressure " + ps_units + "do not match the pressure thickness units" - + delp_units); - - // Get ptop - double ptop = configVariables_.getDouble("air_pressure_at_top_of_atmosphere_model"); - - // Set the array views to manipulate the data - auto delp_view = atlas::array::make_view(delp); - auto ps_view = atlas::array::make_view(ps); - - // Get the grid size - const int gridSize = delp.shape(0); - const int nLevel = delp.shape(1); - - // Set pressure at the surface to surface pressure - for ( size_t jNode = 0; jNode < gridSize ; ++jNode ) { - ps_view(jNode, 0) = ptop; - } - - // Compute pressure from pressure thickness starting at the surface - for (int level = 0; level < delp.shape(1); ++level) { - for ( size_t jNode = 0; jNode < gridSize ; ++jNode ) { - ps_view(jNode, 0) = ps_view(jNode, 0) + delp_view(jNode, level); - } - } - oops::Log::trace() << "SurfaceAirPressure_A::executeNL Done" << std::endl; + oops::Log::trace() << "SurfaceAirPressure_A::executeNL Starting" << std::endl; + + // Get the fields + atlas::Field delp = afieldset.field("air_pressure_thickness"); + atlas::Field ps = afieldset.field("air_pressure_at_surface"); + + // Get the units + std::string delp_units, ps_units; + delp.metadata().get("units", delp_units); + ps.metadata().get("units", ps_units); + + // Assert that the units match + ASSERT_MSG(ps_units == delp_units, "In Vader::SurfaceAirPressure_A::executeNL the units for " + "surface pressure " + ps_units + " do not match the pressure thickness units " + + delp_units); + + // Get ptop + const double ptop = configVariables_.getDouble("air_pressure_at_top_of_atmosphere_model"); + + // Set the array views to manipulate the data + auto delp_view = atlas::array::make_view(delp); + auto ps_view = atlas::array::make_view(ps); + + // Get the grid size + const size_t npoints = delp.shape(0); + const size_t nlevels = delp.shape(1); + + // Set pressure at the surface to ptop initially + for (size_t jn = 0; jn < gridSize; ++jn) { + ps_view(jn, 0) = ptop; + } + + // Compute pressure from pressure thickness starting at the top + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < gridSize; ++jn) { + ps_view(jn, 0) += delp_view(jn, jl); + } + } + oops::Log::trace() << "SurfaceAirPressure_A::executeNL Done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +void SurfaceAirPressure_A::executeTL(atlas::FieldSet & afieldsetTL, + const atlas::FieldSet & /*afieldsetTraj*/) { + oops::Log::trace() << "SurfaceAirPressure_A::executeTL Starting" << std::endl; + + // TL fields + atlas::Field delp_tl = afieldsetTL.field("air_pressure_thickness"); + atlas::Field ps_tl = afieldsetTL.field("air_pressure_at_surface"); + + auto delp_tl_view = atlas::array::make_view(delp_tl); + auto ps_tl_view = atlas::array::make_view(ps_tl); + + const size_t npoints = delp_tl.shape(0); + const size_t nlevels = delp_tl.shape(1); + + // ptop is fixed: start with zero perturbation for surface pressure' + for (size_t jn = 0; jn < gridSize; ++jn) { + ps_tl_view(jn, 0) = 0.0; + } + + // Accumulate delp' into surface pressure' + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < gridSize; ++jn) { + ps_tl_view(jn, 0) += delp_tl_view(jn, jl); + } + } + + oops::Log::trace() << "SurfaceAirPressure_A::executeTL Done" << std::endl; +} + +// ------------------------------------------------------------------------------------------------- + +void SurfaceAirPressure_A::executeAD(atlas::FieldSet & afieldsetAD, + const atlas::FieldSet & /*afieldsetTraj*/) { + oops::Log::trace() << "SurfaceAirPressure_A::executeAD Starting" << std::endl; + + // AD fields + atlas::Field delp_ad = afieldsetAD.field("air_pressure_thickness"); + atlas::Field ps_ad = afieldsetAD.field("air_pressure_at_surface"); + + auto delp_ad_view = atlas::array::make_view(delp_ad); + auto ps_ad_view = atlas::array::make_view(ps_ad); + + const size_t npoints = delp_ad.shape(0); + const size_t nlevels = delp_ad.shape(1); + + // Distribute surface-pressure adjoint into each layer adjoint + for (size_t jl = 0; jl < nlevels; ++jl) { + for (size_t jn = 0; jn < gridSize; ++jn) { + delp_ad_view(jn, jl) += ps_ad_view(jn, 0); + } + } + + // Zero out ps_ad after processing + for (size_t jn = 0; jn < gridSize; ++jn) { + ps_ad_view(jn, 0) = 0.0; + } + + oops::Log::trace() << "SurfaceAirPressure_A::executeAD Done" << std::endl; } } // namespace vader diff --git a/src/vader/recipes/TestRecipes.h b/src/vader/recipes/TestRecipes.h index e132cec..c42182c 100644 --- a/src/vader/recipes/TestRecipes.h +++ b/src/vader/recipes/TestRecipes.h @@ -135,7 +135,7 @@ class Test_VarB_from_E : public RecipeBase { * * \details This instantiation of RecipeBase produces TestVarB * using TestVarA as input. Non-linear recipe only. - * This recipe combines with recipe Test_VarA_from_B to make + * This recipe combines with recipe Test_VarA_from_B to make * sure the Vader algorithm does not end up in an infinte * recursion loop when attempting to create ingredients. */