From 1dd216db8d26fa08f82163876ce1efd921485568 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Fri, 9 Jun 2023 15:53:32 -0400 Subject: [PATCH] oops app for soca incr --- CMakeLists.txt | 5 +- parm/soca/variational/socaincr2mom6.yaml | 27 +++++ .../exgdas_global_marine_analysis_chkpt.sh | 5 +- scripts/exgdas_global_marine_analysis_prep.py | 9 ++ utils/CMakeLists.txt | 10 ++ utils/socaincr2mom6.cc | 8 ++ utils/socaincr2mom6.h | 113 ++++++++++++++++++ 7 files changed, 172 insertions(+), 5 deletions(-) create mode 100644 parm/soca/variational/socaincr2mom6.yaml create mode 100644 utils/CMakeLists.txt create mode 100644 utils/socaincr2mom6.cc create mode 100644 utils/socaincr2mom6.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 03a5be5ed..3358a883d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -94,7 +94,7 @@ if(BUILD_GDASBUNDLE) endif() ecbuild_bundle( PROJECT gsw GIT "https://github.com/jcsda-internal/GSW-Fortran.git" BRANCH develop ) ecbuild_bundle( PROJECT mom6 GIT "https://github.com/jcsda-internal/MOM6.git" BRANCH main-ecbuild RECURSIVE ) - ecbuild_bundle( PROJECT soca GIT "https://github.com/jcsda-internal/soca.git" BRANCH feature/hybrid-B ) + ecbuild_bundle( PROJECT soca GIT "https://github.com/jcsda-internal/soca.git" BRANCH develop ) #feature/hybrid-B ) # Build IODA converters option(BUILD_IODA_CONVERTERS "Build IODA Converters" ON) @@ -106,6 +106,9 @@ if(BUILD_GDASBUNDLE) ecbuild_bundle( PROJECT land-imsproc GIT "https://github.com/NOAA-PSL/land-IMS_proc.git" TAG 6373819 ) ecbuild_bundle( PROJECT land-jediincr GIT "https://github.com/NOAA-PSL/land-apply_jedi_incr.git" TAG ced6576 ) + # Build JEDI/DA or other peripherals + add_subdirectory(utils) + # ioda, ufo, fv3-jedi, and saber test data #--------------------------------- if(CLONE_JCSDADATA) diff --git a/parm/soca/variational/socaincr2mom6.yaml b/parm/soca/variational/socaincr2mom6.yaml new file mode 100644 index 000000000..c6c2a50c3 --- /dev/null +++ b/parm/soca/variational/socaincr2mom6.yaml @@ -0,0 +1,27 @@ +geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + +date: '{{ATM_WINDOW_BEGIN}}' + +layers variable: [hocn] + +increment variables: [tocn, socn, uocn, vocn, ssh] + +vertical geometry: + date: '{{ATM_WINDOW_BEGIN}}' + basename: ./INPUT/ + ocn_filename: MOM.res.nc + read_from_file: 1 + +soca increment: + date: '{{ATM_WINDOW_BEGIN}}' + basename: ./Data/ + ocn_filename: 'ocn.3dvarfgat_pseudo.incr.{{ATM_WINDOW_BEGIN}}.nc' + read_from_file: 1 + +mom6 iau increment: + datadir: ./ + date: '{{ATM_WINDOW_BEGIN}}' + exp: mom6_iau + type: incr diff --git a/scripts/exgdas_global_marine_analysis_chkpt.sh b/scripts/exgdas_global_marine_analysis_chkpt.sh index 6863a29a6..ddc5c5eed 100755 --- a/scripts/exgdas_global_marine_analysis_chkpt.sh +++ b/scripts/exgdas_global_marine_analysis_chkpt.sh @@ -55,10 +55,7 @@ EOF --out "${mom6_iau_incr}" \ --nsst_yaml "nsst.yaml" else - ${HOMEgfs}/sorc/gdas.cd/ush/socaincr2mom6.py --incr "${soca_incr}" \ - --bkg "${DATA}/INPUT/MOM.res.nc" \ - --grid "${DATA}/soca_gridspec.nc" \ - --out "${mom6_iau_incr}" + $APRUN_OCNANAL ${JEDI_BIN}/socaincr2mom6.x socaincr2mom6.yaml fi export err=$? if [ $err -gt 0 ]; then diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index 6525f5ada..a87b906f8 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -502,6 +502,7 @@ def find_clim_ens(input_date): ufsda.yamltools.save_check(varconfig.as_dict(), target=var_yaml, app='var') ################################################################################ +# Prepare the yamls for the "checkpoint" jjob # prepare yaml and CICE restart for soca to cice change of variable logging.info(f"---------------- generate soca to cice yamls") @@ -531,6 +532,14 @@ def find_clim_ens(input_date): yaml.dump(soca2cice_cfg, f, sort_keys=False, default_flow_style=False) ufsda.genYAML.genYAML('tmp.yaml', output=varchgyaml) +# prepare yaml for soca to MOM6 IAU increment +logging.info(f"---------------- generate soca to MOM6 IAU yaml") +socaincr2mom6_yaml = os.path.join(anl_dir, 'socaincr2mom6.yaml') +socaincr2mom6_yaml_template = os.path.join(gdas_home, 'parm', 'soca', 'variational', 'socaincr2mom6.yaml') +s2mconfig = YAMLFile(path=socaincr2mom6_yaml_template) +s2mconfig = Template.substitute_structure(s2mconfig, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) +s2mconfig.save(socaincr2mom6_yaml) + ################################################################################ # Copy initial condition ics_list = [] diff --git a/utils/CMakeLists.txt b/utils/CMakeLists.txt new file mode 100644 index 000000000..b5cea21f1 --- /dev/null +++ b/utils/CMakeLists.txt @@ -0,0 +1,10 @@ +find_package( NetCDF REQUIRED COMPONENTS CXX) +find_package(oops REQUIRED) +find_package(atlas REQUIRED) +find_package(soca REQUIRED) + +ecbuild_add_executable( TARGET socaincr2mom6.x + SOURCES socaincr2mom6.cc ) + +target_compile_features( socaincr2mom6.x PUBLIC cxx_std_17) +target_link_libraries( socaincr2mom6.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) diff --git a/utils/socaincr2mom6.cc b/utils/socaincr2mom6.cc new file mode 100644 index 000000000..d8fa7c788 --- /dev/null +++ b/utils/socaincr2mom6.cc @@ -0,0 +1,8 @@ +#include "socaincr2mom6.h" +#include "oops/runs/Run.h" + +int main(int argc, char ** argv) { + oops::Run run(argc, argv); + gdasapp::SocaIncr2Mom6 socaincr2mom6; + return run.execute(socaincr2mom6); +} diff --git a/utils/socaincr2mom6.h b/utils/socaincr2mom6.h new file mode 100644 index 000000000..01e109988 --- /dev/null +++ b/utils/socaincr2mom6.h @@ -0,0 +1,113 @@ +#include +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "atlas/field.h" + +#include "oops/base/PostProcessor.h" +#include "oops/mpi/mpi.h" +#include "oops/runs/Application.h" +#include "oops/util/DateTime.h" +#include "oops/util/Duration.h" +#include "oops/util/Logger.h" + +#include "soca/State/State.h" +#include "soca/Geometry/Geometry.h" +#include "soca/Increment/Increment.h" + +namespace gdasapp { + + class SocaIncr2Mom6 : public oops::Application { + public: + explicit SocaIncr2Mom6(const eckit::mpi::Comm & comm = oops::mpi::world()) + : Application(comm) {} + static const std::string classname() {return "gdasapp::SocaIncr2Mom6";} + + int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { + + /// Setup the soca geometry + const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); + oops::Log::info() << "geometry: " << std::endl << geomConfig << std::endl; + const soca::Geometry geom(geomConfig, this->getComm()); + + /// Setup the vertical geometry from the background (layer thicknesses) + // get date + std::string strdt; + fullConfig.get("date", strdt); + util::DateTime dt = util::DateTime(strdt); + + // get layer thickness variable name + oops::Variables layerVar(fullConfig, "layers variable"); + ASSERT(layerVar.size() == 1); + + // read layer thicknesses from the relevant background + soca::Increment layerThick(geom, layerVar, dt); + const eckit::LocalConfiguration vertGeomConfig(fullConfig, "vertical geometry"); + layerThick.read(vertGeomConfig); + oops::Log::info() << "layerThick: " << std::endl << layerThick << std::endl; + + // Setup the soca increment + // get the increment variables + oops::Variables socaIncrVar(fullConfig, "increment variables"); + ASSERT(socaIncrVar.size() >= 1); + + // read the soca increment + soca::Increment socaIncr(geom, socaIncrVar, dt); + const eckit::LocalConfiguration socaIncrConfig(fullConfig, "soca increment"); + socaIncr.read(socaIncrConfig); + oops::Log::info() << "socaIncr: " << std::endl << socaIncr << std::endl; + + /// Create the MOM6 IAU increment + // concatenate variables + oops::Variables mom6IauVar(socaIncrVar); + mom6IauVar += layerVar; + oops::Log::info() << "mom6IauVar: " << std::endl << mom6IauVar << std::endl; + + // append layer variable to soca increment + atlas::FieldSet socaIncrFs; + socaIncr.toFieldSet(socaIncrFs); + socaIncr.updateFields(mom6IauVar); + oops::Log::info() << "MOM6 incr: " << std::endl << socaIncr << std::endl; + + // pad layer increment with zeros + atlas::FieldSet socaLayerThickFs; + layerThick.toFieldSet(socaLayerThickFs); + layerThick.updateFields(mom6IauVar); + oops::Log::info() << "h: " << std::endl << layerThick << std::endl; + + // append layer thinckness to increment + socaIncr += layerThick; + oops::Log::info() << "MOM6 IAU increment: " << std::endl << socaIncr << std::endl; + + // Save MOM6 IAU Increment + const eckit::LocalConfiguration mom6IauConfig(fullConfig, "mom6 iau increment"); + socaIncr.write(mom6IauConfig); + + // TODO: the "checkpoint" script expects the ocean increment output to + // be in "inc.nc". Remove what's below, eventually + std::string datadir; + mom6IauConfig.get("datadir", datadir); + std::filesystem::path pathToResolve = datadir; + std::string exp; + mom6IauConfig.get("exp", exp); + std::string outputType; + mom6IauConfig.get("type", outputType); + std::string incrFname = std::filesystem::canonical(pathToResolve); + incrFname += "/ocn." + exp + "." + outputType + "." + dt.toString() + ".nc"; + const char* charPtr = incrFname.c_str(); + oops::Log::info() << "rename: " << incrFname << " to " << "inc.nc" << std::endl; + int result = std::rename(charPtr, "inc.nc"); + + return result; + } + // ----------------------------------------------------------------------------- + private: + std::string appname() const { + return "gdasapp::SocaIncr2Mom6"; + } + // ----------------------------------------------------------------------------- + }; + +} // namespace gdasapp