diff --git a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 index 390db804d..7ba0e20bb 100644 --- a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 +++ b/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 diff --git a/src/soca/LinearVariableChange/LinearVariableChange.cc b/src/soca/LinearVariableChange/LinearVariableChange.cc index f4ef49436..cbf49e8c1 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.cc +++ b/src/soca/LinearVariableChange/LinearVariableChange.cc @@ -5,10 +5,12 @@ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. */ +#include #include #include #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> 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> & + 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: "<{ + "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> & 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( 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{ + "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{ + "sea_water_potential_temperature"}); + oops::Variables varsToDrop(std::vector{ + "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) diff --git a/src/soca/LinearVariableChange/LinearVariableChange.h b/src/soca/LinearVariableChange/LinearVariableChange.h index d1ef99409..a54170f50 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.h +++ b/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> linearVariableChangesWrapper{"linear variable changes", this}; + oops::Parameter vader{"vader", {}, this}; }; // ----------------------------------------------------------------------------- @@ -68,6 +71,10 @@ class LinearVariableChange : public util::Printable { const Geometry & geom_; std::unique_ptr bkg_; LinVarChaVec_ linVarChas_; + + std::unique_ptr vader_; + oops::Variables varsVaderPopulates_; + bool default_; }; // ----------------------------------------------------------------------------- diff --git a/src/soca/VariableChange/VariableChange.cc b/src/soca/VariableChange/VariableChange.cc index 34f899438..1d84139f7 100644 --- a/src/soca/VariableChange/VariableChange.cc +++ b/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); diff --git a/test/testinput/3dvar_godas.yml b/test/testinput/3dvar_godas.yml index 1852d0d24..ebb8baf7b 100644 --- a/test/testinput/3dvar_godas.yml +++ b/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: