Merge branch 'mahf708/eamxx/cosp-expected-behavior' into next (PR #7652)

persist cosp-computed fields when cosp in not run,
properly handle fill-value treatment for night indices,
and update docs to reflect new behavior.

[NBFB]

Author: bartgol

components/eamxx/cime_config/testdefs/testmods_dirs/eamxx/mvkpert/yaml_outs/monthly_average_coarse.yaml

@@ -84,7 +84,6 @@ fields:
     - isccp_ctptau
     - modis_ctptau
     - misr_cthtau
-    - cosp_sunlit
     - isccp_cldtot
 
 output_control:

components/eamxx/cime_config/testdefs/testmods_dirs/eamxx/prod/yaml_outs/eamxx_output.decadal.1hourlyAVG_native.yaml

@@ -11,7 +11,6 @@ fields:
     - isccp_ctptau
     - modis_ctptau
     - misr_cthtau
-    - cosp_sunlit
     - isccp_cldtot
 output_control:
   frequency: 1

components/eamxx/cime_config/testdefs/testmods_dirs/eamxx/prod/yaml_outs/eamxx_output.decadal.hourlyAVG_coarse.yaml

@@ -84,7 +84,6 @@ fields:
     - isccp_ctptau
     - modis_ctptau
     - misr_cthtau
-    - cosp_sunlit
     - isccp_cldtot
 
 output_control:

components/eamxx/docs/user/cosp.md

@@ -55,7 +55,6 @@ fields:
     - isccp_ctptau
     - modis_ctptau
     - misr_cthtau
-    - cosp_sunlit
 max_snapshots_per_file: 1
 filename_prefix: eamxx
 output_control:
@@ -76,19 +75,3 @@ The following output fields are available:
       - MODIS-simulated cloud top pressure/optical depth joint histogram
 - `misr_cthtau`
       - MISR-simulated cloud top height/optical depth joint histogram
-- `cosp_sunlit`
-      - sunlit flag aggregated at COSP frequency for renormalizing daytime averages
-
-ISCCP, MODIS, and MISR outputs are valid only for daytime/sunlit columns
-(to be consistent with available satellite retrievals).
-In order to aggregate only daytime columns in time averages, these outputs are
-multiplied by the sunlit flag (0 or 1) at each COSP calculation time.
-Time averages of these quantities are then aggregated, along with the COSP
-sunlit flag each time COSP is called.
-In order to back out the daytime-only time averages from the outputs,
-one needs to divide the output fields by `cosp_sunlit`.
-E.g.,
-
-```shell
-isccp_ctptau = mean(isccp_ctptau) / mean(cosp_sunlit)
-```

components/eamxx/src/diagnostics/aodvis.cpp

@@ -31,7 +31,7 @@ set_grids(const std::shared_ptr<const GridsManager> grids_manager)
 
   // The fields required for this diagnostic to be computed
   add_field<Required>("aero_tau_sw", vector3d, nondim, grid_name);
-  add_field<Required>("sunlit",      scalar2d, nondim, grid_name);
+  add_field<Required>("sunlit_mask", scalar2d, nondim, grid_name);
 
   // Construct and allocate the aodvis field
   FieldIdentifier fid(name(), scalar2d, nondim, grid_name);
@@ -41,21 +41,7 @@ set_grids(const std::shared_ptr<const GridsManager> grids_manager)
 }
 
 void AODVis::initialize_impl(const RunType /*run_type*/) {
-  // we use initialize_impl to primarily deal with the mask
-  using namespace ekat::units;
-  using namespace ShortFieldTagsNames;
-
-  auto nondim = ekat::units::Units::nondimensional();
-  const auto &grid_name =
-      m_diagnostic_output.get_header().get_identifier().get_grid_name();
-
-  std::string mask_name = name() + " mask";
-  FieldLayout mask_layout({COL}, {m_ncols});
-  FieldIdentifier mask_fid(mask_name, mask_layout, nondim, grid_name);
-  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_field", get_field_in("sunlit_mask"));
 }
 
 void AODVis::compute_diagnostic_impl() {
@@ -67,13 +53,10 @@ void AODVis::compute_diagnostic_impl() {
   constexpr auto fill_value = constants::fill_value<Real>;
 
   const auto aod     = m_diagnostic_output.get_view<Real *>();
-  const auto mask    = m_diagnostic_output.get_header()
-                        .get_extra_data<Field>("mask_field")
-                        .get_view<Real *>();
   const auto tau_vis = get_field_in("aero_tau_sw")
                            .subfield(1, m_vis_bnd)
                            .get_view<const Real **>();
-  const auto sunlit = get_field_in("sunlit").get_view<const Real *>();
+  const auto sunlit = get_field_in("sunlit_mask").get_view<const Real *>();
 
   const auto num_levs = m_nlevs;
   const auto policy   = TPF::get_default_team_policy(m_ncols, m_nlevs);
@@ -82,11 +65,9 @@ void AODVis::compute_diagnostic_impl() {
         const int icol = team.league_rank();
         if(sunlit(icol) == 0.0) {
           aod(icol) = fill_value;
-          Kokkos::single(Kokkos::PerTeam(team), [&] { mask(icol) = 0; });
         } else {
           auto tau_icol = ekat::subview(tau_vis, icol);
           aod(icol)     = RU::view_reduction(team, 0, num_levs, tau_icol);
-          Kokkos::single(Kokkos::PerTeam(team), [&] { mask(icol) = 1; });
         }
       });
 }

