From bab02a5ce1e3f5b14f00ea5fd9e59bf039a32bd0 Mon Sep 17 00:00:00 2001
From: Joao Souza <joao.souza@noaa.gov>
Date: Thu, 23 Apr 2026 11:26:49 -0600
Subject: [PATCH] Refactor ObsADT::Change georef to include only QCed
 observations

Modify the way the reference ADT is calculated to include only the observations that passed the QC steps defined in the IODA yaml.

This avoids contaminating the ADT that is assimilate with any wrog data that is filtered out in the QC.
---
 src/ufo/operators/marine/adt/ObsADT.cc | 61 ++++++++++++++++++--------
 1 file changed, 43 insertions(+), 18 deletions(-)

diff --git a/src/ufo/operators/marine/adt/ObsADT.cc b/src/ufo/operators/marine/adt/ObsADT.cc
index 78019bc3c..7c5921a72 100644
--- a/src/ufo/operators/marine/adt/ObsADT.cc
+++ b/src/ufo/operators/marine/adt/ObsADT.cc
@@ -61,38 +61,63 @@ void ObsADT::simulateObs(const GeoVaLs & geovals, ioda::ObsVector & ovec,
   oops::Log::trace() << "ObsADT::simulateObs start" << std::endl;
 
   const double missing = util::missingValue<double>();
-
-  // get geovals
-  std::vector<double> vec(ovec.nlocs());
-  geovals.getAtLevel(vec, oops::Variable{"sea_surface_height_above_geoid"}, 0);
-
-  // get obs
+  const size_t nlocs = ovec.nlocs();
+
+  // -----------------------
+  // 1. GeoVaLs
+  // -----------------------
+  std::vector<double> vec(nlocs);
+  geovals.getAtLevel(vec,
+                     oops::Variable{"sea_surface_height_above_geoid"}, 0);
+
+  // -----------------------
+  // 2. Observations
+  // -----------------------
   std::vector<double> obs;
-  odb_.get_db("ObsValue", ovec.varnames().variables()[operatorVarIndex_], obs);
+  odb_.get_db("ObsValue",
+              ovec.varnames().variables()[operatorVarIndex_],
+              obs);
+
+  // -----------------------
+  // 3. Compute global offset (QC-filtered)
+  // -----------------------
+  double offset = 0.0;
 
-  // calculate global offsets
-  double offset = 0;
   auto accumVal = odb_.distribution()->createAccumulator<double>();
   auto accumCnt = odb_.distribution()->createAccumulator<int>();
-  for (size_t jloc = 0; jloc < ovec.nlocs(); ++jloc) {
-    // TODO(travis) also remove obs that have failed QC before this point?
-    // Not doing it now to keep answers from changing with this PR.
-    if (obs[jloc] != missing && vec[jloc] != missing) {
+
+  for (size_t jloc = 0; jloc < nlocs; ++jloc) {
+    // ✅ CORRECT QC ACCESS FOR ObsDataVector<int>
+    const int qc = qc_flags[operatorVarIndex_][jloc];
+
+    if (qc == 0 &&
+        obs[jloc] != missing &&
+        vec[jloc] != missing &&
+        std::isfinite(obs[jloc]) &&
+        std::isfinite(vec[jloc])) {
+
       accumVal->addTerm(jloc, vec[jloc] - obs[jloc]);
       accumCnt->addTerm(jloc, 1);
     }
   }
-  auto count = accumCnt->computeResult();
+
+  const int count = accumCnt->computeResult();
+
   if (count > 0) {
     offset = accumVal->computeResult() / count;
   }
-  // oops::Log::debug() << "ObsADT::simulateObs offset: " << offset << std::endl;
 
-  // subtract offset from geoval
-  for (size_t jloc = 0; jloc < ovec.nlocs(); ++jloc) {
+  // -----------------------
+  // 4. Apply offset (no QC filtering here!)
+  // -----------------------
+  for (size_t jloc = 0; jloc < nlocs; ++jloc) {
     const size_t idx = jloc * ovec.nvars() + operatorVarIndex_;
+
     ovec[idx] = vec[jloc];
-    if (ovec[idx] != missing) ovec[idx] -= offset;
+
+    if (ovec[idx] != missing && std::isfinite(ovec[idx])) {
+      ovec[idx] -= offset;
+    }
   }
 
   oops::Log::trace() << "ObsADT::simulateObs done" << std::endl;