Merge Pull Request #2348 from E3SM-Project/scream/aarondonahue/force_sfc_export_from_file_v2

Automatically Merged using E3SM Pull Request AutoTester
PR Title: Add capability to set surface export values via data from file(s)
PR Author: AaronDonahue

Author: E3SM-Bot

components/eamxx/cime_config/namelist_defaults_scream.xml

@@ -170,6 +170,10 @@ be lost if SCREAM_HACK_XML is not enabled.
         <fields type="array(string)">NONE</fields>
         <values type="array(real)">0</values>
       </prescribed_constants>
+      <prescribed_from_file>
+        <fields type="array(string)">NONE</fields>
+        <files type="array(string)">NONE</files>
+      </prescribed_from_file>
     </sc_export>
 
     <!-- Homme dynamics -->

components/eamxx/src/control/atmosphere_driver.hpp

@@ -136,7 +136,6 @@ class AtmosphereDriver
   // NOTE: if already finalized, this is a no-op
   void finalize ();
 
-  field_mgr_ptr get_ref_grid_field_mgr () const;
   field_mgr_ptr get_field_mgr (const std::string& grid_name) const;
 
   // Get atmosphere time stamp

components/eamxx/src/control/atmosphere_surface_coupling_exporter.cpp

@@ -171,6 +171,7 @@ void SurfaceCouplingExporter::initialize_impl (const RunType /* run_type */)
   for (int i=0; i<m_num_scream_exports; ++i) {
 
     std::string fname = m_export_field_names[i];
+    m_export_field_names_vector.push_back(fname);
     EKAT_REQUIRE_MSG(has_helper_field(fname),"Error! Attempting to export "+fname+
                    " which has not been added as a helper field.\n");
     auto& field = m_helper_fields.at(fname);
@@ -209,33 +210,78 @@ void SurfaceCouplingExporter::initialize_impl (const RunType /* run_type */)
   m_export_source_h = Kokkos::create_mirror_view(m_export_source);
   Kokkos::deep_copy(m_export_source_h,FROM_MODEL);  // The default is that all export variables will be derived from the EAMxx state.
   m_num_from_model_exports = m_num_scream_exports;
- 
+
+  // Parse all fields that will be set by interpolating data from file(s):
+  if (m_params.isSublist("prescribed_from_file")) {
+    // Retrieve the parameters for prescribed from file
+    auto export_from_file_params = m_params.sublist("prescribed_from_file");
+    EKAT_REQUIRE_MSG(export_from_file_params.isParameter("fields"),"Error! surface_coupling_exporter::init - prescribed_from_file does not have 'fields' parameter.");
+    auto export_from_file_fields = export_from_file_params.get<vos_type>("fields");
+
+    EKAT_REQUIRE_MSG(export_from_file_params.isParameter("files"),"Error! surface_coupling_exporter::init - prescribed_from_file does not have 'files' parameter.");
+    auto export_from_file_names = export_from_file_params.get<vos_type>("files");
+
+    bool are_fields_present = (std::find(export_from_file_fields.begin(),export_from_file_fields.end(),"NONE") == export_from_file_fields.end()) and (export_from_file_fields.size() > 0);
+    bool are_files_present  = (std::find(export_from_file_names.begin(),export_from_file_names.end(),"NONE") == export_from_file_names.end()) and (export_from_file_names.size() > 0);
+    if (are_fields_present) {
+      EKAT_REQUIRE_MSG(are_files_present,"ERROR!! surface_coupling_exporter::init - prescribed_from_file has fields but no files.");
+    }
+    if (are_files_present) {
+      EKAT_REQUIRE_MSG(are_fields_present,"ERROR!! surface_coupling_exporter::init - prescribed_from_file has files but no fields.");
+    }
+    bool do_export_from_file = are_fields_present and are_files_present;
+    if (do_export_from_file) {
+      // Construct a time interpolation object
+      auto time_interp = util::TimeInterpolation(m_grid,export_from_file_names);
+      for (auto fname : export_from_file_fields) {
+        // Find the index for this field in the list of export fields.
+	auto v_loc = std::find(m_export_field_names_vector.begin(),m_export_field_names_vector.end(),fname);
+	EKAT_REQUIRE_MSG(v_loc != m_export_field_names_vector.end(), "ERROR!! surface_coupling_exporter::init - prescribed_from_file has field with name " << fname << " which can't be found in set of exported fields\n.");
+	auto idx = v_loc - m_export_field_names_vector.begin();
+	EKAT_REQUIRE_MSG(m_export_source_h(idx)==FROM_MODEL,"Error! surface_coupling_exporter::init - attempting to set field " << fname << " export type, which has already been set. Please check namelist options");
+	m_export_source_h(idx) = FROM_FILE;
+        ++m_num_from_file_exports;
+        --m_num_from_model_exports;
+        auto& f_helper = m_helper_fields.at(fname);
+	// We want to add the field as a deep copy so that the helper_fields are automatically updated.
+        time_interp.add_field(f_helper, true);
+        m_export_from_file_field_names.push_back(fname);
+      }
+      m_time_interp = time_interp;
+      m_time_interp.initialize_data_from_files();
+    }
+  }
+
+  // Parse all fields that will be set to a constant value: 
   if (m_params.isSublist("prescribed_constants")) {
     auto export_constant_params = m_params.sublist("prescribed_constants");
     EKAT_REQUIRE_MSG(export_constant_params.isParameter("fields"),"Error! surface_coupling_exporter::init - prescribed_constants does not have 'fields' parameter.");
     EKAT_REQUIRE_MSG(export_constant_params.isParameter("values"),"Error! surface_coupling_exporter::init - prescribed_constants does not have 'values' parameter.");
     auto export_constant_fields = export_constant_params.get<vos_type>("fields");
     auto export_constant_values = export_constant_params.get<vor_type>("values");
     EKAT_REQUIRE_MSG(export_constant_fields.size()==export_constant_values.size(),"Error! surface_coupling_exporter::init - prescribed_constants 'fields' and 'values' are not the same size");
-    if (export_constant_fields.size()>0) {
+    bool are_fields_present = (std::find(export_constant_fields.begin(),export_constant_fields.end(),"NONE") == export_constant_fields.end()) and (export_constant_fields.size() > 0);
+    if (are_fields_present) {
       // Determine which fields need constants
-      for (int i=0; i<m_num_scream_exports; ++i) {  // TODO: This loop would probably be simpler if we just checked which "i" corresponded to each name in the fields list.
-        std::string fname = m_export_field_names[i];
-        auto loc = std::find(export_constant_fields.begin(),export_constant_fields.end(),fname);
-        if (loc != export_constant_fields.end()) {
-          const auto pos = loc-export_constant_fields.begin();
-          m_export_source_h(i) = CONSTANT;
-          ++m_num_const_exports;
-          --m_num_from_model_exports;
-          m_export_constants.emplace(fname,export_constant_values[pos]);
-        }
+      for (int ii=0; ii<export_constant_fields.size(); ii++) {
+        auto fname = export_constant_fields[ii];
+        // Find the index for this field in the list of export fields.
+	auto v_loc = std::find(m_export_field_names_vector.begin(),m_export_field_names_vector.end(),fname);
+	EKAT_REQUIRE_MSG(v_loc != m_export_field_names_vector.end(), "ERROR!! surface_coupling_exporter::init - prescribed_constants has field with name " << fname << " which can't be found in set of exported fields\n.");
+	auto idx = v_loc - m_export_field_names_vector.begin();
+        // This field should not have been set to anything else yet (recall FROM_MODEL is the default)
+        EKAT_REQUIRE_MSG(m_export_source_h(idx)==FROM_MODEL,"Error! surface_coupling_exporter::init - attempting to set field " + fname + " export type, which has already been set.  Please check namelist options");
+        m_export_source_h(idx) = CONSTANT;
+        ++m_num_const_exports;
+        --m_num_from_model_exports;
+        m_export_constants.emplace(fname,export_constant_values[ii]);
       }
     }
   }
   // Copy host view back to device view
   Kokkos::deep_copy(m_export_source,m_export_source_h);
   // Final sanity check
-  EKAT_REQUIRE_MSG(m_num_scream_exports = m_num_const_exports+m_num_from_model_exports,"Error! surface_coupling_exporter - Something went wrong set the type of export for all variables.");
+  EKAT_REQUIRE_MSG(m_num_scream_exports = m_num_from_file_exports+m_num_const_exports+m_num_from_model_exports,"Error! surface_coupling_exporter - Something went wrong set the type of export for all variables.");
   EKAT_REQUIRE_MSG(m_num_from_model_exports>=0,"Error! surface_coupling_exporter - The number of exports derived from EAMxx < 0, something must have gone wrong in assigning the types of exports for all variables.");
 
   // Perform initial export (if any are marked for export during initialization)
@@ -250,8 +296,13 @@ void SurfaceCouplingExporter::run_impl (const double dt)
 void SurfaceCouplingExporter::do_export(const double dt, const bool called_during_initialization)
 {
   if (m_num_const_exports>0) {
-    set_constant_exports(dt,called_during_initialization);
+    set_constant_exports();
   }
+
+  if (m_num_from_file_exports>0) {
+    set_from_file_exports(dt);
+  }
+
   if (m_num_from_model_exports>0) {
     compute_eamxx_exports(dt,called_during_initialization);
   }
@@ -260,17 +311,32 @@ void SurfaceCouplingExporter::do_export(const double dt, const bool called_durin
   do_export_to_cpl(called_during_initialization);
 }
 // =========================================================================================
-void SurfaceCouplingExporter::set_constant_exports(const double dt, const bool called_during_initialization)
+void SurfaceCouplingExporter::set_constant_exports()
 {
   // Cycle through those fields that will be set to a constant value:
+  int num_set = 0; // Checker to make sure we got all the fields we wanted.
   for (int i=0; i<m_num_scream_exports; ++i) {
     if (m_export_source_h(i)==CONSTANT) {
       std::string fname = m_export_field_names[i];
       const auto field_view = m_helper_fields.at(fname).get_view<Real*>();
       Kokkos::deep_copy(field_view,m_export_constants.at(fname));
+      num_set++;
     }
   }
+  // Gotta catch em all
+  EKAT_REQUIRE_MSG(num_set==m_num_const_exports,"ERROR! SurfaceCouplingExporter::set_constant_exports() - Number of fields set to a constant (" + std::to_string(num_set) +") doesn't match the number recorded at initialization (" + std::to_string(m_num_const_exports) +").  Something went wrong.");
 
+}
+// =========================================================================================
+void SurfaceCouplingExporter::set_from_file_exports(const int dt)
+{
+  // Perform interpolation on the data with the latest timestamp
+  auto ts = timestamp();
+  if (dt > 0) {
+    ts += dt;
+  }
+  m_time_interp.perform_time_interpolation(ts);
+
 }
 // =========================================================================================
 // This compute_eamxx_exports routine  handles all export variables that are derived from the EAMxx state.
@@ -479,7 +545,7 @@ void SurfaceCouplingExporter::do_export_to_cpl(const bool called_during_initiali
 // =========================================================================================
 void SurfaceCouplingExporter::finalize_impl()
 {
-
+  // Nothing to do
 }
 // =========================================================================================
 } // namespace scream

components/eamxx/src/control/atmosphere_surface_coupling_exporter.hpp

@@ -4,6 +4,7 @@
 #include "share/atm_process/atmosphere_process.hpp"
 #include "ekat/ekat_parameter_list.hpp"
 #include "share/util/scream_common_physics_functions.hpp"
+#include "share/util/eamxx_time_interpolation.hpp"
 #include "share/atm_process/ATMBufferManager.hpp"
 #include "share/atm_process/SCDataManager.hpp"
 
@@ -79,7 +80,8 @@ class SurfaceCouplingExporter : public AtmosphereProcess
   // which do not have valid entries.
   void do_export(const double dt, const bool called_during_initialization=false);             // Main export routine
   void compute_eamxx_exports(const double dt, const bool called_during_initialization=false); // Export vars are derived from eamxx state
-  void set_constant_exports(const double dt, const bool called_during_initialization=false);  // Export vars are set to a constant
+  void set_constant_exports();                                                                // Export vars are set to a constant
+  void set_from_file_exports(const int dt);                                                   // Export vars are set by interpolation of data from files
   void do_export_to_cpl(const bool called_during_initialization=false);                       // Finish export by copying data to cpl structures.
 
   // Take and store data from SCDataManager
@@ -123,12 +125,16 @@ class SurfaceCouplingExporter : public AtmosphereProcess
   Int m_num_cpl_exports;
 
   // Number of exports from EAMxx and how they will be handled
-  Int                               m_num_scream_exports;
+  int                               m_num_scream_exports;
   view_1d<DefaultDevice,ExportType> m_export_source;
   view_1d<HostDevice,ExportType>    m_export_source_h;
   std::map<std::string,Real>        m_export_constants;
   int                               m_num_from_model_exports=0;
   int                               m_num_const_exports=0;
+  // For exporting from file
+  int                               m_num_from_file_exports=0;
+  util::TimeInterpolation           m_time_interp;
+  std::vector<std::string>          m_export_from_file_field_names;
 
   // Views storing a 2d array with dims (num_cols,num_fields) for cpl export data.
   // The field idx strides faster, since that's what mct does (so we can "view" the
@@ -137,7 +143,8 @@ class SurfaceCouplingExporter : public AtmosphereProcess
   uview_2d<HostDevice,    Real> m_cpl_exports_view_h;
 
   // Array storing the field names for exports
-  name_t* m_export_field_names;
+  name_t*                   m_export_field_names;
+  std::vector<std::string>  m_export_field_names_vector;
 
   // Views storing information for each export
   uview_1d<HostDevice, int>  m_cpl_indices_view;

components/eamxx/src/share/CMakeLists.txt

@@ -34,6 +34,7 @@ set(SHARE_SRC
   util/scream_utils.cpp
   util/eamxx_time_interpolation.cpp
   util/scream_bfbhash.cpp
+  util/eamxx_time_interpolation.cpp
 )
 
 add_library(scream_share ${SHARE_SRC})

components/eamxx/src/share/util/eamxx_time_interpolation.cpp

@@ -27,13 +27,6 @@ TimeInterpolation::TimeInterpolation(
   m_is_data_from_file = true;
 }
 /*-----------------------------------------------------------------------------------------------*/
-// Destructor
-TimeInterpolation::
-~TimeInterpolation ()
-{
-  finalize();
-}
-/*-----------------------------------------------------------------------------------------------*/
 void TimeInterpolation::finalize()
 {
   if (m_is_data_from_file) {
@@ -275,6 +268,10 @@ void TimeInterpolation::set_file_data_triplets(const vos_type& list_of_files) {
         ts_snap += (time_snap*time_mult);
       }
       auto time = ts_snap.seconds_from(ts_ref);
+      // Sanity check that we don't have multiples of the same timesnap
+      EKAT_REQUIRE_MSG(map_of_times_to_vector_idx.count(time)==0,"Error! TimeInterpolation::set_file_data_triplets - The same time step has been encountered more than once in the data files, please check\n"  
+		      << "    TimeStamp: " << ts_snap.to_string() << "\n"  
+		      << "     Filename: " << filename << "\n");
       map_of_times_to_vector_idx.emplace(time,running_idx);
       filenames_tmp.push_back(filename);
       timestamps_tmp.push_back(ts_snap);

components/eamxx/src/share/util/eamxx_time_interpolation.hpp

@@ -23,7 +23,7 @@ class TimeInterpolation {
   TimeInterpolation() = default;
   TimeInterpolation(const grid_ptr_type& grid);
   TimeInterpolation(const grid_ptr_type& grid, const vos_type& list_of_files);
-  ~TimeInterpolation ();
+  ~TimeInterpolation () = default;
 
   // Running the interpolation
   void initialize_timestamps(const TimeStamp& ts_in);