components/eamxx/src/diagnostics/tests/aodvis_test.cpp

@@ -62,7 +62,7 @@ TEST_CASE("aodvis") {
   tau.get_header().get_tracking().update_time_stamp(t0);
   // Input (randomized) sunlit
   FieldLayout scalar2d_layout = grid->get_2d_scalar_layout();
-  FieldIdentifier sunlit_fid("sunlit", scalar2d_layout, nondim, grid->name());
+  FieldIdentifier sunlit_fid("sunlit_mask", scalar2d_layout, nondim, grid->name());
   Field sunlit(sunlit_fid);
   sunlit.allocate_view();
   sunlit.get_header().get_tracking().update_time_stamp(t0);
@@ -112,6 +112,7 @@ TEST_CASE("aodvis") {
 
     const auto tau_h  = tau.get_view<const Real ***, Host>();
     const auto aod_hf = diag->get_diagnostic();
+    const auto aod_mask = aod_hf.get_header().get_extra_data<Field>("mask_field");
 
     Field aod_tf = diag->get_diagnostic().clone();
     auto aod_t = aod_tf.get_view<Real *, Host>();
@@ -133,6 +134,8 @@ TEST_CASE("aodvis") {
     if(SCREAM_BFB_TESTING) {
       REQUIRE(views_are_equal(aod_hf, aod_tf));
     }
+    // Ensure the masks are identical (should be pointing exactly to each other)
+    REQUIRE(views_are_equal(aod_mask, sunlit));
   }
 }
 

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

@@ -1,6 +1,7 @@
 #include "eamxx_cosp.hpp"
 #include "cosp_functions.hpp"
 #include "physics/share/physics_constants.hpp"
+#include "share/util/eamxx_universal_constants.hpp"
 #include "share/property_checks/field_within_interval_check.hpp"
 #include "share/field/field_utils.hpp"
 
@@ -65,7 +66,7 @@ void Cosp::set_grids(const std::shared_ptr<const GridsManager> grids_manager)
   add_field<Required>("surf_radiative_T", scalar2d    , K,      grid_name);
   //add_field<Required>("surfelev",    scalar2d    , m,      grid_name);
   //add_field<Required>("landmask",    scalar2d    , nondim, grid_name);
-  add_field<Required>("sunlit",           scalar2d    , nondim, grid_name);
+  add_field<Required>("sunlit_mask",       scalar2d    , nondim, grid_name);
   add_field<Required>("p_mid",             scalar3d_mid, Pa,     grid_name);
   add_field<Required>("p_int",             scalar3d_int, Pa,     grid_name);
   //add_field<Required>("height_mid",  scalar3d_mid, m,      grid_name);
@@ -88,11 +89,11 @@ void Cosp::set_grids(const std::shared_ptr<const GridsManager> grids_manager)
   add_field<Required>("eff_radius_qc",     scalar3d_mid, micron,      grid_name);
   add_field<Required>("eff_radius_qi",     scalar3d_mid, micron,      grid_name);
   // Set of fields used strictly as output
+  // NOTE we set their corresponding masks in init impl
   add_field<Computed>("isccp_cldtot", scalar2d, percent, grid_name);
   add_field<Computed>("isccp_ctptau", scalar4d_ctptau, percent, grid_name, 1);
   add_field<Computed>("modis_ctptau", scalar4d_ctptau, percent, grid_name, 1);
   add_field<Computed>("misr_cthtau", scalar4d_cthtau, percent, grid_name, 1);
-  add_field<Computed>("cosp_sunlit", scalar2d, nondim, grid_name);
 
   // We can allocate these now
   m_z_mid = Field(FieldIdentifier("z_mid",scalar3d_mid,m,grid_name));
@@ -107,15 +108,11 @@ void Cosp::initialize_impl (const RunType /* run_type */)
   // Set property checks for fields in this process
   CospFunc::initialize(m_num_cols, m_num_subcols, m_num_levs);
 
-
-  // Add note to output files about processing ISCCP fields that are only valid during
-  // daytime. This can go away once I/O can handle masked time averages.
-  using stratts_t = std::map<std::string,std::string>;
+  // Set the mask field for each of the cosp computed fields
   std::list<std::string> vnames = {"isccp_cldtot", "isccp_ctptau", "modis_ctptau", "misr_cthtau"};
-  for (const auto field_name : {"isccp_cldtot", "isccp_ctptau", "modis_ctptau", "misr_cthtau"}) {
-      auto& f = get_field_out(field_name);
-      auto& atts = f.get_header().get_extra_data<stratts_t>("io: string attributes");
-      atts["note"] = "Night values are zero; divide by cosp_sunlit to get daytime mean";
+  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"));
   }
 }
 
