From df0992a90639925f21c10446ab4c031c3a903457 Mon Sep 17 00:00:00 2001
From: Dan Lea <daniel.lea@metoffice.gov.uk>
Date: Fri, 20 Sep 2024 18:21:14 +0100
Subject: [PATCH] Add some additional code to respect fillvalues.

---
 src/orca-jedi/increment/Increment.cc       | 20 +++++++++++++++++---
 src/orca-jedi/interpolator/Interpolator.cc |  3 +++
 src/orca-jedi/state/State.cc               |  9 ++++++++-
 3 files changed, 28 insertions(+), 4 deletions(-)

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<bool>(mv1);
+    atlas::field::MissingValue mv2(field2);
+    bool has_mv2 = static_cast<bool>(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<double, 2>(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<double, 2>(tgt_field);
     atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]);
     bool has_mv = static_cast<bool>(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<bool>(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<bool>(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<double, 2>(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);
+          }
+        }
       }
     }
   }