oops app for soca incr

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)

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

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

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 = []

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)

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);
+}

utils/socaincr2mom6.h

@@ -0,0 +1,113 @@
+#include <netcdf>
+#include <iostream>
+#include <filesystem>
+
+#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