diff --git a/src/orca-jedi/increment/Increment.cc b/src/orca-jedi/increment/Increment.cc index ed124a0..12be537 100644 --- a/src/orca-jedi/increment/Increment.cc +++ b/src/orca-jedi/increment/Increment.cc @@ -212,7 +212,7 @@ Increment & Increment::operator*=(const double & zz) { return *this; } -/// \brief Create increment from the different of two state objects. +/// \brief Create increment from the difference of two state objects. /// \param x1 State object. /// \param x2 State object subtracted. void Increment::diff(const State & x1, const State & x2) { @@ -227,6 +227,13 @@ void Increment::diff(const State & x1, const State & x2) { atlas::Field field2 = x2.getField(i); atlas::Field fieldi = incrementFields_[i]; + atlas::field::MissingValue mv1(field1); + bool has_mv1 = static_cast(mv1); + atlas::field::MissingValue mv2(field2); + bool has_mv2 = static_cast(mv2); + bool has_mv = has_mv1 || has_mv2; + oops::Log::debug() << "DJL Increment::diff mv1 " << mv1 << " mv2 " << mv2 << " has_mv " << has_mv << std::endl; + std::string fieldName1 = field1.name(); std::string fieldName2 = field2.name(); std::string fieldNamei = fieldi.name(); @@ -239,8 +246,15 @@ void Increment::diff(const State & x1, const State & x2) { auto field_viewi = atlas::array::make_view(fieldi); for (atlas::idx_t j = 0; j < field_viewi.shape(0); ++j) { for (atlas::idx_t k = 0; k < field_viewi.shape(1); ++k) { - if (!ghost(j)) { field_viewi(j, k) = field_view1(j, k) - field_view2(j, k); - } else { field_viewi(j, k) = 0; } + if (!ghost(j)) { + if (!has_mv1 || (has_mv1 && !mv1(field_view1(j,k)))) { + if (!has_mv2 || (has_mv2 && !mv2(field_view2(j,k)))) { + field_viewi(j, k) = field_view1(j, k) - field_view2(j, k); + } + } + } else { + field_viewi(j, k) = 0; + } } } } diff --git a/src/orca-jedi/interpolator/Interpolator.cc b/src/orca-jedi/interpolator/Interpolator.cc index 5305409..6ffcbd4 100644 --- a/src/orca-jedi/interpolator/Interpolator.cc +++ b/src/orca-jedi/interpolator/Interpolator.cc @@ -229,6 +229,7 @@ void Interpolator::apply(const oops::Variables& vars, const Increment& inc, auto field_view = atlas::array::make_view(tgt_field); atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]); bool has_mv = static_cast(mv); + oops::Log::debug() << "DJL Interpolator::apply mv " << mv << " has_mv " << has_mv << std::endl; for (std::size_t klev=0; klev < varSizes[jvar]; ++klev) { for (std::size_t iloc=0; iloc < nlocs_; iloc++) { if (has_mv && mv(field_view(iloc, klev))) { @@ -310,6 +311,8 @@ void Interpolator::applyAD(const oops::Variables& vars, Increment& inc, // field_view.assign(0.0); atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]); bool has_mv = static_cast(mv); + oops::Log::debug() << "DJL Interpolator::applyAD mv " << mv << " has_mv " << has_mv << std::endl; + for (std::size_t klev=0; klev < varSizes[jvar]; ++klev) { for (std::size_t iloc=0; iloc < nlocs_; iloc++) { if (has_mv && mv(field_view(iloc, klev))) { diff --git a/src/orca-jedi/state/State.cc b/src/orca-jedi/state/State.cc index 40595cb..f5d72a4 100644 --- a/src/orca-jedi/state/State.cc +++ b/src/orca-jedi/state/State.cc @@ -178,6 +178,9 @@ State & State::operator+=(const Increment & dx) { for (int i = 0; i< stateFields_.size();i++) { atlas::Field field = stateFields_[i]; + atlas::field::MissingValue mv(field); + bool has_mv = static_cast(mv); + atlas::Field fieldi = dx.incrementFields()[i]; std::string fieldName = field.name(); @@ -189,7 +192,11 @@ State & State::operator+=(const Increment & dx) { auto field_viewi = atlas::array::make_view(fieldi); for (atlas::idx_t j = 0; j < field_view.shape(0); ++j) { for (atlas::idx_t k = 0; k < field_view.shape(1); ++k) { - if (!ghost(j)) field_view(j, k) += field_viewi(j, k); + if (!ghost(j)) { + if (!has_mv || (has_mv && !mv(field_view(j,k)))) { + field_view(j, k) += field_viewi(j, k); + } + } } } }