EAMxx: horiz remappers don't deal with fillval well, leading to errors

This leads to errors in things producing fillvals, etc., like cosp; see PR #7992 for more context of potential fix needed, patched below

patch

```patch
From bd476dc25a7b899f9db96af29d9547af3c0582b2 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 15 Jan 2026 16:35:33 +0000
Subject: [PATCH 01/10] Initial plan

From 109f33ec7e4128f8dc0291b371494592c65734e6 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 15 Jan 2026 16:54:14 +0000
Subject: [PATCH 02/10] Add mask_value to COSP outputs and propagate in
remappers

Co-authored-by: mahf708 <122953255+mahf708@users.noreply.github.com>
---
components/eamxx/src/physics/cosp/eamxx_cosp.cpp | 1 +
components/eamxx/src/share/diagnostics/aodvis.cpp | 2 ++
.../src/share/diagnostics/field_at_pressure_level.cpp | 1 +
components/eamxx/src/share/remap/coarsening_remapper.cpp | 9 +++++++++
components/eamxx/src/share/remap/vertical_remapper.cpp | 5 +++++
5 files changed, 18 insertions(+)

diff --git a/components/eamxx/src/physics/cosp/eamxx_cosp.cpp b/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
index 039c23e1165b..5a4a4f4b235b 100644
--- a/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
+++ b/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
@@ -115,6 +115,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_extra_data("mask_value", constants::fill_value);
get_field_out(field_name).get_header().set_may_be_filled(true);
}
}
diff --git a/components/eamxx/src/share/diagnostics/aodvis.cpp b/components/eamxx/src/share/diagnostics/aodvis.cpp
index 0c669909efcf..751d9c7103f7 100644
--- a/components/eamxx/src/share/diagnostics/aodvis.cpp
+++ b/components/eamxx/src/share/diagnostics/aodvis.cpp
@@ -42,6 +42,8 @@ set_grids(const std::shared_ptr grids_manager)

void AODVis::initialize_impl(const RunType /*run_type*/) {
m_diagnostic_output.get_header().set_extra_data("mask_field", get_field_in("sunlit_mask"));
+ m_diagnostic_output.get_header().set_extra_data("mask_value", constants::fill_value);
+ m_diagnostic_output.get_header().set_may_be_filled(true);
}

void AODVis::compute_diagnostic_impl() {
diff --git a/components/eamxx/src/share/diagnostics/field_at_pressure_level.cpp b/components/eamxx/src/share/diagnostics/field_at_pressure_level.cpp
index 9447bf0fce0f..b96476d6eb2b 100644
--- a/components/eamxx/src/share/diagnostics/field_at_pressure_level.cpp
+++ b/components/eamxx/src/share/diagnostics/field_at_pressure_level.cpp
@@ -92,6 +92,7 @@ initialize_impl (const RunType /*run_type*/)
Field diag_mask(mask_fid);
diag_mask.allocate_view();
m_diagnostic_output.get_header().set_extra_data("mask_field",diag_mask);
+ m_diagnostic_output.get_header().set_extra_data("mask_value", constants::fill_value);
m_diagnostic_output.get_header().set_may_be_filled(true);

using stratts_t = std::map;
diff --git a/components/eamxx/src/share/remap/coarsening_remapper.cpp b/components/eamxx/src/share/remap/coarsening_remapper.cpp
index 1938d8de7008..68da6741f9db 100644
--- a/components/eamxx/src/share/remap/coarsening_remapper.cpp
+++ b/components/eamxx/src/share/remap/coarsening_remapper.cpp
@@ -134,6 +134,15 @@ registration_ends_impl ()
mask_idx = masks.size()-1;
}
tgt.get_header().set_extra_data("mask_field",masks[mask_idx].second);
+ // Also propagate mask_value if present
+ if (src.get_header().has_extra_data("mask_value")) {
+ auto src_mask_value = src.get_header().get_extra_data("mask_value");
+ tgt.get_header().set_extra_data("mask_value",src_mask_value);
+ }
+ // Also propagate may_be_filled flag
+ if (src.get_header().may_be_filled()) {
+ tgt.get_header().set_may_be_filled(true);
+ }
}

// Add all masks to the fields to remap
diff --git a/components/eamxx/src/share/remap/vertical_remapper.cpp b/components/eamxx/src/share/remap/vertical_remapper.cpp
index 62c3c32240ff..ed5d98329607 100644
--- a/components/eamxx/src/share/remap/vertical_remapper.cpp
+++ b/components/eamxx/src/share/remap/vertical_remapper.cpp
@@ -208,6 +208,7 @@ registration_ends_impl ()
" - tgt field name: " + tgt.name() + "\n");

tgt.get_header().set_extra_data("mask_field",mask);
+ tgt.get_header().set_extra_data("mask_value", constants::fill_value);

// Since we do mask (at top and/or bot), the tgt field MAY be contain fill_value entries
tgt.get_header().set_may_be_filled(true);
@@ -223,6 +224,10 @@ registration_ends_impl ()
auto src_mask = src.get_header().get_extra_data("mask_field");
tgt.get_header().set_extra_data("mask_field",src_mask);
}
+ if (src.get_header().has_extra_data("mask_value")) {
+ auto src_mask_value = src.get_header().get_extra_data("mask_value");
+ tgt.get_header().set_extra_data("mask_value",src_mask_value);
+ }
if (src.get_header().may_be_filled()) {
tgt.get_header().set_may_be_filled(true);
}

From 745d4b2a4194ef8cac7be914f89957a86f38384a Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 15 Jan 2026 17:11:21 +0000
Subject: [PATCH 03/10] Fix critical bug in coarsening_remapper rescale for
masked columns

This fixes the actual root cause of large COSP values. When rescaling fields with 1D masks (like sunlit_mask), masked columns (mask<=threshold) were not being set to fill_value for rank-2, rank-3, and rank-4 fields. This caused large fill values to remain in the output instead of being properly masked.

The fix adds else branches to set all components to fill_val when mask indicates the column is masked, matching the behavior of the rank-1 case.

Co-authored-by: mahf708 <122953255+mahf708@users.noreply.github.com>
---
.../src/share/remap/coarsening_remapper.cpp | 25 +++++++++++++++++++
1 file changed, 25 insertions(+)

diff --git a/components/eamxx/src/share/remap/coarsening_remapper.cpp b/components/eamxx/src/share/remap/coarsening_remapper.cpp
index 68da6741f9db..ce5b2a0a6e9f 100644
--- a/components/eamxx/src/share/remap/coarsening_remapper.cpp
+++ b/components/eamxx/src/share/remap/coarsening_remapper.cpp
@@ -298,6 +298,12 @@ rescale_masked_fields (const Field& x, const Field& mask) const
[&](const int j){
x_sub(j) /= mask;
});
+ } else {
+ // Set to fill_value when mask indicates this column is masked
+ Kokkos::parallel_for(Kokkos::TeamVectorRange(team,dim1),
+ [&](const int j){
+ x_sub(j) = fill_val;
+ });
}
} else {
auto m_sub = ekat::subview(mask_2d,icol);
@@ -342,6 +348,15 @@ rescale_masked_fields (const Field& x, const Field& mask) const
auto x_sub = ekat::subview(x_view,icol,j);
x_sub(k) /= mask;
});
+ } else {
+ // Set to fill_value when mask indicates this column is masked
+ Kokkos::parallel_for(Kokkos::TeamVectorRange(team,dim1*dim2),
+ [&](const int idx){
+ const int j = idx / dim2;
+ const int k = idx % dim2;
+ auto x_sub = ekat::subview(x_view,icol,j);
+ x_sub(k) = fill_val;
+ });
}
} else {
auto m_sub = ekat::subview(mask_2d,icol);
@@ -392,6 +407,16 @@ rescale_masked_fields (const Field& x, const Field& mask) const
auto x_sub = ekat::subview(x_view,icol,j,k);
x_sub(l) /= mask;
});
+ } else {
+ // Set to fill_value when mask indicates this column is masked
+ Kokkos::parallel_for(Kokkos::TeamVectorRange(team,dim1*dim2*dim3),
+ [&](const int idx){
+ const int j = (idx / dim3) / dim2;
+ const int k = (idx / dim3) % dim2;
+ const int l = idx % dim3;
+ auto x_sub = ekat::subview(x_view,icol,j,k);
+ x_sub(l) = fill_val;
+ });
}
} else {
auto m_sub = ekat::subview(mask_2d,icol);

From ace07b8f5d2227151346b5b3fdc31dd07c25478a Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 15 Jan 2026 17:26:25 +0000
Subject: [PATCH 04/10] Add mask_value propagation to IO fields - critical
missing piece

This was the actual root cause! IO fields were getting may_be_filled but NOT mask_value, causing fill-aware update operations to fail. Without mask_value, the field update methods couldn't identify which values were fill values, leading to large fill values being accumulated and divided during temporal averaging, producing the observed 10^31-10^32 values.

Co-authored-by: mahf708 <122953255+mahf708@users.noreply.github.com>
---
components/eamxx/src/share/io/scorpio_output.cpp | 6 ++++++
1 file changed, 6 insertions(+)

diff --git a/components/eamxx/src/share/io/scorpio_output.cpp b/components/eamxx/src/share/io/scorpio_output.cpp
index 69bb3d1fbfad..a17ccac8ec30 100644
--- a/components/eamxx/src/share/io/scorpio_output.cpp
+++ b/components/eamxx/src/share/io/scorpio_output.cpp
@@ -35,6 +35,12 @@ transfer_extra_data(const scream::Field &src, scream::Field &tgt)
if (src.get_header().may_be_filled()) {
tgt.get_header().set_may_be_filled(true);
}
+
+ // Transfer mask_value if present
+ if (src.get_header().has_extra_data("mask_value")) {
+ auto src_mask_value = src.get_header().get_extra_data("mask_value");
+ tgt.get_header().set_extra_data("mask_value", src_mask_value);
+ }
};

// Note: this is also declared in eamxx_scorpio_interface.cpp. Move it somewhere else?

From f5672ecd22f4ef8518d34281efbbcbfe44cb59f3 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 15 Jan 2026 18:06:18 +0000
Subject: [PATCH 05/10] Fix build error and add COSP dimension coordinates

- Fix build error: qualify Real type with scream:: namespace in scorpio_output.cpp
- Add COSP dimension values and bounds (cosp_tau, cosp_prs, cosp_cth) as geometry data
- Uses standard ISCCP/MODIS/MISR bin definitions for e3sm_diags compatibility
- Tau bins: log-spaced optical depth bins (0.3 to 100000)
- Pressure bins: 7 bins from 50-1000 hPa
- Height bins: 16 MISR bins from 0-21 km

Co-authored-by: mahf708 <122953255+mahf708@users.noreply.github.com>
---
.../eamxx/src/physics/cosp/eamxx_cosp.cpp | 125 ++++++++++++++++++
.../eamxx/src/share/io/scorpio_output.cpp | 2 +-
2 files changed, 126 insertions(+), 1 deletion(-)

diff --git a/components/eamxx/src/physics/cosp/eamxx_cosp.cpp b/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
index 5a4a4f4b235b..4fca050ea250 100644
--- a/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
+++ b/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
@@ -44,6 +44,7 @@ void Cosp::set_grids(const std::shared_ptr grids_manager)
auto micron = micro*m;
auto m2 = pow(m, 2);
auto s2 = pow(s, 2);
+ auto hPa = hecto*Pa; // hectopascal (100 Pa)

m_grid = grids_manager->get_grid("physics");
const auto& grid_name = m_grid->name();
@@ -102,6 +103,130 @@ void Cosp::set_grids(const std::shared_ptr grids_manager)
m_z_int = Field(FieldIdentifier("z_int",scalar3d_int,m,grid_name));
m_z_mid.allocate_view();
m_z_int.allocate_view();
+
+ // Add COSP dimension coordinate values and bounds as geometry data
+ // These are needed for e3sm_diags to produce COSP diagnostics
+ using namespace ShortFieldTagsNames;
+
+ // COSP tau bins (optical depth) - standard ISCCP/MODIS bins
+ // Bin edges: 0.3, 1.3, 3.6, 9.4, 23, 60, 379
+ if (not m_grid->has_geometry_data("cosp_tau_bnds")) {
+ FieldLayout tau_bnds_layout({CMP,CMP},{m_num_tau,2},{"cosp_tau","nbnd"});
+ Field tau_bnds(FieldIdentifier("cosp_tau_bnds", tau_bnds_layout, nondim, grid_name));
+ tau_bnds.allocate_view();
+ auto tau_bnds_h = tau_bnds.get_view();
+
+ // Standard ISCCP/MODIS tau bin boundaries
+ tau_bnds_h(0,0) = 0.3; tau_bnds_h(0,1) = 1.3;
+ tau_bnds_h(1,0) = 1.3; tau_bnds_h(1,1) = 3.6;
+ tau_bnds_h(2,0) = 3.6; tau_bnds_h(2,1) = 9.4;
+ tau_bnds_h(3,0) = 9.4; tau_bnds_h(3,1) = 23.0;
+ tau_bnds_h(4,0) = 23.0; tau_bnds_h(4,1) = 60.0;
+ tau_bnds_h(5,0) = 60.0; tau_bnds_h(5,1) = 379.0;
+ tau_bnds_h(6,0) = 379.0; tau_bnds_h(6,1) = 100000.0; // effectively infinity
+
+ tau_bnds.sync_to_dev();
+ m_grid->set_geometry_data(tau_bnds);
+ }
+
+ // COSP tau centers
+ if (not m_grid->has_geometry_data("cosp_tau")) {
+ auto tau_bnds = m_grid->get_geometry_data("cosp_tau_bnds");
+ auto tau_bnds_h = tau_bnds.get_view();
+
+ FieldLayout tau_layout({CMP},{m_num_tau},{"cosp_tau"});
+ Field tau(FieldIdentifier("cosp_tau", tau_layout, nondim, grid_name));
+ tau.allocate_view();
+ auto tau_h = tau.get_view();
+
+ // Geometric mean of bin edges (log-space midpoint)
+ for (int i=0; iset_geometry_data(tau);
+ }
+
+ // COSP pressure bins (cloud-top pressure in hPa) - standard ISCCP bins
+ // Bin edges: 50, 180, 310, 440, 560, 680, 800, 1000 hPa
+ if (not m_grid->has_geometry_data("cosp_prs_bnds")) {
+ FieldLayout prs_bnds_layout({CMP,CMP},{m_num_ctp,2},{"cosp_prs","nbnd"});
+ Field prs_bnds(FieldIdentifier("cosp_prs_bnds", prs_bnds_layout, hPa, grid_name));
+ prs_bnds.allocate_view();
+ auto prs_bnds_h = prs_bnds.get_view();
+
+ // Standard ISCCP pressure bin boundaries (hPa)
+ prs_bnds_h(0,0) = 50.0; prs_bnds_h(0,1) = 180.0;
+ prs_bnds_h(1,0) = 180.0; prs_bnds_h(1,1) = 310.0;
+ prs_bnds_h(2,0) = 310.0; prs_bnds_h(2,1) = 440.0;
+ prs_bnds_h(3,0) = 440.0; prs_bnds_h(3,1) = 560.0;
+ prs_bnds_h(4,0) = 560.0; prs_bnds_h(4,1) = 680.0;
+ prs_bnds_h(5,0) = 680.0; prs_bnds_h(5,1) = 800.0;
+ prs_bnds_h(6,0) = 800.0; prs_bnds_h(6,1) = 1000.0;
+
+ prs_bnds.sync_to_dev();
+ m_grid->set_geometry_data(prs_bnds);
+ }
+
+ // COSP pressure centers
+ if (not m_grid->has_geometry_data("cosp_prs")) {
+ auto prs_bnds = m_grid->get_geometry_data("cosp_prs_bnds");
+ auto prs_bnds_h = prs_bnds.get_view();
+
+ FieldLayout prs_layout({CMP},{m_num_ctp},{"cosp_prs"});
+ Field prs(FieldIdentifier("cosp_prs", prs_layout, hPa, grid_name));
+ prs.allocate_view();
+ auto prs_h = prs.get_view();
+
+ // Midpoint of pressure bins
+ for (int i=0; iset_geometry_data(prs);
+ }
+
+ // COSP height bins (cloud-top height in m) - standard MISR bins
+ // 16 bins from 0 to 20 km
+ if (not m_grid->has_geometry_data("cosp_cth_bnds")) {
+ FieldLayout cth_bnds_layout({CMP,CMP},{m_num_cth,2},{"cosp_cth","nbnd"});
+ Field cth_bnds(FieldIdentifier("cosp_cth_bnds", cth_bnds_layout, m, grid_name));
+ cth_bnds.allocate_view();
+ auto cth_bnds_h = cth_bnds.get_view();
+
+ // Standard MISR height bin boundaries (m) - 16 bins of varying width
+ // These follow the CFMIP/COSP standard MISR bins
+ const Real cth_edges[17] = {0.0, 500.0, 1000.0, 1500.0, 2000.0, 2500.0, 3000.0, 4000.0,
+ 5000.0, 7000.0, 9000.0, 11000.0, 13000.0, 15000.0, 17000.0, 19000.0, 21000.0};
+ for (int i=0; iset_geometry_data(cth_bnds);
+ }
+
+ // COSP height centers
+ if (not m_grid->has_geometry_data("cosp_cth")) {
+ auto cth_bnds = m_grid->get_geometry_data("cosp_cth_bnds");
+ auto cth_bnds_h = cth_bnds.get_view();
+
+ FieldLayout cth_layout({CMP},{m_num_cth},{"cosp_cth"});
+ Field cth(FieldIdentifier("cosp_cth", cth_layout, m, grid_name));
+ cth.allocate_view();
+ auto cth_h = cth.get_view();
+
+ // Midpoint of height bins
+ for (int i=0; iset_geometry_data(cth);
+ }
}

// =========================================================================================
diff --git a/components/eamxx/src/share/io/scorpio_output.cpp b/components/eamxx/src/share/io/scorpio_output.cpp
index a17ccac8ec30..6db387366645 100644
--- a/components/eamxx/src/share/io/scorpio_output.cpp
+++ b/components/eamxx/src/share/io/scorpio_output.cpp
@@ -38,7 +38,7 @@ transfer_extra_data(const scream::Field &src, scream::Field &tgt)

// Transfer mask_value if present
if (src.get_header().has_extra_data("mask_value")) {
- auto src_mask_value = src.get_header().get_extra_data("mask_value");
+ auto src_mask_value = src.get_header().get_extra_data("mask_value");
tgt.get_header().set_extra_data("mask_value", src_mask_value);
}
};

From beaefe5f6129a004bbcc28e22c6d6ed058138fca Mon Sep 17 00:00:00 2001
From: mahf708
Date: Thu, 15 Jan 2026 17:32:02 -0800
Subject: [PATCH 06/10] properly persist rad cld opt depth

---
.../src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp | 7 +++++--
1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp b/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
index 301f810c5f78..2d13bc4b09e3 100644
--- a/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
+++ b/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
@@ -614,8 +614,6 @@ void RRTMGPRadiation::run_impl (const double dt) {
auto d_dtau105 = get_field_out("dtau105").get_view();
auto d_sunlit = get_field_out("sunlit_mask").get_view();

- Kokkos::deep_copy(d_dtau067,0.0);
- Kokkos::deep_copy(d_dtau105,0.0);
// Outputs for AeroCom cloud-top diagnostics
auto d_T_mid_at_cldtop = get_field_out("T_mid_at_cldtop").get_view();
auto d_p_mid_at_cldtop = get_field_out("p_mid_at_cldtop").get_view();
@@ -643,6 +641,11 @@ void RRTMGPRadiation::run_impl (const double dt) {
auto update_rad = scream::rrtmgp::radiation_do(m_rad_freq_in_steps, ts.get_num_steps());

if (update_rad) {
+ // Init these to zero inside the update_rad conditional
+ // TODO why do we need this anyway???
+ // TODO figure out a way to get rid of this!!!
+ Kokkos::deep_copy(d_dtau067,0.0);
+ Kokkos::deep_copy(d_dtau105,0.0);
// On each chunk, we internally "reset" the GasConcs object to subview the concs 3d array
// with the correct ncol dimension. So let's keep a copy of the original (ref-counted)
// array, to restore at the end inside the m_gast_concs object.

From 25a1e968e9bbed5b2890d0b39c5901e362361389 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 16 Jan 2026 01:49:28 +0000
Subject: [PATCH 07/10] Fix mask threshold in coarsening remapper rescale

Change mask threshold from epsilon (~1e-15) to 0.5 in rescale_masked_fields.
This prevents division by very small mask values which amplify numerical errors.

Root cause: When remapping fields with masks (like COSP sunlit_mask), the
remapped mask represents the weighted fraction of unmasked sources. Using
epsilon threshold meant even target cells with tiny mask values (e.g. 0.01
from 99% masked sources) would divide by that tiny value, amplifying any
stray fill values or numerical errors to produce huge spurious results
(10^31-10^32 range observed).

With 0.5 threshold: Only target cells where >50% of contributing sources
(by weight) are unmasked will be rescaled. Otherwise set to fill_value.
This ensures we only compute averages when we have meaningful data.

Co-authored-by: mahf708 <122953255+mahf708@users.noreply.github.com>
---
components/eamxx/src/share/remap/coarsening_remapper.cpp | 8 +++++++-
1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/components/eamxx/src/share/remap/coarsening_remapper.cpp b/components/eamxx/src/share/remap/coarsening_remapper.cpp
index ce5b2a0a6e9f..e5403d5bddd5 100644
--- a/components/eamxx/src/share/remap/coarsening_remapper.cpp
+++ b/components/eamxx/src/share/remap/coarsening_remapper.cpp
@@ -252,7 +252,13 @@ rescale_masked_fields (const Field& x, const Field& mask) const
const auto& layout = x.get_header().get_identifier().get_layout();
const int rank = layout.rank();
const int ncols = m_tgt_grid->get_num_local_dofs();
- const Real mask_threshold = std::numeric_limits::epsilon(); // TODO: Should we not hardcode the threshold for simply masking out the column.
+ // Use a threshold of 0.5 for the mask: only rescale if more than half of the
+ // contributing sources (by weight) are unmasked. This prevents division by very
+ // small mask values which can amplify numerical errors and cause large spurious values.
+ // For example, if a target cell has contributions from 9 masked and 1 unmasked source
+ // (each weight 0.1), the remapped mask would be 0.1. Dividing by such small values
+ // would amplify any numerical errors or stray fill values, producing huge results.
+ const Real mask_threshold = 0.5;

switch (rank) {
case 1:

From 7c4165b5d2ca8c05080f3339cb213789f97fc88e Mon Sep 17 00:00:00 2001
From: mahf708
Date: Thu, 15 Jan 2026 17:58:24 -0800
Subject: [PATCH 08/10] Fix critical inverted mask logic bug in coarsening
remapper

The rescale_masked_fields function had completely inverted logic for rank-2,
rank-3, and rank-4 fields with 2D masks:
- It was setting VALID (unmasked) values to fill_val
- Instead of setting INVALID (masked) values to fill_val

This caused large fill values (1.7e33) to leak into horizontally remapped COSP
outputs. The bug only manifested when there was a MIX of masked/unmasked cells
in the coarse grid target, explaining why some values were correct while others
showed huge fill values.

The fix:
1. Rename 'masked' -> 'valid' for clarity (mask > threshold means data is valid)
3. Then rescale VALID entries (overwrites those positions)
4. Remove unnecessary conditional check

Also updated radiation code TODOs to document why the dtau067/dtau105 deep_copy
is necessary: these fields must be zeroed before recomputation on radiation
timesteps, but persist between radiation calls for COSP to use.
---
.../rrtmgp/eamxx_rrtmgp_process_interface.cpp | 7 ++--
.../src/share/remap/coarsening_remapper.cpp | 32 +++++++++----------
2 files changed, 19 insertions(+), 20 deletions(-)

diff --git a/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp b/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
index 2d13bc4b09e3..064c110666e7 100644
--- a/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
+++ b/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
@@ -641,9 +641,10 @@ void RRTMGPRadiation::run_impl (const double dt) {
auto update_rad = scream::rrtmgp::radiation_do(m_rad_freq_in_steps, ts.get_num_steps());

if (update_rad) {
- // Init these to zero inside the update_rad conditional
- // TODO why do we need this anyway???
- // TODO figure out a way to get rid of this!!!
+ // Init optical depths to zero before recomputing.
+ // These fields are only updated on radiation timesteps (controlled by update_rad),
+ // and they must be zeroed before accumulation to avoid carrying over stale values.
+ // Between radiation calls, these fields persist in the FieldManager so COSP can use them.
Kokkos::deep_copy(d_dtau067,0.0);
Kokkos::deep_copy(d_dtau105,0.0);
// On each chunk, we internally "reset" the GasConcs object to subview the concs 3d array
diff --git a/components/eamxx/src/share/remap/coarsening_remapper.cpp b/components/eamxx/src/share/remap/coarsening_remapper.cpp
index e5403d5bddd5..cd9ce935d56c 100644
--- a/components/eamxx/src/share/remap/coarsening_remapper.cpp
+++ b/components/eamxx/src/share/remap/coarsening_remapper.cpp
@@ -315,11 +315,11 @@ rescale_masked_fields (const Field& x, const Field& mask) const
auto m_sub = ekat::subview(mask_2d,icol);
Kokkos::parallel_for(Kokkos::TeamVectorRange(team,dim1),
[&](const int j){
- auto masked = m_sub(j) > mask_threshold;
- if (masked.any()) {
- x_sub(j).set(masked,x_sub(j)/m_sub(j));
- }
- x_sub(j).set(!masked,fill_val);
+ auto valid = m_sub(j) > mask_threshold;
+ // First set invalid values to fill_val
+ x_sub(j).set(!valid,fill_val);
+ // Then rescale valid values (overwrites the valid entries that were just set)
+ x_sub(j).set(valid,x_sub(j)/m_sub(j));
});
}
});
@@ -371,12 +371,11 @@ rescale_masked_fields (const Field& x, const Field& mask) const
const int j = idx / dim2;
const int k = idx % dim2;
auto x_sub = ekat::subview(x_view,icol,j);
- auto masked = m_sub(k) > mask_threshold;
-
- if (masked.any()) {
- x_sub(k).set(masked,x_sub(k)/m_sub(k));
- }
- x_sub(k).set(!masked,fill_val);
+ auto valid = m_sub(k) > mask_threshold;
+ // First set invalid values to fill_val
+ x_sub(k).set(!valid,fill_val);
+ // Then rescale valid values (overwrites the valid entries that were just set)
+ x_sub(k).set(valid,x_sub(k)/m_sub(k));
});
}
});
@@ -432,12 +431,11 @@ rescale_masked_fields (const Field& x, const Field& mask) const
const int k = (idx / dim3) % dim2;
const int l = idx % dim3;
auto x_sub = ekat::subview(x_view,icol,j,k);
- auto masked = m_sub(l) > mask_threshold;
-
- if (masked.any()) {
- x_sub(l).set(masked,x_sub(l)/m_sub(l));
- }
- x_sub(l).set(!masked,fill_val);
+ auto valid = m_sub(l) > mask_threshold;
+ // First set invalid values to fill_val
+ x_sub(l).set(!valid,fill_val);
+ // Then rescale valid values (overwrites the valid entries that were just set)
+ x_sub(l).set(valid,x_sub(l)/m_sub(l));
});
}
});

From 36c310a857e0d8b833009d53c8310ddb98e0137c Mon Sep 17 00:00:00 2001
From: mahf708
Date: Thu, 15 Jan 2026 18:04:42 -0800
Subject: [PATCH 09/10] Fix critical fill_val leak in local_mat_vec

The actual root cause of large COSP values: when source fields contain
fill_val at masked locations, local_mat_vec was blindly multiplying them
by the mask and including them in the weighted sum. Even with mask=0,
floating point arithmetic can cause fill_val * tiny_mask to produce huge
spurious values.

The fix: Check if source value equals fill_val BEFORE multiplying by mask.
If it's fill_val, treat as zero contribution (skip it entirely).

Applied to all ranks (1-4) in local_mat_vec function.

Also reverted mask_threshold back to epsilon - the 0.5 threshold was too
aggressive and discarded valid data. The threshold of 0.5 was a band-aid
for the inverted logic bug we fixed earlier.
---
.../src/share/remap/coarsening_remapper.cpp | 67 +++++++++++++------
1 file changed, 46 insertions(+), 21 deletions(-)

diff --git a/components/eamxx/src/share/remap/coarsening_remapper.cpp b/components/eamxx/src/share/remap/coarsening_remapper.cpp
index cd9ce935d56c..78777839962d 100644
--- a/components/eamxx/src/share/remap/coarsening_remapper.cpp
+++ b/components/eamxx/src/share/remap/coarsening_remapper.cpp
@@ -252,13 +252,11 @@ rescale_masked_fields (const Field& x, const Field& mask) const
const auto& layout = x.get_header().get_identifier().get_layout();
const int rank = layout.rank();
const int ncols = m_tgt_grid->get_num_local_dofs();
- // Use a threshold of 0.5 for the mask: only rescale if more than half of the
- // contributing sources (by weight) are unmasked. This prevents division by very
- // small mask values which can amplify numerical errors and cause large spurious values.
- // For example, if a target cell has contributions from 9 masked and 1 unmasked source
- // (each weight 0.1), the remapped mask would be 0.1. Dividing by such small values
- // would amplify any numerical errors or stray fill values, producing huge results.
- const Real mask_threshold = 0.5;
+ // Use epsilon as threshold: rescale if there are any valid (unmasked) contributions.
+ // The mask value represents the weighted fraction of valid source cells.
+ // Values near zero indicate almost entirely masked sources -> set to fill_val.
+ // Any non-negligible mask value -> valid data to preserve by rescaling.
+ const Real mask_threshold = std::numeric_limits::epsilon();

switch (rank) {
case 1:
@@ -466,6 +464,7 @@ local_mat_vec (const Field& x, const Field& y, const Field& mask) const
// loop to zero out y before the mat-vec.
case 1:
{
+ constexpr auto fill_val = constants::fill_value;
// Unlike get_view, get_strided_view returns a LayoutStride view,
// therefore allowing the 1d field to be a subfield of a 2d field
// along the 2nd dimension.
@@ -476,15 +475,21 @@ local_mat_vec (const Field& x, const Field& y, const Field& mask) const
KOKKOS_LAMBDA(const int& row) {
const auto beg = row_offsets(row);
const auto end = row_offsets(row+1);
- y_view(row) = weights(beg)*x_view(col_lids(beg))*mask_view(col_lids(beg));
+ // Check if source value is fill_val - if so, don't include it (treat as 0)
+ auto x_val = x_view(col_lids(beg));
+ y_view(row) = (x_val == fill_val) ? 0.0 : weights(beg)*x_val*mask_view(col_lids(beg));
for (int icol=beg+1; icol;
auto x_view = x.get_view();
auto y_view = y.get_view< Pack**>();
view_1d mask_1d;
@@ -507,11 +512,17 @@ local_mat_vec (const Field& x, const Field& y, const Field& mask) const
const auto end = row_offsets(row+1);
Kokkos::parallel_for(Kokkos::TeamVectorRange(team,dim1),
[&](const int j){
- y_view(row,j) = weights(beg)*x_view(col_lids(beg),j) *
- (mask1d ? mask_1d (col_lids(beg)) : mask_2d(col_lids(beg),j));
+ // Check if source value is fill_val - if so, don't include it (treat as 0)
+ auto x_val = x_view(col_lids(beg),j);
+ auto is_fill = (x_val == fill_val);
+ y_view(row,j) = is_fill.any() ? Pack(0.0) :
+ weights(beg)*x_val*(mask1d ? mask_1d(col_lids(beg)) : mask_2d(col_lids(beg),j));
for (int icol=beg+1; icol;
auto x_view = x.get_view();
auto y_view = y.get_view< Pack***>();
// Note, the mask is still assumed to be defined on COLxLEV so still only 2D for case 3.
@@ -545,11 +557,17 @@ local_mat_vec (const Field& x, const Field& y, const Field& mask) const
[&](const int idx){
const int j = idx / dim2;
const int k = idx % dim2;
- y_view(row,j,k) = weights(beg)*x_view(col_lids(beg),j,k) *
- (mask1d ? mask_1d (col_lids(beg)) : mask_2d(col_lids(beg),k));
+ // Check if source value is fill_val - if so, don't include it (treat as 0)
+ auto x_val = x_view(col_lids(beg),j,k);
+ auto is_fill = (x_val == fill_val);
+ y_view(row,j,k) = is_fill.any() ? Pack(0.0) :
+ weights(beg)*x_val*(mask1d ? mask_1d(col_lids(beg)) : mask_2d(col_lids(beg),k));
for (int icol=beg+1; icol;
auto x_view = x.get_view();
auto y_view = y.get_view< Pack****>();
// Note, the mask is still assumed to be defined on COLxLEV so still only 2D for case 3.
@@ -585,11 +604,17 @@ local_mat_vec (const Field& x, const Field& y, const Field& mask) const
const int j = (idx / dim3) / dim2;
const int k = (idx / dim3) % dim2;
const int l = idx % dim3;
- y_view(row,j,k,l) = weights(beg)*x_view(col_lids(beg),j,k,l) *
- (mask1d ? mask_1d (col_lids(beg)) : mask_2d(col_lids(beg),l));
+ // Check if source value is fill_val - if so, don't include it (treat as 0)
+ auto x_val = x_view(col_lids(beg),j,k,l);
+ auto is_fill = (x_val == fill_val);
+ y_view(row,j,k,l) = is_fill.any() ? Pack(0.0) :
+ weights(beg)*x_val*(mask1d ? mask_1d(col_lids(beg)) : mask_2d(col_lids(beg),l));
for (int icol=beg+1; icol
Date: Sat, 7 Feb 2026 04:44:40 +0000
Subject: [PATCH 10/10] Address PR review feedback: fix COSP dims, use
create_geometry_data, prevent FPE, remove unneeded deep_copy

Key changes per review feedback:
- Fix COSP dimension values to match exact values in cosp_config.F90 (brhillman)
* tau bins/centers now match ISCCP standard (0.15, 0.80, 2.45, 6.5, 16.2, 41.5, 100.0)
* pressure bins now in Pa with correct values (90000, 74000, 62000, 50000, 37500, 24500, 9000)
* height bins include -99000 to 0 "no retrieval" bin per MISR standard
- Use create_geometry_data API instead of manual field creation (bartgol)
- Reinstate if checks in rescale_masked_fields to prevent FPE from division by zero (bartgol)
- Remove unnecessary deep_copy in radiation code (mahf708/bartgol/brhillman)

Co-authored-by: mahf708 <122953255+mahf708@users.noreply.github.com>
---
.../eamxx/src/physics/cosp/eamxx_cosp.cpp | 157 ++++++++++--------
.../rrtmgp/eamxx_rrtmgp_process_interface.cpp | 6 -
.../src/share/remap/coarsening_remapper.cpp | 21 ++-
3 files changed, 97 insertions(+), 87 deletions(-)

diff --git a/components/eamxx/src/physics/cosp/eamxx_cosp.cpp b/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
index 4fca050ea250..084582fe737d 100644
--- a/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
+++ b/components/eamxx/src/physics/cosp/eamxx_cosp.cpp
@@ -44,7 +44,6 @@ void Cosp::set_grids(const std::shared_ptr grids_manager)
auto micron = micro*m;
auto m2 = pow(m, 2);
auto s2 = pow(s, 2);
- auto hPa = hecto*Pa; // hectopascal (100 Pa)

m_grid = grids_manager->get_grid("physics");
const auto& grid_name = m_grid->name();
@@ -106,126 +105,140 @@ void Cosp::set_grids(const std::shared_ptr grids_manager)

// Add COSP dimension coordinate values and bounds as geometry data
// These are needed for e3sm_diags to produce COSP diagnostics
+ // Values match those hard-coded in components/eam/src/physics/cosp2/local/cosp_config.F90
using namespace ShortFieldTagsNames;

- // COSP tau bins (optical depth) - standard ISCCP/MODIS bins
- // Bin edges: 0.3, 1.3, 3.6, 9.4, 23, 60, 379
+ // COSP tau bins (optical depth) - standard ISCCP bins
+ // From COSP: tau_binEdges = reshape((/0.0, 0.3, 0.3, 1.3, 1.3, 3.6, 3.6, 9.4, 9.4, 23.0, 23.0, 60.0, 60.0, 100000.0/))
if (not m_grid->has_geometry_data("cosp_tau_bnds")) {
FieldLayout tau_bnds_layout({CMP,CMP},{m_num_tau,2},{"cosp_tau","nbnd"});
- Field tau_bnds(FieldIdentifier("cosp_tau_bnds", tau_bnds_layout, nondim, grid_name));
- tau_bnds.allocate_view();
+ auto& tau_bnds = m_grid->create_geometry_data("cosp_tau_bnds", tau_bnds_layout, nondim);
auto tau_bnds_h = tau_bnds.get_view();

- // Standard ISCCP/MODIS tau bin boundaries
- tau_bnds_h(0,0) = 0.3; tau_bnds_h(0,1) = 1.3;
- tau_bnds_h(1,0) = 1.3; tau_bnds_h(1,1) = 3.6;
- tau_bnds_h(2,0) = 3.6; tau_bnds_h(2,1) = 9.4;
- tau_bnds_h(3,0) = 9.4; tau_bnds_h(3,1) = 23.0;
- tau_bnds_h(4,0) = 23.0; tau_bnds_h(4,1) = 60.0;
- tau_bnds_h(5,0) = 60.0; tau_bnds_h(5,1) = 379.0;
- tau_bnds_h(6,0) = 379.0; tau_bnds_h(6,1) = 100000.0; // effectively infinity
+ // Standard ISCCP tau bin edges from COSP Fortran source
+ tau_bnds_h(0,0) = 0.0; tau_bnds_h(0,1) = 0.3;
+ tau_bnds_h(1,0) = 0.3; tau_bnds_h(1,1) = 1.3;
+ tau_bnds_h(2,0) = 1.3; tau_bnds_h(2,1) = 3.6;
+ tau_bnds_h(3,0) = 3.6; tau_bnds_h(3,1) = 9.4;
+ tau_bnds_h(4,0) = 9.4; tau_bnds_h(4,1) = 23.0;
+ tau_bnds_h(5,0) = 23.0; tau_bnds_h(5,1) = 60.0;
+ tau_bnds_h(6,0) = 60.0; tau_bnds_h(6,1) = 100000.0;

tau_bnds.sync_to_dev();
- m_grid->set_geometry_data(tau_bnds);
}

- // COSP tau centers
+ // COSP tau centers - use hard-coded values from COSP
+ // From COSP: tau_binCenters = (/0.15, 0.80, 2.45, 6.5, 16.2, 41.5, 100.0/)
if (not m_grid->has_geometry_data("cosp_tau")) {
- auto tau_bnds = m_grid->get_geometry_data("cosp_tau_bnds");
- auto tau_bnds_h = tau_bnds.get_view();
-
FieldLayout tau_layout({CMP},{m_num_tau},{"cosp_tau"});
- Field tau(FieldIdentifier("cosp_tau", tau_layout, nondim, grid_name));
- tau.allocate_view();
+ auto& tau = m_grid->create_geometry_data("cosp_tau", tau_layout, nondim);
auto tau_h = tau.get_view();

- // Geometric mean of bin edges (log-space midpoint)
- for (int i=0; iset_geometry_data(tau);
}

- // COSP pressure bins (cloud-top pressure in hPa) - standard ISCCP bins
- // Bin edges: 50, 180, 310, 440, 560, 680, 800, 1000 hPa
+ // COSP pressure bins (cloud-top pressure in Pa) - standard ISCCP bins
+ // From COSP: pres_binEdges = reshape((/100000, 80000, 80000, 68000, 68000, 56000, 56000, 44000, 44000, 31000, 31000, 18000, 18000, 0/))
if (not m_grid->has_geometry_data("cosp_prs_bnds")) {
FieldLayout prs_bnds_layout({CMP,CMP},{m_num_ctp,2},{"cosp_prs","nbnd"});
- Field prs_bnds(FieldIdentifier("cosp_prs_bnds", prs_bnds_layout, hPa, grid_name));
- prs_bnds.allocate_view();
+ auto& prs_bnds = m_grid->create_geometry_data("cosp_prs_bnds", prs_bnds_layout, Pa);
auto prs_bnds_h = prs_bnds.get_view();

- // Standard ISCCP pressure bin boundaries (hPa)
- prs_bnds_h(0,0) = 50.0; prs_bnds_h(0,1) = 180.0;
- prs_bnds_h(1,0) = 180.0; prs_bnds_h(1,1) = 310.0;
- prs_bnds_h(2,0) = 310.0; prs_bnds_h(2,1) = 440.0;
- prs_bnds_h(3,0) = 440.0; prs_bnds_h(3,1) = 560.0;
- prs_bnds_h(4,0) = 560.0; prs_bnds_h(4,1) = 680.0;
- prs_bnds_h(5,0) = 680.0; prs_bnds_h(5,1) = 800.0;
- prs_bnds_h(6,0) = 800.0; prs_bnds_h(6,1) = 1000.0;
+ // Standard ISCCP pressure bin edges from COSP Fortran source (in Pa)
+ prs_bnds_h(0,0) = 100000.0; prs_bnds_h(0,1) = 80000.0;
+ prs_bnds_h(1,0) = 80000.0; prs_bnds_h(1,1) = 68000.0;
+ prs_bnds_h(2,0) = 68000.0; prs_bnds_h(2,1) = 56000.0;
+ prs_bnds_h(3,0) = 56000.0; prs_bnds_h(3,1) = 44000.0;
+ prs_bnds_h(4,0) = 44000.0; prs_bnds_h(4,1) = 31000.0;
+ prs_bnds_h(5,0) = 31000.0; prs_bnds_h(5,1) = 18000.0;
+ prs_bnds_h(6,0) = 18000.0; prs_bnds_h(6,1) = 0.0;

prs_bnds.sync_to_dev();
- m_grid->set_geometry_data(prs_bnds);
}

- // COSP pressure centers
+ // COSP pressure centers - use hard-coded values from COSP
+ // From COSP: pres_binCenters = (/90000., 74000., 62000., 50000., 37500., 24500., 9000./)
if (not m_grid->has_geometry_data("cosp_prs")) {
- auto prs_bnds = m_grid->get_geometry_data("cosp_prs_bnds");
- auto prs_bnds_h = prs_bnds.get_view();
-
FieldLayout prs_layout({CMP},{m_num_ctp},{"cosp_prs"});
- Field prs(FieldIdentifier("cosp_prs", prs_layout, hPa, grid_name));
- prs.allocate_view();
+ auto& prs = m_grid->create_geometry_data("cosp_prs", prs_layout, Pa);
auto prs_h = prs.get_view();

- // Midpoint of pressure bins
- for (int i=0; iset_geometry_data(prs);
}

// COSP height bins (cloud-top height in m) - standard MISR bins
- // 16 bins from 0 to 20 km
+ // From COSP: hgt_binEdges = 1000*reshape((/-99.0, 0.0, 0.0, 0.5, 0.5, 1.0, ..., 17.0, 99.0/))
+ // Note: first bin (-99 to 0) is "no retrieval" bin for cases with no cloud top height retrieval
if (not m_grid->has_geometry_data("cosp_cth_bnds")) {
FieldLayout cth_bnds_layout({CMP,CMP},{m_num_cth,2},{"cosp_cth","nbnd"});
- Field cth_bnds(FieldIdentifier("cosp_cth_bnds", cth_bnds_layout, m, grid_name));
- cth_bnds.allocate_view();
+ auto& cth_bnds = m_grid->create_geometry_data("cosp_cth_bnds", cth_bnds_layout, m);
auto cth_bnds_h = cth_bnds.get_view();

- // Standard MISR height bin boundaries (m) - 16 bins of varying width
- // These follow the CFMIP/COSP standard MISR bins
- const Real cth_edges[17] = {0.0, 500.0, 1000.0, 1500.0, 2000.0, 2500.0, 3000.0, 4000.0,
- 5000.0, 7000.0, 9000.0, 11000.0, 13000.0, 15000.0, 17000.0, 19000.0, 21000.0};
- for (int i=0; iset_geometry_data(cth_bnds);
}

- // COSP height centers
+ // COSP height centers - use hard-coded values from COSP
+ // From COSP: hgt_binCenters = 1000*(/0., 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.5, 4.5, 6., 8., 10., 12., 14.5, 16., 18./)
if (not m_grid->has_geometry_data("cosp_cth")) {
- auto cth_bnds = m_grid->get_geometry_data("cosp_cth_bnds");
- auto cth_bnds_h = cth_bnds.get_view();
-
FieldLayout cth_layout({CMP},{m_num_cth},{"cosp_cth"});
- Field cth(FieldIdentifier("cosp_cth", cth_layout, m, grid_name));
- cth.allocate_view();
+ auto& cth = m_grid->create_geometry_data("cosp_cth", cth_layout, m);
auto cth_h = cth.get_view();

- // Midpoint of height bins
- for (int i=0; iset_geometry_data(cth);
}
}

diff --git a/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp b/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
index 064c110666e7..d6de68ba3e8c 100644
--- a/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
+++ b/components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
@@ -641,12 +641,6 @@ void RRTMGPRadiation::run_impl (const double dt) {
auto update_rad = scream::rrtmgp::radiation_do(m_rad_freq_in_steps, ts.get_num_steps());

if (update_rad) {
- // Init optical depths to zero before recomputing.
- // These fields are only updated on radiation timesteps (controlled by update_rad),
- // and they must be zeroed before accumulation to avoid carrying over stale values.
- // Between radiation calls, these fields persist in the FieldManager so COSP can use them.
- Kokkos::deep_copy(d_dtau067,0.0);
- Kokkos::deep_copy(d_dtau105,0.0);
// On each chunk, we internally "reset" the GasConcs object to subview the concs 3d array
// with the correct ncol dimension. So let's keep a copy of the original (ref-counted)
// array, to restore at the end inside the m_gast_concs object.
diff --git a/components/eamxx/src/share/remap/coarsening_remapper.cpp b/components/eamxx/src/share/remap/coarsening_remapper.cpp
index 78777839962d..1e200f137c69 100644
--- a/components/eamxx/src/share/remap/coarsening_remapper.cpp
+++ b/components/eamxx/src/share/remap/coarsening_remapper.cpp
@@ -314,10 +314,11 @@ rescale_masked_fields (const Field& x, const Field& mask) const
Kokkos::parallel_for(Kokkos::TeamVectorRange(team,dim1),
[&](const int j){
auto valid = m_sub(j) > mask_threshold;
- // First set invalid values to fill_val
+ // Avoid division by zero FPE: only divide if there are valid entries
+ if (valid.any()) {
+ x_sub(j).set(valid,x_sub(j)/m_sub(j));
+ }
x_sub(j).set(!valid,fill_val);
- // Then rescale valid values (overwrites the valid entries that were just set)
- x_sub(j).set(valid,x_sub(j)/m_sub(j));
});
}
});
@@ -370,10 +371,11 @@ rescale_masked_fields (const Field& x, const Field& mask) const
const int k = idx % dim2;
auto x_sub = ekat::subview(x_view,icol,j);
auto valid = m_sub(k) > mask_threshold;
- // First set invalid values to fill_val
+ // Avoid division by zero FPE: only divide if there are valid entries
+ if (valid.any()) {
+ x_sub(k).set(valid,x_sub(k)/m_sub(k));
+ }
x_sub(k).set(!valid,fill_val);
- // Then rescale valid values (overwrites the valid entries that were just set)
- x_sub(k).set(valid,x_sub(k)/m_sub(k));
});
}
});
@@ -430,10 +432,11 @@ rescale_masked_fields (const Field& x, const Field& mask) const
const int l = idx % dim3;
auto x_sub = ekat::subview(x_view,icol,j,k);
auto valid = m_sub(l) > mask_threshold;
- // First set invalid values to fill_val
+ // Avoid division by zero FPE: only divide if there are valid entries
+ if (valid.any()) {
+ x_sub(l).set(valid,x_sub(l)/m_sub(l));
+ }
x_sub(l).set(!valid,fill_val);
- // Then rescale valid values (overwrites the valid entries that were just set)
- x_sub(l).set(valid,x_sub(l)/m_sub(l));
});
}
});
```