@@ -150,7 +147,7 @@ void Cosp::run_impl (const double dt)
     get_field_in("qv").sync_to_host();
     get_field_in("qc").sync_to_host();
     get_field_in("qi").sync_to_host();
-    get_field_in("sunlit").sync_to_host();
+    get_field_in("sunlit_mask").sync_to_host();
     get_field_in("surf_radiative_T").sync_to_host();
     get_field_in("T_mid").sync_to_host();
     get_field_in("p_mid").sync_to_host();
@@ -211,7 +208,7 @@ void Cosp::run_impl (const double dt)
     const auto p_mid_h   = get_field_in("p_mid").get_view<const Real**,Host>();
     const auto qc_h      = get_field_in("qc").get_view<const Real**, Host>();
     const auto qi_h      = get_field_in("qi").get_view<const Real**, Host>();
-    const auto sunlit_h  = get_field_in("sunlit").get_view<const Real*, Host>();
+    const auto sunlit_h  = get_field_in("sunlit_mask").get_view<const Real*, Host>();
     const auto skt_h     = get_field_in("surf_radiative_T").get_view<const Real*, Host>();
     const auto p_int_h   = get_field_in("p_int").get_view<const Real**, Host>();
     const auto cldfrac_h = get_field_in("cldfrac_rad").get_view<const Real**, Host>();
@@ -223,29 +220,28 @@ void Cosp::run_impl (const double dt)
     auto isccp_cldtot_h = get_field_out("isccp_cldtot").get_view<Real*, Host>();
     auto isccp_ctptau_h = get_field_out("isccp_ctptau").get_view<Real***, Host>();
     auto modis_ctptau_h = get_field_out("modis_ctptau").get_view<Real***, Host>();
-    auto misr_cthtau_h  = get_field_out("misr_cthtau").get_view<Real***, Host>();
-    auto cosp_sunlit_h  = get_field_out("cosp_sunlit").get_view<Real*, Host>();  // Copy of sunlit flag with COSP frequency for proper averaging
+    auto misr_cthtau_h  = get_field_out("misr_cthtau"). get_view<Real***, Host>();
 
     Real emsfc_lw = 0.99;
-    Kokkos::deep_copy(cosp_sunlit_h, sunlit_h);
     CospFunc::main(
             m_num_cols, m_num_subcols, m_num_levs, m_num_tau, m_num_ctp, m_num_cth, emsfc_lw,
             sunlit_h, skt_h, T_mid_h, p_mid_h, p_int_h, z_mid_h, qv_h, qc_h, qi_h,
             cldfrac_h, reff_qc_h, reff_qi_h, dtau067_h, dtau105_h,
             isccp_cldtot_h, isccp_ctptau_h, modis_ctptau_h, misr_cthtau_h
     );
