diff --git a/components/eamxx/src/dynamics/homme/homme_grids_manager.cpp b/components/eamxx/src/dynamics/homme/homme_grids_manager.cpp
index 87009c7d074..df5de6827f6 100644
--- a/components/eamxx/src/dynamics/homme/homme_grids_manager.cpp
+++ b/components/eamxx/src/dynamics/homme/homme_grids_manager.cpp
@@ -271,6 +271,7 @@ build_physics_grid (const ci_string& type, const ci_string& rebalance) {
     auto hyam = phys_grid->create_geometry_data("hyam",layout_mid,nondim);
     auto hybm = phys_grid->create_geometry_data("hybm",layout_mid,nondim);
     auto lev  = phys_grid->create_geometry_data("lev", layout_mid,mbar);
+    auto ilev = phys_grid->create_geometry_data("ilev",layout_int,mbar);
 
     for (auto f : {hyai, hybi, hyam, hybm}) {
       auto f_d = get_grid("Dynamics")->get_geometry_data(f.name());
@@ -281,13 +282,20 @@ build_physics_grid (const ci_string& type, const ci_string& rebalance) {
     // Build lev from hyam and hybm
     const Real ps0        = 100000.0;
   
-    auto hya_v = hyam.get_view<const Real*,Host>();
-    auto hyb_v = hybm.get_view<const Real*,Host>();
-    auto lev_v = lev.get_view<Real*,Host>();
-    for (int ii=0;ii<phys_grid->get_num_vertical_levels();ii++) {
-      lev_v(ii) = 0.01*ps0*(hya_v(ii)+hyb_v(ii));
+    auto hyam_v = hyam.get_view<const Real*,Host>();
+    auto hybm_v = hybm.get_view<const Real*,Host>();
+    auto hyai_v = hyai.get_view<const Real*,Host>();
+    auto hybi_v = hybi.get_view<const Real*,Host>();
+    auto lev_v  = lev.get_view<Real*,Host>();
+    auto ilev_v = ilev.get_view<Real*,Host>();
+    auto num_v_levs = phys_grid->get_num_vertical_levels();
+    for (int ii=0;ii<num_v_levs;ii++) {
+      lev_v(ii)  = 0.01*ps0*(hyam_v(ii)+hybm_v(ii));
+      ilev_v(ii) = 0.01*ps0*(hyai_v(ii)+hybi_v(ii));
     }
+    ilev_v(num_v_levs) = 0.01*ps0*(hyai_v(num_v_levs)+hybi_v(num_v_levs));
     lev.sync_to_dev();
+    ilev.sync_to_dev();
   }
 
   if (is_planar_geometry_f90()) {
diff --git a/components/eamxx/src/share/grid/mesh_free_grids_manager.cpp b/components/eamxx/src/share/grid/mesh_free_grids_manager.cpp
index c930013f839..b083dfcfa48 100644
--- a/components/eamxx/src/share/grid/mesh_free_grids_manager.cpp
+++ b/components/eamxx/src/share/grid/mesh_free_grids_manager.cpp
@@ -147,24 +147,30 @@ add_geo_data (const nonconstgrid_ptr_type& grid) const
   if (geo_data_source=="CREATE_EMPTY_DATA") {
     using namespace ShortFieldTagsNames;
     FieldLayout layout_mid ({LEV},{grid->get_num_vertical_levels()});
+    FieldLayout layout_int ({ILEV},{grid->get_num_vertical_levels()+1});
     const auto units = ekat::units::Units::nondimensional();
 
     auto lat  = grid->create_geometry_data("lat" ,  grid->get_2d_scalar_layout(), units);
     auto lon  = grid->create_geometry_data("lon" ,  grid->get_2d_scalar_layout(), units);
     auto hyam = grid->create_geometry_data("hyam" , layout_mid, units);
     auto hybm = grid->create_geometry_data("hybm" , layout_mid, units);
+    auto hyai = grid->create_geometry_data("hyai" , layout_int, units);
+    auto hybi = grid->create_geometry_data("hybi" , layout_int, units);
     auto lev  = grid->create_geometry_data("lev" ,  layout_mid, units);
+    auto ilev = grid->create_geometry_data("ilev" , layout_int, units);
 
     lat.deep_copy(ekat::ScalarTraits<Real>::invalid());
     lon.deep_copy(ekat::ScalarTraits<Real>::invalid());
     hyam.deep_copy(ekat::ScalarTraits<Real>::invalid());
     hybm.deep_copy(ekat::ScalarTraits<Real>::invalid());
     lev.deep_copy(ekat::ScalarTraits<Real>::invalid());
+    ilev.deep_copy(ekat::ScalarTraits<Real>::invalid());
     lat.sync_to_dev();
     lon.sync_to_dev();
     hyam.sync_to_dev();
     hybm.sync_to_dev();
     lev.sync_to_dev();
+    ilev.sync_to_dev();
   } else if (geo_data_source=="IC_FILE"){
     const auto& filename = m_params.get<std::string>("ic_filename");
     if (scorpio::has_var(filename,"lat") &&
@@ -173,7 +179,9 @@ add_geo_data (const nonconstgrid_ptr_type& grid) const
     }
 
     if (scorpio::has_var(filename,"hyam") &&
-        scorpio::has_var(filename,"hybm")) {
+        scorpio::has_var(filename,"hybm") &&
+        scorpio::has_var(filename,"hyai") &&
+        scorpio::has_var(filename,"hybi") ) {
       load_vertical_coordinates(grid,filename);
     }
   }
@@ -234,53 +242,71 @@ load_vertical_coordinates (const nonconstgrid_ptr_type& grid, const std::string&
   using namespace ekat::units;
 
   FieldLayout layout_mid ({LEV},{grid->get_num_vertical_levels()});
+  FieldLayout layout_int ({ILEV},{grid->get_num_vertical_levels()+1});
   Units nondim = Units::nondimensional();
   Units mbar (100*Pa,"mb");
 
   auto hyam = grid->create_geometry_data("hyam", layout_mid, nondim);
   auto hybm = grid->create_geometry_data("hybm", layout_mid, nondim);
+  auto hyai = grid->create_geometry_data("hyai", layout_int, nondim);
+  auto hybi = grid->create_geometry_data("hybi", layout_int, nondim);
   auto lev  = grid->create_geometry_data("lev",  layout_mid, mbar);
+  auto ilev = grid->create_geometry_data("ilev", layout_int, mbar);
 
   // Create host mirrors for reading in data
   std::map<std::string,geo_view_host> host_views = {
     { "hyam", hyam.get_view<Real*,Host>() },
-    { "hybm", hybm.get_view<Real*,Host>() }
+    { "hybm", hybm.get_view<Real*,Host>() },
+    { "hyai", hyai.get_view<Real*,Host>() },
+    { "hybi", hybi.get_view<Real*,Host>() }
   };
 
   // Store view layouts
   using namespace ShortFieldTagsNames;
   std::map<std::string,FieldLayout> layouts = {
     { "hyam", hyam.get_header().get_identifier().get_layout() },
-    { "hybm", hybm.get_header().get_identifier().get_layout() }
+    { "hybm", hybm.get_header().get_identifier().get_layout() },
+    { "hyai", hyai.get_header().get_identifier().get_layout() },
+    { "hybi", hybi.get_header().get_identifier().get_layout() }
   };
 
   // Read hyam/hybm into host views
   ekat::ParameterList vcoord_reader_pl;
   vcoord_reader_pl.set("Filename",filename);
-  vcoord_reader_pl.set<std::vector<std::string>>("Field Names",{"hyam","hybm"});
+  vcoord_reader_pl.set<std::vector<std::string>>("Field Names",{"hyam","hybm","hyai","hybi"});
 
   AtmosphereInput vcoord_reader(vcoord_reader_pl,grid, host_views, layouts);
   vcoord_reader.read_variables();
   vcoord_reader.finalize();
 
-  // Build lev from hyam and hybm
+  // Build lev and ilev from hyam and hybm, and ilev from hyai and hybi
   using PC             = scream::physics::Constants<Real>;
   const Real ps0        = PC::P0;
 
-  auto hya_v = hyam.get_view<const Real*,Host>();
-  auto hyb_v = hybm.get_view<const Real*,Host>();
-  auto lev_v = lev.get_view<Real*,Host>();
-  for (int ii=0;ii<grid->get_num_vertical_levels();ii++) {
-    lev_v(ii) = 0.01*ps0*(hya_v(ii)+hyb_v(ii));
+  auto hyam_v  = hyam.get_view<const Real*,Host>();
+  auto hybm_v  = hybm.get_view<const Real*,Host>();
+  auto lev_v   = lev.get_view<Real*,Host>();
+  auto hyai_v  = hyai.get_view<const Real*,Host>();
+  auto hybi_v  = hybi.get_view<const Real*,Host>();
+  auto ilev_v  = ilev.get_view<Real*,Host>();
+  auto num_lev = grid->get_num_vertical_levels();
+  for (int ii=0;ii<num_lev;ii++) {
+    lev_v(ii)  = 0.01*ps0*(hyam_v(ii)+hybm_v(ii));
+    ilev_v(ii) = 0.01*ps0*(hyai_v(ii)+hybi_v(ii));
   }
+  // Note, ilev is just 1 more level than the number of midpoint levs
+  ilev_v(num_lev) = 0.01*ps0*(hyai_v(num_lev)+hybi_v(num_lev));
 
   // Sync to dev
   hyam.sync_to_dev();
   hybm.sync_to_dev();
+  hyai.sync_to_dev();
+  hybi.sync_to_dev();
   lev.sync_to_dev();
+  ilev.sync_to_dev();
 
 #ifndef NDEBUG
-  for (auto f : {hyam, hybm}) {
+  for (auto f : {hyam, hybm, hyai, hybi}) {
     auto nan_check = std::make_shared<FieldNaNCheck>(f,grid)->check();
     EKAT_REQUIRE_MSG (nan_check.result==CheckResult::Pass,
         "ERROR! NaN values detected in " + f.name() + " field.\n" + nan_check.msg);
diff --git a/components/eamxx/src/share/io/scream_io_utils.hpp b/components/eamxx/src/share/io/scream_io_utils.hpp
index 01bc46e1e60..d6c36ab3758 100644
--- a/components/eamxx/src/share/io/scream_io_utils.hpp
+++ b/components/eamxx/src/share/io/scream_io_utils.hpp
@@ -85,8 +85,9 @@ struct LongNames {
 
   // Create map of longnames, can be added to as developers see fit.
   std::map<std::string,std::string> name_2_longname = {
-	  {"lev","hybrid level at midpoints (1000*(A+B))"},
-	  {"hyai","hybrid A coefficient at layer interfaces"},
+    {"lev","hybrid level at midpoints (1000*(A+B))"},
+    {"ilev","hybrid level at interfaces (1000*(A+B))"},
+    {"hyai","hybrid A coefficient at layer interfaces"},
     {"hybi","hybrid B coefficient at layer interfaces"},
     {"hyam","hybrid A coefficient at layer midpoints"},
     {"hybm","hybrid B coefficient at layer midpoints"}
diff --git a/components/eamxx/src/share/io/scream_output_manager.cpp b/components/eamxx/src/share/io/scream_output_manager.cpp
index fd05ca193b9..050a55936b4 100644
--- a/components/eamxx/src/share/io/scream_output_manager.cpp
+++ b/components/eamxx/src/share/io/scream_output_manager.cpp
@@ -402,7 +402,7 @@ void OutputManager::run(const util::TimeStamp& timestamp)
     // Check if we need to open a new file
     if (not filespecs.is_open) {
       filespecs.filename = compute_filename (filespecs,timestamp);
-      // Register all dims/vars, write geometry data (e.g. lat/lon/hyam/hybm)
+      // Register all dims/vars, write geometry data (e.g. lat/lon/hyam/hybm/hyai/hybi)
       setup_file(filespecs,control);
     }