EAMxx: fix accumul for non-diag fields

Author: mahf708

components/eamxx/src/physics/cosp/eamxx_cosp.cpp

@@ -113,6 +113,7 @@ void Cosp::initialize_impl (const RunType /* run_type */)
   for (const auto& field_name : vnames) {
     // the mask here is just the sunlit mask, so set it
     get_field_out(field_name).get_header().set_extra_data("mask_field", get_field_in("sunlit_mask"));
+    get_field_out(field_name).get_header().set_may_be_filled(true);
   }
 }
 

components/eamxx/src/share/io/scorpio_output.cpp

@@ -146,7 +146,7 @@ AtmosphereOutput (const ekat::Comm& comm, const ekat::ParameterList& params,
   EKAT_REQUIRE_MSG (not has_duplicates(m_fields_names),
       "[AtmosphereOutput] Error! One of the output yaml files has duplicate field entries.\n"
       " - yaml file: " + params.name() + "\n"
-      " - fields names; " + ekat::join(m_fields_names,",") + "\n");
+      " - fields names: " + ekat::join(m_fields_names,",") + "\n");
 
   // Check if remapping and if so create the appropriate remapper
   // Note: We currently support three remappers
@@ -174,6 +174,7 @@ AtmosphereOutput (const ekat::Comm& comm, const ekat::ParameterList& params,
   }
 
   // ... then add diagnostic fields
+  // ... and we also use this to init track avg cnt for regular fields
   init_diagnostics ();
 
   // Avg count only makes sense if we have
@@ -508,6 +509,16 @@ run (const std::string& filename,
     const auto& f_in  = fm_after_hr->get_field(field_name);
           auto& f_out = fm_scorpio->get_field(field_name);
 
+    // Safety check: if a field may contain fill values and we are computing an Average,
+    // we must have created an avg-count tracking field; otherwise division by the raw
+    // number of steps would bias the result wherever fill values occurred.
+    if (m_avg_type==OutputAvgType::Average && f_in.get_header().may_be_filled()) {
+      EKAT_REQUIRE_MSG(m_field_to_avg_count.count(field_name),
+        "[AtmosphereOutput::run] Error! Averaging a fill-aware field without avg-count tracking.\n"
+        " - field name : " + field_name + "\n"
+        "This indicates the field was marked may_be_filled after output initialization or tracking logic missed it." );
+    }
+
     switch (m_avg_type) {
       case OutputAvgType::Instant:
         f_out.deep_copy(f_in);  break; // Note: if f_in aliases f_out, this is a no-op
@@ -1003,6 +1014,26 @@ init_diagnostics ()
   for (const auto& fname : m_fields_names) {
     if (not m_field_mgrs[FromModel]->has_field(fname)) {
       create_diag(fname);
+    } else {
+      // This is a regular field, not a diagnostic.
+      // Still, we might need to do some extra setup, like for avg_count.
+      const auto& f = m_field_mgrs[FromModel]->get_field(fname);
+      // We need avg-count tracking for any averaged (non-instant) field that:
+      //  - supplies explicit mask info (mask_data or mask_field), OR
+      //  - is marked as potentially containing fill values (may_be_filled()).
+      // Without this, fill-aware updates skip fill_value during accumulation (good)
+      // but we would still divide by the raw nsteps, biasing the result low.
+      if (m_avg_type!=OutputAvgType::Instant) {
+        const bool has_mask = f.get_header().has_extra_data("mask_data") || f.get_header().has_extra_data("mask_field");
+        const bool may_be_filled = f.get_header().may_be_filled();
+        if (has_mask || may_be_filled) {
+          m_track_avg_cnt = true;
+          // Avoid duplicate insertion if already present (e.g., mask + filled both true)
+          if (m_field_to_avg_cnt_suffix.count(f.name())==0) {
+            m_field_to_avg_cnt_suffix.emplace(f.name(), "_" + f.name());
+          }
+        }
+      }
     }
   }
 }