From 30f3634d2301d07f7dea0e23084a53398b938ed4 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 27 Feb 2024 12:49:12 -0500 Subject: [PATCH 01/21] wip --- test/soca/testinput/gdas_diagb.yaml | 21 ++++ utils/soca/gdas_diagb.cc | 8 ++ utils/soca/gdas_diagb.h | 156 ++++++++++++++++++++++++++++ 3 files changed, 185 insertions(+) create mode 100644 test/soca/testinput/gdas_diagb.yaml create mode 100644 utils/soca/gdas_diagb.cc create mode 100644 utils/soca/gdas_diagb.h diff --git a/test/soca/testinput/gdas_diagb.yaml b/test/soca/testinput/gdas_diagb.yaml new file mode 100644 index 000000000..234f03aa9 --- /dev/null +++ b/test/soca/testinput/gdas_diagb.yaml @@ -0,0 +1,21 @@ +geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + +date: 2018-04-15T09:00:00Z + +background: + date: 2018-04-15T09:00:00Z + basename: ./INPUT/ + ocn_filename: MOM.res.nc + ice_filename: cice.res.nc + read_from_file: 1 + +background error: + datadir: ./ + date: '2018-04-15T09:00:00Z' + exp: bkgerr + type: incr + +variables: + name: [tocn, socn, hocn, ssh, cicen, hicen, hsnon] diff --git a/utils/soca/gdas_diagb.cc b/utils/soca/gdas_diagb.cc new file mode 100644 index 000000000..09332cd2e --- /dev/null +++ b/utils/soca/gdas_diagb.cc @@ -0,0 +1,8 @@ +#include "gdas_diagb.h" +#include "oops/runs/Run.h" + +int main(int argc, char ** argv) { + oops::Run run(argc, argv); + gdasapp::DiagB diagb; + return run.execute(diagb); +} diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h new file mode 100644 index 000000000..02f49acf7 --- /dev/null +++ b/utils/soca/gdas_diagb.h @@ -0,0 +1,156 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "atlas/field.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/Logger.h" + +#include "soca/Geometry/Geometry.h" +#include "soca/Increment/Increment.h" +#include "soca/State/State.h" + +namespace gdasapp { + + class DiagB : public oops::Application { + public: + explicit DiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) + : Application(comm) {} + static const std::string classname() {return "gdasapp::DiagB";} + + int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { + /// Setup the soca geometry + oops::Log::info() << "====================== geometry" << std::endl; + const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); + const soca::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 socaVars(fullConfig, "variables.name"); + + /// Read the background + oops::Log::info() << "====================== read bkg" << std::endl; + soca::State xb(geom, socaVars, cycleDate); + const eckit::LocalConfiguration bkgConfig(fullConfig, "background"); + xb.read(bkgConfig); + std::cout << xb << std::endl; + atlas::FieldSet xbFs; + xb.toFieldSet(xbFs); + + // Create the GeometryData object + oops::Log::info() << "====================== geometryData" << std::endl; + oops::GeometryData geometryData(geom.functionSpace(), xbFs["tocn"], true, this->getComm()); + + // Compute local tree (for now) + oops::Log::info() << "====================== kd tree" << std::endl; + const auto lonlat = + atlas::array::make_view(geometryData.functionSpace().lonlat()); + size_t nNodes = xbFs["tocn"].shape(0); + std::vector lons(nNodes); + std::vector lats(nNodes); + for (int jnode = 0; jnode < nNodes; ++jnode) { + lons[jnode] = lonlat(jnode, 0); + lats[jnode] = lonlat(jnode, 1); + } + geometryData.setLocalTree(lats, lons); + + /// Compute local std. dev. as a proxy of the bkg error + oops::Log::info() << "====================== std dev " << xbFs["tocn"].shape(0) << std::endl; + int nb = 4; // Number of closest point (horizontal) + int nbz = 2; // Number of closest point (vertical) + const auto ghostView = + atlas::array::make_view(geom.functionSpace().ghost()); + + // Create the background error fieldset + soca::Increment bkgErr(geom, socaVars, cycleDate); + bkgErr.zero(); + atlas::FieldSet bkgErrFs; + bkgErr.toFieldSet(bkgErrFs); + + // TODO(G): Need to loop through variables + auto stdDevBkg = atlas::array::make_view(bkgErrFs["tocn"]); + auto bkg = atlas::array::make_view(xbFs["tocn"]); + auto h = atlas::array::make_view(xbFs["hocn"]); + for (atlas::idx_t level = nbz; level < xbFs["tocn"].shape(1) - nbz; ++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 (ghostView(jnode) > 0 || abs(h(jnode, 0)) <= 0.1) { + continue; + } + + // Ocean or ice node, do something + auto jn = geometryData.closestPoints(lonlat(jnode, 1), lonlat(jnode, 0), nb); + std::vector local; + std::vector amIG; + for (int nn = 0; nn < nb; ++nn) { + auto nbNode = jn[nn].payload(); + for (int ll = level - nbz; ll < level + nbz; ++ll) { + if (ghostView(nbNode) > 0 || abs(h(nbNode, ll)) <= 0.1 ) { + continue; + } + local.push_back(bkg(nbNode, ll)); + amIG.push_back(ghostView(nbNode)); + } + } + + if (local.size() > 3) { + // Mean + double mean(0); + mean = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); + + // Standard deviation + double stdDev(0); + for (int nn = 0; nn < nb; ++nn) { + stdDev += std::pow(local[nn] - mean, 2.0); + } + stdDevBkg(jnode, level) = std::sqrt(stdDev / local.size()); + + if (stdDevBkg(jnode, level) > 5.0 ) { + std::cout << " ------------------------- " << std::endl; + std::cout << "mean : " << mean << std::endl; + std::cout << "stdDev : " << stdDevBkg(jnode, level) + << " " << ghostView(jnode) << " " << h(jnode, level) << std::endl; + std::cout << "local : " << local << std::endl; + std::cout << "am I G : " << amIG << std::endl; + } + } + } + } + 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::DiagB"; + } + // ----------------------------------------------------------------------------- + }; + +} // namespace gdasapp From a651fa437269070442e8fb182322f6268c7a05f4 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 27 Feb 2024 12:59:35 -0500 Subject: [PATCH 02/21] wip --- utils/soca/CMakeLists.txt | 7 +++++++ utils/soca/gdas_diagb.h | 13 +++++++------ 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/utils/soca/CMakeLists.txt b/utils/soca/CMakeLists.txt index 9e639ccb3..3a53b6f31 100644 --- a/utils/soca/CMakeLists.txt +++ b/utils/soca/CMakeLists.txt @@ -15,3 +15,10 @@ ecbuild_add_executable( TARGET gdas_socahybridweights.x SOURCES gdas_socahybridweights.cc ) target_compile_features( gdas_socahybridweights.x PUBLIC cxx_std_17) target_link_libraries( gdas_socahybridweights.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) + + +# Diag B +ecbuild_add_executable( TARGET gdas_diagb.x + SOURCES gdas_diagb.cc ) +target_compile_features( gdas_diagb.x PUBLIC cxx_std_17) +target_link_libraries( gdas_diagb.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index 02f49acf7..d6f6be035 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -59,7 +59,8 @@ namespace gdasapp { atlas::FieldSet xbFs; xb.toFieldSet(xbFs); - // Create the GeometryData object + /// Create the GeometryData object + // This is used to initialize a local KD-Tree oops::Log::info() << "====================== geometryData" << std::endl; oops::GeometryData geometryData(geom.functionSpace(), xbFs["tocn"], true, this->getComm()); @@ -78,8 +79,8 @@ namespace gdasapp { /// Compute local std. dev. as a proxy of the bkg error oops::Log::info() << "====================== std dev " << xbFs["tocn"].shape(0) << std::endl; - int nb = 4; // Number of closest point (horizontal) - int nbz = 2; // Number of closest point (vertical) + int nbh = 8; // Number of closest point (horizontal) + int nbz = 1; // Number of closest point (vertical) const auto ghostView = atlas::array::make_view(geom.functionSpace().ghost()); @@ -101,10 +102,10 @@ namespace gdasapp { } // Ocean or ice node, do something - auto jn = geometryData.closestPoints(lonlat(jnode, 1), lonlat(jnode, 0), nb); + auto jn = geometryData.closestPoints(lonlat(jnode, 1), lonlat(jnode, 0), nbh); std::vector local; std::vector amIG; - for (int nn = 0; nn < nb; ++nn) { + for (int nn = 0; nn < nbh; ++nn) { auto nbNode = jn[nn].payload(); for (int ll = level - nbz; ll < level + nbz; ++ll) { if (ghostView(nbNode) > 0 || abs(h(nbNode, ll)) <= 0.1 ) { @@ -122,7 +123,7 @@ namespace gdasapp { // Standard deviation double stdDev(0); - for (int nn = 0; nn < nb; ++nn) { + for (int nn = 0; nn < nbh; ++nn) { stdDev += std::pow(local[nn] - mean, 2.0); } stdDevBkg(jnode, level) = std::sqrt(stdDev / local.size()); From 27a408ba6bc72882448c4eb408e24b33a798c216 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Thu, 29 Feb 2024 13:53:57 -0500 Subject: [PATCH 03/21] using FH's lambda, but still wip --- utils/soca/gdas_diagb.h | 131 ++++++++++++++++++++++++---------------- 1 file changed, 79 insertions(+), 52 deletions(-) diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index d6f6be035..fd5810c5a 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -9,12 +9,15 @@ #include "eckit/config/LocalConfiguration.h" #include "atlas/field.h" +#include "atlas/mesh/actions/BuildEdges.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/base/GeometryData.h" #include "oops/mpi/mpi.h" #include "oops/runs/Application.h" #include "oops/util/DateTime.h" @@ -62,24 +65,45 @@ namespace gdasapp { /// Create the GeometryData object // This is used to initialize a local KD-Tree oops::Log::info() << "====================== geometryData" << std::endl; - oops::GeometryData geometryData(geom.functionSpace(), xbFs["tocn"], true, this->getComm()); - - // Compute local tree (for now) - oops::Log::info() << "====================== kd tree" << std::endl; - const auto lonlat = - atlas::array::make_view(geometryData.functionSpace().lonlat()); - size_t nNodes = xbFs["tocn"].shape(0); - std::vector lons(nNodes); - std::vector lats(nNodes); - for (int jnode = 0; jnode < nNodes; ++jnode) { - lons[jnode] = lonlat(jnode, 0); - lats[jnode] = lonlat(jnode, 1); - } - geometryData.setLocalTree(lats, lons); + + /// Create the mesh connectivity (Copy/past of Francois's stuff) + // Build edges, then connections between nodes and edges + atlas::functionspace::NodeColumns test = geom.functionSpace(); + atlas::Mesh mesh = test.mesh(); + atlas::mesh::actions::build_edges(mesh); + atlas::mesh::actions::build_node_to_edge_connectivity(mesh); + const auto & node2edge = mesh.nodes().edge_connectivity(); + const auto & edge2node = mesh.edges().node_connectivity(); + + oops::Log::info() << "====================== get neighbors" << std::endl; + // Lambda function to get the neighbors of a node + const auto get_neighbors_of_node = [&](const int node) { + std::vector neighbors{}; + neighbors.reserve(6); + if (node >= mesh.nodes().size()) { + std::cout << "Node " << node << " is out of bounds!" << std::endl; + 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 oops::Log::info() << "====================== std dev " << xbFs["tocn"].shape(0) << std::endl; - int nbh = 8; // Number of closest point (horizontal) + //int nbh = 8; // Number of closest point (horizontal) int nbz = 1; // Number of closest point (vertical) const auto ghostView = atlas::array::make_view(geom.functionSpace().ghost()); @@ -95,57 +119,60 @@ namespace gdasapp { auto bkg = atlas::array::make_view(xbFs["tocn"]); auto h = atlas::array::make_view(xbFs["hocn"]); for (atlas::idx_t level = nbz; level < xbFs["tocn"].shape(1) - nbz; ++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 (ghostView(jnode) > 0 || abs(h(jnode, 0)) <= 0.1) { - continue; - } + for (atlas::idx_t jnode = 1; jnode < xbFs["tocn"].shape(0)-1; ++jnode) { + // Early exit if thickness is 0 or on a ghost cell + if (ghostView(jnode) > 0 || abs(h(jnode, 0)) <= 0.1) { + continue; + } - // Ocean or ice node, do something - auto jn = geometryData.closestPoints(lonlat(jnode, 1), lonlat(jnode, 0), nbh); - std::vector local; - std::vector amIG; - for (int nn = 0; nn < nbh; ++nn) { - auto nbNode = jn[nn].payload(); - for (int ll = level - nbz; ll < level + nbz; ++ll) { - if (ghostView(nbNode) > 0 || abs(h(nbNode, ll)) <= 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]; + for (int ll = level - nbz; ll < level + nbz; ++ll) { + if (ghostView(nbNode) > 0 || abs(h(nbNode, ll)) <= 0.1 ) { + continue; + } + local.push_back(bkg(nbNode, ll)); } - local.push_back(bkg(nbNode, ll)); - amIG.push_back(ghostView(nbNode)); } - } - if (local.size() > 3) { - // Mean - double mean(0); - mean = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); + if (local.size() > 3) { + // Mean + double mean = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); - // Standard deviation - double stdDev(0); - for (int nn = 0; nn < nbh; ++nn) { - stdDev += std::pow(local[nn] - mean, 2.0); - } - stdDevBkg(jnode, level) = std::sqrt(stdDev / local.size()); - - if (stdDevBkg(jnode, level) > 5.0 ) { - std::cout << " ------------------------- " << std::endl; - std::cout << "mean : " << mean << std::endl; - std::cout << "stdDev : " << stdDevBkg(jnode, level) - << " " << ghostView(jnode) << " " << h(jnode, level) << std::endl; - std::cout << "local : " << local << std::endl; - std::cout << "am I G : " << amIG << std::endl; + // 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() != 0) { + std::cout << "stddev: " << std::sqrt(stdDev / local.size()) << std::endl; + std::cout << "h: " << h(jnode, 0) << std::endl; + stdDevBkg(jnode, level) = stdDev; //std::sqrt(stdDev / local.size()); + } + + if (stdDevBkg(jnode, level) > 5.0 ) { + std::cout << " ------------------------- " << std::endl; + std::cout << "mean : " << mean << std::endl; + std::cout << "stdDev : " << stdDevBkg(jnode, level) + << " " << ghostView(jnode) << " " << h(jnode, level) << std::endl; + std::cout << "local : " << local << std::endl; + } } } } - } bkgErr.fromFieldSet(bkgErrFs); // Save the background error const eckit::LocalConfiguration bkgErrorConfig(fullConfig, "background error"); bkgErr.write(bkgErrorConfig); + return 0; } + // ----------------------------------------------------------------------------- private: std::string appname() const { From 4dc9d4f23d5680deb7fff87c2deeb127a207f2c9 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Thu, 29 Feb 2024 15:34:36 -0600 Subject: [PATCH 04/21] wip --- utils/soca/gdas_diagb.h | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index fd5810c5a..9258bf4b5 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -79,7 +79,7 @@ namespace gdasapp { // Lambda function to get the neighbors of a node const auto get_neighbors_of_node = [&](const int node) { std::vector neighbors{}; - neighbors.reserve(6); + neighbors.reserve(8); if (node >= mesh.nodes().size()) { std::cout << "Node " << node << " is out of bounds!" << std::endl; return neighbors; @@ -132,7 +132,7 @@ namespace gdasapp { for (int nn = 0; nn < neighbors.size(); ++nn) { int nbNode = neighbors[nn]; for (int ll = level - nbz; ll < level + nbz; ++ll) { - if (ghostView(nbNode) > 0 || abs(h(nbNode, ll)) <= 0.1 ) { + if ( abs(h(nbNode, ll)) <= 0.1 ) { continue; } local.push_back(bkg(nbNode, ll)); @@ -149,11 +149,17 @@ namespace gdasapp { stdDev += std::pow(local[nn] - mean, 2.0); } if (stdDev > 0.0 || local.size() != 0) { - std::cout << "stddev: " << std::sqrt(stdDev / local.size()) << std::endl; - std::cout << "h: " << h(jnode, 0) << std::endl; - stdDevBkg(jnode, level) = stdDev; //std::sqrt(stdDev / local.size()); + //std::cout << "stddev: " << std::sqrt(stdDev / local.size()) << std::endl; + //std::cout << "h: " << h(jnode, 0) << std::endl; + stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)); } + // Extrapolate upper levels + for (int ll = 0; ll < nbz; ++ll) { + stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); + } + + /* if (stdDevBkg(jnode, level) > 5.0 ) { std::cout << " ------------------------- " << std::endl; std::cout << "mean : " << mean << std::endl; @@ -161,9 +167,11 @@ namespace gdasapp { << " " << ghostView(jnode) << " " << h(jnode, level) << std::endl; std::cout << "local : " << local << std::endl; } + */ } } } + bkgErr.fromFieldSet(bkgErrFs); // Save the background error From 89c00e2540ad9062da3b5d20afc59e7f9192b8e7 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Mon, 4 Mar 2024 18:52:58 -0600 Subject: [PATCH 05/21] produces credible results but needs tidy --- .gitmodules | 4 + bundle/CMakeLists.txt | 3 + sorc/fv3-jedi | 2 +- sorc/ioda | 2 +- sorc/iodaconv | 2 +- sorc/land-imsproc | 2 +- sorc/land-jediincr | 2 +- sorc/oops | 2 +- sorc/saber | 2 +- sorc/soca | 2 +- sorc/ufo | 2 +- sorc/vader | 2 +- utils/CMakeLists.txt | 2 +- utils/soca/gdas_diagb.h | 162 +++++++++++++++++++++++++++------------- 14 files changed, 129 insertions(+), 62 deletions(-) diff --git a/.gitmodules b/.gitmodules index f8f34c6fa..40b595645 100644 --- a/.gitmodules +++ b/.gitmodules @@ -77,3 +77,7 @@ path = sorc/fms url = https://github.com/jcsda/fms.git branch = release-stable +[submodule "sorc/daml"] + path = sorc/daml + url = https://github.com/NOAA-EMC/daml.git + branch = develop diff --git a/bundle/CMakeLists.txt b/bundle/CMakeLists.txt index 9493ac0c8..328d049f1 100644 --- a/bundle/CMakeLists.txt +++ b/bundle/CMakeLists.txt @@ -94,6 +94,9 @@ if(BUILD_GDASBUNDLE) option(ENABLE_FV3_JEDI_DATA "Obtain fv3-jedi test data from fv3-jedi-data repository (vs tarball)" ON) ecbuild_bundle( PROJECT fv3-jedi SOURCE "../sorc/fv3-jedi" ) +# DAML + ecbuild_bundle( PROJECT daml SOURCE "../sorc/daml" ) + # SOCA associated repositories # TODO: Move the Icepack fork to EMC github set(BUILD_ICEPACK "ON" CACHE STRING "Build the icepack library") diff --git a/sorc/fv3-jedi b/sorc/fv3-jedi index 4a2dd7ead..6f5fe3f87 160000 --- a/sorc/fv3-jedi +++ b/sorc/fv3-jedi @@ -1 +1 @@ -Subproject commit 4a2dd7ead78798a18f0678572314a9fea48ba8a3 +Subproject commit 6f5fe3f87ed4ca8202bc89e11c6c37ad73154b84 diff --git a/sorc/ioda b/sorc/ioda index 79f96592d..35f2dac24 160000 --- a/sorc/ioda +++ b/sorc/ioda @@ -1 +1 @@ -Subproject commit 79f96592d6cdb5fa5b7696735d7ff4296d583d88 +Subproject commit 35f2dac24beff87406cd5cd5732814bc43496218 diff --git a/sorc/iodaconv b/sorc/iodaconv index 18094c03f..a2f420d73 160000 --- a/sorc/iodaconv +++ b/sorc/iodaconv @@ -1 +1 @@ -Subproject commit 18094c03f0ed2d2468b6fc974bfe4a6798aa1025 +Subproject commit a2f420d73d7cd0e69dfff833eefc74e71e71ee2d diff --git a/sorc/land-imsproc b/sorc/land-imsproc index 2ac7c0d90..6373819ca 160000 --- a/sorc/land-imsproc +++ b/sorc/land-imsproc @@ -1 +1 @@ -Subproject commit 2ac7c0d90d3ebd449ad11bce155a71d381b4534a +Subproject commit 6373819ca034d66523cb7852cd7b1b66f3f8ae07 diff --git a/sorc/land-jediincr b/sorc/land-jediincr index 79576b79d..2923344b3 160000 --- a/sorc/land-jediincr +++ b/sorc/land-jediincr @@ -1 +1 @@ -Subproject commit 79576b79df6a3f1ab8e7dde4425821f2915006b3 +Subproject commit 2923344b33511d80c685c30261ad37896eb50768 diff --git a/sorc/oops b/sorc/oops index 5c6deff5a..201239805 160000 --- a/sorc/oops +++ b/sorc/oops @@ -1 +1 @@ -Subproject commit 5c6deff5ad8399b67b1c73ef9cf9ad58b941c76d +Subproject commit 2012398058043bda46db2b9d80a909810df04cea diff --git a/sorc/saber b/sorc/saber index 6d610aa19..ace19a0ee 160000 --- a/sorc/saber +++ b/sorc/saber @@ -1 +1 @@ -Subproject commit 6d610aa195ff994ba07c85ba4585bfb4055d6264 +Subproject commit ace19a0ee9789f507e10ca6951a5099981a8f706 diff --git a/sorc/soca b/sorc/soca index 2ce6ad95f..729e9b73c 160000 --- a/sorc/soca +++ b/sorc/soca @@ -1 +1 @@ -Subproject commit 2ce6ad95f757b192178e9984a9ad24a206a5047b +Subproject commit 729e9b73c635f450fc817c3b2d1de24d3581f09a diff --git a/sorc/ufo b/sorc/ufo index f8e4a3509..f3734dab1 160000 --- a/sorc/ufo +++ b/sorc/ufo @@ -1 +1 @@ -Subproject commit f8e4a3509b314503623480f7390f6fd4955d747e +Subproject commit f3734dab113dcc3cfa2ac693191da24c3f93022e diff --git a/sorc/vader b/sorc/vader index 2baeb8dfa..b8404a832 160000 --- a/sorc/vader +++ b/sorc/vader @@ -1 +1 @@ -Subproject commit 2baeb8dfa8781a67bcf386bf152f02619b748298 +Subproject commit b8404a83208d46bbfb5add30b6561c54da9c6d82 diff --git a/utils/CMakeLists.txt b/utils/CMakeLists.txt index 49290c0dc..095c387aa 100644 --- a/utils/CMakeLists.txt +++ b/utils/CMakeLists.txt @@ -16,6 +16,6 @@ find_package(atlas REQUIRED) find_package(soca REQUIRED) add_subdirectory(soca) -add_subdirectory(ioda_example) +#add_subdirectory(ioda_example) add_subdirectory(test) add_subdirectory(obsproc) diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index 9258bf4b5..f211c219c 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -17,7 +17,7 @@ #include "atlas/functionspace/NodeColumns.h" #include "oops/base/FieldSet3D.h" -//#include "oops/base/GeometryData.h" +#include "oops/base/GeometryData.h" #include "oops/mpi/mpi.h" #include "oops/runs/Application.h" #include "oops/util/DateTime.h" @@ -28,6 +28,8 @@ #include "soca/Geometry/Geometry.h" #include "soca/Increment/Increment.h" #include "soca/State/State.h" +#include "soca/ExplicitDiffusion/ExplicitDiffusion.h" +#include "soca/ExplicitDiffusion/ExplicitDiffusionParameters.h" namespace gdasapp { @@ -53,6 +55,10 @@ namespace gdasapp { oops::Log::info() << "====================== variables" << std::endl; oops::Variables socaVars(fullConfig, "variables.name"); + /// Read rescaling factor + double rescale; + fullConfig.get("rescale", rescale); + /// Read the background oops::Log::info() << "====================== read bkg" << std::endl; soca::State xb(geom, socaVars, cycleDate); @@ -62,11 +68,7 @@ namespace gdasapp { atlas::FieldSet xbFs; xb.toFieldSet(xbFs); - /// Create the GeometryData object - // This is used to initialize a local KD-Tree - oops::Log::info() << "====================== geometryData" << std::endl; - - /// Create the mesh connectivity (Copy/past of Francois's stuff) + /// Create the mesh connectivity (Copy/paste of Francois's stuff) // Build edges, then connections between nodes and edges atlas::functionspace::NodeColumns test = geom.functionSpace(); atlas::Mesh mesh = test.mesh(); @@ -102,9 +104,8 @@ namespace gdasapp { }; /// Compute local std. dev. as a proxy of the bkg error - oops::Log::info() << "====================== std dev " << xbFs["tocn"].shape(0) << std::endl; //int nbh = 8; // Number of closest point (horizontal) - int nbz = 1; // Number of closest point (vertical) + const auto ghostView = atlas::array::make_view(geom.functionSpace().ghost()); @@ -114,14 +115,82 @@ namespace gdasapp { atlas::FieldSet bkgErrFs; bkgErr.toFieldSet(bkgErrFs); - // TODO(G): Need to loop through variables - auto stdDevBkg = atlas::array::make_view(bkgErrFs["tocn"]); - auto bkg = atlas::array::make_view(xbFs["tocn"]); + // Loop through variables auto h = atlas::array::make_view(xbFs["hocn"]); - for (atlas::idx_t level = nbz; level < xbFs["tocn"].shape(1) - nbz; ++level) { - for (atlas::idx_t jnode = 1; jnode < xbFs["tocn"].shape(0)-1; ++jnode) { + + for (auto & var : socaVars.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]); + + int nbz = 1; // Number of closest point in the vertical, above and below + nbz = std::min(nbz, xbFs[var].shape(1) - 1); + for (atlas::idx_t level = nbz; level < xbFs[var].shape(1) - nbz; ++level) { + oops::Log::info() << " level: " << level << std::endl; + 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; + } + + // 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]; + for (int ll = level - nbz; ll <= level + nbz; ++ll) { + if ( abs(h(nbNode, ll)) <= 0.1 ) { + continue; + } + local.push_back(bkg(nbNode, ll)); + } + } + + if (local.size() > 3) { + // 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() != 0) { + //std::cout << "stddev: " << std::sqrt(stdDev / local.size()) << std::endl; + //std::cout << "h: " << h(jnode, 0) << std::endl; + stdDevBkg(jnode, level) = rescale * std::sqrt(stdDev / (local.size() - 1)); + } + + // Extrapolate upper levels + for (int ll = 0; ll < nbz; ++ll) { + stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); + } + + /* + if (stdDevBkg(jnode, level) > 5.0 ) { + std::cout << " ------------------------- " << std::endl; + std::cout << "mean : " << mean << std::endl; + std::cout << "stdDev : " << stdDevBkg(jnode, level) + << " " << ghostView(jnode) << " " << h(jnode, level) << std::endl; + std::cout << "local : " << local << std::endl; + } + */ + } + } + } + } + + // TODO(G): Assume that the steric balance explains 97% of ssh ... or do it properly ... maybe + + /* + bkgErrFs["tocn"].haloExchange(); + + /// Smooth the fields + for (atlas::idx_t level = 0; level < xbFs["tocn"].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 (ghostView(jnode) > 0 || abs(h(jnode, 0)) <= 0.1) { + if (abs(h(jnode, 0)) <= 0.1) { continue; } @@ -131,47 +200,38 @@ namespace gdasapp { int nbh = neighbors.size(); for (int nn = 0; nn < neighbors.size(); ++nn) { int nbNode = neighbors[nn]; - for (int ll = level - nbz; ll < level + nbz; ++ll) { - if ( abs(h(nbNode, ll)) <= 0.1 ) { + if ( abs(h(nbNode, level)) <= 0.1 ) { continue; - } - local.push_back(bkg(nbNode, ll)); - } - } - - if (local.size() > 3) { - // 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() != 0) { - //std::cout << "stddev: " << std::sqrt(stdDev / local.size()) << std::endl; - //std::cout << "h: " << h(jnode, 0) << std::endl; - stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)); - } - - // Extrapolate upper levels - for (int ll = 0; ll < nbz; ++ll) { - stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); } + local.push_back(stdDevBkg(nbNode, level)); + } - /* - if (stdDevBkg(jnode, level) > 5.0 ) { - std::cout << " ------------------------- " << std::endl; - std::cout << "mean : " << mean << std::endl; - std::cout << "stdDev : " << stdDevBkg(jnode, level) - << " " << ghostView(jnode) << " " << h(jnode, level) << std::endl; - std::cout << "local : " << local << std::endl; - } - */ - } - } + if (local.size() > 2) { + stdDevBkg(jnode, level) = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); + } + } } - + */ + /// Smooth the background error + // That doesn't seem to work, the output is in [0, ~1000] + // Initialize the diffusion central block + if (fullConfig.has("diffusion")) { + const eckit::LocalConfiguration diffusionConfig(fullConfig, "diffusion"); + soca::ExplicitDiffusionParameters params; + params.deserialize(diffusionConfig); + oops::GeometryData geometryData(geom.functionSpace(), bkgErrFs["tocn"], true, this->getComm()); + const oops::FieldSet3D dumyXb(cycleDate, this->getComm()); + soca::ExplicitDiffusion diffuse(geometryData, socaVars, diffusionConfig, params, dumyXb, dumyXb); + diffuse.read(); + + // Smooth the field + oops::FieldSet3D dx(cycleDate, this->getComm()); + dx.deepCopy(bkgErrFs); + diffuse.multiply(dx); + bkgErrFs = dx.fieldSet(); + } + + // We want to write with soca, not atlas: Syncronize with soca Increment bkgErr.fromFieldSet(bkgErrFs); // Save the background error From 37288e7bec03f417ed99bc8c60a5417447283b2e Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 5 Mar 2024 14:09:25 -0500 Subject: [PATCH 06/21] ctests ok --- parm/soca/berror/saber_blocks.yaml | 138 +++++---------- parm/soca/berror/soca_diagb.yaml | 23 +++ scripts/exgdas_global_marine_analysis_bmat.sh | 23 ++- scripts/exgdas_global_marine_analysis_post.py | 12 +- scripts/exgdas_global_marine_analysis_prep.py | 38 ++-- test/soca/testinput/gdas_diagb.yaml | 2 +- utils/soca/gdas_diagb.h | 165 ++++++++---------- 7 files changed, 188 insertions(+), 213 deletions(-) create mode 100644 parm/soca/berror/soca_diagb.yaml diff --git a/parm/soca/berror/saber_blocks.yaml b/parm/soca/berror/saber_blocks.yaml index 7957d3b50..dc2666786 100644 --- a/parm/soca/berror/saber_blocks.yaml +++ b/parm/soca/berror/saber_blocks.yaml @@ -1,101 +1,43 @@ -covariance model: hybrid -components: -- covariance: - # This will setup B = K D C_v C_h C_h^{T} C_v^{T} D K^{T} - covariance model: SABER - saber central block: - saber block name: EXPLICIT_DIFFUSION - active variables: [tocn, socn, ssh] - geometry: - mom6_input_nml: mom_input.nml - fields metadata: ./fields_metadata.yaml - group mapping: - - name: ocean - variables: [tocn, socn, socn, ssh] - - name: ice - variables: [cicen] - read: - groups: - - name: ocean - horizontal: - filename: hz_ocean.nc - vertical: - filename: vt_ocean.nc - - name: ice - horizontal: - filename: hz_ice.nc +covariance model: SABER +saber central block: + saber block name: EXPLICIT_DIFFUSION + active variables: [tocn, socn, ssh, cicen] + geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + group mapping: + - name: ocean + variables: + - tocn + - socn + - ssh + - name: ice + variables: + - cicen + read: + groups: + - name: ocean + horizontal: + filename: hz_ocean.nc + vertical: + filename: vt_ocean.nc + - name: ice + horizontal: + filename: hz_ice.nc - linear variable change: - input variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] - output variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] - linear variable changes: - - linear variable change name: BkgErrFILT - ocean_depth_min: 500 # [m] - rescale_bkgerr: 1.0 - efold_z: 1500.0 # [m] +saber outer blocks: +- saber block name: StdDev + read: + model file: + date: '{{ATM_WINDOW_MIDDLE}}' + basename: ./ + ocn_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' + ice_filename: 'ice.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' + read_from_file: 3 - - linear variable change name: BkgErrSOCA - read_from_file: 3 - basename: ./static_ens/ - ocn_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc' - ice_filename: 'ice.bkgerr_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc' - date: '{{ATM_WINDOW_MIDDLE}}' - t_min: 0.1 - t_max: 5.0 - s_min: 0.1 - s_max: 1.0 - ssh_min: 0.0 # std ssh=0 => ssh balance applied as - ssh_max: 1.0 # strong constraint - cicen_min: 0.1 - cicen_max: 0.5 - hicen_min: 0.0 - hicen_max: 0.0 - standard deviation: true - - - linear variable change name: BalanceSOCA - ksshts: - nlayers: 10 - - weight: - value: 0.25 -#- covariance: -# covariance model: ensemble -# members from template: -# template: -# read_from_file: 1 -# date: '{{ATM_WINDOW_MIDDLE}}' -# basename: ./static_ens/ -# ocn_filename: 'ocn.pert.steric.%mem%.nc' -# ice_filename: 'ice.%mem%.nc' -# state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon] -# pattern: '%mem%' -# nmembers: ${ENS_SIZE} -# localization: -# localization method: SABER -# saber central block: -# saber block name: EXPLICIT_DIFFUSION -# active variables: [tocn, socn, ssh] -# geometry: -# mom6_input_nml: mom_input.nml -# fields metadata: ./fields_metadata.yaml -# group mapping: -# - name: ocean -# variables: [tocn, socn, ssh, uocn, vocn] -# - name: ice -# variables: [cicen, hicen, hsnon] -# read: -# groups: -# - name: ocean -# horizontal: -# filename: hz_ocean.nc -# - name: ice -# horizontal: -# filename: hz_ice.nc -# -# weight: -# read_from_file: 3 -# basename: ./ -# ocn_filename: 'ocn.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc' -# ice_filename: 'ice.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc' -# date: '{{ATM_WINDOW_MIDDLE}}' +linear variable change: + input variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] + output variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] + linear variable changes: + - linear variable change name: BalanceSOCA diff --git a/parm/soca/berror/soca_diagb.yaml b/parm/soca/berror/soca_diagb.yaml new file mode 100644 index 000000000..f8c70b545 --- /dev/null +++ b/parm/soca/berror/soca_diagb.yaml @@ -0,0 +1,23 @@ +geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + +date: '{{ATM_WINDOW_MIDDLE}}' + +background: + date: '{{ATM_WINDOW_MIDDLE}}' + basename: ./bkg/ + ocn_filename: '{{RUN}}.ocean.t{{gcyc}}z.inst.f009.nc' + ice_filename: '{{RUN}}.agg_ice.t{{gcyc}}z.inst.f009.nc' + read_from_file: 1 + +background error: + datadir: ./ + date: '{{ATM_WINDOW_MIDDLE}}' + exp: bkgerr_stddev + type: incr + +variables: + name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon] + +rescale: 2.0 diff --git a/scripts/exgdas_global_marine_analysis_bmat.sh b/scripts/exgdas_global_marine_analysis_bmat.sh index 1b7f7c215..ec4facc16 100755 --- a/scripts/exgdas_global_marine_analysis_bmat.sh +++ b/scripts/exgdas_global_marine_analysis_bmat.sh @@ -58,11 +58,11 @@ fi ################################################################################ # Write ensemble weights for the hybrid envar -$APRUN_OCNANAL $JEDI_BIN/gdas_socahybridweights.x soca_ensweights.yaml -export err=$?; err_chk -if [ $err -gt 0 ]; then - exit $err -fi +#$APRUN_OCNANAL $JEDI_BIN/gdas_socahybridweights.x soca_ensweights.yaml +#export err=$?; err_chk +#if [ $err -gt 0 ]; then +# exit $err +#fi ################################################################################ # Ensemble perturbations for the EnVAR and diagonal of static B @@ -79,8 +79,17 @@ fi # - moments # Process static ensemble -clean_yaml soca_clim_ens_moments.yaml -$APRUN_OCNANAL $JEDI_BIN/gdas_ens_handler.x soca_ensb.yaml +#clean_yaml soca_clim_ens_moments.yaml +#$APRUN_OCNANAL $JEDI_BIN/gdas_ens_handler.x soca_ensb.yaml +#export err=$?; err_chk +#if [ $err -gt 0 ]; then +# exit $err +#fi + +################################################################################ +# Compute the parametric diag of B +echo "++++++++++++++++++++++++++++++ $APRUN_OCNANAL" +$APRUN_OCNANAL $JEDI_BIN/gdas_diagb.x soca_diagb.yaml export err=$?; err_chk if [ $err -gt 0 ]; then exit $err diff --git a/scripts/exgdas_global_marine_analysis_post.py b/scripts/exgdas_global_marine_analysis_post.py index 6e30f7cad..bb2f8ed45 100755 --- a/scripts/exgdas_global_marine_analysis_post.py +++ b/scripts/exgdas_global_marine_analysis_post.py @@ -61,12 +61,12 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): domains = ['ocn', 'ice'] for domain in domains: # Copy of the diagonal of the background error for the cycle - post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.bkgerr_stddev.incr.{bdate}.nc'), + post_file_list.append([os.path.join(anl_dir, f'{domain}.bkgerr_stddev.incr.{mdate}.nc'), os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.bkgerr_stddev.nc')]) # Copy the recentering error - post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.ssh_recentering_error.incr.{bdate}.nc'), - os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.recentering_error.nc')]) + #post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.ssh_recentering_error.incr.{bdate}.nc'), + # os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.recentering_error.nc')]) # Copy the ice and ocean increments post_file_list.append([os.path.join(anl_dir, 'Data', f'{domain}.3dvarfgat_pseudo.incr.{mdate}.nc'), @@ -77,9 +77,9 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}ana.nc')]) # Copy of the ssh diagnostics -for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']: - post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'), - os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')]) +#for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']: +# post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'), +# os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')]) # Copy DA grid (computed for the start of the window) post_file_list.append([os.path.join(anl_dir, 'soca_gridspec.nc'), diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index e3e52acea..4fd87d905 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -101,7 +101,9 @@ def find_clim_ens(input_date): envconfig = {'window_begin': f"{window_begin.strftime('%Y-%m-%dT%H:%M:%SZ')}", 'ATM_WINDOW_BEGIN': window_begin_iso, 'ATM_WINDOW_MIDDLE': window_middle_iso, - 'ATM_WINDOW_LENGTH': f"PT{os.getenv('assim_freq')}H"} + 'ATM_WINDOW_LENGTH': f"PT{os.getenv('assim_freq')}H", + 'gcyc': gcyc, + 'RUN': RUN} stage_cfg = YAMLFile(path=os.path.join(gdas_home, 'parm', 'templates', 'stage.yaml')) stage_cfg = Template.substitute_structure(stage_cfg, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) stage_cfg = Template.substitute_structure(stage_cfg, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get) @@ -158,18 +160,20 @@ def find_clim_ens(input_date): for mem in range(1, nmem_ens+1): cice_fname = os.path.realpath(os.path.join(static_ens, "ice."+str(mem)+".nc")) bkg_utils.cice_hist2fms(cice_fname, cice_fname) -else: - logger.info("---------------- Stage offline ensemble members") - ens_member_list = [] - clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None)) - nmem_ens = len(glob.glob(os.path.abspath(os.path.join(clim_ens_dir, 'ocn.*.nc')))) - for domain in ['ocn', 'ice']: - for mem in range(1, nmem_ens+1): - fname = domain+"."+str(mem)+".nc" - fname_in = os.path.abspath(os.path.join(clim_ens_dir, fname)) - fname_out = os.path.abspath(os.path.join(static_ens, fname)) - ens_member_list.append([fname_in, fname_out]) - FileHandler({'copy': ens_member_list}).sync() + +# Commented out while testing the parametric diagb +#else: +# logger.info("---------------- Stage offline ensemble members") +# ens_member_list = [] +# clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None)) +# nmem_ens = len(glob.glob(os.path.abspath(os.path.join(clim_ens_dir, 'ocn.*.nc')))) +# for domain in ['ocn', 'ice']: +# for mem in range(1, nmem_ens+1): +# fname = domain+"."+str(mem)+".nc" +# fname_in = os.path.abspath(os.path.join(clim_ens_dir, fname)) +# fname_out = os.path.abspath(os.path.join(static_ens, fname)) +# ens_member_list.append([fname_in, fname_out]) +# FileHandler({'copy': ens_member_list}).sync() os.environ['ENS_SIZE'] = str(nmem_ens) ################################################################################ @@ -189,6 +193,14 @@ def find_clim_ens(input_date): # generate the YAML file for the post processing of the clim. ens. B berror_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror') +logger.info(f"---------------- generate soca_diagb.yaml") +diagb_yaml = os.path.join(anl_dir, 'soca_diagb.yaml') +diagb_yaml_template = os.path.join(berror_yaml_dir, 'soca_diagb.yaml') +config = YAMLFile(path=diagb_yaml_template) +config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) +config = Template.substitute_structure(config, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get) +config.save(diagb_yaml) + logger.info(f"---------------- generate soca_ensb.yaml") berr_yaml = os.path.join(anl_dir, 'soca_ensb.yaml') berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_ensb.yaml') diff --git a/test/soca/testinput/gdas_diagb.yaml b/test/soca/testinput/gdas_diagb.yaml index 234f03aa9..e865db768 100644 --- a/test/soca/testinput/gdas_diagb.yaml +++ b/test/soca/testinput/gdas_diagb.yaml @@ -18,4 +18,4 @@ background error: type: incr variables: - name: [tocn, socn, hocn, ssh, cicen, hicen, hsnon] + name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon] diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index f211c219c..ec8b9de66 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -64,7 +64,6 @@ namespace gdasapp { soca::State xb(geom, socaVars, cycleDate); const eckit::LocalConfiguration bkgConfig(fullConfig, "background"); xb.read(bkgConfig); - std::cout << xb << std::endl; atlas::FieldSet xbFs; xb.toFieldSet(xbFs); @@ -83,7 +82,6 @@ namespace gdasapp { std::vector neighbors{}; neighbors.reserve(8); if (node >= mesh.nodes().size()) { - std::cout << "Node " << node << " is out of bounds!" << std::endl; return neighbors; } // Use node2edge and edge2node connectivities to find neighbors of node @@ -104,8 +102,6 @@ namespace gdasapp { }; /// Compute local std. dev. as a proxy of the bkg error - //int nbh = 8; // Number of closest point (horizontal) - const auto ghostView = atlas::array::make_view(geom.functionSpace().ghost()); @@ -118,75 +114,67 @@ namespace gdasapp { // Loop through variables auto h = atlas::array::make_view(xbFs["hocn"]); - for (auto & var : socaVars.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]); - - int nbz = 1; // Number of closest point in the vertical, above and below - nbz = std::min(nbz, xbFs[var].shape(1) - 1); - for (atlas::idx_t level = nbz; level < xbFs[var].shape(1) - nbz; ++level) { - oops::Log::info() << " level: " << level << std::endl; - 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; - } - - // 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]; - for (int ll = level - nbz; ll <= level + nbz; ++ll) { - if ( abs(h(nbNode, ll)) <= 0.1 ) { - continue; - } - local.push_back(bkg(nbNode, ll)); - } - } - - if (local.size() > 3) { - // 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() != 0) { - //std::cout << "stddev: " << std::sqrt(stdDev / local.size()) << std::endl; - //std::cout << "h: " << h(jnode, 0) << std::endl; - stdDevBkg(jnode, level) = rescale * std::sqrt(stdDev / (local.size() - 1)); - } - - // Extrapolate upper levels - for (int ll = 0; ll < nbz; ++ll) { - stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); - } - - /* - if (stdDevBkg(jnode, level) > 5.0 ) { - std::cout << " ------------------------- " << std::endl; - std::cout << "mean : " << mean << std::endl; - std::cout << "stdDev : " << stdDevBkg(jnode, level) - << " " << ghostView(jnode) << " " << h(jnode, level) << std::endl; - std::cout << "local : " << local << std::endl; - } - */ - } - } - } - } + for (auto & var : socaVars.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]); + + int nbz = 1; // Number of closest point in the vertical, above and below + nbz = std::min(nbz, xbFs[var].shape(1) - 1); + for (atlas::idx_t level = nbz; level < xbFs[var].shape(1) - nbz; ++level) { + oops::Log::info() << " level: " << level << std::endl; + 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; + } + + // 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]; + for (int ll = level - nbz; ll <= level + nbz; ++ll) { + if ( abs(h(nbNode, ll)) <= 0.1 ) { + continue; + } + local.push_back(bkg(nbNode, ll)); + } + } + + if (local.size() > 3) { + // 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() != 0) { + stdDevBkg(jnode, level) = rescale * std::sqrt(stdDev / (local.size() - 1)); + } + + // Extrapolate upper levels + for (int ll = 0; ll < nbz; ++ll) { + stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); + } + } + } + } + } // TODO(G): Assume that the steric balance explains 97% of ssh ... or do it properly ... maybe - + /* bkgErrFs["tocn"].haloExchange(); - - /// Smooth the fields + + /// Smooth the fields for (atlas::idx_t level = 0; level < xbFs["tocn"].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 @@ -200,35 +188,36 @@ 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 ) { + if ( abs(h(nbNode, level)) <= 0.1 ) { continue; - } - local.push_back(stdDevBkg(nbNode, level)); - } + } + local.push_back(stdDevBkg(nbNode, level)); + } if (local.size() > 2) { stdDevBkg(jnode, level) = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); - } - } + } + } } */ + /// Smooth the background error // That doesn't seem to work, the output is in [0, ~1000] // Initialize the diffusion central block if (fullConfig.has("diffusion")) { - const eckit::LocalConfiguration diffusionConfig(fullConfig, "diffusion"); - soca::ExplicitDiffusionParameters params; - params.deserialize(diffusionConfig); - oops::GeometryData geometryData(geom.functionSpace(), bkgErrFs["tocn"], true, this->getComm()); - const oops::FieldSet3D dumyXb(cycleDate, this->getComm()); - soca::ExplicitDiffusion diffuse(geometryData, socaVars, diffusionConfig, params, dumyXb, dumyXb); - diffuse.read(); - - // Smooth the field - oops::FieldSet3D dx(cycleDate, this->getComm()); - dx.deepCopy(bkgErrFs); - diffuse.multiply(dx); - bkgErrFs = dx.fieldSet(); + const eckit::LocalConfiguration diffusionConfig(fullConfig, "diffusion"); + soca::ExplicitDiffusionParameters params; + params.deserialize(diffusionConfig); + oops::GeometryData geometryData(geom.functionSpace(), bkgErrFs["tocn"], true, this->getComm()); + const oops::FieldSet3D dumyXb(cycleDate, this->getComm()); + soca::ExplicitDiffusion diffuse(geometryData, socaVars, diffusionConfig, params, dumyXb, dumyXb); + diffuse.read(); + + // Smooth the field + oops::FieldSet3D dx(cycleDate, this->getComm()); + dx.deepCopy(bkgErrFs); + diffuse.multiply(dx); + bkgErrFs = dx.fieldSet(); } // We want to write with soca, not atlas: Syncronize with soca Increment From f8a77a3ef768510d88f6ac995d84e52049e1699d Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 5 Mar 2024 14:23:37 -0500 Subject: [PATCH 07/21] norms --- bundle/CMakeLists.txt | 3 --- scripts/exgdas_global_marine_analysis_post.py | 10 +++++----- scripts/exgdas_global_marine_analysis_prep.py | 3 ++- 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/bundle/CMakeLists.txt b/bundle/CMakeLists.txt index 328d049f1..9493ac0c8 100644 --- a/bundle/CMakeLists.txt +++ b/bundle/CMakeLists.txt @@ -94,9 +94,6 @@ if(BUILD_GDASBUNDLE) option(ENABLE_FV3_JEDI_DATA "Obtain fv3-jedi test data from fv3-jedi-data repository (vs tarball)" ON) ecbuild_bundle( PROJECT fv3-jedi SOURCE "../sorc/fv3-jedi" ) -# DAML - ecbuild_bundle( PROJECT daml SOURCE "../sorc/daml" ) - # SOCA associated repositories # TODO: Move the Icepack fork to EMC github set(BUILD_ICEPACK "ON" CACHE STRING "Build the icepack library") diff --git a/scripts/exgdas_global_marine_analysis_post.py b/scripts/exgdas_global_marine_analysis_post.py index bb2f8ed45..d58ed1a80 100755 --- a/scripts/exgdas_global_marine_analysis_post.py +++ b/scripts/exgdas_global_marine_analysis_post.py @@ -65,8 +65,8 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.bkgerr_stddev.nc')]) # Copy the recentering error - #post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.ssh_recentering_error.incr.{bdate}.nc'), - # os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.recentering_error.nc')]) + # post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.ssh_recentering_error.incr.{bdate}.nc'), + # os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.recentering_error.nc')]) # Copy the ice and ocean increments post_file_list.append([os.path.join(anl_dir, 'Data', f'{domain}.3dvarfgat_pseudo.incr.{mdate}.nc'), @@ -77,9 +77,9 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}ana.nc')]) # Copy of the ssh diagnostics -#for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']: -# post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'), -# os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')]) +# for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']: +# post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'), +# os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')]) # Copy DA grid (computed for the start of the window) post_file_list.append([os.path.join(anl_dir, 'soca_gridspec.nc'), diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index 9bdbe4948..1997940e3 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -137,6 +137,7 @@ def find_clim_ens(input_date): ################################################################################ # stage ensemble members +nmem_ens = 0 if dohybvar: logger.info("---------------- Stage ensemble members") nmem_ens = int(os.getenv('NMEM_ENS')) @@ -162,7 +163,7 @@ def find_clim_ens(input_date): bkg_utils.cice_hist2fms(cice_fname, cice_fname) # Commented out while testing the parametric diagb -#else: +# else: # logger.info("---------------- Stage offline ensemble members") # ens_member_list = [] # clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None)) From cfdcc20844276bf31a66907071f36ba5b015c310 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Fri, 8 Mar 2024 16:00:00 +0000 Subject: [PATCH 08/21] wip --- parm/soca/obs/config/sst_abi_g16_l3c.yaml | 10 +++++ parm/soca/obs/config/sst_abi_g17_l3c.yaml | 10 +++++ parm/soca/obs/config/sst_ahi_h08_l3c.yaml | 10 +++++ parm/soca/obs/config/sst_avhrr_ma_l3u.yaml | 10 +++++ parm/soca/obs/config/sst_avhrr_mb_l3u.yaml | 10 +++++ parm/soca/obs/config/sst_avhrr_mc_l3u.yaml | 10 +++++ parm/soca/obs/config/sst_viirs_n20_l3u.yaml | 10 +++++ parm/soca/obs/config/sst_viirs_npp_l3u.yaml | 10 +++++ parm/soca/obsprep/obsprep_config.yaml | 1 + sorc/ioda | 2 +- sorc/land-imsproc | 2 +- sorc/land-jediincr | 2 +- sorc/ufo | 2 +- sorc/vader | 2 +- utils/soca/gdas_diagb.h | 46 ++++++++++++++++++++- 15 files changed, 131 insertions(+), 6 deletions(-) diff --git a/parm/soca/obs/config/sst_abi_g16_l3c.yaml b/parm/soca/obs/config/sst_abi_g16_l3c.yaml index 689207b14..f104f2351 100644 --- a/parm/soca/obs/config/sst_abi_g16_l3c.yaml +++ b/parm/soca/obs/config/sst_abi_g16_l3c.yaml @@ -44,3 +44,13 @@ obs filters: where: - variable: {name: GeoVaLs/distance_from_coast} minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_abi_g17_l3c.yaml b/parm/soca/obs/config/sst_abi_g17_l3c.yaml index b2c3cb692..ece84669a 100644 --- a/parm/soca/obs/config/sst_abi_g17_l3c.yaml +++ b/parm/soca/obs/config/sst_abi_g17_l3c.yaml @@ -44,3 +44,13 @@ obs filters: where: - variable: {name: GeoVaLs/distance_from_coast} minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_ahi_h08_l3c.yaml b/parm/soca/obs/config/sst_ahi_h08_l3c.yaml index 5597d53c5..b29049e8d 100644 --- a/parm/soca/obs/config/sst_ahi_h08_l3c.yaml +++ b/parm/soca/obs/config/sst_ahi_h08_l3c.yaml @@ -44,3 +44,13 @@ obs filters: where: - variable: {name: GeoVaLs/distance_from_coast} minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml b/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml index a3fbfe6ca..1c558c139 100644 --- a/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml @@ -44,3 +44,13 @@ obs filters: where: - variable: {name: GeoVaLs/distance_from_coast} minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml b/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml index 9c9f9cfbf..0d29c70c5 100644 --- a/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml @@ -44,3 +44,13 @@ obs filters: where: - variable: {name: GeoVaLs/distance_from_coast} minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml b/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml index 57e759017..2f37b3ef2 100644 --- a/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml @@ -44,3 +44,13 @@ obs filters: where: - variable: {name: GeoVaLs/distance_from_coast} minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_viirs_n20_l3u.yaml b/parm/soca/obs/config/sst_viirs_n20_l3u.yaml index 36f155576..d533ed37b 100644 --- a/parm/soca/obs/config/sst_viirs_n20_l3u.yaml +++ b/parm/soca/obs/config/sst_viirs_n20_l3u.yaml @@ -44,3 +44,13 @@ obs filters: where: - variable: {name: GeoVaLs/distance_from_coast} minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_viirs_npp_l3u.yaml b/parm/soca/obs/config/sst_viirs_npp_l3u.yaml index c1647ada8..2f6eb1476 100644 --- a/parm/soca/obs/config/sst_viirs_npp_l3u.yaml +++ b/parm/soca/obs/config/sst_viirs_npp_l3u.yaml @@ -44,3 +44,13 @@ obs filters: where: - variable: {name: GeoVaLs/distance_from_coast} minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obsprep/obsprep_config.yaml b/parm/soca/obsprep/obsprep_config.yaml index 38ca66759..6a720861a 100644 --- a/parm/soca/obsprep/obsprep_config.yaml +++ b/parm/soca/obsprep/obsprep_config.yaml @@ -19,6 +19,7 @@ observations: dmpdir regex: rads_adt_*.nc provider: RADS output file: adt_rads_all.nc4 + error ratio: 0.00001157407407407407 # [m/s] correspond to 1.0m/day window: back: 8 # look back 8 six-hourly obs dumps forward: 1 # look forward 1 six-hourly bin diff --git a/sorc/ioda b/sorc/ioda index 35f2dac24..59a03c59c 160000 --- a/sorc/ioda +++ b/sorc/ioda @@ -1 +1 @@ -Subproject commit 35f2dac24beff87406cd5cd5732814bc43496218 +Subproject commit 59a03c59c9b8cc4865c1ff104c87b01f18ed82d3 diff --git a/sorc/land-imsproc b/sorc/land-imsproc index 6373819ca..2ac7c0d90 160000 --- a/sorc/land-imsproc +++ b/sorc/land-imsproc @@ -1 +1 @@ -Subproject commit 6373819ca034d66523cb7852cd7b1b66f3f8ae07 +Subproject commit 2ac7c0d90d3ebd449ad11bce155a71d381b4534a diff --git a/sorc/land-jediincr b/sorc/land-jediincr index 2923344b3..ab0b1c655 160000 --- a/sorc/land-jediincr +++ b/sorc/land-jediincr @@ -1 +1 @@ -Subproject commit 2923344b33511d80c685c30261ad37896eb50768 +Subproject commit ab0b1c655acc7fcef4894f864a800f1785a85a37 diff --git a/sorc/ufo b/sorc/ufo index f3734dab1..f0337cdd9 160000 --- a/sorc/ufo +++ b/sorc/ufo @@ -1 +1 @@ -Subproject commit f3734dab113dcc3cfa2ac693191da24c3f93022e +Subproject commit f0337cdd96ea8447131c9ab0686430ae605ab6d8 diff --git a/sorc/vader b/sorc/vader index b8404a832..82a2c4509 160000 --- a/sorc/vader +++ b/sorc/vader @@ -1 +1 @@ -Subproject commit b8404a83208d46bbfb5add30b6561c54da9c6d82 +Subproject commit 82a2c4509e9fbf08b2c9e1cec71e8901f371bc42 diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index ec8b9de66..691d0725d 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -147,7 +147,12 @@ namespace gdasapp { } } - if (local.size() > 3) { + //Set the minimum number of points + int minn = 12; /// probably should be passed through the config + if (xbFs[var].shape(1) == 1) { + minn = 4; + } + if (local.size() >= minn) { // Mean double mean = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); @@ -230,6 +235,45 @@ namespace gdasapp { return 0; } +/* + // ----------------------------------------------------------------------------- + std::vector getHorizNeighbors(const soca::Geometry) { + + /// Create the mesh connectivity (Copy/paste of Francois's stuff) + // Build edges, then connections between nodes and edges + atlas::functionspace::NodeColumns test = geom.functionSpace(); + atlas::Mesh mesh = test.mesh(); + atlas::mesh::actions::build_edges(mesh); + atlas::mesh::actions::build_node_to_edge_connectivity(mesh); + const auto & node2edge = mesh.nodes().edge_connectivity(); + const auto & edge2node = mesh.edges().node_connectivity(); + + oops::Log::info() << "====================== get neighbors" << std::endl; + // Lambda function to get the neighbors of a node + const auto get_neighbors_of_node = [&](const int node) { + std::vector neighbors{}; + neighbors.reserve(8); + 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; + }; +*/ + // ----------------------------------------------------------------------------- private: std::string appname() const { From 6c67e1f0699e59df4274f3f2e10a3f0620cf55a1 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Fri, 8 Mar 2024 21:37:22 +0000 Subject: [PATCH 09/21] wip --- parm/soca/obs/config/adt_rads_all.yaml | 2 +- parm/soca/obsprep/obsprep_config.yaml | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/parm/soca/obs/config/adt_rads_all.yaml b/parm/soca/obs/config/adt_rads_all.yaml index fe41a50da..74469902b 100644 --- a/parm/soca/obs/config/adt_rads_all.yaml +++ b/parm/soca/obs/config/adt_rads_all.yaml @@ -41,7 +41,7 @@ obs filters: variables: [GeoVaLs/mesoscale_representation_error, ObsError/absoluteDynamicTopography] coefs: [0.1, - 0.5] + 1.0] - filter: Domain Check where: - variable: { name: GeoVaLs/sea_ice_area_fraction} diff --git a/parm/soca/obsprep/obsprep_config.yaml b/parm/soca/obsprep/obsprep_config.yaml index 6a720861a..456ac17cc 100644 --- a/parm/soca/obsprep/obsprep_config.yaml +++ b/parm/soca/obsprep/obsprep_config.yaml @@ -21,8 +21,8 @@ observations: output file: adt_rads_all.nc4 error ratio: 0.00001157407407407407 # [m/s] correspond to 1.0m/day window: - back: 8 # look back 8 six-hourly obs dumps - forward: 1 # look forward 1 six-hourly bin + back: 0 # look back 8 six-hourly obs dumps + forward: 0 # look forward 1 six-hourly bin # Ice concentration - obs space: From 7fb6d22c66c5c55d6452b24ef6509e0700f77a41 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Mon, 11 Mar 2024 17:33:38 +0000 Subject: [PATCH 10/21] wip --- utils/soca/gdas_diagb.h | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index 691d0725d..7eb9c6892 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -123,9 +123,10 @@ namespace gdasapp { auto bkg = atlas::array::make_view(xbFs[var]); auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); - int nbz = 1; // Number of closest point in the vertical, above and below + int nbz = 2; // Number of closest point in the vertical, above and below + int nzMld = 10; // Vertical index for the mixed layer depth nbz = std::min(nbz, xbFs[var].shape(1) - 1); - for (atlas::idx_t level = nbz; level < xbFs[var].shape(1) - nbz; ++level) { + for (atlas::idx_t level = 0; level < xbFs[var].shape(1) - nbz; ++level) { oops::Log::info() << " level: " << level << std::endl; for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { // Early exit if thickness is 0 or on a ghost cell @@ -139,7 +140,13 @@ namespace gdasapp { int nbh = neighbors.size(); for (int nn = 0; nn < neighbors.size(); ++nn) { int nbNode = neighbors[nn]; - for (int ll = level - nbz; ll <= level + nbz; ++ll) { + int levelMin = level - nbz; + int levelMax = level + nbz; + if (level < nzMld) { + int levelMin = 0; + int levelMax = nzMld; + } + for (int ll = levelMin; ll <= levelMax; ++ll) { if ( abs(h(nbNode, ll)) <= 0.1 ) { continue; } @@ -166,8 +173,10 @@ namespace gdasapp { } // Extrapolate upper levels - for (int ll = 0; ll < nbz; ++ll) { - stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); + double meanMld = std::accumulate(local.begin(), local.begin() + nzMld, 0.0) / nzMld; + for (int ll = 0; ll < nzMld; ++ll) { + //stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); + stdDevBkg(jnode, ll) = meanMld; } } } From 54914eb6d120382211c04a22c9ad59b6082ea190 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 12 Mar 2024 19:29:13 -0400 Subject: [PATCH 11/21] working halo update --- parm/soca/berror/soca_diagb.yaml | 3 ++ utils/soca/gdas_diagb.h | 81 ++++++++++++++++++++------------ 2 files changed, 53 insertions(+), 31 deletions(-) diff --git a/parm/soca/berror/soca_diagb.yaml b/parm/soca/berror/soca_diagb.yaml index f8c70b545..c4baafcb3 100644 --- a/parm/soca/berror/soca_diagb.yaml +++ b/parm/soca/berror/soca_diagb.yaml @@ -21,3 +21,6 @@ variables: name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon] rescale: 2.0 +smoothing iterations: 1 +number of halo points: 5 +number of neighbors: 16 diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index 7eb9c6892..25b4119da 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -9,7 +9,9 @@ #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" @@ -69,18 +71,22 @@ namespace gdasapp { /// Create the mesh connectivity (Copy/paste of Francois's stuff) // Build edges, then connections between nodes and edges - atlas::functionspace::NodeColumns test = geom.functionSpace(); - atlas::Mesh mesh = test.mesh(); + int nbHalos(2); + fullConfig.get("number of halo points", nbHalos); + int nbNeighbors(4); + fullConfig.get("number of halo points", nbHalos); + 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(); - oops::Log::info() << "====================== get neighbors" << std::endl; - // 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{}; - neighbors.reserve(8); + neighbors.reserve(nbNeighbors); if (node >= mesh.nodes().size()) { return neighbors; } @@ -175,8 +181,8 @@ namespace gdasapp { // Extrapolate upper levels double meanMld = std::accumulate(local.begin(), local.begin() + nzMld, 0.0) / nzMld; for (int ll = 0; ll < nzMld; ++ll) { - //stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); - stdDevBkg(jnode, ll) = meanMld; + stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); + //stdDevBkg(jnode, ll) = meanMld; } } } @@ -185,35 +191,48 @@ namespace gdasapp { // TODO(G): Assume that the steric balance explains 97% of ssh ... or do it properly ... maybe - /* - bkgErrFs["tocn"].haloExchange(); - /// Smooth the fields - for (atlas::idx_t level = 0; level < xbFs["tocn"].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, 0)) <= 0.1) { - continue; - } + int niter(0); + fullConfig.get("smoothing iterations", niter); + for (auto & var : socaVars.variables()) { + // Skip the layer thickness variable + if (var == "hocn") { + 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 ) { + 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, 0)) <= 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(); + // 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(); + } + } } } } - */ /// Smooth the background error // That doesn't seem to work, the output is in [0, ~1000] @@ -244,7 +263,7 @@ namespace gdasapp { return 0; } -/* + /* // ----------------------------------------------------------------------------- std::vector getHorizNeighbors(const soca::Geometry) { @@ -281,7 +300,7 @@ namespace gdasapp { } return neighbors; }; -*/ + */ // ----------------------------------------------------------------------------- private: From d8007eeab14eb4d62f8b2a423ac83e79e7a76dc9 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Wed, 20 Mar 2024 14:26:35 +0000 Subject: [PATCH 12/21] more options --- parm/soca/berror/soca_diagb.yaml | 10 +- utils/soca/gdas_diagb.h | 228 ++++++++++++++++++++----------- 2 files changed, 156 insertions(+), 82 deletions(-) diff --git a/parm/soca/berror/soca_diagb.yaml b/parm/soca/berror/soca_diagb.yaml index c4baafcb3..f33bcf05d 100644 --- a/parm/soca/berror/soca_diagb.yaml +++ b/parm/soca/berror/soca_diagb.yaml @@ -20,7 +20,11 @@ background error: variables: name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon] -rescale: 2.0 -smoothing iterations: 1 -number of halo points: 5 +rescale: 2.0 # rescales the filtered std. dev. by "rescale" +min sst: 0.25 # Added to sst bkg. err. +max ssh: 0.0 # Limits the amplitude of the unbalanced bkg err +number of halo points: 8 number of neighbors: 16 +simple smoothing: + horizontal iterations: 5 + vertical iterations: 2 diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index 25b4119da..c3579a0ad 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -25,6 +25,7 @@ #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 "soca/Geometry/Geometry.h" @@ -49,19 +50,18 @@ namespace gdasapp { 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 socaVars(fullConfig, "variables.name"); - /// Read rescaling factor - double rescale; - fullConfig.get("rescale", rescale); - /// Read the background + // -------------------- oops::Log::info() << "====================== read bkg" << std::endl; soca::State xb(geom, socaVars, cycleDate); const eckit::LocalConfiguration bkgConfig(fullConfig, "background"); @@ -70,11 +70,12 @@ namespace gdasapp { xb.toFieldSet(xbFs); /// 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 halo points", nbHalos); + fullConfig.get("number of neighbors", nbNeighbors); atlas::functionspace::NodeColumns nodeColumns = geom.functionSpace(); atlas::Mesh mesh = nodeColumns.mesh(); atlas::mesh::actions::build_edges(mesh); @@ -108,6 +109,8 @@ namespace gdasapp { }; /// Compute local std. dev. as a proxy of the bkg error + // ---------------------------------------------------- + // Identify halo nodes const auto ghostView = atlas::array::make_view(geom.functionSpace().ghost()); @@ -117,9 +120,17 @@ namespace gdasapp { atlas::FieldSet bkgErrFs; bkgErr.toFieldSet(bkgErrFs); - // Loop through variables + // 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); + auto h = atlas::array::make_view(xbFs["hocn"]); + // Loop through variables for (auto & var : socaVars.variables()) { // Skip the layer thickness variable if (var == "hocn") { @@ -129,113 +140,167 @@ namespace gdasapp { auto bkg = atlas::array::make_view(xbFs[var]); auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); - int nbz = 2; // Number of closest point in the vertical, above and below - int nzMld = 10; // Vertical index for the mixed layer depth - nbz = std::min(nbz, xbFs[var].shape(1) - 1); - for (atlas::idx_t level = 0; level < xbFs[var].shape(1) - nbz; ++level) { - oops::Log::info() << " level: " << level << std::endl; - 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; - } + // 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(); - // Ocean or ice node, do something + // 2D case + if (xbFs[var].shape(1) == 1) { 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]; - int levelMin = level - nbz; - int levelMax = level + nbz; - if (level < nzMld) { - int levelMin = 0; - int levelMax = nzMld; - } - for (int ll = levelMin; ll <= levelMax; ++ll) { - if ( abs(h(nbNode, ll)) <= 0.1 ) { - continue; - } - local.push_back(bkg(nbNode, ll)); + if ( abs(h(neighbors[nn], 0)) <= 0.1 ) { + continue; } + local.push_back(bkg(neighbors[nn], 0)); } - //Set the minimum number of points - int minn = 12; /// probably should be passed through the config - if (xbFs[var].shape(1) == 1) { - minn = 4; - } + 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() != 0) { - stdDevBkg(jnode, level) = rescale * std::sqrt(stdDev / (local.size() - 1)); - } - // Extrapolate upper levels - double meanMld = std::accumulate(local.begin(), local.begin() + nzMld, 0.0) / nzMld; - for (int ll = 0; ll < nzMld; ++ll) { - stdDevBkg(jnode, ll) = stdDevBkg(jnode, nbz); - //stdDevBkg(jnode, ll) = meanMld; + 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 = 10; // Vertical index for the mixed layer depth, hard-wired to the MOM6 UFS config + 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 water column + levelMin = 0; + levelMax = 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 = 12; /// 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 - int niter(0); - fullConfig.get("smoothing iterations", niter); - for (auto & var : socaVars.variables()) { - // Skip the layer thickness variable - if (var == "hocn") { - continue; - } - for (int iter = 0; iter < niter; ++iter) { + /// Smooth the fields + // ------------------ + if (fullConfig.has("simple smoothing")) { + int niter(0); + fullConfig.get("simple smoothing.horizontal iterations", niter); + for (auto & var : socaVars.variables()) { + // Skip the layer thickness variable + if (var == "hocn") { + continue; + } - // Update the halo points - nodeColumns.haloExchange(bkgErrFs[var]); - auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); + // Horizontal averaging + for (int iter = 0; iter < niter; ++iter) { - // 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, 0)) <= 0.1) { - continue; - } + // Update the halo points + nodeColumns.haloExchange(bkgErrFs[var]); + auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); - // 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 ) { + // 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; } - local.push_back(stdDevBkg(nbNode, level)); + + // 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(); + } } + } + } - 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); } } } } } - /// Smooth the background error - // That doesn't seem to work, the output is in [0, ~1000] + /// 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"); @@ -253,6 +318,11 @@ namespace gdasapp { bkgErrFs = dx.fieldSet(); } + // Rescale + double rescale; + fullConfig.get("rescale", rescale); + util::multiplyFieldSet(bkgErrFs, rescale); + // We want to write with soca, not atlas: Syncronize with soca Increment bkgErr.fromFieldSet(bkgErrFs); From 92a119774d6abbd8bfe1c6f914dabfa497954b37 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Fri, 29 Mar 2024 13:19:51 +0000 Subject: [PATCH 13/21] wip, unrealistic near-shore analysis --- parm/soca/berror/soca_diagb.yaml | 10 +-- .../berror/soca_parameters_diffusion_vt.yaml | 6 +- parm/soca/berror/soca_vtscales.yaml | 13 +++ parm/soca/obs/config/sst_abi_g16_l3c.yaml | 8 +- parm/soca/obs/config/sst_abi_g17_l3c.yaml | 8 +- parm/soca/obs/config/sst_ahi_h08_l3c.yaml | 8 +- parm/soca/obs/config/sst_avhrr_ma_l3u.yaml | 8 +- parm/soca/obs/config/sst_avhrr_mb_l3u.yaml | 8 +- parm/soca/obs/config/sst_avhrr_mc_l3u.yaml | 8 +- parm/soca/obs/config/sst_viirs_n20_l3u.yaml | 8 +- parm/soca/obs/config/sst_viirs_npp_l3u.yaml | 8 +- scripts/exgdas_global_marine_analysis_bmat.sh | 9 +- scripts/exgdas_global_marine_analysis_prep.py | 8 ++ utils/soca/gdas_diagb.h | 85 +++++++++---------- 14 files changed, 112 insertions(+), 83 deletions(-) create mode 100644 parm/soca/berror/soca_vtscales.yaml diff --git a/parm/soca/berror/soca_diagb.yaml b/parm/soca/berror/soca_diagb.yaml index f33bcf05d..27ce039c5 100644 --- a/parm/soca/berror/soca_diagb.yaml +++ b/parm/soca/berror/soca_diagb.yaml @@ -18,13 +18,13 @@ background error: type: incr variables: - name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon] + name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon, mom6_mld] rescale: 2.0 # rescales the filtered std. dev. by "rescale" -min sst: 0.25 # Added to sst bkg. err. +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: 8 +number of halo points: 4 number of neighbors: 16 simple smoothing: - horizontal iterations: 5 - vertical iterations: 2 + horizontal iterations: 2 + vertical iterations: 1 diff --git a/parm/soca/berror/soca_parameters_diffusion_vt.yaml b/parm/soca/berror/soca_parameters_diffusion_vt.yaml index 8597bfe15..0c63938c9 100644 --- a/parm/soca/berror/soca_parameters_diffusion_vt.yaml +++ b/parm/soca/berror/soca_parameters_diffusion_vt.yaml @@ -25,6 +25,10 @@ background error: - name: vt_ocean vertical: as gaussian: true - fixed value: 5.0 + from file: + filename: vt_scales.nc + variable name: vt + #as gaussian: true + #fixed value: 5.0 write: filename: vt_ocean.nc diff --git a/parm/soca/berror/soca_vtscales.yaml b/parm/soca/berror/soca_vtscales.yaml new file mode 100644 index 000000000..aa8975c89 --- /dev/null +++ b/parm/soca/berror/soca_vtscales.yaml @@ -0,0 +1,13 @@ +gridspec_filename: soca_gridspec.nc +restart_filename: ./INPUT/MOM.res.nc +mld_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' +output_filename: ./vt_scales.nc +output_variable_vt: vt +output_variable_hz: hz + +VT_MIN: 5 +VT_MAX: 15 + +HZ_ROSSBY_MULT: 1.0 +HZ_MAX: 200e3 +HZ_MIN_GRID_MULT: 2.0 diff --git a/parm/soca/obs/config/sst_abi_g16_l3c.yaml b/parm/soca/obs/config/sst_abi_g16_l3c.yaml index f104f2351..b9d999843 100644 --- a/parm/soca/obs/config/sst_abi_g16_l3c.yaml +++ b/parm/soca/obs/config/sst_abi_g16_l3c.yaml @@ -40,10 +40,10 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 - filter: Perform Action action: name: assign error diff --git a/parm/soca/obs/config/sst_abi_g17_l3c.yaml b/parm/soca/obs/config/sst_abi_g17_l3c.yaml index ece84669a..a7ada5774 100644 --- a/parm/soca/obs/config/sst_abi_g17_l3c.yaml +++ b/parm/soca/obs/config/sst_abi_g17_l3c.yaml @@ -40,10 +40,10 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 - filter: Perform Action action: name: assign error diff --git a/parm/soca/obs/config/sst_ahi_h08_l3c.yaml b/parm/soca/obs/config/sst_ahi_h08_l3c.yaml index b29049e8d..ed0efd721 100644 --- a/parm/soca/obs/config/sst_ahi_h08_l3c.yaml +++ b/parm/soca/obs/config/sst_ahi_h08_l3c.yaml @@ -40,10 +40,10 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 - filter: Perform Action action: name: assign error diff --git a/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml b/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml index 1c558c139..73c048891 100644 --- a/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml @@ -40,10 +40,10 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 - filter: Perform Action action: name: assign error diff --git a/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml b/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml index 0d29c70c5..ccac3cc4e 100644 --- a/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml @@ -40,10 +40,10 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 - filter: Perform Action action: name: assign error diff --git a/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml b/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml index 2f37b3ef2..9c34bcd56 100644 --- a/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml @@ -40,10 +40,10 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 - filter: Perform Action action: name: assign error diff --git a/parm/soca/obs/config/sst_viirs_n20_l3u.yaml b/parm/soca/obs/config/sst_viirs_n20_l3u.yaml index d533ed37b..04021235f 100644 --- a/parm/soca/obs/config/sst_viirs_n20_l3u.yaml +++ b/parm/soca/obs/config/sst_viirs_n20_l3u.yaml @@ -40,10 +40,10 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 - filter: Perform Action action: name: assign error diff --git a/parm/soca/obs/config/sst_viirs_npp_l3u.yaml b/parm/soca/obs/config/sst_viirs_npp_l3u.yaml index 2f6eb1476..6df13028a 100644 --- a/parm/soca/obs/config/sst_viirs_npp_l3u.yaml +++ b/parm/soca/obs/config/sst_viirs_npp_l3u.yaml @@ -40,10 +40,10 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 - filter: Perform Action action: name: assign error diff --git a/scripts/exgdas_global_marine_analysis_bmat.sh b/scripts/exgdas_global_marine_analysis_bmat.sh index ec4facc16..ce8eec6da 100755 --- a/scripts/exgdas_global_marine_analysis_bmat.sh +++ b/scripts/exgdas_global_marine_analysis_bmat.sh @@ -88,7 +88,6 @@ fi ################################################################################ # Compute the parametric diag of B -echo "++++++++++++++++++++++++++++++ $APRUN_OCNANAL" $APRUN_OCNANAL $JEDI_BIN/gdas_diagb.x soca_diagb.yaml export err=$?; err_chk if [ $err -gt 0 ]; then @@ -104,7 +103,7 @@ if [ $err -gt 0 ]; then fi ################################################################################ -# Initialize diffusion blocks +# Initialize the horizontal diffusion blocks clean_yaml soca_parameters_diffusion_hz.yaml $APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_hz.yaml export err=$?; err_chk @@ -112,7 +111,13 @@ if [ $err -gt 0 ]; then exit $err fi +################################################################################ +# Initialize the vertical diffusion blocks +# Compute the MLD dependent vertical scales +python ${HOMEgfs}/sorc/gdas.cd/sorc/soca/tools/calc_scales.py soca_vtscales.yaml clean_yaml soca_parameters_diffusion_vt.yaml + +# Vertical fiffusion and normalization $APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_vt.yaml export err=$?; err_chk if [ $err -gt 0 ]; then diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index 1997940e3..2dafbbddd 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -244,6 +244,14 @@ def find_clim_ens(input_date): config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) config.save(diffu_hz_yaml) +logger.info(f"---------------- generate soca_vtscales.yaml") +vtscales_yaml = os.path.join(anl_dir, 'soca_vtscales.yaml') +vtscales_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror') +vtscales_yaml_template = os.path.join(berror_yaml_dir, 'soca_vtscales.yaml') +config = YAMLFile(path=vtscales_yaml_template) +config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) +config.save(vtscales_yaml) + logger.info(f"---------------- generate soca_parameters_diffusion_vt.yaml") diffu_vt_yaml = os.path.join(anl_dir, 'soca_parameters_diffusion_vt.yaml') diffu_vt_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror') diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index c3579a0ad..b3f5dc9c8 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -35,6 +35,15 @@ #include "soca/ExplicitDiffusion/ExplicitDiffusionParameters.h" namespace gdasapp { + /** + * DiagB Class Implementation + * + * Implements variance partitioning within the GDAS for the ocean and sea-ice. It is used as a proxy for the + * diagonal of B. + * This class is designed to partition the variance of ocean and sea-ice fields by analyzing the variance within + * predefined geographical bins. + * This approach allows for an ensemble-free estimate of a flow-dependent background error. + */ class DiagB : public oops::Application { public: @@ -68,6 +77,8 @@ namespace gdasapp { 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) // -------------------------------------------------------------- @@ -128,7 +139,34 @@ 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 : socaVars.variables()) { @@ -183,16 +221,16 @@ namespace gdasapp { } else { // 3D case int nbz = 1; // Number of closest point in the vertical, above and below - int nzMld = 10; // Vertical index for the mixed layer depth, hard-wired to the MOM6 UFS config + 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 water column + // If in the MLD, compute the std. dev. throughout the MLD levelMin = 0; - levelMax = nzMld; + levelMax = 1; //nzMld; } for (int ll = levelMin; ll <= levelMax; ++ll) { if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { @@ -202,7 +240,7 @@ namespace gdasapp { } } //Set the minimum number of points - int minn = 12; /// 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(); @@ -333,45 +371,6 @@ namespace gdasapp { return 0; } - /* - // ----------------------------------------------------------------------------- - std::vector getHorizNeighbors(const soca::Geometry) { - - /// Create the mesh connectivity (Copy/paste of Francois's stuff) - // Build edges, then connections between nodes and edges - atlas::functionspace::NodeColumns test = geom.functionSpace(); - atlas::Mesh mesh = test.mesh(); - atlas::mesh::actions::build_edges(mesh); - atlas::mesh::actions::build_node_to_edge_connectivity(mesh); - const auto & node2edge = mesh.nodes().edge_connectivity(); - const auto & edge2node = mesh.edges().node_connectivity(); - - oops::Log::info() << "====================== get neighbors" << std::endl; - // Lambda function to get the neighbors of a node - const auto get_neighbors_of_node = [&](const int node) { - std::vector neighbors{}; - neighbors.reserve(8); - 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; - }; - */ - // ----------------------------------------------------------------------------- private: std::string appname() const { From 284930a624ef832961a98d2179d21f8c18c1bc09 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Mon, 1 Apr 2024 13:37:58 +0000 Subject: [PATCH 14/21] wip --- utils/soca/gdas_diagb.h | 54 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 51 insertions(+), 3 deletions(-) diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index b3f5dc9c8..2a7cdf55b 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -45,6 +45,50 @@ namespace gdasapp { * This approach allows for an ensemble-free estimate of a flow-dependent background error. */ + void stdDevFilt(const int level, + const int nbz, + const std::vector neighbors, + const int nzMld, + double stdDev) { + + 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; + } + } + } + class DiagB : public oops::Application { public: explicit DiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) @@ -167,7 +211,6 @@ namespace gdasapp { viewMldindex(jnode, 0) = std::distance(testMld.begin(), minVal); } - // Loop through variables for (auto & var : socaVars.variables()) { // Skip the layer thickness variable @@ -288,8 +331,8 @@ 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) { - // Early exit if thickness is 0 or on a ghost cell - if (abs(h(jnode, level)) <= 0.1) { + // Early exit if on a ghost cell + if (ghostView(jnode) > 0) { continue; } @@ -308,6 +351,11 @@ namespace gdasapp { if (local.size() > 2) { stdDevBkg(jnode, level) = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); } + + // Reset to 0 over land + if (abs(h(jnode, level)) <= 0.1) { + stdDevBkg(jnode, level) = 0.0; + } } } } From 80d9e48d6d24379431a0fdc2b6e34da2a311f870 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 2 Apr 2024 01:58:42 +0000 Subject: [PATCH 15/21] std dev fctn --- utils/soca/gdas_diagb.h | 170 ++++++++++++++++++++-------------------- 1 file changed, 83 insertions(+), 87 deletions(-) diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index 2a7cdf55b..e7efb24b3 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -45,56 +45,66 @@ namespace gdasapp { * This approach allows for an ensemble-free estimate of a flow-dependent background error. */ - void stdDevFilt(const int level, - const int nbz, - const std::vector neighbors, - const int nzMld, - double stdDev) { - - 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; - } - } - } - + // ----------------------------------------------------------------------------- class DiagB : public oops::Application { public: explicit DiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) : Application(comm) {} static const std::string classname() {return "gdasapp::DiagB";} + // ----------------------------------------------------------------------------- + void stdDevFilt(const int jnode, + const int level, + const int nbz, + const double depthMin, + const std::vector neighbors, + const int nzMld, + const atlas::array::ArrayView h, + const atlas::array::ArrayView bkg, + const atlas::array::ArrayView bathy, + atlas::array::ArrayView& stdDevBkg) const { + + // Early exit if too shallow + if ( bathy(jnode, 0) < depthMin ) { + stdDevBkg(jnode, level) = 0.0; + return; + } + + int nbh = neighbors.size(); + 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); + } + if (stdDev > 0.0 || local.size() > 2) { + stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)); + } + } + } + + // ----------------------------------------------------------------------------- int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { /// Setup the soca geometry oops::Log::info() << "====================== geometry" << std::endl; @@ -127,6 +137,7 @@ namespace gdasapp { /// Create the mesh connectivity (Copy/paste of Francois's stuff) // -------------------------------------------------------------- // Build edges, then connections between nodes and edges + oops::Log::info() << "====================== build mesh connectivity" << std::endl; int nbHalos(2); fullConfig.get("number of halo points", nbHalos); int nbNeighbors(4); @@ -165,11 +176,13 @@ namespace gdasapp { /// Compute local std. dev. as a proxy of the bkg error // ---------------------------------------------------- + oops::Log::info() << "====================== start variance partitioning" << std::endl; // Identify halo nodes const auto ghostView = atlas::array::make_view(geom.functionSpace().ghost()); // Create the background error fieldset + oops::Log::info() << "====================== allocate std. dev. field set" << std::endl; soca::Increment bkgErr(geom, socaVars, cycleDate); bkgErr.zero(); atlas::FieldSet bkgErrFs; @@ -183,19 +196,35 @@ namespace gdasapp { double sshMax(0.0); fullConfig.get("max ssh", sshMax); - // Get the layer thicknesses and convert to depth + // Get the minimum depth for which to partition the 3D field's std. dev. + double depthMin(0.0); + fullConfig.get("min depth", depthMin); + + // Get the layer thicknesses and convert to layer depth + oops::Log::info() << "====================== calculate layer depth" << std::endl; 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) { + oops::Log::info() << depth.shape() << std::endl; + oops::Log::info() << viewDepth.shape(0) << "-" << viewDepth.shape(1) << std::endl; + for (atlas::idx_t jnode = 0; jnode < depth.shape(0); ++jnode) { viewDepth(jnode, 0) = 0.5 * h(jnode, 0); - for (atlas::idx_t level = 1; level < h.shape(1); ++level) { + for (atlas::idx_t level = 1; level < depth.shape(1); ++level) { viewDepth(jnode, level) = viewDepth(jnode, level-1) + 0.5 * (h(jnode, level-1 ) + h(jnode, level)); } } + // Compute the bathymetry + oops::Log::info() << "====================== calculate bathymetry" << std::endl; + atlas::array::ArrayT bathy(h.shape(0), 1); + auto viewBathy = atlas::array::make_view(bathy); + for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { + for (atlas::idx_t level = 0; level < h.shape(1); ++level) { + viewBathy(jnode, 0) += 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)); @@ -266,45 +295,12 @@ namespace gdasapp { 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 + // Std. dev. of the partition + stdDevFilt(jnode, level, nbz, depthMin, neighbors, nzMld, h, bkg, viewBathy, stdDevBkg); + } } // end 3D case } // end jnode + this->getComm().barrier(); } // end var // TODO(G): Assume that the steric balance explains 97% of ssh ... or do it properly ... maybe @@ -403,6 +399,7 @@ namespace gdasapp { diffuse.multiply(dx); bkgErrFs = dx.fieldSet(); } + this->getComm().barrier(); // Rescale double rescale; @@ -420,11 +417,10 @@ namespace gdasapp { } // ----------------------------------------------------------------------------- - private: + private: std::string appname() const { return "gdasapp::DiagB"; } // ----------------------------------------------------------------------------- }; - } // namespace gdasapp From 6b421127c8e72a79da114053fc1b5f04648a1ebd Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 2 Apr 2024 15:45:41 -0500 Subject: [PATCH 16/21] fixed pe with 0 variance --- utils/soca/gdas_diagb.h | 49 +++++++++++++++++++++++------------------ 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index e7efb24b3..aa8818d3e 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -62,10 +62,12 @@ namespace gdasapp { const atlas::array::ArrayView h, const atlas::array::ArrayView bkg, const atlas::array::ArrayView bathy, - atlas::array::ArrayView& stdDevBkg) const { - + atlas::array::ArrayView& stdDevBkg, + bool doBathy = true) const { + // Early exit if too shallow - if ( bathy(jnode, 0) < depthMin ) { + + if ( doBathy & bathy(jnode, 0) < depthMin ) { stdDevBkg(jnode, level) = 0.0; return; } @@ -78,7 +80,7 @@ namespace gdasapp { if (level < nzMld) { // If in the MLD, compute the std. dev. throughout the MLD levelMin = 0; - levelMax = 1; //nzMld; + levelMax = 1; //nzMld; Projecting the surface variance down the water column for now } for (int ll = levelMin; ll <= levelMax; ++ll) { if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { @@ -202,26 +204,25 @@ namespace gdasapp { // Get the layer thicknesses and convert to layer depth oops::Log::info() << "====================== calculate layer depth" << std::endl; - auto h = atlas::array::make_view(xbFs["hocn"]); - atlas::array::ArrayT depth(h.shape(0), h.shape(1)); + auto viewHocn = atlas::array::make_view(xbFs["hocn"]); + atlas::array::ArrayT depth(viewHocn.shape(0), viewHocn.shape(1)); auto viewDepth = atlas::array::make_view(depth); - oops::Log::info() << depth.shape() << std::endl; - oops::Log::info() << viewDepth.shape(0) << "-" << viewDepth.shape(1) << std::endl; for (atlas::idx_t jnode = 0; jnode < depth.shape(0); ++jnode) { - viewDepth(jnode, 0) = 0.5 * h(jnode, 0); + viewDepth(jnode, 0) = 0.5 * viewHocn(jnode, 0); for (atlas::idx_t level = 1; level < depth.shape(1); ++level) { viewDepth(jnode, level) = viewDepth(jnode, level-1) + - 0.5 * (h(jnode, level-1 ) + h(jnode, level)); + 0.5 * (viewHocn(jnode, level-1 ) + viewHocn(jnode, level)); } } // Compute the bathymetry oops::Log::info() << "====================== calculate bathymetry" << std::endl; - atlas::array::ArrayT bathy(h.shape(0), 1); + atlas::array::ArrayT bathy(viewHocn.shape(0), 1); auto viewBathy = atlas::array::make_view(bathy); - for (atlas::idx_t jnode = 0; jnode < h.shape(0); ++jnode) { - for (atlas::idx_t level = 0; level < h.shape(1); ++level) { - viewBathy(jnode, 0) += h(jnode, level); + for (atlas::idx_t jnode = 0; jnode < viewHocn.shape(0); ++jnode) { + viewBathy(jnode, 0) = 0.0; + for (atlas::idx_t level = 0; level < viewHocn.shape(1); ++level) { + viewBathy(jnode, 0) += viewHocn(jnode, level); } } @@ -230,7 +231,7 @@ namespace gdasapp { 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) { + for (atlas::idx_t jnode = 0; jnode < viewHocn.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))); @@ -240,8 +241,14 @@ namespace gdasapp { viewMldindex(jnode, 0) = std::distance(testMld.begin(), minVal); } + // Update the layer thickness halo + nodeColumns.haloExchange(xbFs["hocn"]); + // Loop through variables for (auto & var : socaVars.variables()) { + // Update the halo + nodeColumns.haloExchange(xbFs[var]); + // Skip the layer thickness variable if (var == "hocn") { continue; @@ -253,7 +260,7 @@ 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) { + if (ghostView(jnode) > 0 || abs(viewHocn(jnode, 0)) <= 0.1) { continue; } @@ -265,7 +272,7 @@ 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 ) { + if ( abs(viewHocn(neighbors[nn], 0)) <= 0.1 ) { continue; } local.push_back(bkg(neighbors[nn], 0)); @@ -296,11 +303,10 @@ namespace gdasapp { int nzMld = std::max(10, viewMldindex(jnode, 0)); for (atlas::idx_t level = 0; level < xbFs[var].shape(1) - nbz; ++level) { // Std. dev. of the partition - stdDevFilt(jnode, level, nbz, depthMin, neighbors, nzMld, h, bkg, viewBathy, stdDevBkg); + stdDevFilt(jnode, level, nbz, depthMin, neighbors, nzMld, viewHocn, bkg, viewBathy, stdDevBkg); } } // end 3D case } // end jnode - this->getComm().barrier(); } // end var // TODO(G): Assume that the steric balance explains 97% of ssh ... or do it properly ... maybe @@ -338,7 +344,7 @@ 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 ) { + if ( abs(viewHocn(nbNode, level)) <= 0.1 ) { continue; } local.push_back(stdDevBkg(nbNode, level)); @@ -349,7 +355,7 @@ namespace gdasapp { } // Reset to 0 over land - if (abs(h(jnode, level)) <= 0.1) { + if (abs(viewHocn(jnode, level)) <= 0.1) { stdDevBkg(jnode, level) = 0.0; } } @@ -399,7 +405,6 @@ namespace gdasapp { diffuse.multiply(dx); bkgErrFs = dx.fieldSet(); } - this->getComm().barrier(); // Rescale double rescale; From 8a188118c1c81e302f1edf878a6494f7eac91360 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Wed, 3 Apr 2024 09:10:43 -0500 Subject: [PATCH 17/21] code tidy --- test/soca/testinput/gdas_diagb.yaml | 40 ++++++ utils/soca/gdas_diagb.h | 189 ++++++++++++++++------------ 2 files changed, 146 insertions(+), 83 deletions(-) diff --git a/test/soca/testinput/gdas_diagb.yaml b/test/soca/testinput/gdas_diagb.yaml index e865db768..24c9bb30c 100644 --- a/test/soca/testinput/gdas_diagb.yaml +++ b/test/soca/testinput/gdas_diagb.yaml @@ -19,3 +19,43 @@ background error: variables: name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon] + +rescale: 2.0 +min sst: 0.0 +max ssh: 0.0 +min depth: 300.0 +number of halo points: 4 +number of neighbors: 16 + +simple smoothing: + horizontal iterations: 1 + vertical iterations: 1 + +#diffusion: +# saber block name: EXPLICIT_DIFFUSION +# active variables: [tocn, socn, ssh, cicen, hicen, hsnon] +# geometry: +# mom6_input_nml: mom_input.nml +# fields metadata: ./fields_metadata.yaml +# group mapping: +# - name: ocean +# variables: +# - tocn +# - socn +# - ssh +# - name: ice +# variables: +# - cicen +# - hicen +# - hsnon +# read: +# groups: +# - name: ocean +# horizontal: +# filename: hz_ocean.nc +# vertical: +# filename: vt_ocean.nc +# - name: ice +# horizontal: +# filename: hz_ice.nc + diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_diagb.h index aa8818d3e..cce9989fb 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_diagb.h @@ -45,11 +45,77 @@ namespace gdasapp { * This approach allows for an ensemble-free estimate of a flow-dependent background error. */ + // ----------------------------------------------------------------------------- + // Since we can't have useful data members in DiagB because execute is "const" + // and the constructor does not know about the configuration, + // we're going around the issue with the use of a simple data structure + struct DiagBConfig { + util::DateTime cycleDate; + oops::Variables socaVars; + double sshMax; // poor man's alternative to limiting the unbalanced ssh variance + double depthMin; // set the bkg error to 0 for depth < depthMin + bool diffusion; // apply explicit diffusion to the std. dev. fields + bool simpleSmoothing; // Simple spatial averaging + int niter; // number of iterations for the horizontal averaging + int niterVert; // number of iterations for the vertical averaging + double rescale; // inflation/deflation of the variance after filtering + }; + // ----------------------------------------------------------------------------- class DiagB : public oops::Application { - public: + public: + // ----------------------------------------------------------------------------- explicit DiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) - : Application(comm) {} + : Application(comm) { + } + + // ----------------------------------------------------------------------------- + DiagBConfig setup(const eckit::Configuration & fullConfig) const { + + DiagBConfig diagBConfig; + + /// Get the date + // ------------- + oops::Log::info() << "====================== date" << std::endl; + std::string strdt; + fullConfig.get("date", strdt); + diagBConfig.cycleDate = util::DateTime(strdt); + + /// Get the list of variables + // -------------------------- + oops::Log::info() << "====================== variables" << std::endl; + oops::Variables vars(fullConfig, "variables.name"); + diagBConfig.socaVars = vars; + + // Get the max std dev for ssh + diagBConfig.sshMax = 0.0; + fullConfig.get("max ssh", diagBConfig.sshMax); + + // Get the minimum depth for which to partition the 3D field's std. dev. + diagBConfig.depthMin = 0.0; + fullConfig.get("min depth", diagBConfig.depthMin); + + // Explicit diffusion + diagBConfig.diffusion = false; + if (fullConfig.has("diffusion")) { + diagBConfig.diffusion = true; + } + + // Simple smoothing parameters + diagBConfig.simpleSmoothing = false; + if (fullConfig.has("simple smoothing")) { + diagBConfig.simpleSmoothing = true; + fullConfig.get("simple smoothing.horizontal iterations", diagBConfig.niter); + fullConfig.get("simple smoothing.vertical iterations", diagBConfig.niterVert); + } + + // Variance rescaling + fullConfig.get("rescale", diagBConfig.rescale); + + return diagBConfig; + } + + // ----------------------------------------------------------------------------- static const std::string classname() {return "gdasapp::DiagB";} // ----------------------------------------------------------------------------- @@ -63,10 +129,10 @@ namespace gdasapp { const atlas::array::ArrayView bkg, const atlas::array::ArrayView bathy, atlas::array::ArrayView& stdDevBkg, - bool doBathy = true) const { - - // Early exit if too shallow + bool doBathy = true, + int minn = 6) const { + // Early exit if too shallow if ( doBathy & bathy(jnode, 0) < depthMin ) { stdDevBkg(jnode, level) = 0.0; return; @@ -77,11 +143,17 @@ namespace gdasapp { for (int nn = 0; nn < neighbors.size(); ++nn) { int levelMin = std::max(0, level - nbz); int levelMax = level + nbz; + // Do weird things withing the MLD if (level < nzMld) { // If in the MLD, compute the std. dev. throughout the MLD levelMin = 0; levelMax = 1; //nzMld; Projecting the surface variance down the water column for now } + // 2D case + if (bkg.shape(1) == 1) { + levelMin = 0; + levelMax = 0; //nzMld; Projecting the surface variance down the water column for now + } for (int ll = levelMin; ll <= levelMax; ++ll) { if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { continue; @@ -89,8 +161,6 @@ namespace gdasapp { 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(); @@ -108,27 +178,19 @@ namespace gdasapp { // ----------------------------------------------------------------------------- int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { + + /// Initialize the paramters for D + DiagBConfig configD = setup(fullConfig); + /// Setup the soca geometry oops::Log::info() << "====================== geometry" << std::endl; const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); const soca::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 socaVars(fullConfig, "variables.name"); - /// Read the background // -------------------- oops::Log::info() << "====================== read bkg" << std::endl; - soca::State xb(geom, socaVars, cycleDate); + soca::State xb(geom, configD.socaVars, configD.cycleDate); const eckit::LocalConfiguration bkgConfig(fullConfig, "background"); xb.read(bkgConfig); atlas::FieldSet xbFs; @@ -185,23 +247,11 @@ namespace gdasapp { // Create the background error fieldset oops::Log::info() << "====================== allocate std. dev. field set" << std::endl; - soca::Increment bkgErr(geom, socaVars, cycleDate); + soca::Increment bkgErr(geom, configD.socaVars, configD.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 minimum depth for which to partition the 3D field's std. dev. - double depthMin(0.0); - fullConfig.get("min depth", depthMin); - // Get the layer thicknesses and convert to layer depth oops::Log::info() << "====================== calculate layer depth" << std::endl; auto viewHocn = atlas::array::make_view(xbFs["hocn"]); @@ -245,7 +295,7 @@ namespace gdasapp { nodeColumns.haloExchange(xbFs["hocn"]); // Loop through variables - for (auto & var : socaVars.variables()) { + for (auto & var : configD.socaVars.variables()) { // Update the halo nodeColumns.haloExchange(xbFs[var]); @@ -264,67 +314,42 @@ namespace gdasapp { continue; } - // get indices of neighbor cells + // get indices of the cell's neighbors auto neighbors = get_neighbors_of_node(jnode); - int nbh = neighbors.size(); - // 2D case + // 2D case if (xbFs[var].shape(1) == 1) { - std::vector local; - for (int nn = 0; nn < neighbors.size(); ++nn) { - if ( abs(viewHocn(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 + // Std. dev. of the partition + stdDevFilt(jnode, 0, 0, + configD.depthMin, neighbors, 0, viewHocn, bkg, viewBathy, stdDevBkg, false, 4); + if (var == "ssh") { + // TODO(G): Extract the unbalanced ssh variance, in the mean time, do this: + stdDevBkg(jnode, 0) = std::min(configD.sshMax, stdDevBkg(jnode, 0)); + } } 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. dev. of the partition - stdDevFilt(jnode, level, nbz, depthMin, neighbors, nzMld, viewHocn, bkg, viewBathy, stdDevBkg); + stdDevFilt(jnode, level, nbz, + configD.depthMin, neighbors, nzMld, viewHocn, bkg, viewBathy, stdDevBkg, true, 6); } } // 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 : socaVars.variables()) { + if (configD.simpleSmoothing) { + for (auto & var : configD.socaVars.variables()) { // Skip the layer thickness variable if (var == "hocn") { continue; } // Horizontal averaging - for (int iter = 0; iter < niter; ++iter) { + for (int iter = 0; iter < configD.niter; ++iter) { // Update the halo points nodeColumns.haloExchange(bkgErrFs[var]); @@ -366,12 +391,9 @@ namespace gdasapp { 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 (int iter = 0; iter < configD.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) + @@ -395,21 +417,19 @@ namespace gdasapp { soca::ExplicitDiffusionParameters params; params.deserialize(diffusionConfig); oops::GeometryData geometryData(geom.functionSpace(), bkgErrFs["tocn"], true, this->getComm()); - const oops::FieldSet3D dumyXb(cycleDate, this->getComm()); - soca::ExplicitDiffusion diffuse(geometryData, socaVars, diffusionConfig, params, dumyXb, dumyXb); + const oops::FieldSet3D dumyXb(configD.cycleDate, this->getComm()); + soca::ExplicitDiffusion diffuse(geometryData, configD.socaVars, diffusionConfig, params, dumyXb, dumyXb); diffuse.read(); // Smooth the field - oops::FieldSet3D dx(cycleDate, this->getComm()); + oops::FieldSet3D dx(configD.cycleDate, this->getComm()); dx.deepCopy(bkgErrFs); diffuse.multiply(dx); bkgErrFs = dx.fieldSet(); } // Rescale - double rescale; - fullConfig.get("rescale", rescale); - util::multiplyFieldSet(bkgErrFs, rescale); + util::multiplyFieldSet(bkgErrFs, configD.rescale); // We want to write with soca, not atlas: Syncronize with soca Increment bkgErr.fromFieldSet(bkgErrFs); @@ -422,10 +442,13 @@ namespace gdasapp { } // ----------------------------------------------------------------------------- + private: std::string appname() const { return "gdasapp::DiagB"; } + // ----------------------------------------------------------------------------- + }; } // namespace gdasapp From 522d7d244b92f438899f64250199834a6e3f406d Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 16 Apr 2024 10:57:50 -0400 Subject: [PATCH 18/21] wip --- .gitmodules | 4 ---- sorc/land-jediincr | 2 +- utils/CMakeLists.txt | 2 +- utils/soca/CMakeLists.txt | 1 - 4 files changed, 2 insertions(+), 7 deletions(-) diff --git a/.gitmodules b/.gitmodules index 40b595645..f8f34c6fa 100644 --- a/.gitmodules +++ b/.gitmodules @@ -77,7 +77,3 @@ path = sorc/fms url = https://github.com/jcsda/fms.git branch = release-stable -[submodule "sorc/daml"] - path = sorc/daml - url = https://github.com/NOAA-EMC/daml.git - branch = develop diff --git a/sorc/land-jediincr b/sorc/land-jediincr index ab0b1c655..79576b79d 160000 --- a/sorc/land-jediincr +++ b/sorc/land-jediincr @@ -1 +1 @@ -Subproject commit ab0b1c655acc7fcef4894f864a800f1785a85a37 +Subproject commit 79576b79df6a3f1ab8e7dde4425821f2915006b3 diff --git a/utils/CMakeLists.txt b/utils/CMakeLists.txt index 767026d01..0195c15dc 100644 --- a/utils/CMakeLists.txt +++ b/utils/CMakeLists.txt @@ -17,7 +17,7 @@ find_package(soca REQUIRED) find_package(fv3jedi REQUIRED) add_subdirectory(soca) -#add_subdirectory(ioda_example) +add_subdirectory(ioda_example) add_subdirectory(test) add_subdirectory(obsproc) add_subdirectory(fv3jedi) diff --git a/utils/soca/CMakeLists.txt b/utils/soca/CMakeLists.txt index 963b8a6e0..78e9781f2 100644 --- a/utils/soca/CMakeLists.txt +++ b/utils/soca/CMakeLists.txt @@ -18,7 +18,6 @@ ecbuild_add_executable( TARGET gdas_socahybridweights.x target_compile_features( gdas_socahybridweights.x PUBLIC cxx_std_17) target_link_libraries( gdas_socahybridweights.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) - # Diag B ecbuild_add_executable( TARGET gdas_diagb.x SOURCES gdas_diagb.cc ) From fd6a57195422ddc7d0e4580f64a9799902c5e05d Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 16 Apr 2024 15:11:41 -0400 Subject: [PATCH 19/21] code tidy --- parm/soca/berror/soca_diagb.yaml | 30 --- .../soca/berror/soca_diagb.yaml.j2 | 28 ++- .../berror/soca_parameters_diffusion_vt.yaml | 2 - ...ca_vtscales.yaml => soca_vtscales.yaml.j2} | 2 +- parm/soca/fields_metadata.yaml | 180 ++++++++++++++++++ parm/soca/obsprep/obsprep_config.yaml | 9 +- scripts/exgdas_global_marine_analysis_bmat.sh | 3 +- scripts/exgdas_global_marine_analysis_prep.py | 21 +- test/soca/gw/static.sh | 4 +- utils/soca/CMakeLists.txt | 12 +- .../{gdas_diagb.cc => gdas_soca_diagb.cc} | 4 +- .../soca/{gdas_diagb.h => gdas_soca_diagb.h} | 133 +++++++------ utils/test/testref/ghrsst2ioda.test | 6 +- 13 files changed, 294 insertions(+), 140 deletions(-) delete mode 100644 parm/soca/berror/soca_diagb.yaml rename test/soca/testinput/gdas_diagb.yaml => parm/soca/berror/soca_diagb.yaml.j2 (59%) rename parm/soca/berror/{soca_vtscales.yaml => soca_vtscales.yaml.j2} (77%) create mode 100644 parm/soca/fields_metadata.yaml rename utils/soca/{gdas_diagb.cc => gdas_soca_diagb.cc} (68%) rename utils/soca/{gdas_diagb.h => gdas_soca_diagb.h} (82%) diff --git a/parm/soca/berror/soca_diagb.yaml b/parm/soca/berror/soca_diagb.yaml deleted file mode 100644 index 27ce039c5..000000000 --- a/parm/soca/berror/soca_diagb.yaml +++ /dev/null @@ -1,30 +0,0 @@ -geometry: - mom6_input_nml: mom_input.nml - fields metadata: ./fields_metadata.yaml - -date: '{{ATM_WINDOW_MIDDLE}}' - -background: - date: '{{ATM_WINDOW_MIDDLE}}' - basename: ./bkg/ - ocn_filename: '{{RUN}}.ocean.t{{gcyc}}z.inst.f009.nc' - ice_filename: '{{RUN}}.agg_ice.t{{gcyc}}z.inst.f009.nc' - read_from_file: 1 - -background error: - datadir: ./ - date: '{{ATM_WINDOW_MIDDLE}}' - exp: bkgerr_stddev - type: incr - -variables: - name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon, mom6_mld] - -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/test/soca/testinput/gdas_diagb.yaml b/parm/soca/berror/soca_diagb.yaml.j2 similarity index 59% rename from test/soca/testinput/gdas_diagb.yaml rename to parm/soca/berror/soca_diagb.yaml.j2 index 24c9bb30c..c57a4d82f 100644 --- a/test/soca/testinput/gdas_diagb.yaml +++ b/parm/soca/berror/soca_diagb.yaml.j2 @@ -2,35 +2,34 @@ geometry: mom6_input_nml: mom_input.nml fields metadata: ./fields_metadata.yaml -date: 2018-04-15T09:00:00Z +date: '{{ ATM_WINDOW_MIDDLE }}' background: - date: 2018-04-15T09:00:00Z - basename: ./INPUT/ - ocn_filename: MOM.res.nc - ice_filename: cice.res.nc + date: '{{ ATM_WINDOW_MIDDLE }}' + basename: ./bkg/ + ocn_filename: '{{ RUN }}.ocean.t{{ gcyc }}z.inst.f009.nc' + ice_filename: '{{ RUN }}.agg_ice.t{{ gcyc }}z.inst.f009.nc' read_from_file: 1 background error: datadir: ./ - date: '2018-04-15T09:00:00Z' - exp: bkgerr + date: '{{ ATM_WINDOW_MIDDLE }}' + exp: bkgerr_stddev type: incr variables: - name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon] + name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon, mom6_mld] -rescale: 2.0 -min sst: 0.0 -max ssh: 0.0 -min depth: 300.0 +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: 1 + horizontal iterations: 2 vertical iterations: 1 +# TODO(G): Start using when the normalization is optional #diffusion: # saber block name: EXPLICIT_DIFFUSION # active variables: [tocn, socn, ssh, cicen, hicen, hsnon] @@ -58,4 +57,3 @@ simple smoothing: # - name: ice # horizontal: # filename: hz_ice.nc - diff --git a/parm/soca/berror/soca_parameters_diffusion_vt.yaml b/parm/soca/berror/soca_parameters_diffusion_vt.yaml index 0c63938c9..98370be7d 100644 --- a/parm/soca/berror/soca_parameters_diffusion_vt.yaml +++ b/parm/soca/berror/soca_parameters_diffusion_vt.yaml @@ -28,7 +28,5 @@ background error: from file: filename: vt_scales.nc variable name: vt - #as gaussian: true - #fixed value: 5.0 write: filename: vt_ocean.nc diff --git a/parm/soca/berror/soca_vtscales.yaml b/parm/soca/berror/soca_vtscales.yaml.j2 similarity index 77% rename from parm/soca/berror/soca_vtscales.yaml rename to parm/soca/berror/soca_vtscales.yaml.j2 index aa8975c89..bad7d623f 100644 --- a/parm/soca/berror/soca_vtscales.yaml +++ b/parm/soca/berror/soca_vtscales.yaml.j2 @@ -1,6 +1,6 @@ gridspec_filename: soca_gridspec.nc restart_filename: ./INPUT/MOM.res.nc -mld_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' +mld_filename: 'ocn.bkgerr_stddev.incr.{{ ATM_WINDOW_MIDDLE }}.nc' output_filename: ./vt_scales.nc output_variable_vt: vt output_variable_hz: hz diff --git a/parm/soca/fields_metadata.yaml b/parm/soca/fields_metadata.yaml new file mode 100644 index 000000000..d239c4000 --- /dev/null +++ b/parm/soca/fields_metadata.yaml @@ -0,0 +1,180 @@ +# -------------------------------------------------------------------------------------------------- +# Field metadata for SOCA. Each field can contain the following information: +# +# name: Internal name used by soca code and config files +# grid: "h", "u", or "v" (Default: h) +# masked: use land mask if true (Default: true) +# levels: "1" or "full_ocn" (Default: 1) +# getval_name: variable name expected by GetValues (Default: ) +# getval_name_surface: GetValues variable name for 2D surface of a 3D field (Default: ) +# io_file: The restart file domain "ocn", "sfc", or "ice" (Default: ) +# io_name: The variable name used in the restart IO (Default: ) +# -------------------------------------------------------------------------------------------------- + + +# -------------------------------------------------------------------------------------------------- +# Ocean state variables +# -------------------------------------------------------------------------------------------------- +- name: tocn + levels: full_ocn + getval name: sea_water_potential_temperature + getval name surface: sea_surface_temperature + io file: ocn + io name: Temp + fill value: 0.0 + +- name: socn + levels: full_ocn + getval name: sea_water_salinity + getval name surface: sea_surface_salinity + io file: ocn + io name: Salt + property: positive_definite + fill value: 0.0 + +- name: uocn + grid: u + levels: full_ocn + getval name: eastward_sea_water_velocity + getval name surface: surface_eastward_sea_water_velocity + io file: ocn + io name: u + fill value: 0.0 + +- name: vocn + grid: v + levels: full_ocn + getval name: northward_sea_water_velocity + getval name surface: surface_northward_sea_water_velocity + io file: ocn + io name: v + fill value: 0.0 + +- name: hocn + levels: full_ocn + getval name: sea_water_cell_thickness + io file: ocn + io name: h + fill value: 0.001 + vert interp: false + +- name: ssh + getval name: sea_surface_height_above_geoid + io file: ocn + io name: ave_ssh + fill value: 0.0 + +- name: mom6_mld + io file: ocn + io name: MLD + fill value: 0.0 + +# -------------------------------------------------------------------------------------------------- +# ice state variables +# -------------------------------------------------------------------------------------------------- +- name: hicen + getval name: sea_ice_category_thickness + io file: ice + io name: hicen + property: positive_definite + fill value: 0.0 + +- name: cicen + getval name: sea_ice_category_area_fraction + getval name surface: sea_ice_area_fraction # note: not accurate, should be "sum" not "surface" + io file: ice + io name: aicen + fill value: 0.0 + +- name: hsnon + getval name: sea_ice_category_snow_thickness + io file: ice + io name: hsnon + property: positive_definite + fill value: 0.0 + +# -------------------------------------------------------------------------------------------------- +# wave state variables +# -------------------------------------------------------------------------------------------------- +- name: swh + getval name: sea_surface_wave_significant_height + io file: wav + io name: hs + property: positive_definite + +# -------------------------------------------------------------------------------------------------- +# sea surface variables +# -------------------------------------------------------------------------------------------------- +- name: sw + masked: false + getval name: net_downwelling_shortwave_radiation + io file: sfc + io name: sw_rad + +- name: lw + masked: false + getval name: net_downwelling_longwave_radiation + io file: sfc + io name: lw_rad + +- name: lhf + masked: false + getval name: upward_latent_heat_flux_in_air + io file: sfc + io name: latent_heat + +- name: shf + masked: false + getval name: upward_sensible_heat_flux_in_air + io file: sfc + io name: sens_heat + +- name: us + masked: false + getval name: friction_velocity_over_water + io file: sfc + io name: fric_vel + + +# -------------------------------------------------------------------------------------------------- +# BGC +# -------------------------------------------------------------------------------------------------- +- name: chl + levels: full_ocn + getval name: mass_concentration_of_chlorophyll_in_sea_water + getval name surface: sea_surface_chlorophyll + io file: ocn + io name: chl + property: positive_definite + +- name: biop + levels: full_ocn + getval name: molar_concentration_of_biomass_in_sea_water_in_p_units + getval name surface: sea_surface_biomass_in_p_units + io file: ocn + io name: biomass_p + property: positive_definite + +# -------------------------------------------------------------------------------------------------- +# built-in derived variables that you don't need to worry about +# -------------------------------------------------------------------------------------------------- +- name: distance_from_coast + masked: false + +- name: layer_depth + levels: full_ocn + vert interp: false + +- name: mesoscale_representation_error + +- name: mld + +- name: sea_floor_depth_below_sea_surface + +- name: sea_area_fraction + masked: false + +- name: surface_temperature_where_sea + +- name: sea_water_depth + levels: full_ocn diff --git a/parm/soca/obsprep/obsprep_config.yaml b/parm/soca/obsprep/obsprep_config.yaml index 982873644..cb7e4e9dc 100644 --- a/parm/soca/obsprep/obsprep_config.yaml +++ b/parm/soca/obsprep/obsprep_config.yaml @@ -19,11 +19,12 @@ observations: name: adt_rads_all dmpdir subdir: ocean/adt provider: RADS - output file: adt_rads_all.nc4 - error ratio: 1.0 # [m/s] correspond to 1.0m/day + type: nc + dmpdir regex: 'rads_adt_*.nc' window: - back: 0 # look back 8 six-hourly obs dumps - forward: 0 # look forward 1 six-hourly bin + back: 0 + forward: 0 + error ratio: 0.4 # meters per day # Ice concentration - obs space: diff --git a/scripts/exgdas_global_marine_analysis_bmat.sh b/scripts/exgdas_global_marine_analysis_bmat.sh index ce8eec6da..4548dd13c 100755 --- a/scripts/exgdas_global_marine_analysis_bmat.sh +++ b/scripts/exgdas_global_marine_analysis_bmat.sh @@ -88,7 +88,7 @@ fi ################################################################################ # Compute the parametric diag of B -$APRUN_OCNANAL $JEDI_BIN/gdas_diagb.x soca_diagb.yaml +$APRUN_OCNANAL $JEDI_BIN/gdas_soca_diagb.x soca_diagb.yaml export err=$?; err_chk if [ $err -gt 0 ]; then exit $err @@ -114,6 +114,7 @@ fi ################################################################################ # Initialize the vertical diffusion blocks # Compute the MLD dependent vertical scales +# TODO(G): Probably not the right way to execute this script python ${HOMEgfs}/sorc/gdas.cd/sorc/soca/tools/calc_scales.py soca_vtscales.yaml clean_yaml soca_parameters_diffusion_vt.yaml diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index 5bce33e0d..85b611cb9 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -10,7 +10,8 @@ from soca import bkg_utils from datetime import datetime, timedelta import pytz -from wxflow import (Logger, Template, TemplateConstants, YAMLFile, FileHandler) +from wxflow import (Logger, Template, TemplateConstants, + YAMLFile, FileHandler, AttrDict, parse_j2yaml) logger = Logger() @@ -196,12 +197,9 @@ def find_clim_ens(input_date): berror_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror') logger.info(f"---------------- generate soca_diagb.yaml") -diagb_yaml = os.path.join(anl_dir, 'soca_diagb.yaml') -diagb_yaml_template = os.path.join(berror_yaml_dir, 'soca_diagb.yaml') -config = YAMLFile(path=diagb_yaml_template) -config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) -config = Template.substitute_structure(config, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get) -config.save(diagb_yaml) +conf = parse_j2yaml(path=os.path.join(berror_yaml_dir, 'soca_diagb.yaml.j2'), + data=envconfig) +conf.save(os.path.join(anl_dir, 'soca_diagb.yaml')) logger.info(f"---------------- generate soca_ensb.yaml") berr_yaml = os.path.join(anl_dir, 'soca_ensb.yaml') @@ -245,12 +243,9 @@ def find_clim_ens(input_date): config.save(diffu_hz_yaml) logger.info(f"---------------- generate soca_vtscales.yaml") -vtscales_yaml = os.path.join(anl_dir, 'soca_vtscales.yaml') -vtscales_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror') -vtscales_yaml_template = os.path.join(berror_yaml_dir, 'soca_vtscales.yaml') -config = YAMLFile(path=vtscales_yaml_template) -config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) -config.save(vtscales_yaml) +conf = parse_j2yaml(path=os.path.join(gdas_home, 'parm', 'soca', 'berror', 'soca_vtscales.yaml.j2'), + data=envconfig) +conf.save(os.path.join(anl_dir, 'soca_vtscales.yaml')) logger.info(f"---------------- generate soca_parameters_diffusion_vt.yaml") diffu_vt_yaml = os.path.join(anl_dir, 'soca_parameters_diffusion_vt.yaml') diff --git a/test/soca/gw/static.sh b/test/soca/gw/static.sh index 05407ea92..30e6b25c9 100755 --- a/test/soca/gw/static.sh +++ b/test/soca/gw/static.sh @@ -14,8 +14,8 @@ lowres=${project_source_dir}/sorc/soca/test/Data cp -L ${lowres}/workdir/{diag_table,field_table} ${soca_static} cp -L ${project_source_dir}/test/soca/fix/MOM_input ${soca_static} -cp -L ${lowres}/{fields_metadata.yml,godas_sst_bgerr.nc,rossrad.dat} ${soca_static} -mv ${soca_static}/fields_metadata.yml ${soca_static}/fields_metadata.yaml +cp -L ${lowres}/rossrad.dat ${soca_static} +cp -L ${project_source_dir}/parm/soca/fields_metadata.yaml ${soca_static} cp -L ${project_source_dir}/sorc/soca/test/testinput/obsop_name_map.yml ${soca_static}/obsop_name_map.yaml cp -L ${lowres}/72x35x25/input.nml ${soca_static}/inputnml cp -L ${lowres}/72x35x25/INPUT/{hycom1_25.nc,ocean_mosaic.nc,grid_spec.nc,layer_coord25.nc,ocean_hgrid.nc,ocean_topog.nc} ${soca_static}/INPUT diff --git a/utils/soca/CMakeLists.txt b/utils/soca/CMakeLists.txt index a1fa005dd..62996f2af 100644 --- a/utils/soca/CMakeLists.txt +++ b/utils/soca/CMakeLists.txt @@ -18,14 +18,14 @@ ecbuild_add_executable( TARGET gdas_socahybridweights.x target_compile_features( gdas_socahybridweights.x PUBLIC cxx_std_17) target_link_libraries( gdas_socahybridweights.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) -# Diag B -ecbuild_add_executable( TARGET gdas_diagb.x - SOURCES gdas_diagb.cc ) -target_compile_features( gdas_diagb.x PUBLIC cxx_std_17) -target_link_libraries( gdas_diagb.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) - # Obs Stats ecbuild_add_executable( TARGET gdassoca_obsstats.x SOURCES gdassoca_obsstats.cc ) target_compile_features( gdassoca_obsstats.x PUBLIC cxx_std_17) target_link_libraries( gdassoca_obsstats.x PUBLIC NetCDF::NetCDF_CXX oops ioda) + +# Diag B +ecbuild_add_executable( TARGET gdas_soca_diagb.x + SOURCES gdas_soca_diagb.cc ) +target_compile_features( gdas_soca_diagb.x PUBLIC cxx_std_17) +target_link_libraries( gdas_soca_diagb.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) diff --git a/utils/soca/gdas_diagb.cc b/utils/soca/gdas_soca_diagb.cc similarity index 68% rename from utils/soca/gdas_diagb.cc rename to utils/soca/gdas_soca_diagb.cc index 09332cd2e..7273f2e75 100644 --- a/utils/soca/gdas_diagb.cc +++ b/utils/soca/gdas_soca_diagb.cc @@ -1,8 +1,8 @@ -#include "gdas_diagb.h" +#include "gdas_soca_diagb.h" #include "oops/runs/Run.h" int main(int argc, char ** argv) { oops::Run run(argc, argv); - gdasapp::DiagB diagb; + gdasapp::SocaDiagB diagb; return run.execute(diagb); } diff --git a/utils/soca/gdas_diagb.h b/utils/soca/gdas_soca_diagb.h similarity index 82% rename from utils/soca/gdas_diagb.h rename to utils/soca/gdas_soca_diagb.h index cce9989fb..62cff84f7 100644 --- a/utils/soca/gdas_diagb.h +++ b/utils/soca/gdas_soca_diagb.h @@ -1,6 +1,6 @@ #pragma once -#include +#include #include #include #include @@ -9,6 +9,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" @@ -16,7 +17,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" @@ -28,15 +28,15 @@ #include "oops/util/FieldSetOperations.h" #include "oops/util/Logger.h" +#include "soca/ExplicitDiffusion/ExplicitDiffusion.h" +#include "soca/ExplicitDiffusion/ExplicitDiffusionParameters.h" #include "soca/Geometry/Geometry.h" #include "soca/Increment/Increment.h" #include "soca/State/State.h" -#include "soca/ExplicitDiffusion/ExplicitDiffusion.h" -#include "soca/ExplicitDiffusion/ExplicitDiffusionParameters.h" namespace gdasapp { /** - * DiagB Class Implementation + * SocaDiagB Class Implementation * * Implements variance partitioning within the GDAS for the ocean and sea-ice. It is used as a proxy for the * diagonal of B. @@ -46,33 +46,36 @@ namespace gdasapp { */ // ----------------------------------------------------------------------------- - // Since we can't have useful data members in DiagB because execute is "const" - // and the constructor does not know about the configuration, - // we're going around the issue with the use of a simple data structure - struct DiagBConfig { + // Since we can't have useful data members in SocaDiagB because execute is "const" + // and the constructor does not know about the configuration, + // we're going around the issue with the use of a simple data structure + struct SocaDiagBConfig { util::DateTime cycleDate; oops::Variables socaVars; double sshMax; // poor man's alternative to limiting the unbalanced ssh variance double depthMin; // set the bkg error to 0 for depth < depthMin bool diffusion; // apply explicit diffusion to the std. dev. fields bool simpleSmoothing; // Simple spatial averaging - int niter; // number of iterations for the horizontal averaging - int niterVert; // number of iterations for the vertical averaging + int niterHoriz; // number of iterations for the horizontal averaging + int niterVert; // number of iterations for the vertical averaging double rescale; // inflation/deflation of the variance after filtering }; // ----------------------------------------------------------------------------- - class DiagB : public oops::Application { - public: + class SocaDiagB : public oops::Application { + public: // ----------------------------------------------------------------------------- - explicit DiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) + explicit SocaDiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) : Application(comm) { } // ----------------------------------------------------------------------------- - DiagBConfig setup(const eckit::Configuration & fullConfig) const { - - DiagBConfig diagBConfig; + /** + * Sets up and returns a SocaDiagBConfig object configured according to the provided + * eckit::Configuration + */ + SocaDiagBConfig setup(const eckit::Configuration & fullConfig) const { + SocaDiagBConfig diagBConfig; /// Get the date // ------------- @@ -98,15 +101,15 @@ namespace gdasapp { // Explicit diffusion diagBConfig.diffusion = false; if (fullConfig.has("diffusion")) { - diagBConfig.diffusion = true; + diagBConfig.diffusion = true; } // Simple smoothing parameters diagBConfig.simpleSmoothing = false; if (fullConfig.has("simple smoothing")) { - diagBConfig.simpleSmoothing = true; - fullConfig.get("simple smoothing.horizontal iterations", diagBConfig.niter); - fullConfig.get("simple smoothing.vertical iterations", diagBConfig.niterVert); + diagBConfig.simpleSmoothing = true; + fullConfig.get("simple smoothing.horizontal iterations", diagBConfig.niterHoriz); + fullConfig.get("simple smoothing.vertical iterations", diagBConfig.niterVert); } // Variance rescaling @@ -116,9 +119,12 @@ namespace gdasapp { } // ----------------------------------------------------------------------------- - static const std::string classname() {return "gdasapp::DiagB";} + static const std::string classname() {return "gdasapp::SocaDiagB";} // ----------------------------------------------------------------------------- + /** + * Variance partitioning at jnode and level + */ void stdDevFilt(const int jnode, const int level, const int nbz, @@ -129,10 +135,9 @@ namespace gdasapp { const atlas::array::ArrayView bkg, const atlas::array::ArrayView bathy, atlas::array::ArrayView& stdDevBkg, - bool doBathy = true, - int minn = 6) const { - - // Early exit if too shallow + bool doBathy = true, + int minn = 6) const { + // Early exit if too shallow if ( doBathy & bathy(jnode, 0) < depthMin ) { stdDevBkg(jnode, level) = 0.0; return; @@ -143,17 +148,17 @@ namespace gdasapp { for (int nn = 0; nn < neighbors.size(); ++nn) { int levelMin = std::max(0, level - nbz); int levelMax = level + nbz; - // Do weird things withing the MLD + // Do weird things withing the MLD if (level < nzMld) { // If in the MLD, compute the std. dev. throughout the MLD levelMin = 0; - levelMax = 1; //nzMld; Projecting the surface variance down the water column for now + levelMax = 1; // nzMld; Projecting the surface variance down the water column for now } - // 2D case + // 2D case if (bkg.shape(1) == 1) { levelMin = 0; - levelMax = 0; //nzMld; Projecting the surface variance down the water column for now - } + levelMax = 0; // nzMld; Projecting the surface variance down the water column for now + } for (int ll = levelMin; ll <= levelMax; ++ll) { if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { continue; @@ -177,10 +182,12 @@ namespace gdasapp { } // ----------------------------------------------------------------------------- + /** + * Implementation of the virtual execute method from the Application parent class + */ int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { - - /// Initialize the paramters for D - DiagBConfig configD = setup(fullConfig); + /// Initialize the paramters for D + SocaDiagBConfig configD = setup(fullConfig); /// Setup the soca geometry oops::Log::info() << "====================== geometry" << std::endl; @@ -214,7 +221,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); @@ -261,7 +269,7 @@ namespace gdasapp { viewDepth(jnode, 0) = 0.5 * viewHocn(jnode, 0); for (atlas::idx_t level = 1; level < depth.shape(1); ++level) { viewDepth(jnode, level) = viewDepth(jnode, level-1) + - 0.5 * (viewHocn(jnode, level-1 ) + viewHocn(jnode, level)); + 0.5 * (viewHocn(jnode, level-1) + viewHocn(jnode, level)); } } @@ -270,7 +278,7 @@ namespace gdasapp { atlas::array::ArrayT bathy(viewHocn.shape(0), 1); auto viewBathy = atlas::array::make_view(bathy); for (atlas::idx_t jnode = 0; jnode < viewHocn.shape(0); ++jnode) { - viewBathy(jnode, 0) = 0.0; + viewBathy(jnode, 0) = 0.0; for (atlas::idx_t level = 0; level < viewHocn.shape(1); ++level) { viewBathy(jnode, 0) += viewHocn(jnode, level); } @@ -296,8 +304,8 @@ namespace gdasapp { // Loop through variables for (auto & var : configD.socaVars.variables()) { - // Update the halo - nodeColumns.haloExchange(xbFs[var]); + // Update the halo + nodeColumns.haloExchange(xbFs[var]); // Skip the layer thickness variable if (var == "hocn") { @@ -317,23 +325,25 @@ namespace gdasapp { // get indices of the cell's neighbors auto neighbors = get_neighbors_of_node(jnode); - // 2D case + // 2D case if (xbFs[var].shape(1) == 1) { - // Std. dev. of the partition - stdDevFilt(jnode, 0, 0, - configD.depthMin, neighbors, 0, viewHocn, bkg, viewBathy, stdDevBkg, false, 4); - if (var == "ssh") { - // TODO(G): Extract the unbalanced ssh variance, in the mean time, do this: - stdDevBkg(jnode, 0) = std::min(configD.sshMax, stdDevBkg(jnode, 0)); - } + // Std. dev. of the partition + stdDevFilt(jnode, 0, 0, + configD.depthMin, neighbors, 0, + viewHocn, bkg, viewBathy, stdDevBkg, false, 4); + if (var == "ssh") { + // TODO(G): Extract the unbalanced ssh variance, in the mean time, do this: + stdDevBkg(jnode, 0) = std::min(configD.sshMax, stdDevBkg(jnode, 0)); + } } 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. dev. of the partition - stdDevFilt(jnode, level, nbz, - configD.depthMin, neighbors, nzMld, viewHocn, bkg, viewBathy, stdDevBkg, true, 6); + stdDevFilt(jnode, level, nbz, + configD.depthMin, neighbors, nzMld, + viewHocn, bkg, viewBathy, stdDevBkg, true, 6); } } // end 3D case } // end jnode @@ -349,8 +359,7 @@ namespace gdasapp { } // Horizontal averaging - for (int iter = 0; iter < configD.niter; ++iter) { - + for (int iter = 0; iter < configD.niterHoriz; ++iter) { // Update the halo points nodeColumns.haloExchange(bkgErrFs[var]); auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); @@ -376,7 +385,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(); } // Reset to 0 over land @@ -396,9 +406,9 @@ namespace gdasapp { for (int iter = 0; iter < configD.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, level) = (tmpArray(jnode, level-1) + + tmpArray(jnode, level) + + tmpArray(jnode, level+1)) / 3.0; } stdDevBkg(jnode, 0) = stdDevBkg(jnode, 1); } @@ -409,16 +419,18 @@ namespace gdasapp { /// Use explicit diffusion to smooth the background error // ------------------------------------------------------ - // TODO: Use this once Travis adds the option to skip the normalization. + // TODO(G): 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"); soca::ExplicitDiffusionParameters params; params.deserialize(diffusionConfig); - oops::GeometryData geometryData(geom.functionSpace(), bkgErrFs["tocn"], true, this->getComm()); + oops::GeometryData geometryData(geom.functionSpace(), + bkgErrFs["tocn"], true, this->getComm()); const oops::FieldSet3D dumyXb(configD.cycleDate, this->getComm()); - soca::ExplicitDiffusion diffuse(geometryData, configD.socaVars, diffusionConfig, params, dumyXb, dumyXb); + soca::ExplicitDiffusion diffuse(geometryData, configD.socaVars, + diffusionConfig, params, dumyXb, dumyXb); diffuse.read(); // Smooth the field @@ -443,12 +455,11 @@ namespace gdasapp { // ----------------------------------------------------------------------------- - private: + private: std::string appname() const { - return "gdasapp::DiagB"; + return "gdasapp::SocaDiagB"; } // ----------------------------------------------------------------------------- - }; } // namespace gdasapp diff --git a/utils/test/testref/ghrsst2ioda.test b/utils/test/testref/ghrsst2ioda.test index 1413f9d04..ab0a0b561 100644 --- a/utils/test/testref/ghrsst2ioda.test +++ b/utils/test/testref/ghrsst2ioda.test @@ -21,6 +21,6 @@ latitude: Max: 77.95 Sum: 623.423 datetime: - Min: 1269445504 - Max: 1269445504 - Sum: 10155564032 + Min: 1269445446 + Max: 1269445446 + Sum: 10155563568 From 5ee0fe097bef725709fff2e42539fbd972b4e5be Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 16 Apr 2024 20:04:56 -0400 Subject: [PATCH 20/21] Update exgdas_global_marine_analysis_bmat.sh --- scripts/exgdas_global_marine_analysis_bmat.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/exgdas_global_marine_analysis_bmat.sh b/scripts/exgdas_global_marine_analysis_bmat.sh index 4548dd13c..db6efa87f 100755 --- a/scripts/exgdas_global_marine_analysis_bmat.sh +++ b/scripts/exgdas_global_marine_analysis_bmat.sh @@ -118,7 +118,7 @@ fi python ${HOMEgfs}/sorc/gdas.cd/sorc/soca/tools/calc_scales.py soca_vtscales.yaml clean_yaml soca_parameters_diffusion_vt.yaml -# Vertical fiffusion and normalization +# Vertical diffusion and normalization $APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_vt.yaml export err=$?; err_chk if [ $err -gt 0 ]; then From f455a697fb973e5d30c033e20215770d972ffacd Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Wed, 17 Apr 2024 14:23:28 -0400 Subject: [PATCH 21/21] added back hybrid --- parm/soca/berror/soca_hybrid_bmat.yaml | 92 +++++++++++++++++++ ...aber_blocks.yaml => soca_static_bmat.yaml} | 0 scripts/exgdas_global_marine_analysis_bmat.sh | 90 +++++++++--------- scripts/exgdas_global_marine_analysis_post.py | 13 ++- scripts/exgdas_global_marine_analysis_prep.py | 46 +++++----- scripts/exgdas_global_marine_analysis_run.sh | 2 +- 6 files changed, 172 insertions(+), 71 deletions(-) create mode 100644 parm/soca/berror/soca_hybrid_bmat.yaml rename parm/soca/berror/{saber_blocks.yaml => soca_static_bmat.yaml} (100%) diff --git a/parm/soca/berror/soca_hybrid_bmat.yaml b/parm/soca/berror/soca_hybrid_bmat.yaml new file mode 100644 index 000000000..6e3d35351 --- /dev/null +++ b/parm/soca/berror/soca_hybrid_bmat.yaml @@ -0,0 +1,92 @@ +covariance model: hybrid +components: +- covariance: + covariance model: SABER + saber central block: + saber block name: EXPLICIT_DIFFUSION + active variables: [tocn, socn, ssh, cicen] + geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + group mapping: + - name: ocean + variables: + - tocn + - socn + - ssh + - name: ice + variables: + - cicen + read: + groups: + - name: ocean + horizontal: + filename: hz_ocean.nc + vertical: + filename: vt_ocean.nc + - name: ice + horizontal: + filename: hz_ice.nc + + saber outer blocks: + - saber block name: StdDev + read: + model file: + date: '{{ATM_WINDOW_MIDDLE}}' + basename: ./ + ocn_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' + ice_filename: 'ice.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' + read_from_file: 3 + + linear variable change: + input variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] + output variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] + linear variable changes: + - linear variable change name: BalanceSOCA + + weight: + value: 1.00 + +- covariance: + covariance model: ensemble + members from template: + template: + read_from_file: 1 + date: '{{ATM_WINDOW_MIDDLE}}' + basename: ./static_ens/ + ocn_filename: 'ocn.pert.steric.%mem%.nc' + ice_filename: 'ice.%mem%.nc' + state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon] + pattern: '%mem%' + nmembers: ${ENS_SIZE} + localization: + localization method: SABER + saber central block: + saber block name: EXPLICIT_DIFFUSION + active variables: [tocn, socn, ssh] + geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + group mapping: + - name: ocean + variables: [tocn, socn, ssh, uocn, vocn] + - name: ice + variables: [cicen, hicen, hsnon] + read: + groups: + - name: ocean + multivariate strategy: duplicated + horizontal: + filename: hz_ocean.nc + vertical: + strategy: duplicated + - name: ice + horizontal: + filename: hz_ice.nc + + weight: + read_from_file: 3 + basename: ./ + ocn_filename: 'ocn.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc' + ice_filename: 'ice.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc' + date: '{{ATM_WINDOW_MIDDLE}}' diff --git a/parm/soca/berror/saber_blocks.yaml b/parm/soca/berror/soca_static_bmat.yaml similarity index 100% rename from parm/soca/berror/saber_blocks.yaml rename to parm/soca/berror/soca_static_bmat.yaml diff --git a/scripts/exgdas_global_marine_analysis_bmat.sh b/scripts/exgdas_global_marine_analysis_bmat.sh index db6efa87f..9ae2b3646 100755 --- a/scripts/exgdas_global_marine_analysis_bmat.sh +++ b/scripts/exgdas_global_marine_analysis_bmat.sh @@ -56,35 +56,38 @@ else fi fi -################################################################################ -# Write ensemble weights for the hybrid envar -#$APRUN_OCNANAL $JEDI_BIN/gdas_socahybridweights.x soca_ensweights.yaml -#export err=$?; err_chk -#if [ $err -gt 0 ]; then -# exit $err -#fi ################################################################################ -# Ensemble perturbations for the EnVAR and diagonal of static B -# Static B: -# - Compute moments of original ensemble with the balanced ssh removed -# from the statistics -# -# Ensemble of perturbations: -# - apply the linear steric height balance to each members, using the -# deterministic for trajectory. -# - add the unblanced ssh to the steric ssh field -# Diagnostics: -# - variance explained by the steric heigh -# - moments +# Ensemble processing +if [ "${NMEM_ENS}" -gt 3 ]; then + ################################################################################ + # Write ensemble weights for the hybrid envar + $APRUN_OCNANAL $JEDI_BIN/gdas_socahybridweights.x soca_ensweights.yaml + export err=$?; err_chk + if [ $err -gt 0 ]; then + exit $err + fi -# Process static ensemble -#clean_yaml soca_clim_ens_moments.yaml -#$APRUN_OCNANAL $JEDI_BIN/gdas_ens_handler.x soca_ensb.yaml -#export err=$?; err_chk -#if [ $err -gt 0 ]; then -# exit $err -#fi + ################################################################################ + # Ensemble perturbations for the EnVAR + # Ensemble diagnostics: + # - Compute moments of original ensemble with the balanced ssh removed + # from the statistics + # + # Ensemble of perturbations: + # - apply the linear steric height balance to each members, using the + # deterministic for trajectory. + # - add the unblanced ssh to the steric ssh field + # Diagnostics: + # - variance explained by the steric heigh + # - moments + clean_yaml soca_clim_ens_moments.yaml + $APRUN_OCNANAL $JEDI_BIN/gdas_ens_handler.x soca_ensb.yaml + export err=$?; err_chk + if [ $err -gt 0 ]; then + exit $err + fi +fi ################################################################################ # Compute the parametric diag of B @@ -95,30 +98,33 @@ if [ $err -gt 0 ]; then fi ################################################################################ -# Set decorrelation scales for the static B -$APRUN_OCNANAL $JEDI_BIN/soca_setcorscales.x soca_setcorscales.yaml -export err=$?; err_chk -if [ $err -gt 0 ]; then - exit $err -fi - -################################################################################ -# Initialize the horizontal diffusion blocks -clean_yaml soca_parameters_diffusion_hz.yaml -$APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_hz.yaml -export err=$?; err_chk -if [ $err -gt 0 ]; then - exit $err +# Horizontal diffusion +if [ ! -f "ocn.cor_rh.incr.0001-01-01T00:00:00Z.nc" ]; then + # Set decorrelation scales for the static B + $APRUN_OCNANAL $JEDI_BIN/soca_setcorscales.x soca_setcorscales.yaml + export err=$?; err_chk + if [ $err -gt 0 ]; then + exit $err + fi + # Initialize the horizontal diffusion block and normalize + clean_yaml soca_parameters_diffusion_hz.yaml + $APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_hz.yaml + export err=$?; err_chk + if [ $err -gt 0 ]; then + exit $err + fi +else + echo "Diffusion weights already calculated, skipping the init of the horizontal diffusion" fi ################################################################################ -# Initialize the vertical diffusion blocks +# Initialize the vertical diffusion block # Compute the MLD dependent vertical scales # TODO(G): Probably not the right way to execute this script python ${HOMEgfs}/sorc/gdas.cd/sorc/soca/tools/calc_scales.py soca_vtscales.yaml clean_yaml soca_parameters_diffusion_vt.yaml -# Vertical diffusion and normalization +# Initialize the vertical diffusion block and normalize $APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_vt.yaml export err=$?; err_chk if [ $err -gt 0 ]; then diff --git a/scripts/exgdas_global_marine_analysis_post.py b/scripts/exgdas_global_marine_analysis_post.py index 14fa8a8f3..54cce2363 100755 --- a/scripts/exgdas_global_marine_analysis_post.py +++ b/scripts/exgdas_global_marine_analysis_post.py @@ -53,6 +53,7 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): bdatedt = datetime.strptime(cdate, '%Y%m%d%H') - timedelta(hours=3) bdate = datetime.strftime(bdatedt, '%Y-%m-%dT%H:00:00Z') mdate = datetime.strftime(datetime.strptime(cdate, '%Y%m%d%H'), '%Y-%m-%dT%H:00:00Z') +nmem_ens = int(os.getenv('NMEM_ENS')) logger.info(f"---------------- Copy from RUNDIR to COMOUT") @@ -69,8 +70,9 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.bkgerr_stddev.nc')]) # Copy the recentering error - # post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.ssh_recentering_error.incr.{bdate}.nc'), - # os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.recentering_error.nc')]) + if nmem_ens > 2: + post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.ssh_recentering_error.incr.{bdate}.nc'), + os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.recentering_error.nc')]) # Copy the ice and ocean increments post_file_list.append([os.path.join(anl_dir, 'Data', f'{domain}.3dvarfgat_pseudo.incr.{mdate}.nc'), @@ -81,9 +83,10 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}ana.nc')]) # Copy of the ssh diagnostics -# for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']: -# post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'), -# os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')]) +if nmem_ens > 2: + for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']: + post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'), + os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')]) # Copy DA grid (computed for the start of the window) post_file_list.append([os.path.join(anl_dir, 'soca_gridspec.nc'), diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index 85b611cb9..d80bfc834 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -53,6 +53,7 @@ def find_clim_ens(input_date): """ Find the clim. ens. that is the closest to the DA window """ + logger.info(f"$$$$$$$$$$$$$$$$$$$$$$$ {os.getenv('SOCA_INPUT_FIX_DIR')}") ens_clim_dir = os.path.join(os.getenv('SOCA_INPUT_FIX_DIR'), 'bkgerr', 'ens') dirs = glob.glob(os.path.join(ens_clim_dir, '*')) @@ -67,9 +68,10 @@ def find_clim_ens(input_date): comin_obs = os.getenv('COMIN_OBS') anl_dir = os.getenv('DATA') staticsoca_dir = os.getenv('SOCA_INPUT_FIX_DIR') +nmem_ens = 0 +nmem_ens = int(os.getenv('NMEM_ENS')) if os.getenv('DOHYBVAR') == "YES": dohybvar = True - nmem_ens = int(os.getenv('NMEM_ENS')) else: dohybvar = False @@ -138,10 +140,8 @@ def find_clim_ens(input_date): ################################################################################ # stage ensemble members -nmem_ens = 0 if dohybvar: logger.info("---------------- Stage ensemble members") - nmem_ens = int(os.getenv('NMEM_ENS')) ens_member_list = [] for mem in range(1, nmem_ens+1): for domain in ['ocean', 'ice']: @@ -162,20 +162,20 @@ def find_clim_ens(input_date): for mem in range(1, nmem_ens+1): cice_fname = os.path.realpath(os.path.join(static_ens, "ice."+str(mem)+".nc")) bkg_utils.cice_hist2fms(cice_fname, cice_fname) - -# Commented out while testing the parametric diagb -# else: -# logger.info("---------------- Stage offline ensemble members") -# ens_member_list = [] -# clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None)) -# nmem_ens = len(glob.glob(os.path.abspath(os.path.join(clim_ens_dir, 'ocn.*.nc')))) -# for domain in ['ocn', 'ice']: -# for mem in range(1, nmem_ens+1): -# fname = domain+"."+str(mem)+".nc" -# fname_in = os.path.abspath(os.path.join(clim_ens_dir, fname)) -# fname_out = os.path.abspath(os.path.join(static_ens, fname)) -# ens_member_list.append([fname_in, fname_out]) -# FileHandler({'copy': ens_member_list}).sync() +else: + if nmem_ens >= 3: + logger.info("---------------- Stage offline ensemble members") + ens_member_list = [] + clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None)) + list_of_members = glob.glob(os.path.join(clim_ens_dir, 'ocean.*.nc')) + nmem_ens = min(nmem_ens, len(list_of_members)) + for domain in ['ocean', 'ice']: + for mem in range(1, nmem_ens+1): + fname = domain+"."+str(mem)+".nc" + fname_in = os.path.join(clim_ens_dir, fname) + fname_out = os.path.join(static_ens, fname) + ens_member_list.append([fname_in, fname_out]) + FileHandler({'copy': ens_member_list}).sync() os.environ['ENS_SIZE'] = str(nmem_ens) @@ -271,13 +271,13 @@ def find_clim_ens(input_date): yaml_name='bkg_list.yaml') os.environ['BKG_LIST'] = 'bkg_list.yaml' -# select the SABER BLOCKS to use -if 'SABER_BLOCKS_YAML' in os.environ and os.environ['SABER_BLOCKS_YAML']: - saber_blocks_yaml = os.getenv('SABER_BLOCKS_YAML') - logger.info(f"using non-default SABER blocks yaml: {saber_blocks_yaml}") +# select the B-matrix to use +if nmem_ens > 3: + # Use a hybrid static-ensemble B + os.environ['SABER_BLOCKS_YAML'] = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'soca_hybrid_bmat.yaml') else: - logger.info(f"using default SABER blocks yaml") - os.environ['SABER_BLOCKS_YAML'] = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'saber_blocks.yaml') + # Use a static B + os.environ['SABER_BLOCKS_YAML'] = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'soca_static_bmat.yaml') # substitute templated variables in the var config logger.info(f"{config}") diff --git a/scripts/exgdas_global_marine_analysis_run.sh b/scripts/exgdas_global_marine_analysis_run.sh index 3e938df37..7ee52cd5f 100755 --- a/scripts/exgdas_global_marine_analysis_run.sh +++ b/scripts/exgdas_global_marine_analysis_run.sh @@ -38,7 +38,7 @@ function clean_yaml() } ################################################################################ -# run 3DVAR FGAT +# run the variational application cp var.yaml var_original.yaml clean_yaml var.yaml $APRUN_OCNANAL $JEDI_BIN/soca_var.x var.yaml