-    // Remask night values to ZERO since our I/O does not know how to handle masked/missing values
-    // in temporal averages; this is all host data, so we can just use host loops like its the 1980s
+    // Mask night values
+    constexpr auto fill_value = constants::fill_value<Real>;
     for (int i = 0; i < m_num_cols; i++) {
       if (sunlit_h(i) == 0) {
-        isccp_cldtot_h(i) = 0;
+        // if night, set to fill val
+        isccp_cldtot_h(i) = fill_value;
         for (int j = 0; j < m_num_tau; j++) {
           for (int k = 0; k < m_num_ctp; k++) {
-            isccp_ctptau_h(i,j,k) = 0;
-            modis_ctptau_h(i,j,k) = 0;
+            isccp_ctptau_h(i,j,k) = fill_value;
+            modis_ctptau_h(i,j,k) = fill_value;
           }
           for (int k = 0; k < m_num_cth; k++) {
-            misr_cthtau_h (i,j,k) = 0;
+            misr_cthtau_h (i,j,k) = fill_value;
           }
         }
       }
@@ -256,22 +252,6 @@ void Cosp::run_impl (const double dt)
     get_field_out("isccp_ctptau").sync_to_dev();
     get_field_out("modis_ctptau").sync_to_dev();
     get_field_out("misr_cthtau").sync_to_dev();
-    get_field_out("cosp_sunlit").sync_to_dev();
-  } else {
-    // If not updating COSP statistics, set these to ZERO; this essentially weights
-    // the ISCCP cloud properties by the sunlit mask. What will be output for time-averages
-    // then is the time-average mask-weighted statistics; to get true averages, we need to
-    // divide by the time-average of the mask. I.e., if M is the sunlit mask, and X is the ISCCP
-    // statistic, then
-    //
-    //     avg(X) = sum(M * X) / sum(M) = (sum(M * X)/N) / (sum(M)/N) = avg(M * X) / avg(M)
-    //
-    // TODO: mask this when/if the AD ever supports masked averages
-    get_field_out("isccp_cldtot").deep_copy(0);
-    get_field_out("isccp_ctptau").deep_copy(0);
-    get_field_out("modis_ctptau").deep_copy(0);
-    get_field_out("misr_cthtau").deep_copy(0);
-    get_field_out("cosp_sunlit").deep_copy(0);
   }
 }
 

components/eamxx/src/physics/rrtmgp/eamxx_rrtmgp_process_interface.cpp

@@ -209,7 +209,7 @@ void RRTMGPRadiation::set_grids(const std::shared_ptr<const GridsManager> grids_
   // 0.67 micron and 10.5 micron optical depth (needed for COSP)
   add_field<Computed>("dtau067"       , scalar3d_mid, nondim, grid_name);
   add_field<Computed>("dtau105"       , scalar3d_mid, nondim, grid_name);
-  add_field<Computed>("sunlit"        , scalar2d    , nondim, grid_name);
+  add_field<Computed>("sunlit_mask"   , scalar2d    , nondim, grid_name);
   add_field<Computed>("cldfrac_rad"   , scalar3d_mid, nondim, grid_name);
   // Cloud-top diagnostics following AeroCom recommendation
   add_field<Computed>("T_mid_at_cldtop", scalar2d, K, grid_name);
@@ -612,7 +612,7 @@ void RRTMGPRadiation::run_impl (const double dt) {
   // Outputs for COSP
   auto d_dtau067 = get_field_out("dtau067").get_view<Real**>();
   auto d_dtau105 = get_field_out("dtau105").get_view<Real**>();
-  auto d_sunlit = get_field_out("sunlit").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);

components/eamxx/src/share/grid/remap/coarsening_remapper.cpp

@@ -432,9 +432,9 @@ local_mat_vec (const Field& x, const Field& y, const Field& mask) const
       // 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.
-      auto x_view = x.get_strided_view<const Real*>();
-      auto y_view = y.get_strided_view<      Real*>();
-      auto mask_view = mask.get_strided_view<Real*>();
+      auto x_view    =    x.get_strided_view<const Real*>();
+      auto y_view    =    y.get_strided_view<      Real*>();
+      auto mask_view = mask.get_strided_view<const Real*>();
       Kokkos::parallel_for(RangePolicy(0,nrows),
                            KOKKOS_LAMBDA(const int& row) {
         const auto beg = row_offsets(row);

components/eamxx/tests/single-process/cosp/input.yaml

@@ -34,7 +34,7 @@ initial_conditions:
   cldfrac_rad: 0.5
   eff_radius_qc: 10.0
   eff_radius_qi: 10.0
-  sunlit: 1.0
+  sunlit_mask: 1.0
   surf_radiative_T: 288.0
   pseudo_density: 1.0