Use VADER for insitu/potential TL/AD

src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90

@@ -46,18 +46,32 @@ subroutine soca_model2geovals_linear_changevar_f90(c_key_geom, c_key_dxin, c_key
 
   ! identity operators
   do i=1, size(dxout%fields)
-    call dxin%get(dxout%fields(i)%metadata%name, field)
-    if (dxout%fields(i)%name == field%metadata%name .or. &
-        dxout%fields(i)%name == field%metadata%getval_name) then
-      dxout%fields(i)%val(:,:,:) =  field%val(:,:,:)  !< full field
-    elseif (field%metadata%getval_name_surface == dxout%fields(i)%name) then
-      dxout%fields(i)%val(:,:,1) = field%val(:,:,1) !< surface only of a 3D field
-
-    else
-      call abor1_ftn( 'error in soca_model2geovals_linear_changevar_f90 processing ' &
-                       // dxout%fields(i)%name )
-    endif
-
+    select case (dxout%fields(i)%name)
+
+    ! TODO: These should NOT be needed in an increment. Need to make some changes to vader...
+    case ( 'sea_area_fraction' )
+      dxout%fields(i)%val(:,:,1) = 0.0    
+    case ( 'latitude' )
+      dxout%fields(i)%val(:,:,1) = 0.0
+    case ( 'longitude' )
+      dxout%fields(i)%val(:,:,1) = 0.0
+    case ( 'sea_water_depth')
+      dxout%fields(i)%val(:,:,1) = 0.0
+      ! call geom%thickness2depth(geom%h, dxout%fields(i)%val)
+
+    case default      
+      call dxin%get(dxout%fields(i)%metadata%name, field)
+      if (dxout%fields(i)%name == field%metadata%name .or. &
+          dxout%fields(i)%name == field%metadata%getval_name) then
+        dxout%fields(i)%val(:,:,:) =  field%val(:,:,:)  !< full field
+      elseif (field%metadata%getval_name_surface == dxout%fields(i)%name) then
+        dxout%fields(i)%val(:,:,1) = field%val(:,:,1) !< surface only of a 3D field
+
+      else
+        call abor1_ftn( 'error in soca_model2geovals_linear_changevar_f90 processing ' &
+                        // dxout%fields(i)%name )
+      endif
+    end select
   end do
 end subroutine
 

src/soca/LinearVariableChange/LinearVariableChange.cc

@@ -5,10 +5,12 @@
  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
  */
 
+#include <map>
 #include <ostream>
 #include <string>
 
 #include "soca/LinearVariableChange/LinearVariableChange.h"
+#include "soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.h"
 
 #include "soca/Geometry/Geometry.h"
 #include "soca/Increment/Increment.h"
@@ -22,10 +24,26 @@ namespace soca {
 
 // -----------------------------------------------------------------------------
 
+std::map<std::string, std::vector<std::string>> SocaLinVaderCookbook {
+  {"sea_water_temperature", {"SeaWaterTemperature_A"}},
+};
+
+// -----------------------------------------------------------------------------
+
 LinearVariableChange::LinearVariableChange(const Geometry & geom,
                                            const eckit::Configuration & config)
-  : geom_(geom), params_(), linVarChas_() {
+  : geom_(geom), params_(), linVarChas_(), vader_(), default_(false) {
     params_.deserialize(config);
+
+  // setup vader
+  eckit::LocalConfiguration vaderConfig, vaderCookbookConfig;
+  for (auto kv : SocaLinVaderCookbook) vaderCookbookConfig.set(kv.first, kv.second);
+  vaderConfig.set(vader::configCookbookKey, vaderCookbookConfig);
+  vader_.reset(new vader::Vader(params_.vader, vaderConfig));
+
+  const boost::optional<std::vector<LinearVariableChangeParametersWrapper>> &
+    linVarChgs = params_.linearVariableChangesWrapper;
+  default_ = linVarChgs == boost::none;
 }
 
 // -----------------------------------------------------------------------------
@@ -36,14 +54,56 @@ LinearVariableChange::~LinearVariableChange() {}
 
 void LinearVariableChange::changeVarTraj(const State & xfg,
                                          const oops::Variables & vars) {
-  Log::trace() << "LinearVariableChange::setTrajectory starting" << std::endl;
+  Log::trace() << "LinearVariableChange::changeVarTraj starting, default: "<<default_ << std::endl;
 
-  // TODO(travis): do something with vars?
+  Log::debug() << "LinearVariableChange::changeVarTraj vars in: "
+               << xfg.variables() << std::endl;
+  Log::debug() << "LinearVariableChange::changeVarTraj vars out: "
+               << vars << std::endl;
+
+  // The following is TEMPORARY.
+  // ----------------------------------------------------------------------------
+  // We need to do some variable renaming BEFORE we run VADER.
+  // Eventually, we will internally rename these variables when they are
+  // first loaded in so that we won't have to worry about it here.
+  State xfg2(xfg);
+  if (default_ && vars.has("sea_water_temperature")) {
+    Log::debug() << "LinearVariableChange::changeVarTraj Pre-VADER variable changes. " << std::endl;
+    oops::Variables preVaderVars(std::vector<std::string>{
+      "latitude",
+      "longitude",
+      "sea_water_potential_temperature",
+      "sea_water_salinity",
+      "sea_water_depth",
+      "sea_area_fraction",
+      });
+    preVaderVars += xfg.variables();
+    xfg2.updateFields(preVaderVars);
+    Model2GeoVaLs m2gv(xfg.geometry(), params_.toConfiguration());
+    m2gv.changeVar(xfg, xfg2);
+    Log::debug() << "LinearVariableChange::changeVarTraj variables after var change: "
+        << xfg2.variables() << std::endl;
+  }
+
+  // call Vader
+  // ----------------------------------------------------------------------------
+  Log::debug() << "LinearVariableChange::changeVarTraj VADER variable changes. " << std::endl;
+  oops::Variables varsFilled = xfg2.variables();
+  oops::Variables varsVader = vars;
+  varsVader -= varsFilled; // pass only the needed variables
+  atlas::FieldSet xfs;
+  xfg2.toFieldSet(xfs);
+  varsVaderPopulates_ = vader_->changeVarTraj(xfs, varsVader);
+  varsFilled += varsVaderPopulates_;
+  xfg2.updateFields(varsFilled);
+  xfg2.fromFieldSet(xfs);
+  Log::debug() << "LinearVariableChange::changeVarTraj variables after var change: "
+    << xfg2.variables() << std::endl;
 
   // TODO(travis): this is not ideal. We are saving the first trajectory and
   // assuming it is the background. This should all get ripped out when the
   // variable changes that rely on the background are dealt with properly.
-  if (!bkg_) { bkg_.reset(new State(xfg)); }
+  if (!bkg_) { bkg_.reset(new State(xfg2)); }
 
   const boost::optional<std::vector<LinearVariableChangeParametersWrapper>> &
     linVarChgs = params_.linearVariableChangesWrapper;
@@ -59,13 +119,13 @@ void LinearVariableChange::changeVarTraj(const State & xfg,
             linVarChaParWra.linearVariableChangeParameters;
       // Add linear variable change to vector
       linVarChas_.push_back(
-        LinearVariableChangeFactory::create(*bkg_, xfg, geom_, linVarChaPar));
+        LinearVariableChangeFactory::create(*bkg_, xfg2, geom_, linVarChaPar));
     }
-  } else {
+  } else {  
     // No variable changes were specified, use the default (LinearModel2GeoVaLs)
     eckit::LocalConfiguration conf;
     conf.set("linear variable change name", "default");
-    linVarChas_.push_back(LinearVariableChangeFactory::create(*bkg_, xfg, geom_,
+    linVarChas_.push_back(LinearVariableChangeFactory::create(*bkg_, xfg2, geom_,
       oops::validateAndDeserialize<GenericLinearVariableChangeParameters>(
         conf)));
   }
@@ -76,16 +136,54 @@ void LinearVariableChange::changeVarTraj(const State & xfg,
 
 void LinearVariableChange::changeVarTL(Increment & dx,
                                        const oops::Variables & vars) const {
-  Log::trace() << "LinearVariableChange::multiply starting" << std::endl;
+  Log::trace() << "LinearVariableChange::changeVarTL starting, default: " << default_ << std::endl;
+
+  Log::debug() << "LinearVariableChange::changeVarTL vars in: "
+               << dx.variables() << std::endl;
+  Log::debug() << "LinearVariableChange::changeVarTL vars out: "
+               << vars << std::endl;
 
-  // If all variables already in incoming state just remove the no longer
-  // needed fields
-  // if (hasAllFields) {
-  //   dx.updateFields(vars);
-  //   oops::Log::trace() << "VariableChange::changeVar done (identity)"
-  //                      << std::endl;
-  //   return;
-  // }
+
+  // The following is TEMPORARY.
+  // ----------------------------------------------------------------------------
+  // We need to do some variable renaming BEFORE we run VADER.
+  // Eventually, we will internally rename these variables when they are
+  // first loaded in so that we won't have to worry about it here.
+  if (default_ && vars.has("sea_water_temperature")) {
+    Log::debug() << "LinearVariableChange::changeVarTL Pre-VADER variable changes. " << std::endl;
+    oops::Variables preVaderVars(std::vector<std::string>{
+      "latitude",
+      "longitude",
+      "sea_water_potential_temperature",
+      "sea_water_salinity",
+      "sea_water_depth",
+      "sea_area_fraction"});
+    preVaderVars += dx.variables();    
+    Increment preVader(dx.geometry(), preVaderVars, dx.time());
+
+    for (icst_ it = linVarChas_.begin(); it != linVarChas_.end(); ++it) {
+      it->multiply(dx, preVader);
+      dx.updateFields(preVaderVars);
+      dx = preVader;
+    }
+    Log::debug() << "LinearVariableChange::changeVarTL variables after var change: "
+                << dx.variables() << std::endl;  
+  }
+
+
+  // call Vader
+  // ----------------------------------------------------------------------------
+  Log::debug() << "LinearVariableChange::changeVarTL VADER variable changes. " << std::endl;
+  atlas::FieldSet xfs;
+  dx.toFieldSet(xfs);
+  oops::Variables varsFilled = dx.variables();
+  oops::Variables varsVader = vars;
+  varsVader -= varsFilled;
+  varsFilled += vader_->changeVarTL(xfs, varsVader);
+  dx.updateFields(varsFilled);
+  dx.fromFieldSet(xfs);
+  Log::debug() << "LinearVariableChange::changeVarTL variables after var change: "
+               << dx.variables() << std::endl;
 
   // Create output state
   Increment dxout(dx.geometry(), vars, dx.time());
@@ -98,12 +196,6 @@ void LinearVariableChange::changeVarTL(Increment & dx,
      dx = dxout;
   }
 
-  // Allocate any extra fields and remove fields no longer needed
-  // dx.updateFields(vars);
-
-  // Copy data from temporary state
-  // dx = dxout;
-
   Log::trace() << "LinearVariableChange::multiply done" << std::endl;
 }
 
@@ -132,7 +224,42 @@ void LinearVariableChange::changeVarInverseTL(Increment & dx,
 
 void LinearVariableChange::changeVarAD(Increment & dx,
                                        const oops::Variables & vars) const {
-  Log::trace() << "LinearVariableChange::multiplyAD starting" << std::endl;
+  Log::trace() << "LinearVariableChange::changeVarAD starting" << std::endl;
+
+  Log::debug() << "LinearVariableChange::changeVarAD vars in: "
+               << dx.variables() << std::endl;
+  Log::debug() << "LinearVariableChange::changeVarAD vars out: "
+               << vars << std::endl;
+
+  // NOTE: the IF is temporary, we need to do some variable renaming afterward
+  if (default_ && dx.variables().has("sea_water_temperature")) {    
+    // call vader    
+    Log::debug() << "LinearVariableChange::changeVarAD VADER variable changes. " << std::endl;    
+    oops::Variables varsToAdj(std::vector<std::string>{
+      "sea_water_potential_temperature"});
+    oops::Variables varsToDrop(std::vector<std::string>{
+      "sea_water_temperature"});
+
+    // run Vader    
+    oops::Variables vaderVars = dx.variables();
+    vaderVars += varsToAdj;
+    dx.updateFields(vaderVars);
+
+    atlas::FieldSet xfs;
+    dx.toFieldSet(xfs);
+
+    oops::Variables varsVaderWillAdjoint = varsVaderPopulates_;
+    vader_->changeVarAD(xfs, varsVaderWillAdjoint);
+    ASSERT(varsVaderWillAdjoint.size() == 0);
+
+    vaderVars -= varsToDrop;
+    dx.updateFields(vaderVars);
+    dx.fromFieldSet(xfs);
+    Log::debug() << "LinearVariableChange::changeVarTL variables after var change: "
+                 << dx.variables() << std::endl;  
+  }
+
+
   Increment dxout(dx.geometry(), vars, dx.time());
 
   // Call variable change(s)

src/soca/LinearVariableChange/LinearVariableChange.h

@@ -21,6 +21,8 @@
 #include "oops/util/parameters/RequiredParameter.h"
 #include "oops/util/Printable.h"
 
+#include "vader/vader.h"
+
 #include "soca/LinearVariableChange/Base/LinearVariableChangeBase.h"
 
 namespace soca {
@@ -37,6 +39,7 @@ class LinearVariableChangeParameters :
  public:
   oops::OptionalParameter<std::vector<LinearVariableChangeParametersWrapper>>
          linearVariableChangesWrapper{"linear variable changes", this};
+  oops::Parameter<vader::VaderParameters> vader{"vader", {}, this};         
 };
 
 // -----------------------------------------------------------------------------
@@ -68,6 +71,10 @@ class LinearVariableChange : public util::Printable {
   const Geometry & geom_;
   std::unique_ptr<State> bkg_;
   LinVarChaVec_ linVarChas_;
+
+  std::unique_ptr<vader::Vader> vader_;
+  oops::Variables varsVaderPopulates_;
+  bool default_;
 };
 
 // -----------------------------------------------------------------------------

src/soca/VariableChange/VariableChange.cc

@@ -76,7 +76,8 @@ void VariableChange::changeVar(State & x, const oops::Variables & vars) const {
       "longitude",
       "sea_water_potential_temperature",
       "sea_water_salinity",
-      "sea_water_depth"});
+      "sea_water_depth",
+      "sea_area_fraction"});
     preVaderVars += x.variables();
     State preVader(x.geometry(), preVaderVars, x.time());
     variableChange_->changeVar(x, preVader);

test/testinput/3dvar_godas.yml

@@ -192,7 +192,13 @@ cost function:
             obsfile: data_static/obs/prof.nc
         simulated variables: [waterTemperature]
       obs operator:
-        name: InsituTemperature
+        name: VertInterp
+        observation alias file: testinput/obsop_name_map.yml
+        variables:
+          - waterTemperature
+        vertical coordinate: sea_water_depth
+        observation vertical coordinate: depth
+        interpolation method: linear          
       obs error:
         covariance model: diagonal
       obs filters: