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(); - - // get geovals - std::vector 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 vec(nlocs); + geovals.getAtLevel(vec, + oops::Variable{"sea_surface_height_above_geoid"}, 0); + + // ----------------------- + // 2. Observations + // ----------------------- std::vector 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(); auto accumCnt = odb_.distribution()->createAccumulator(); - 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 + 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;