JEDI-based ensemble recentering and analysis calculation

mains/gdas.cc

@@ -14,6 +14,7 @@
 #include "ufo/instantiateObsFilterFactory.h"
 #include "ufo/ObsTraits.h"
 
+#include "oops/runs/AddIncrement.h"
 #include "oops/runs/ConvertToStructuredGrid.h"
 #include "oops/runs/ConvertState.h"
 #include "oops/runs/EnsMeanAndVariance.h"
@@ -54,6 +55,9 @@ int runApp(int argc, char** argv, const std::string traits, const std::string ap
   // Define a map from app names to lambda functions that create unique_ptr to Applications
   std::map<std::string, std::function<std::unique_ptr<oops::Application>()>> apps;
 
+  apps["addincrement"] = []() {
+      return std::make_unique<oops::AddIncrement<Traits>>();
+  };
   apps["converttostructuredgrid"] = []() {
       return std::make_unique<oops::ConvertToStructuredGrid<Traits>>();
   };
@@ -104,6 +108,7 @@ int main(int argc,  char ** argv) {
   // Check that the application is recognized
   // ----------------------------------------
   const std::set<std::string> validApps = {
+    "addincrement",
     "converttostructuredgrid",
     "convertstate",
     "ensmean",

parm/aero/jcb-base.yaml.j2

@@ -50,6 +50,8 @@ aero_npx_clim_b: {{ npx_clim_b | default(49, true) }}
 aero_npy_clim_b: {{ npy_clim_b | default(49, true) }}
 aero_npz_clim_b: {{ npz_anl | default(127, true) }}
 
+aero_variational_history_prefix: "{{ GPREFIX }}"
+
 aero_fv3jedi_files_path: ./fv3jedi  # Ideally this would be {{DATA}}/fv3jedi but FMS
 
 # Background

parm/atm/jcb-base.yaml.j2

@@ -12,6 +12,7 @@ app_path_observation_chronicle: "{{PARMgfs}}/gdas/jcb-gdas/observation_chronicle
 #final_increment_file: final_increment_gaussian
 final_increment_file: atmosphere_final_increment_cubed_sphere
 output_ensemble_increments_file: atmosphere_output_ensemble_increments_cubed_sphere
+posterior_output_file: atmosphere_posterior_output_cubed_sphere
 model_file: atmosphere_model_pseudo
 initial_condition_file: atmosphere_background  # Initial conditions for 4D apps is background
 background_error_file: "{{BERROR_YAML}}"
@@ -55,14 +56,18 @@ atmosphere_number_ensemble_members: {{ NMEM_ENS | default(1, true) }}
 atmosphere_layout_gsib_x: {{ layout_gsib_x | default(1, true) }}
 atmosphere_layout_gsib_y: {{ layout_gsib_y | default(1, true) }}
 
-
 # Forecasting
 atmosphere_forecast_timestep: "{{ BKG_TSTEP }}"
 
 # Write final increment on Guassian grid in variational
 atmosphere_final_increment_prefix: "./anl/atminc."
 
 atmosphere_variational_history_prefix: "{{GPREFIX}}"
+atmosphere_variational_analysis_prefix: "{{APREFIX}}"
+
+# IAU hours
+atmosphere_iau_hours: {{ IAUFHRS }}
+atmosphere_iau_times_iso: {{ iau_times_iso | default([current_cycle | to_isotime], true) }}
 
 # Observation things
 # ------------------
@@ -125,7 +130,7 @@ inflation_mult: 1.0
 
 # Driver
 driver_update_obs_config_with_geometry_info: true
-driver_save_posterior_mean: false
+driver_save_posterior_mean: true
 driver_save_posterior_ensemble: false
 driver_save_prior_mean: false
 driver_save_posterior_mean_increment: false
@@ -138,4 +143,5 @@ distribution_type: Halo
 atmosphere_ensemble_increment_prefix: "./anl/mem%{member}%/atminc."
 atmosphere_posterior_output_gaussian: "./mem%{member}%/atmanl."
 
-atmosphere_ensemble_history_prefix: "enkf{{GPREFIX}}"
+atmosphere_ensemble_history_prefix: "{{ GPREFIX_ENS }}"
+atmosphere_ensemble_analysis_prefix: "{{ APREFIX_ENS }}"

parm/snow/jcb-base.yaml.j2

@@ -50,6 +50,8 @@ snow_npx_anl: {{ npx_ges | default(49, true) }}
 snow_npy_anl: {{ npy_ges | default(49, true) }}
 snow_npz_anl: {{ npz_ges | default(127, true) }}
 
+snow_variational_history_prefix: "{{ GPREFIX }}"
+
 snow_fv3jedi_files_path: ./fv3jedi  # Ideally this would be {{DATA}}/fv3jedi but FMS
 snow_orog_files_path: "{{ FIXgfs }}/orog/{{ CASE }}"
 snow_orog_prefix: "{{ CASE }}.mx{{ OCNRES }}"

utils/fv3jedi/CMakeLists.txt

@@ -3,3 +3,15 @@ ecbuild_add_executable( TARGET fv3jedi_fv3inc.x
                         SOURCES fv3jedi_fv3inc.cc fv3jedi_fv3inc.h)
 target_compile_features( fv3jedi_fv3inc.x PUBLIC cxx_std_17)
 target_link_libraries( fv3jedi_fv3inc.x PUBLIC NetCDF::NetCDF_CXX oops fv3jedi)
+
+# Correction increment calculation (for ensemble recentering)
+ecbuild_add_executable( TARGET fv3jedi_correction_increment.x
+                        SOURCES fv3jedi_correction_increment.cc fv3jedi_correction_increment.h)
+target_compile_features( fv3jedi_correction_increment.x PUBLIC cxx_std_17)
+target_link_libraries( fv3jedi_correction_increment.x PUBLIC NetCDF::NetCDF_CXX oops fv3jedi)
+
+# Ensemble AddIncrement
+ecbuild_add_executable( TARGET fv3jedi_ensemble_add_increment.x
+                        SOURCES fv3jedi_ensemble_add_increment.cc)
+target_compile_features( fv3jedi_ensemble_add_increment.x PUBLIC cxx_std_17)
+target_link_libraries( fv3jedi_ensemble_add_increment.x PUBLIC NetCDF::NetCDF_CXX oops fv3jedi)

utils/fv3jedi/fv3jedi_correction_increment.cc

@@ -0,0 +1,8 @@
+#include "fv3jedi_correction_increment.h"
+#include "oops/runs/Run.h"
+
+int main(int argc,  char ** argv) {
+  oops::Run run(argc, argv);
+  gdasapp::CorrectionIncrement CorrectionIncrement;
+  return run.execute(CorrectionIncrement);
+}

utils/fv3jedi/fv3jedi_correction_increment.h

@@ -0,0 +1,120 @@
+#pragma once
+
+#include <iostream>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "eckit/config/LocalConfiguration.h"
+
+#include "fv3jedi/Geometry/Geometry.h"
+#include "fv3jedi/Increment/Increment.h"
+#include "fv3jedi/State/State.h"
+
+#include "oops/mpi/mpi.h"
+#include "oops/runs/Application.h"
+#include "oops/util/ConfigFunctions.h"
+#include "oops/util/DateTime.h"
+#include "oops/util/Duration.h"
+#include "oops/util/Logger.h"
+#include "oops/util/parameters/OptionalParameter.h"
+#include "oops/util/parameters/Parameter.h"
+#include "oops/util/parameters/Parameters.h"
+#include "oops/util/parameters/RequiredParameter.h"
+
+class ForecastHourParameters : public oops::Parameters {
+  OOPS_CONCRETE_PARAMETERS(ForecastHourParameters, Parameters)
+ public:
+  oops::RequiredParameter<std::string> datetimeStr{"datetime", \
+                                                   this};
+  oops::RequiredParameter<eckit::LocalConfiguration> detBkgConfig{"deterministic background", \
+                                                                  this};
+  oops::RequiredParameter<eckit::LocalConfiguration> varIncrConfig{"variational increment", \
+                                                                   this};
+  oops::RequiredParameter<eckit::LocalConfiguration> ensMeanAnlConfig{"ensemble mean analysis", \
+                                                                      this};
+  oops::RequiredParameter<eckit::LocalConfiguration> corIncrConfig{"correction increment", \
+                                                                   this};
+};
+
+class CorrectionIncrementParameters : public oops::Parameters {
+  OOPS_CONCRETE_PARAMETERS(CorrectionIncrementParameters, Parameters)
+ public:
+  oops::RequiredParameter<oops::Variables> \
+    incrVars{"increment variables", this};
+  oops::RequiredParameter<eckit::LocalConfiguration> \
+    varIncrGeomConfig{"variational increment geometry", this};
+  oops::RequiredParameter<eckit::LocalConfiguration> \
+    detBkgGeomConfig{"deterministic background geometry", this};
+  oops::RequiredParameter<eckit::LocalConfiguration> \
+    ensMeanAnlGeomConfig{"ensemble mean analysis geometry", this};
+  oops::RequiredParameter<eckit::LocalConfiguration> \
+    corIncrGeomConfig{"correction increment geometry", this};
+  oops::RequiredParameter<std::vector<ForecastHourParameters>> \
+    fcstHourParams{"forecast hours", this};
+};
+
+namespace gdasapp {
+
+  // Main application class
+  class CorrectionIncrement : public oops::Application {
+   public:
+    explicit CorrectionIncrement(const eckit::mpi::Comm & comm = oops::mpi::world())
+      : Application(comm) {}
+    static const std::string classname() {return "gdasapp::CorrectionIncrement";}
+
+    int execute(const eckit::Configuration & fullConfig) const {
+      // Deserialize parameters
+      CorrectionIncrementParameters params;
+      params.deserialize(fullConfig);
+
+      // Setup geometries
+      const fv3jedi::Geometry varIncrGeom(params.varIncrGeomConfig.value(), this->getComm());
+      const fv3jedi::Geometry detBkgGeom(params.detBkgGeomConfig.value(), this->getComm());
+      const fv3jedi::Geometry ensMeanAnlGeom(params.ensMeanAnlGeomConfig.value(), this->getComm());
+      const fv3jedi::Geometry corIncrGeom(params.corIncrGeomConfig.value(), this->getComm());
+
+      // Loop through forecast hours ("recenterings")
+      const int nhours = params.fcstHourParams.value().size();
+      for ( int ihour = 0; ihour < nhours; ihour++ ) {
+        const ForecastHourParameters fcstHourParams = params.fcstHourParams.value()[ihour];
+
+        // Get forecast time
+        const util::DateTime datetime(fcstHourParams.datetimeStr.value());
+
+        // Initialize background
+        fv3jedi::State xxBkgDet(detBkgGeom, params.incrVars.value(), datetime);
+        xxBkgDet.read(fcstHourParams.detBkgConfig.value());
+
+        // Initialize increment
+        fv3jedi::Increment dxVar(varIncrGeom, params.incrVars.value(), datetime);
+        dxVar.read(fcstHourParams.varIncrConfig.value());
+
+        // Initialize ensemble mean analysis
+        fv3jedi::State xxAnlEnsMean(ensMeanAnlGeom, params.incrVars.value(), datetime);
+        xxAnlEnsMean.read(fcstHourParams.ensMeanAnlConfig.value());
+
+        // Compute analysis
+        fv3jedi::State xxAnlVar(detBkgGeom, xxBkgDet);
+        xxAnlVar += dxVar;
+
+        // Interpolate full resolution analysis to ensemble resolution
+        fv3jedi::State xxAnlVarEnsRes(corIncrGeom, xxAnlVar);
+
+        // Compute correction increment
+        fv3jedi::Increment dxCor(corIncrGeom, params.incrVars.value(), datetime);
+        dxCor.diff(xxAnlVarEnsRes, xxAnlEnsMean);
+
+        // Write correction increment
+        dxCor.write(fcstHourParams.corIncrConfig.value());
+      }
+
+      return 0;
+    }
+
+   private:
+    std::string appname() const {
+      return "gdasapp::CorrectionIncrement";
+    }
+  };
+}  // namespace gdasapp

utils/fv3jedi/fv3jedi_ensemble_add_increment.cc

@@ -0,0 +1,11 @@
+#include "fv3jedi/Utilities/Traits.h"
+
+#include "oops/runs/AddIncrement.h"
+#include "oops/runs/EnsembleApplication.h"
+#include "oops/runs/Run.h"
+
+int main(int argc,  char ** argv) {
+  oops::Run run(argc, argv);
+  oops::EnsembleApplication<oops::AddIncrement<fv3jedi::Traits>> addIncrement;
+  return run.execute(addIncrement);
+}