diff --git a/src/ufo/filters/CMakeLists.txt b/src/ufo/filters/CMakeLists.txt index 044728843..c91ebde6f 100644 --- a/src/ufo/filters/CMakeLists.txt +++ b/src/ufo/filters/CMakeLists.txt @@ -323,6 +323,8 @@ set ( filters_files obsfunctions/ObsErrorFactorSfcPressure.h obsfunctions/ObsErrorFactorConventional.cc obsfunctions/ObsErrorFactorConventional.h + obsfunctions/ObsErrorFactorSnowElevationDiff.cc + obsfunctions/ObsErrorFactorSnowElevationDiff.h obsfunctions/PotentialTemperatureFromTemperature.cc obsfunctions/PotentialTemperatureFromTemperature.h obsfunctions/NearSSTRetCheckIR.cc diff --git a/src/ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.cc b/src/ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.cc new file mode 100644 index 000000000..cc414858b --- /dev/null +++ b/src/ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.cc @@ -0,0 +1,113 @@ +/* + * (C) Copyright 2020 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 "ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.h" + +#include + +#include "eckit/exception/Exceptions.h" +#include "ioda/ObsDataVector.h" +#include "oops/util/Logger.h" +#include "oops/util/missingValues.h" +#include "ufo/filters/ObsFilterData.h" + +namespace ufo { + +static ObsFunctionMaker + makerObsErrorFactorSnowElevationDiff_("ObsErrorFactorSnowElevationDiff"); + +// ----------------------------------------------------------------------------- + +ObsErrorFactorSnowElevationDiff::ObsErrorFactorSnowElevationDiff( + const eckit::Configuration &config) + : invars_() { + oops::Log::trace() << "ObsErrorFactorSnowElevationDiff constructor" << std::endl; + oops::Log::debug() << "ObsErrorFactorSnowElevationDiff: config = " << config << std::endl; + + // Initialize options + options_.reset(new ObsErrorFactorSnowElevationDiffParameters()); + options_->deserialize(config); + + // Include observation elevation variable + const std::string obs_elev_var = options_->obs_elevation_var.value(); + invars_ += Variable(obs_elev_var); + + // Include model elevation variable + const std::string model_elev_var = options_->model_elevation_var.value(); + invars_ += Variable(model_elev_var); +} + +// ----------------------------------------------------------------------------- + +ObsErrorFactorSnowElevationDiff::~ObsErrorFactorSnowElevationDiff() { + oops::Log::trace() << "ObsErrorFactorSnowElevationDiff destructor" << std::endl; +} + +// ----------------------------------------------------------------------------- + +void ObsErrorFactorSnowElevationDiff::compute( + const ObsFilterData & data, + ioda::ObsDataVector & obserr) const { + oops::Log::trace() << "ObsErrorFactorSnowElevationDiff compute start" << std::endl; + + const float missing = util::missingValue(); + const float h_scale = options_->elevation_scale_m.value(); + + // If no observations on this processor then nothing to do + if (data.nlocs() == 0) return; + + // Ensure that only one output variable is expected + ASSERT(obserr.nvars() == 1); + + // Get dimensions + size_t nlocs = data.nlocs(); + + // Get observation elevation variable + std::vector ob_elevation(nlocs); + const std::string obs_elev_var = options_->obs_elevation_var.value(); + data.get(Variable(obs_elev_var), ob_elevation); + + // Get model surface elevation variable + std::vector model_elevation(nlocs); + const std::string model_elev_var = options_->model_elevation_var.value(); + data.get(Variable(model_elev_var), model_elevation); + + // Compute inflation factor for each observation + float dz, inflation_factor_1, inflation_factor; + int iv = 0; + + for (size_t iloc = 0; iloc < nlocs; ++iloc) { + // If missing observation or model elevation, set obserror to missing + if (ob_elevation[iloc] == missing || model_elevation[iloc] == missing) { + obserr[iv][iloc] = missing; + } else { + // Compute elevation difference + dz = std::abs(model_elevation[iloc] - ob_elevation[iloc]); + + // Compute inflation_factor_1 = exp(-1 * dz^2 / (h^2)) + inflation_factor_1 = std::exp(-1.0f * dz * dz / (h_scale * h_scale)); + + // Compute output inflation_factor = 1 / inflation_factor_1 + inflation_factor = 1.0f / inflation_factor_1; + + // Inflate error by the inflation factor + obserr[iv][iloc] *= inflation_factor; + } + } + + oops::Log::trace() << "ObsErrorFactorSnowElevationDiff compute complete" << std::endl; +} + +// ----------------------------------------------------------------------------- + +const ufo::Variables & ObsErrorFactorSnowElevationDiff::requiredVariables() const { + return invars_; +} + +// ----------------------------------------------------------------------------- + +} // namespace ufo diff --git a/src/ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.h b/src/ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.h new file mode 100644 index 000000000..e2f38bede --- /dev/null +++ b/src/ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.h @@ -0,0 +1,78 @@ +/* + * (C) Copyright 2020 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. + */ + +#ifndef UFO_FILTERS_OBSFUNCTIONS_OBSERRORFACTORSNOWELEVATIONDIFF_H_ +#define UFO_FILTERS_OBSFUNCTIONS_OBSERRORFACTORSNOWELEVATIONDIFF_H_ + +#include +#include +#include + +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/Parameters.h" +#include "oops/util/parameters/RequiredParameter.h" + +#include "ufo/filters/obsfunctions/ObsFunctionBase.h" +#include "ufo/filters/Variables.h" + +namespace ufo { + +/// \brief Options controlling ObsErrorFactorSnowElevationDiff ObsFunction +class ObsErrorFactorSnowElevationDiffParameters : public oops::Parameters { + OOPS_CONCRETE_PARAMETERS(ObsErrorFactorSnowElevationDiffParameters, Parameters) + + public: + oops::Parameter obs_elevation_var{"observation_elevation", "MetaData/stationElevation", this}; + oops::Parameter model_elevation_var{"model_elevation", "GeoVaLs/filtered_orography", this}; + oops::Parameter elevation_scale_m{"elevation_scale_m", 800., this}; +}; + +// ----------------------------------------------------------------------------- + +/// \brief Inflate observation error based on elevation difference between model and observation. +/// +/// This routine computes an observation error inflation factor based on the elevation difference +/// between the model surface elevation and the observed station elevation. +/// The inflation factor is computed as: 1 / exp(-1 * dz^2 / (h^2)) +/// where dz = |model_elevation - obs_elevation| +/// and h is the elevation_scale parameter in m (Elevation difference at which the obs-error will be inflated by a factor of exp(1)). +/// +/// Authors: Tseganeh Z Gichamo and Gihub Copilot. All the ML generated code sections have been revised and tested by Tseganeh Z. Gichamo +/// +/// ### example configurations for application of this filter: ### +/// +/// - filter: Perform Action +/// filter variables: +/// - name: totalSnowDepth +/// action: +/// name: inflate error +/// inflation variable: +/// name: ObsFunction/ObsErrorFactorSnowElevationDiff +/// options: +/// observation_elevation: MetaData/stationElevation +/// model_elevation: GeoVaLs/filtered_orography +/// elevation_scale_m: 800.0 // (Unit m) +/// +class ObsErrorFactorSnowElevationDiff : public ObsFunctionBase { + public: + static const std::string classname() {return "ObsErrorFactorSnowElevationDiff";} + + explicit ObsErrorFactorSnowElevationDiff(const eckit::Configuration &config); + ~ObsErrorFactorSnowElevationDiff(); + + void compute(const ObsFilterData &, ioda::ObsDataVector &) const; + const ufo::Variables & requiredVariables() const; + private: + ufo::Variables invars_; + std::unique_ptr options_; +}; + +// ----------------------------------------------------------------------------- + +} // namespace ufo + +#endif // UFO_FILTERS_OBSFUNCTIONS_OBSERRORFACTORSNOWELEVATIONDIFF_H_