diff --git a/parm/aero/berror/aero_diagb.yaml.j2 b/parm/aero/berror/aero_diagb.yaml.j2 new file mode 100644 index 000000000..b39eae148 --- /dev/null +++ b/parm/aero/berror/aero_diagb.yaml.j2 @@ -0,0 +1,58 @@ +geometry: + fms initialization: + namelist filename: ./fv3jedi/fmsmpp.nml + field table filename: ./fv3jedi/field_table + akbk: ./fv3jedi/akbk.nc4 + layout: + - {{ layout_x }} + - {{ layout_y }} + npx: {{ npx_ges }} + npy: {{ npy_ges }} + npz: {{ npz_ges }} + field metadata override: ./fv3jedi/fv3jedi_fieldmetadata_restart.yaml +background: + datetime: '{{ current_cycle | to_isotime }}' + filetype: fms restart + state variables: [mass_fraction_of_sulfate_in_air, + mass_fraction_of_hydrophobic_black_carbon_in_air, + mass_fraction_of_hydrophilic_black_carbon_in_air, + mass_fraction_of_hydrophobic_organic_carbon_in_air, + mass_fraction_of_hydrophilic_organic_carbon_in_air, + mass_fraction_of_dust001_in_air, mass_fraction_of_dust002_in_air, + mass_fraction_of_dust003_in_air, mass_fraction_of_dust004_in_air, + mass_fraction_of_dust005_in_air, mass_fraction_of_sea_salt001_in_air, + mass_fraction_of_sea_salt002_in_air, mass_fraction_of_sea_salt003_in_air, + mass_fraction_of_sea_salt004_in_air] + datapath: ./bkg + filename_core: '{{ current_cycle | to_fv3time }}.fv_core.res.nc' + filename_trcr: '{{ current_cycle | to_fv3time }}.fv_tracer.res.nc' + filename_cplr: '{{ current_cycle | to_fv3time }}.coupler.res' + +background error: + datadir: ./ + date: '{{ current_cycle | to_isotime }}' + exp: bkgerr_stddev + type: incr + +variables: + - mass_fraction_of_sulfate_in_air + - mass_fraction_of_hydrophobic_black_carbon_in_air + - mass_fraction_of_hydrophilic_black_carbon_in_air + - mass_fraction_of_hydrophobic_organic_carbon_in_air + - mass_fraction_of_hydrophilic_organic_carbon_in_air + - mass_fraction_of_dust001_in_air + - mass_fraction_of_dust002_in_air + - mass_fraction_of_dust003_in_air + - mass_fraction_of_dust004_in_air + - mass_fraction_of_dust005_in_air + - mass_fraction_of_sea_salt001_in_air + - mass_fraction_of_sea_salt002_in_air + - mass_fraction_of_sea_salt003_in_air + - mass_fraction_of_sea_salt004_in_air + +rescale: 2.0 # rescales the filtered std. dev. by "rescale" +number of halo points: 4 +number of neighbors: 16 +simple smoothing: + horizontal iterations: 2 + vertical iterations: 1 diff --git a/parm/aero/berror/staticb_bump.yaml.j2 b/parm/aero/berror/staticb_bump.yaml.j2 index 4ce32187d..34b19ad06 100644 --- a/parm/aero/berror/staticb_bump.yaml.j2 +++ b/parm/aero/berror/staticb_bump.yaml.j2 @@ -14,12 +14,12 @@ saber outer blocks: - saber block name: StdDev read: model file: - datetime: '2021-02-27T00:00:00Z' + datetime: '{{ current_cycle | to_isotime }}' set datetime on read: true filetype: fms restart psinfile: true datapath: "{{ DATA }}/berror" - filename_core: 20210227.000000.stddev.fv_core.res.nc - filename_trcr: 20210227.000000.stddev.fv_tracer.res.nc - filename_cplr: 20210227.000000.stddev.coupler.res - date: '2021-02-27T18:00:00Z' + filename_core: '{{ current_cycle | to_fv3time }}.stddev.fv_core.res.nc' + filename_trcr: '{{ current_cycle | to_fv3time }}.stddev.fv_tracer.res.nc' + filename_cplr: '{{ current_cycle | to_fv3time }}.stddev.coupler.res' + date: '{{ current_cycle | to_isotime }}' diff --git a/utils/CMakeLists.txt b/utils/CMakeLists.txt index 0195c15dc..4f376c26a 100644 --- a/utils/CMakeLists.txt +++ b/utils/CMakeLists.txt @@ -21,3 +21,4 @@ add_subdirectory(ioda_example) add_subdirectory(test) add_subdirectory(obsproc) add_subdirectory(fv3jedi) +add_subdirectory(chem) diff --git a/utils/chem/CMakeLists.txt b/utils/chem/CMakeLists.txt new file mode 100644 index 000000000..0a0b064b8 --- /dev/null +++ b/utils/chem/CMakeLists.txt @@ -0,0 +1,5 @@ +# Diag B +ecbuild_add_executable( TARGET gdasapp_chem_diagb.x + SOURCES chem_diagb.cc ) +target_compile_features( gdasapp_chem_diagb.x PUBLIC cxx_std_17) +target_link_libraries( gdasapp_chem_diagb.x PUBLIC NetCDF::NetCDF_CXX oops atlas fv3jedi) diff --git a/utils/chem/chem_diagb.cc b/utils/chem/chem_diagb.cc new file mode 100644 index 000000000..aec91d104 --- /dev/null +++ b/utils/chem/chem_diagb.cc @@ -0,0 +1,8 @@ +#include "chem_diagb.h" +#include "oops/runs/Run.h" + +int main(int argc, char ** argv) { + oops::Run run(argc, argv); + gdasapp::FV3ChemDiagB fv3chemdiagb; + return run.execute(fv3chemdiagb); +} diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h new file mode 100644 index 000000000..b8bfcdae9 --- /dev/null +++ b/utils/chem/chem_diagb.h @@ -0,0 +1,255 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "atlas/field.h" +#include "atlas/functionspace/NodeColumns.h" +#include "atlas/mesh.h" +#include "atlas/mesh/actions/BuildEdges.h" +#include "atlas/mesh/actions/BuildHalo.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/util/Earth.h" +#include "atlas/util/Geometry.h" +#include "atlas/util/Point.h" + +#include "oops/base/FieldSet3D.h" +#include "oops/base/GeometryData.h" +#include "oops/mpi/mpi.h" +#include "oops/runs/Application.h" +#include "oops/util/DateTime.h" +#include "oops/util/Duration.h" +#include "oops/util/FieldSetHelpers.h" +#include "oops/util/FieldSetOperations.h" +#include "oops/util/Logger.h" + +#include "fv3jedi/Geometry/Geometry.h" +#include "fv3jedi/Increment/Increment.h" +#include "fv3jedi/State/State.h" + +namespace gdasapp { + /** + * FV3ChemDiagB Class Implementation + * + * Implements variance partitioning within the GDAS for the aerosols and tracers. It is used as a proxy for the + * diagonal of B. + * This class is designed to partition the variance of aerosol and tracer fields by analyzing the variance within + * predefined geographical bins. + * This approach allows for an ensemble-free estimate of a flow-dependent background error. + */ + + class FV3ChemDiagB : public oops::Application { + public: + explicit FV3ChemDiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) + : Application(comm) {} + static const std::string classname() {return "gdasapp::FV3ChemDiagB";} + + int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { + /// Setup the fv3jedi geometry + oops::Log::info() << "====================== geometry" << std::endl; + const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); + const fv3jedi::Geometry geom(geomConfig, this->getComm()); + + oops::Log::info() << "====================== date" << std::endl; + /// Get the date + // ------------- + std::string strdt; + fullConfig.get("date", strdt); + util::DateTime cycleDate = util::DateTime(strdt); + + /// Get the list of variables + // -------------------------- + oops::Log::info() << "====================== variables" << std::endl; + oops::Variables chemVars(fullConfig, "variables.name"); + + /// Read the background + // -------------------- + oops::Log::info() << "====================== read bkg" << std::endl; + fv3jedi::State xb(geom, chemVars, cycleDate); + const eckit::LocalConfiguration bkgConfig(fullConfig, "background"); + xb.read(bkgConfig); + atlas::FieldSet xbFs; + xb.toFieldSet(xbFs); + oops::Log::info() << "Background:" << std::endl; + oops::Log::info() << xb << std::endl; + + /// Create the mesh connectivity (Copy/paste of Francois's stuff) + // -------------------------------------------------------------- + // Build edges, then connections between nodes and edges + int nbHalos(2); + fullConfig.get("number of halo points", nbHalos); + int nbNeighbors(4); + fullConfig.get("number of neighbors", nbNeighbors); + atlas::functionspace::NodeColumns nodeColumns = geom.functionSpace(); + atlas::Mesh mesh = nodeColumns.mesh(); + atlas::mesh::actions::build_edges(mesh); + atlas::mesh::actions::build_node_to_edge_connectivity(mesh); + atlas::mesh::actions::build_halo(mesh, nbHalos); + const auto & node2edge = mesh.nodes().edge_connectivity(); + const auto & edge2node = mesh.edges().node_connectivity(); + + // Lambda function to get the neighbors of a node + // (Copy/paste from Francois's un-merged oops branch) + const auto get_neighbors_of_node = [&](const int node) { + std::vector neighbors{}; + neighbors.reserve(nbNeighbors); + if (node >= mesh.nodes().size()) { + return neighbors; + } + // Use node2edge and edge2node connectivities to find neighbors of node + const int nb_edges = node2edge.cols(node); + for (size_t ie = 0; ie < nb_edges; ++ie) { + const int edge = node2edge(node, ie); + ASSERT(edge2node.cols(edge) == 2); // sanity check, maybe not in production/release build + const int node0 = edge2node(edge, 0); + const int node1 = edge2node(edge, 1); + ASSERT(node == node0 || node == node1); // sanity check edge contains initial node + if (node != node0) { + neighbors.push_back(node0); + } else { + neighbors.push_back(node1); + } + } + return neighbors; + }; + + /// Compute local std. dev. as a proxy of the bkg error + // ---------------------------------------------------- + // Identify halo nodes + const auto ghostView = + atlas::array::make_view(geom.functionSpace().ghost()); + + // Create the background error fieldset + fv3jedi::Increment bkgErr(geom, chemVars, cycleDate); + bkgErr.zero(); + atlas::FieldSet bkgErrFs; + bkgErr.toFieldSet(bkgErrFs); + + // Loop through variables + for (auto & var : chemVars.variables()) { + oops::Log::info() << "====================== std dev for " << var << std::endl; + auto bkg = atlas::array::make_view(xbFs[var]); + auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); + + // Loop through nodes + for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { + // get indices of neighbor cells + auto neighbors = get_neighbors_of_node(jnode); + int nbh = neighbors.size(); + + int nbz = 1; // Number of closest point in the vertical, above and below + for (atlas::idx_t level = 0; level <= xbFs[var].shape(1) - nbz; ++level) { + std::vector local; + for (int nn = 0; nn < neighbors.size(); ++nn) { + int levelMin = std::max(0, level - nbz); + int levelMax = level + nbz; + for (int ll = levelMin; ll <= levelMax; ++ll) { + local.push_back(bkg(neighbors[nn], ll)); + } + } + // Set the minimum number of points + int minn = 6; /// probably should be passed through the config + if (local.size() >= minn) { + // Mean + double mean = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); + + // Standard deviation + double stdDev(0.0); + for (int nn = 0; nn < nbh; ++nn) { + stdDev += std::pow(local[nn] - mean, 2.0); + } + // Setup the additive variance (only used ofr sst) + double additiveStdDev(0.0); + if (stdDev > 0.0 || local.size() > 2) { + stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)); + } + } + } // end level + } // end jnode + } // end var + + + + /// Smooth the fields + // ------------------ + if (fullConfig.has("simple smoothing")) { + int niter(0); + fullConfig.get("simple smoothing.horizontal iterations", niter); + for (auto & var : chemVars.variables()) { + // Horizontal averaging + for (int iter = 0; iter < niter; ++iter) { + // Update the halo points + nodeColumns.haloExchange(bkgErrFs[var]); + auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); + + // Loops through nodes and levels + for (atlas::idx_t level = 0; level <= xbFs[var].shape(1); ++level) { + for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { + std::vector local; + auto neighbors = get_neighbors_of_node(jnode); + int nbh = neighbors.size(); + for (int nn = 0; nn < neighbors.size(); ++nn) { + int nbNode = neighbors[nn]; + local.push_back(stdDevBkg(nbNode, level)); + } + + if (local.size() > 2) { + stdDevBkg(jnode, level) = std::accumulate(local.begin(), local.end(), 0.0) / + local.size(); + } + } + } + } + + // Vertical averaging + if (xbFs[var].shape(1) == 1) { + oops::Log::info() << "skipping vertical smoothing for " << var << std::endl; + } else { + int niterVert(0); + fullConfig.get("simple smoothing.vertical iterations", niterVert); + + auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); + auto tmpArray(stdDevBkg); + for (int iter = 0; iter < niterVert; ++iter) { + for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { + for (atlas::idx_t level = 1; level < xbFs[var].shape(1)-1; ++level) { + stdDevBkg(jnode, level) = (tmpArray(jnode, level-1) + + tmpArray(jnode, level) + + tmpArray(jnode, level+1)) / 3.0; + } + stdDevBkg(jnode, 0) = stdDevBkg(jnode, 1); + } + } + } + } + } + + // Rescale + double rescale; + fullConfig.get("rescale", rescale); + util::multiplyFieldSet(bkgErrFs, rescale); + + // We want to write with fv3jedi, not atlas: Syncronize with fv3jedi Increment + bkgErr.fromFieldSet(bkgErrFs); + + // Save the background error + const eckit::LocalConfiguration bkgErrorConfig(fullConfig, "background error"); + bkgErr.write(bkgErrorConfig); + + return 0; + } + + // ----------------------------------------------------------------------------- + private: + std::string appname() const { + return "gdasapp::FV3ChemDiagB"; + } + // ----------------------------------------------------------------------------- + }; + +} // namespace gdasapp