Skip to content

EAMxx: fix COSP fill value leak and add dimension coordinates for cosp diagnostics #7992

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Copilot wants to merge 9 commits into
base: master
Choose a base branch
from
+241 −34
Open

EAMxx: fix COSP fill value leak and add dimension coordinates for cosp diagnostics #7992

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
bd476dc
Initial plan
Copilot Jan 15, 2026
109f33e
Add mask_value to COSP outputs and propagate in remappers
Copilot Jan 15, 2026
745d4b2
Fix critical bug in coarsening_remapper rescale for masked columns
Copilot Jan 15, 2026
ace07b8
Add mask_value propagation to IO fields - critical missing piece
Copilot Jan 15, 2026
f5672ec
Fix build error and add COSP dimension coordinates
Copilot Jan 15, 2026
beaefe5
properly persist rad cld opt depth
mahf708 Jan 16, 2026
25a1e96
Fix mask threshold in coarsening remapper rescale
Copilot Jan 16, 2026
7c4165b
Fix critical inverted mask logic bug in coarsening remapper
mahf708 Jan 16, 2026
36c310a
Fix critical fill_val leak in local_mat_vec
mahf708 Jan 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 126 additions & 0 deletions components/eamxx/src/physics/cosp/eamxx_cosp.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ void Cosp::set_grids(const std::shared_ptr<const GridsManager> 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();
Expand Down Expand Up @@ -102,6 +103,130 @@ void Cosp::set_grids(const std::shared_ptr<const GridsManager> 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
Copy link
Contributor

@mahf708 mahf708 Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no idea where the bot came up with these values --- good job bot for either finding the real ones or making up convincing looking ones!

// 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));
Copy link
Contributor

@bartgol bartgol Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could do

auto& tau_bnds = m_grid->create_geometry_data("cosp_tau_bnds", tau_bnds_layout,nondim);

and avoid to have to create/alloc the field manually, and call set_geometry_data later, saving some lines and avoid using lower-level interfaces.

Copy link
Contributor

@bartgol bartgol Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similarly for all the other geo data.

tau_bnds.allocate_view();
auto tau_bnds_h = tau_bnds.get_view<Real**,Host>();

// 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;
Copy link
Contributor

@brhillman brhillman Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Close but not quite right. These are defined in the COSP source (components/eam/src/physics/cosp2/external/src/cosp_config.F90):

    ! Optical depth bin axis
    integer,parameter :: &
         ntau=7  
    real(wp),parameter,dimension(ntau+1) :: &
       tau_binBounds = (/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60., 10000./)
    real(wp),parameter,dimension(ntau) :: &
         tau_binCenters = (/0.15, 0.80, 2.45, 6.5, 16.2, 41.5, 100.0/)
    real(wp),parameter,dimension(2,ntau) :: &
         tau_binEdges = reshape(source=(/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/),   &
                                         shape=(/2,ntau/))

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<const Real**,Host>();

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<Real*,Host>();

// Geometric mean of bin edges (log-space midpoint)
Copy link
Contributor

@brhillman brhillman Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would be better to do here, if it were consistent with what COSP does under the hood, but I think the midpoints are generally expected to have the sort of standard ISCCP-defined midpoints, which are specified explicitly in the above code block I copied. Also note that, unfortunately, the MODIS simulator uses a different set (although it's the same thing just with the first index dropped):

    ! Optical depth bin axis
    integer,parameter :: &
         ntau=7  
    real(wp),parameter,dimension(ntau+1) :: &
       tau_binBounds = (/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60., 10000./)
    real(wp),parameter,dimension(ntau) :: &
         tau_binCenters = (/0.15, 0.80, 2.45, 6.5, 16.2, 41.5, 100.0/)
    real(wp),parameter,dimension(2,ntau) :: &
         tau_binEdges = reshape(source=(/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/),   &
                                         shape=(/2,ntau/)) 

    ! Optical depth bin axes (ONLY USED BY MODIS SIMULATOR IN v1.4)
    integer :: l,k
    integer,parameter :: &
         ntauV1p4 = 6
    real(wp),parameter,dimension(ntauV1p4+1) :: &
         tau_binBoundsV1p4 = (/0.3, 1.3, 3.6, 9.4, 23., 60., 10000./)
    real(wp),parameter,dimension(2,ntauV1p4) :: &
         tau_binEdgesV1p4 = reshape(source =(/tau_binBoundsV1p4(1),((tau_binBoundsV1p4(k),l=1,2),   &
                                             k=2,ntauV1p4),100000._wp/),shape = (/2,ntauV1p4/)) 
    real(wp),parameter,dimension(ntauV1p4) :: &
         tau_binCentersV1p4 = (tau_binEdgesV1p4(1,:)+tau_binEdgesV1p4(2,:))/2._wp

for (int i=0; i<m_num_tau; ++i) {
tau_h(i) = std::sqrt(tau_bnds_h(i,0) * tau_bnds_h(i,1));
}

tau.sync_to_dev();
m_grid->set_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<Real**,Host>();

// Standard ISCCP pressure bin boundaries (hPa)
prs_bnds_h(0,0) = 50.0; prs_bnds_h(0,1) = 180.0;
Copy link
Contributor

@brhillman brhillman Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

close but again differs slightly from what is hard-coded in cosp:

    ! Cloud-top height pressure bin axis
    integer,parameter :: &
         npres = 7     
    real(wp),parameter,dimension(npres+1) :: &
         pres_binBounds = (/0., 180., 310., 440., 560., 680., 800., 10000./)
    real(wp),parameter,dimension(npres) :: &
         pres_binCenters = (/90000., 74000., 62000., 50000., 37500., 24500., 9000./)   
    real(wp),parameter,dimension(2,npres) :: &
         pres_binEdges = reshape(source=(/100000.0, 80000.0, 80000.0, 68000.0, 68000.0,    &
                                           56000.0, 56000.0, 44000.0, 44000.0, 31000.0,    &
                                           31000.0, 18000.0, 18000.0,     0.0/),           &
                                           shape=(/2,npres/))

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
Copy link
Contributor

@brhillman brhillman Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that pressure centers are also defined explicitly (hard-coded in COSP)

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<const Real**,Host>();

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<Real*,Host>();

// Midpoint of pressure bins
Copy link
Contributor

@brhillman brhillman Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think these will come out the same, but they are defined explicitly in the COSP code, we may as well use those here.

for (int i=0; i<m_num_ctp; ++i) {
prs_h(i) = (prs_bnds_h(i,0) + prs_bnds_h(i,1)) / 2.0;
}

prs.sync_to_dev();
m_grid->set_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<Real**,Host>();

// Standard MISR height bin boundaries (m) - 16 bins of varying width
Copy link
Contributor

@brhillman brhillman Jan 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again there are defined MISR CTH bins, and they differ from what is here:

    real(wp),parameter,dimension(nhgt+1) :: &
         hgt_binBounds = (/-.99,0.,0.5,1.,1.5,2.,2.5,3.,4.,5.,7.,9.,11.,13.,15.,17.,99./)
    real(wp),parameter,dimension(nhgt) :: &
         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./)
    real(wp),parameter,dimension(2,nhgt) :: &
         hgt_binEdges = 1000.0*reshape(source=(/-99.0, 0.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.5,  &
                                                  1.5, 2.0, 2.0, 2.5, 2.5, 3.0, 3.0, 4.0,  &
                                                  4.0, 5.0, 5.0, 7.0, 7.0, 9.0, 9.0,11.0,  &
                                                  11.0,13.0,13.0,15.0,15.0,17.0,17.0,99.0/),&
                                                  shape=(/2,nhgt/))

The first, -99 - 0 bin, is a "no retrieval" bin, indicating cases expected to have no cloud top height retrieval (usually because optical depth is too low I think).

// 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; i<m_num_cth; ++i) {
cth_bnds_h(i,0) = cth_edges[i];
cth_bnds_h(i,1) = cth_edges[i+1];
}

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

// COSP height centers
Copy link
Contributor

@brhillman brhillman Jan 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Centers hard-coded in COSP code. See above comments.

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<const Real**,Host>();

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<Real*,Host>();

// Midpoint of height bins
for (int i=0; i<m_num_cth; ++i) {
cth_h(i) = (cth_bnds_h(i,0) + cth_bnds_h(i,1)) / 2.0;
}

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

// =========================================================================================
Expand All @@ -115,6 +240,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<Real>);
get_field_out(field_name).get_header().set_may_be_filled(true);
}
}
Expand Down
8 changes: 6 additions & 2 deletions components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -614,8 +614,6 @@ void RRTMGPRadiation::run_impl (const double dt) {
auto d_dtau105 = get_field_out("dtau105").get_view<Real**>();
auto d_sunlit = get_field_out("sunlit_mask").get_view<Real*>();

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<Real *>();
auto d_p_mid_at_cldtop = get_field_out("p_mid_at_cldtop").get_view<Real *>();
Expand Down Expand Up @@ -643,6 +641,12 @@ 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.
Copy link
Contributor

@mahf708 mahf708 Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

meh, i don't believe this, but the bot kept me it is the case

Copy link
Contributor

@bartgol bartgol Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we don't accum into the fields, and the "fix" is not needed, let's remove it before we merge.

Copy link
Contributor

@brhillman brhillman Jan 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I'm not sure this is right. Fundamentally, we should NOT be calling COSP between radiation calls anyways, and if we do, it looks like zero-ing them out just gives us nonsense anyways.

// 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.
Expand Down
2 changes: 2 additions & 0 deletions components/eamxx/src/share/diagnostics/aodvis.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,8 @@ set_grids(const std::shared_ptr<const GridsManager> 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<Real>);
m_diagnostic_output.get_header().set_may_be_filled(true);
}

void AODVis::compute_diagnostic_impl() {
Expand Down
1 change: 1 addition & 0 deletions components/eamxx/src/share/diagnostics/field_at_pressure_level.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real>);
Copy link
Contributor

@bartgol bartgol Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The extra data "mask_value" is no longer used in IO. We only guard against fill value, no other arbitrary "special" value is considered.

Copy link
Contributor

@bartgol bartgol Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And we need to clean up the mask_xyz metadata. We have some "mask_data" and some "mask_field", both referring to the same thing. Not in this PR though.

Copy link
Contributor

@bartgol bartgol Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Uh, I do see the extra data "mask_value" being used in some diags though. I will bump up the priority of cleaning up the mask_xyz ambiguities.

m_diagnostic_output.get_header().set_may_be_filled(true);

using stratts_t = std::map<std::string,std::string>;
Expand Down
6 changes: 6 additions & 0 deletions components/eamxx/src/share/io/scorpio_output.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")) {
Copy link
Contributor

@bartgol bartgol Jan 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It may be pointless, since IO does not use mask_value.

auto src_mask_value = src.get_header().get_extra_data<scream::Real>("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?
Expand Down
Loading
Loading