From 86370cfd365631e1a5e8915be98205ec0bc65fd7 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Tue, 2 Apr 2024 14:06:41 -0500 Subject: [PATCH 01/14] Added code to compute diagb file using template from soca. --- utils/CMakeLists.txt | 1 + utils/chem/CMakeLists.txt | 5 + utils/chem/chem_diagb.cc | 8 + utils/chem/chem_diagb.h | 382 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 396 insertions(+) create mode 100644 utils/chem/CMakeLists.txt create mode 100644 utils/chem/chem_diagb.cc create mode 100644 utils/chem/chem_diagb.h 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..5cb698209 --- /dev/null +++ b/utils/chem/CMakeLists.txt @@ -0,0 +1,5 @@ +# Diag B +ecbuild_add_executable( TARGET chem_diagb.x + SOURCES chem_diagb.cc ) +target_compile_features( chem_diagb.x PUBLIC cxx_std_17) +target_link_libraries( 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..6c2b2b5e5 --- /dev/null +++ b/utils/chem/chem_diagb.h @@ -0,0 +1,382 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "atlas/field.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 "atlas/functionspace/NodeColumns.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" +//#include "fv3jedi/ExplicitDiffusion/ExplicitDiffusion.h" +//#include "fv3jedi/ExplicitDiffusion/ExplicitDiffusionParameters.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; + chem::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 + chem::Increment bkgErr(geom, chemVars, cycleDate); + bkgErr.zero(); + atlas::FieldSet bkgErrFs; + bkgErr.toFieldSet(bkgErrFs); + + // Get the std dev to add to sst + double sstBkgErrMin(0.0); + fullConfig.get("min sst", sstBkgErrMin); + + // Get the max std dev for ssh + double sshMax(0.0); + fullConfig.get("max ssh", sshMax); + + // Get the layer thicknesses and convert to depth + auto h = atlas::array::make_view(xbFs["hocn"]); + atlas::array::ArrayT depth(h.shape(0), h.shape(1)); + auto viewDepth = atlas::array::make_view(depth); + + for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { + viewDepth(jnode, 0) = 0.5 * h(jnode, 0); + for (atlas::idx_t level = 1; level < h.shape(1); ++level) { + viewDepth(jnode, level) = viewDepth(jnode, level-1) + + 0.5 * (h(jnode, level-1 ) + h(jnode, level)); + } + } + + // Get the mixed layer depth + auto mld = atlas::array::make_view(xbFs["mom6_mld"]); + atlas::array::ArrayT mldindex(mld.shape(0), mld.shape(1)); + auto viewMldindex = atlas::array::make_view(mldindex); + + for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { + std::vector testMld; + for (atlas::idx_t level = 0; level < viewDepth.shape(1); ++level) { + testMld.push_back(std::abs(viewDepth(jnode, level) - mld(jnode, 0))); + } + auto minVal = std::min_element(testMld.begin(), testMld.end()); + + viewMldindex(jnode, 0) = std::distance(testMld.begin(), minVal); + } + + + // Loop through variables + for (auto & var : chemVars.variables()) { + // Skip the layer thickness variable + if (var == "hocn") { + continue; + } + 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) { + // Early exit if thickness is 0 or on a ghost cell + if (ghostView(jnode) > 0 || abs(h(jnode, 0)) <= 0.1) { + continue; + } + + // get indices of neighbor cells + auto neighbors = get_neighbors_of_node(jnode); + int nbh = neighbors.size(); + + // 2D case + if (xbFs[var].shape(1) == 1) { + std::vector local; + for (int nn = 0; nn < neighbors.size(); ++nn) { + if ( abs(h(neighbors[nn], 0)) <= 0.1 ) { + continue; + } + local.push_back(bkg(neighbors[nn], 0)); + } + //Set the minimum number of points + int minn = 4; /// 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); + } + + if (stdDev > 0.0 || local.size() > 2) { + stdDevBkg(jnode, 0) = std::sqrt(stdDev / (local.size() - 1)); + } + // Filter out ssh std. dev. + if (var == "ssh") { + stdDevBkg(jnode, 0) = std::min(stdDevBkg(jnode, 0), sshMax); + } + } + continue; // early exit, no need to loop through levels + } else { + // 3D case + int nbz = 1; // Number of closest point in the vertical, above and below + int nzMld = std::max(10, viewMldindex(jnode, 0)); + 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; + if (level < nzMld) { + // If in the MLD, compute the std. dev. throughout the MLD + levelMin = 0; + levelMax = 1; //nzMld; + } + for (int ll = levelMin; ll <= levelMax; ++ll) { + if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { + continue; + } + 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 (var == "tocn" & level == 0) { + additiveStdDev = sstBkgErrMin; + } + if (stdDev > 0.0 || local.size() > 2) { + stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)) + additiveStdDev; + } + } + } // end level + } // end 3D case + } // end jnode + } // end var + + // TODO(G): Assume that the steric balance explains 97% of ssh ... or do it properly ... maybe + + + /// Smooth the fields + // ------------------ + if (fullConfig.has("simple smoothing")) { + int niter(0); + fullConfig.get("simple smoothing.horizontal iterations", niter); + for (auto & var : chemVars.variables()) { + // Skip the layer thickness variable + if (var == "hocn") { + continue; + } + + // 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["tocn"].shape(0); ++jnode) { + // Early exit if thickness is 0 or on a ghost cell + if (abs(h(jnode, level)) <= 0.1) { + continue; + } + + // Ocean or ice node, do something + 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]; + if ( abs(h(nbNode, level)) <= 0.1 ) { + continue; + } + 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["tocn"].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); + } + } + } + } + } + + /// Use explicit diffusion to smooth the background error + // ------------------------------------------------------ + // TODO: Use this once Travis adds the option to skip the normalization. + // The output is currently in [0, ~1000] + // Initialize the diffusion central block +// if (fullConfig.has("diffusion")) { +// const eckit::LocalConfiguration diffusionConfig(fullConfig, "diffusion"); +// chem::ExplicitDiffusionParameters params; +// params.deserialize(diffusionConfig); +// oops::GeometryData geometryData(geom.functionSpace(), bkgErrFs["tocn"], true, this->getComm()); +// const oops::FieldSet3D dumyXb(cycleDate, this->getComm()); +// chem::ExplicitDiffusion diffuse(geometryData, chemVars, diffusionConfig, params, dumyXb, dumyXb); +// diffuse.read(); +// +// // Smooth the field +// oops::FieldSet3D dx(cycleDate, this->getComm()); +// dx.deepCopy(bkgErrFs); +// diffuse.multiply(dx); +// bkgErrFs = dx.fieldSet(); +// } + + // 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 From 5e8315027c6df6c5896ee4edc6a2559615708170 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Fri, 5 Apr 2024 12:49:52 -0500 Subject: [PATCH 02/14] Changes to chem_diagb files --- utils/chem/chem_diagb.cc | 1 + utils/chem/chem_diagb.h | 119 +++++++++++++++++++++------------------ 2 files changed, 65 insertions(+), 55 deletions(-) diff --git a/utils/chem/chem_diagb.cc b/utils/chem/chem_diagb.cc index aec91d104..218f6820d 100644 --- a/utils/chem/chem_diagb.cc +++ b/utils/chem/chem_diagb.cc @@ -3,6 +3,7 @@ int main(int argc, char ** argv) { oops::Run run(argc, argv); + std::cout << "argv = " << argv ; gdasapp::FV3ChemDiagB fv3chemdiagb; return run.execute(fv3chemdiagb); } diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index 6c2b2b5e5..2dfe7f966 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -47,12 +47,13 @@ namespace gdasapp { class FV3ChemDiagB : public oops::Application { public: - explicit Fv3ChemDiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) + 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 + std::cout << "int execute"; oops::Log::info() << "====================== geometry" << std::endl; const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); const fv3jedi::Geometry geom(geomConfig, this->getComm()); @@ -71,8 +72,9 @@ namespace gdasapp { /// Read the background // -------------------- + std::cout << "Read the background"; oops::Log::info() << "====================== read bkg" << std::endl; - chem::State xb(geom, chemVars, cycleDate); + fv3jedi::State xb(geom, chemVars, cycleDate); const eckit::LocalConfiguration bkgConfig(fullConfig, "background"); xb.read(bkgConfig); atlas::FieldSet xbFs; @@ -83,6 +85,7 @@ namespace gdasapp { /// Create the mesh connectivity (Copy/paste of Francois's stuff) // -------------------------------------------------------------- // Build edges, then connections between nodes and edges + std::cout << "Create the mesh connectivity"; int nbHalos(2); fullConfig.get("number of halo points", nbHalos); int nbNeighbors(4); @@ -126,7 +129,7 @@ namespace gdasapp { atlas::array::make_view(geom.functionSpace().ghost()); // Create the background error fieldset - chem::Increment bkgErr(geom, chemVars, cycleDate); + fv3jedi::Increment bkgErr(geom, chemVars, cycleDate); bkgErr.zero(); atlas::FieldSet bkgErrFs; bkgErr.toFieldSet(bkgErrFs); @@ -140,32 +143,32 @@ namespace gdasapp { fullConfig.get("max ssh", sshMax); // Get the layer thicknesses and convert to depth - auto h = atlas::array::make_view(xbFs["hocn"]); - atlas::array::ArrayT depth(h.shape(0), h.shape(1)); - auto viewDepth = atlas::array::make_view(depth); - - for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { - viewDepth(jnode, 0) = 0.5 * h(jnode, 0); - for (atlas::idx_t level = 1; level < h.shape(1); ++level) { - viewDepth(jnode, level) = viewDepth(jnode, level-1) + - 0.5 * (h(jnode, level-1 ) + h(jnode, level)); - } - } +// auto h = atlas::array::make_view(xbFs["hocn"]); +// atlas::array::ArrayT depth(h.shape(0), h.shape(1)); +// auto viewDepth = atlas::array::make_view(depth); + +// for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { +// viewDepth(jnode, 0) = 0.5 * h(jnode, 0); +// for (atlas::idx_t level = 1; level < h.shape(1); ++level) { +// viewDepth(jnode, level) = viewDepth(jnode, level-1) + +// 0.5 * (h(jnode, level-1 ) + h(jnode, level)); +// } +// } // Get the mixed layer depth - auto mld = atlas::array::make_view(xbFs["mom6_mld"]); - atlas::array::ArrayT mldindex(mld.shape(0), mld.shape(1)); - auto viewMldindex = atlas::array::make_view(mldindex); - - for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { - std::vector testMld; - for (atlas::idx_t level = 0; level < viewDepth.shape(1); ++level) { - testMld.push_back(std::abs(viewDepth(jnode, level) - mld(jnode, 0))); - } - auto minVal = std::min_element(testMld.begin(), testMld.end()); - - viewMldindex(jnode, 0) = std::distance(testMld.begin(), minVal); - } +// auto mld = atlas::array::make_view(xbFs["mom6_mld"]); +// atlas::array::ArrayT mldindex(mld.shape(0), mld.shape(1)); +// auto viewMldindex = atlas::array::make_view(mldindex); + +// for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { +// std::vector testMld; +// for (atlas::idx_t level = 0; level < viewDepth.shape(1); ++level) { +// testMld.push_back(std::abs(viewDepth(jnode, level) - mld(jnode, 0))); +// } +// auto minVal = std::min_element(testMld.begin(), testMld.end()); + +// viewMldindex(jnode, 0) = std::distance(testMld.begin(), minVal); +// } // Loop through variables @@ -181,9 +184,9 @@ namespace gdasapp { // Loop through nodes for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { // Early exit if thickness is 0 or on a ghost cell - if (ghostView(jnode) > 0 || abs(h(jnode, 0)) <= 0.1) { - continue; - } +// if (ghostView(jnode) > 0 || abs(h(jnode, 0)) <= 0.1) { +// continue; +// } // get indices of neighbor cells auto neighbors = get_neighbors_of_node(jnode); @@ -193,9 +196,9 @@ namespace gdasapp { if (xbFs[var].shape(1) == 1) { std::vector local; for (int nn = 0; nn < neighbors.size(); ++nn) { - if ( abs(h(neighbors[nn], 0)) <= 0.1 ) { - continue; - } +// if ( abs(h(neighbors[nn], 0)) <= 0.1 ) { +// continue; +// } local.push_back(bkg(neighbors[nn], 0)); } //Set the minimum number of points @@ -221,21 +224,21 @@ namespace gdasapp { } else { // 3D case int nbz = 1; // Number of closest point in the vertical, above and below - int nzMld = std::max(10, viewMldindex(jnode, 0)); +// int nzMld = std::max(10, viewMldindex(jnode, 0)); 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; - if (level < nzMld) { +// if (level < nzMld) { // If in the MLD, compute the std. dev. throughout the MLD levelMin = 0; levelMax = 1; //nzMld; - } +// } for (int ll = levelMin; ll <= levelMax; ++ll) { - if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { - continue; - } +// if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { +// continue; +// } local.push_back(bkg(neighbors[nn], ll)); } } @@ -252,9 +255,9 @@ namespace gdasapp { } // Setup the additive variance (only used ofr sst) double additiveStdDev(0.0); - if (var == "tocn" & level == 0) { - additiveStdDev = sstBkgErrMin; - } +// if (var == "tocn" & level == 0) { +// additiveStdDev = sstBkgErrMin; +// } if (stdDev > 0.0 || local.size() > 2) { stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)) + additiveStdDev; } @@ -274,9 +277,9 @@ namespace gdasapp { fullConfig.get("simple smoothing.horizontal iterations", niter); for (auto & var : chemVars.variables()) { // Skip the layer thickness variable - if (var == "hocn") { - continue; - } +// if (var == "hocn") { +// continue; +// } // Horizontal averaging for (int iter = 0; iter < niter; ++iter) { @@ -287,11 +290,12 @@ namespace gdasapp { // 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["tocn"].shape(0); ++jnode) { +// for (atlas::idx_t jnode = 0; jnode < xbFs["tocn"].shape(0); ++jnode) { + for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { // Early exit if thickness is 0 or on a ghost cell - if (abs(h(jnode, level)) <= 0.1) { - continue; - } +// if (abs(h(jnode, level)) <= 0.1) { +// continue; +// } // Ocean or ice node, do something std::vector local; @@ -299,9 +303,9 @@ namespace gdasapp { int nbh = neighbors.size(); for (int nn = 0; nn < neighbors.size(); ++nn) { int nbNode = neighbors[nn]; - if ( abs(h(nbNode, level)) <= 0.1 ) { - continue; - } +// if ( abs(h(nbNode, level)) <= 0.1 ) { +// continue; +// } local.push_back(stdDevBkg(nbNode, level)); } @@ -322,7 +326,8 @@ namespace gdasapp { 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["tocn"].shape(0); ++jnode) { +// for (atlas::idx_t jnode = 0; jnode < xbFs["tocn"].shape(0); ++jnode) { + 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) + @@ -342,11 +347,11 @@ namespace gdasapp { // Initialize the diffusion central block // if (fullConfig.has("diffusion")) { // const eckit::LocalConfiguration diffusionConfig(fullConfig, "diffusion"); -// chem::ExplicitDiffusionParameters params; +// fv3jedi::ExplicitDiffusionParameters params; // params.deserialize(diffusionConfig); // oops::GeometryData geometryData(geom.functionSpace(), bkgErrFs["tocn"], true, this->getComm()); // const oops::FieldSet3D dumyXb(cycleDate, this->getComm()); -// chem::ExplicitDiffusion diffuse(geometryData, chemVars, diffusionConfig, params, dumyXb, dumyXb); +// fv3jedi::ExplicitDiffusion diffuse(geometryData, chemVars, diffusionConfig, params, dumyXb, dumyXb); // diffuse.read(); // // // Smooth the field @@ -357,17 +362,21 @@ namespace gdasapp { // } // Rescale + std::cout << "Rescale"; double rescale; fullConfig.get("rescale", rescale); util::multiplyFieldSet(bkgErrFs, rescale); // We want to write with fv3jedi, not atlas: Syncronize with fv3jedi Increment + std::cout << "write with fv3jedi"; bkgErr.fromFieldSet(bkgErrFs); // Save the background error + std::cout << "save the background error"; const eckit::LocalConfiguration bkgErrorConfig(fullConfig, "background error"); bkgErr.write(bkgErrorConfig); + std::cout << "after save background"; return 0; } From 407e24497b4153b66a52f2d975366ccba28aceb1 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Fri, 5 Apr 2024 13:25:05 -0500 Subject: [PATCH 03/14] Cleaned up chem_diagb.h --- utils/chem/chem_diagb.h | 143 +++++----------------------------------- 1 file changed, 16 insertions(+), 127 deletions(-) diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index 2dfe7f966..11b79b86e 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -142,128 +142,51 @@ namespace gdasapp { double sshMax(0.0); fullConfig.get("max ssh", sshMax); - // Get the layer thicknesses and convert to depth -// auto h = atlas::array::make_view(xbFs["hocn"]); -// atlas::array::ArrayT depth(h.shape(0), h.shape(1)); -// auto viewDepth = atlas::array::make_view(depth); - -// for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { -// viewDepth(jnode, 0) = 0.5 * h(jnode, 0); -// for (atlas::idx_t level = 1; level < h.shape(1); ++level) { -// viewDepth(jnode, level) = viewDepth(jnode, level-1) + -// 0.5 * (h(jnode, level-1 ) + h(jnode, level)); -// } -// } - - // Get the mixed layer depth -// auto mld = atlas::array::make_view(xbFs["mom6_mld"]); -// atlas::array::ArrayT mldindex(mld.shape(0), mld.shape(1)); -// auto viewMldindex = atlas::array::make_view(mldindex); - -// for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { -// std::vector testMld; -// for (atlas::idx_t level = 0; level < viewDepth.shape(1); ++level) { -// testMld.push_back(std::abs(viewDepth(jnode, level) - mld(jnode, 0))); -// } -// auto minVal = std::min_element(testMld.begin(), testMld.end()); - -// viewMldindex(jnode, 0) = std::distance(testMld.begin(), minVal); -// } - // Loop through variables for (auto & var : chemVars.variables()) { - // Skip the layer thickness variable - if (var == "hocn") { - continue; - } 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) { - // Early exit if thickness is 0 or on a ghost cell -// if (ghostView(jnode) > 0 || abs(h(jnode, 0)) <= 0.1) { -// continue; -// } // get indices of neighbor cells auto neighbors = get_neighbors_of_node(jnode); int nbh = neighbors.size(); - // 2D case - if (xbFs[var].shape(1) == 1) { + // 3D case + 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) { -// if ( abs(h(neighbors[nn], 0)) <= 0.1 ) { -// continue; -// } - local.push_back(bkg(neighbors[nn], 0)); + int levelMin = std::max(0, level - nbz); + int levelMax = level + nbz; + levelMin = 0; + levelMax = 1; //nzMld; + for (int ll = levelMin; ll <= levelMax; ++ll) { + local.push_back(bkg(neighbors[nn], ll)); + } } //Set the minimum number of points - int minn = 4; /// probably should be passed through the config + 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, 0) = std::sqrt(stdDev / (local.size() - 1)); - } - // Filter out ssh std. dev. - if (var == "ssh") { - stdDevBkg(jnode, 0) = std::min(stdDevBkg(jnode, 0), sshMax); + stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)) + additiveStdDev; } } - continue; // early exit, no need to loop through levels - } else { - // 3D case - int nbz = 1; // Number of closest point in the vertical, above and below -// int nzMld = std::max(10, viewMldindex(jnode, 0)); - 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; -// if (level < nzMld) { - // If in the MLD, compute the std. dev. throughout the MLD - levelMin = 0; - levelMax = 1; //nzMld; -// } - for (int ll = levelMin; ll <= levelMax; ++ll) { -// if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { -// continue; -// } - 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 (var == "tocn" & level == 0) { -// additiveStdDev = sstBkgErrMin; -// } - if (stdDev > 0.0 || local.size() > 2) { - stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)) + additiveStdDev; - } - } - } // end level - } // end 3D case + } // end level } // end jnode } // end var @@ -276,10 +199,6 @@ namespace gdasapp { int niter(0); fullConfig.get("simple smoothing.horizontal iterations", niter); for (auto & var : chemVars.variables()) { - // Skip the layer thickness variable -// if (var == "hocn") { -// continue; -// } // Horizontal averaging for (int iter = 0; iter < niter; ++iter) { @@ -290,12 +209,7 @@ namespace gdasapp { // 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["tocn"].shape(0); ++jnode) { for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { - // Early exit if thickness is 0 or on a ghost cell -// if (abs(h(jnode, level)) <= 0.1) { -// continue; -// } // Ocean or ice node, do something std::vector local; @@ -303,9 +217,6 @@ namespace gdasapp { int nbh = neighbors.size(); for (int nn = 0; nn < neighbors.size(); ++nn) { int nbNode = neighbors[nn]; -// if ( abs(h(nbNode, level)) <= 0.1 ) { -// continue; -// } local.push_back(stdDevBkg(nbNode, level)); } @@ -326,7 +237,6 @@ namespace gdasapp { 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["tocn"].shape(0); ++jnode) { 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) + @@ -340,27 +250,6 @@ namespace gdasapp { } } - /// Use explicit diffusion to smooth the background error - // ------------------------------------------------------ - // TODO: Use this once Travis adds the option to skip the normalization. - // The output is currently in [0, ~1000] - // Initialize the diffusion central block -// if (fullConfig.has("diffusion")) { -// const eckit::LocalConfiguration diffusionConfig(fullConfig, "diffusion"); -// fv3jedi::ExplicitDiffusionParameters params; -// params.deserialize(diffusionConfig); -// oops::GeometryData geometryData(geom.functionSpace(), bkgErrFs["tocn"], true, this->getComm()); -// const oops::FieldSet3D dumyXb(cycleDate, this->getComm()); -// fv3jedi::ExplicitDiffusion diffuse(geometryData, chemVars, diffusionConfig, params, dumyXb, dumyXb); -// diffuse.read(); -// -// // Smooth the field -// oops::FieldSet3D dx(cycleDate, this->getComm()); -// dx.deepCopy(bkgErrFs); -// diffuse.multiply(dx); -// bkgErrFs = dx.fieldSet(); -// } - // Rescale std::cout << "Rescale"; double rescale; From a66bd2c1ad6385253bc41ff99a0ad8bb04da3686 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Thu, 11 Apr 2024 13:43:39 -0500 Subject: [PATCH 04/14] Added yaml file and revised based on github comments --- parm/aero/berror/aero_diagb.yaml | 60 ++++++++++++++++++++++++++++++++ utils/chem/chem_diagb.cc | 1 - utils/chem/chem_diagb.h | 18 ++-------- 3 files changed, 63 insertions(+), 16 deletions(-) create mode 100644 parm/aero/berror/aero_diagb.yaml diff --git a/parm/aero/berror/aero_diagb.yaml b/parm/aero/berror/aero_diagb.yaml new file mode 100644 index 000000000..fc4a37323 --- /dev/null +++ b/parm/aero/berror/aero_diagb.yaml @@ -0,0 +1,60 @@ +geometry: + fms initialization: + namelist filename: Data/fv3files/fmsmpp.nml + field table filename: Data/fv3files/field_table_gfdl + akbk: Data/fv3files/akbk64.nc4 + layout: [1,2] + npx: 13 + npy: 13 + npz: 64 + field metadata override: Data/fieldmetadata/gfs-aerosol.yaml +background: + datetime: 2018-04-15T00:00:00Z + 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: Data/inputs/gfs_aero_c12/mem001 + filename_core: 20180415.000000.fv_core.res.nc + filename_trcr: 20180415.000000.fv_tracer.res.nc + filename_sfcd: 20180415.000000.sfc_data.nc + filename_sfcw: 20180415.000000.fv_srf_wnd.res.nc + filename_cplr: 20180415.000000.coupler.res + +background error: + datadir: ./ + date: 2018-04-15T00:00:00Z + 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" +min sst: 0.0 # Added to sst bkg. err. +max ssh: 0.0 # Limits the amplitude of the unbalanced bkg err +number of halo points: 4 +number of neighbors: 16 +simple smoothing: + horizontal iterations: 2 + vertical iterations: 1 diff --git a/utils/chem/chem_diagb.cc b/utils/chem/chem_diagb.cc index 218f6820d..aec91d104 100644 --- a/utils/chem/chem_diagb.cc +++ b/utils/chem/chem_diagb.cc @@ -3,7 +3,6 @@ int main(int argc, char ** argv) { oops::Run run(argc, argv); - std::cout << "argv = " << argv ; gdasapp::FV3ChemDiagB fv3chemdiagb; return run.execute(fv3chemdiagb); } diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index 11b79b86e..9b347c4df 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -1,6 +1,5 @@ #pragma once -#include #include #include #include @@ -53,7 +52,6 @@ namespace gdasapp { int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { /// Setup the fv3jedi geometry - std::cout << "int execute"; oops::Log::info() << "====================== geometry" << std::endl; const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); const fv3jedi::Geometry geom(geomConfig, this->getComm()); @@ -72,7 +70,6 @@ namespace gdasapp { /// Read the background // -------------------- - std::cout << "Read the background"; oops::Log::info() << "====================== read bkg" << std::endl; fv3jedi::State xb(geom, chemVars, cycleDate); const eckit::LocalConfiguration bkgConfig(fullConfig, "background"); @@ -85,7 +82,6 @@ namespace gdasapp { /// Create the mesh connectivity (Copy/paste of Francois's stuff) // -------------------------------------------------------------- // Build edges, then connections between nodes and edges - std::cout << "Create the mesh connectivity"; int nbHalos(2); fullConfig.get("number of halo points", nbHalos); int nbNeighbors(4); @@ -158,13 +154,11 @@ namespace gdasapp { // 3D case 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) { + 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; - levelMin = 0; - levelMax = 1; //nzMld; for (int ll = levelMin; ll <= levelMax; ++ll) { local.push_back(bkg(neighbors[nn], ll)); } @@ -183,14 +177,13 @@ namespace gdasapp { // 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)) + additiveStdDev; + stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)); } } } // end level } // end jnode } // end var - // TODO(G): Assume that the steric balance explains 97% of ssh ... or do it properly ... maybe /// Smooth the fields @@ -208,10 +201,9 @@ namespace gdasapp { 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 level = 0; level <= xbFs[var].shape(1); ++level) { for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { - // Ocean or ice node, do something std::vector local; auto neighbors = get_neighbors_of_node(jnode); int nbh = neighbors.size(); @@ -251,21 +243,17 @@ namespace gdasapp { } // Rescale - std::cout << "Rescale"; double rescale; fullConfig.get("rescale", rescale); util::multiplyFieldSet(bkgErrFs, rescale); // We want to write with fv3jedi, not atlas: Syncronize with fv3jedi Increment - std::cout << "write with fv3jedi"; bkgErr.fromFieldSet(bkgErrFs); // Save the background error - std::cout << "save the background error"; const eckit::LocalConfiguration bkgErrorConfig(fullConfig, "background error"); bkgErr.write(bkgErrorConfig); - std::cout << "after save background"; return 0; } From 2a74d81ffa0befdb0c89acbf0eb9ac477c43bf4e Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Fri, 12 Apr 2024 11:31:49 -0500 Subject: [PATCH 05/14] Converted aero_diagb.yaml to template named aero_diagb.yaml.j2 --- .../{aero_diagb.yaml => aero_diagb.yaml.j2} | 34 +++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) rename parm/aero/berror/{aero_diagb.yaml => aero_diagb.yaml.j2} (73%) diff --git a/parm/aero/berror/aero_diagb.yaml b/parm/aero/berror/aero_diagb.yaml.j2 similarity index 73% rename from parm/aero/berror/aero_diagb.yaml rename to parm/aero/berror/aero_diagb.yaml.j2 index fc4a37323..a54a0dd94 100644 --- a/parm/aero/berror/aero_diagb.yaml +++ b/parm/aero/berror/aero_diagb.yaml.j2 @@ -1,15 +1,17 @@ geometry: - fms initialization: - namelist filename: Data/fv3files/fmsmpp.nml - field table filename: Data/fv3files/field_table_gfdl - akbk: Data/fv3files/akbk64.nc4 - layout: [1,2] - npx: 13 - npy: 13 - npz: 64 - field metadata override: Data/fieldmetadata/gfs-aerosol.yaml + 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: 2018-04-15T00:00:00Z + datetime: '{{ current_cycle | to_isotime }}' filetype: fms restart state variables: [mass_fraction_of_sulfate_in_air, mass_fraction_of_hydrophobic_black_carbon_in_air, @@ -21,16 +23,14 @@ background: 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: Data/inputs/gfs_aero_c12/mem001 - filename_core: 20180415.000000.fv_core.res.nc - filename_trcr: 20180415.000000.fv_tracer.res.nc - filename_sfcd: 20180415.000000.sfc_data.nc - filename_sfcw: 20180415.000000.fv_srf_wnd.res.nc - filename_cplr: 20180415.000000.coupler.res + 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: 2018-04-15T00:00:00Z + date: '{{ current_cycle | to_isotime }}' exp: bkgerr_stddev type: incr From b021a79c49b0bdb3d8c8c7e1e81aeeb198faee30 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Fri, 12 Apr 2024 12:53:30 -0500 Subject: [PATCH 06/14] Changes to fix coding norm errors. --- utils/chem/chem_diagb.h | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index 9b347c4df..e3c97ec12 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -8,6 +8,7 @@ #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" @@ -15,7 +16,6 @@ #include "atlas/util/Earth.h" #include "atlas/util/Geometry.h" #include "atlas/util/Point.h" -#include "atlas/functionspace/NodeColumns.h" #include "oops/base/FieldSet3D.h" #include "oops/base/GeometryData.h" @@ -30,8 +30,8 @@ #include "fv3jedi/Geometry/Geometry.h" #include "fv3jedi/Increment/Increment.h" #include "fv3jedi/State/State.h" -//#include "fv3jedi/ExplicitDiffusion/ExplicitDiffusion.h" -//#include "fv3jedi/ExplicitDiffusion/ExplicitDiffusionParameters.h" +// #include "fv3jedi/ExplicitDiffusion/ExplicitDiffusion.h" +// #include "fv3jedi/ExplicitDiffusion/ExplicitDiffusionParameters.h" namespace gdasapp { /** @@ -94,7 +94,8 @@ namespace gdasapp { 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) + // 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); @@ -147,7 +148,6 @@ namespace gdasapp { // 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(); @@ -163,7 +163,7 @@ namespace gdasapp { local.push_back(bkg(neighbors[nn], ll)); } } - //Set the minimum number of points + // Set the minimum number of points int minn = 6; /// probably should be passed through the config if (local.size() >= minn) { // Mean @@ -192,10 +192,8 @@ namespace gdasapp { 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]); @@ -203,7 +201,6 @@ namespace gdasapp { // 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(); @@ -213,7 +210,8 @@ namespace gdasapp { } if (local.size() > 2) { - stdDevBkg(jnode, level) = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); + stdDevBkg(jnode, level) = std::accumulate(local.begin(), local.end(), 0.0) / + local.size(); } } } @@ -231,9 +229,9 @@ namespace gdasapp { 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) + + stdDevBkg(jnode, level) = (tmpArray(jnode, level-1) + tmpArray(jnode, level) + - tmpArray(jnode, level+1) ) / 3.0; + tmpArray(jnode, level+1)) / 3.0; } stdDevBkg(jnode, 0) = stdDevBkg(jnode, 1); } From d0c6f24a8000a1201ddf4728a6e8697237ab67b8 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Fri, 12 Apr 2024 13:06:29 -0500 Subject: [PATCH 07/14] More coding norm fixes. --- utils/chem/chem_diagb.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index e3c97ec12..a01c56dd7 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include @@ -94,7 +95,7 @@ namespace gdasapp { const auto & node2edge = mesh.nodes().edge_connectivity(); const auto & edge2node = mesh.edges().node_connectivity(); - // Lambda function to get the neighbors of a node + // 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{}; @@ -210,7 +211,7 @@ namespace gdasapp { } if (local.size() > 2) { - stdDevBkg(jnode, level) = std::accumulate(local.begin(), local.end(), 0.0) / + stdDevBkg(jnode, level) = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); } } From 66b9bd056fc22ca27c2310ee8d55268d7776e4ee Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Mon, 15 Apr 2024 08:47:03 -0500 Subject: [PATCH 08/14] Removed 2 unneeded lines from chem_diagb.h --- utils/chem/chem_diagb.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index a01c56dd7..4aabe0b1b 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -31,8 +31,6 @@ #include "fv3jedi/Geometry/Geometry.h" #include "fv3jedi/Increment/Increment.h" #include "fv3jedi/State/State.h" -// #include "fv3jedi/ExplicitDiffusion/ExplicitDiffusion.h" -// #include "fv3jedi/ExplicitDiffusion/ExplicitDiffusionParameters.h" namespace gdasapp { /** From a3891f72a4e8c14111412e33bf415353ddfc8451 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Mon, 15 Apr 2024 09:30:44 -0500 Subject: [PATCH 09/14] Remove ssh, sst items from aero_diag.yaml.j2 and chem_digb.h --- parm/aero/berror/aero_diagb.yaml.j2 | 2 -- utils/chem/chem_diagb.h | 1 - 2 files changed, 3 deletions(-) diff --git a/parm/aero/berror/aero_diagb.yaml.j2 b/parm/aero/berror/aero_diagb.yaml.j2 index a54a0dd94..b39eae148 100644 --- a/parm/aero/berror/aero_diagb.yaml.j2 +++ b/parm/aero/berror/aero_diagb.yaml.j2 @@ -51,8 +51,6 @@ variables: - mass_fraction_of_sea_salt004_in_air rescale: 2.0 # rescales the filtered std. dev. by "rescale" -min sst: 0.0 # Added to sst bkg. err. -max ssh: 0.0 # Limits the amplitude of the unbalanced bkg err number of halo points: 4 number of neighbors: 16 simple smoothing: diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index 4aabe0b1b..55c5e1ade 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -151,7 +151,6 @@ namespace gdasapp { auto neighbors = get_neighbors_of_node(jnode); int nbh = neighbors.size(); - // 3D case 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; From bf78cffba0c2143da78b6a9e420fb1b0cbcf8012 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Mon, 15 Apr 2024 09:54:28 -0500 Subject: [PATCH 10/14] Removed further ssh and sst lines from chem_diagb.h --- utils/chem/chem_diagb.h | 6 ------ 1 file changed, 6 deletions(-) diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index 55c5e1ade..208a36fb0 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -130,13 +130,7 @@ namespace gdasapp { atlas::FieldSet bkgErrFs; bkgErr.toFieldSet(bkgErrFs); - // Get the std dev to add to sst - double sstBkgErrMin(0.0); - fullConfig.get("min sst", sstBkgErrMin); - // Get the max std dev for ssh - double sshMax(0.0); - fullConfig.get("max ssh", sshMax); // Loop through variables From 33a7bb20c7273e505fbdc5b2afb463bd2c8cfe3a Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Mon, 15 Apr 2024 10:51:28 -0500 Subject: [PATCH 11/14] updated executable name to gdasapp_chem_diagb.x --- utils/chem/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/chem/CMakeLists.txt b/utils/chem/CMakeLists.txt index 5cb698209..ad4b544f4 100644 --- a/utils/chem/CMakeLists.txt +++ b/utils/chem/CMakeLists.txt @@ -2,4 +2,4 @@ ecbuild_add_executable( TARGET chem_diagb.x SOURCES chem_diagb.cc ) target_compile_features( chem_diagb.x PUBLIC cxx_std_17) -target_link_libraries( chem_diagb.x PUBLIC NetCDF::NetCDF_CXX oops atlas fv3jedi) +target_link_libraries( gdasapp_chem_diagb.x PUBLIC NetCDF::NetCDF_CXX oops atlas fv3jedi) From 5e33713e25a00d5519b0835265acbcb149af5e6f Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Mon, 15 Apr 2024 11:54:41 -0500 Subject: [PATCH 12/14] Removing blank space. --- utils/chem/chem_diagb.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/utils/chem/chem_diagb.h b/utils/chem/chem_diagb.h index 208a36fb0..b8bfcdae9 100644 --- a/utils/chem/chem_diagb.h +++ b/utils/chem/chem_diagb.h @@ -130,9 +130,6 @@ namespace gdasapp { atlas::FieldSet bkgErrFs; bkgErr.toFieldSet(bkgErrFs); - - - // Loop through variables for (auto & var : chemVars.variables()) { oops::Log::info() << "====================== std dev for " << var << std::endl; From c955d986bb84e45d8560c38555be95ac67a2ff15 Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Tue, 16 Apr 2024 07:38:01 -0500 Subject: [PATCH 13/14] Changed staticb_bump.yaml.j2 to have templated date and datetime. --- parm/aero/berror/staticb_bump.yaml.j2 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) 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 }}' From 1b469aebb263e03955001e54ce69e1c97ad4f4cf Mon Sep 17 00:00:00 2001 From: Andrew Tangborn Date: Wed, 17 Apr 2024 13:16:27 -0500 Subject: [PATCH 14/14] changed executable name in CMakeLists.txt --- utils/chem/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/chem/CMakeLists.txt b/utils/chem/CMakeLists.txt index ad4b544f4..0a0b064b8 100644 --- a/utils/chem/CMakeLists.txt +++ b/utils/chem/CMakeLists.txt @@ -1,5 +1,5 @@ # Diag B -ecbuild_add_executable( TARGET chem_diagb.x +ecbuild_add_executable( TARGET gdasapp_chem_diagb.x SOURCES chem_diagb.cc ) -target_compile_features( chem_diagb.x PUBLIC cxx_std_17) +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)