modify theta calculation for coupling exporter

Author: whannah1

components/eamxx/src/control/atmosphere_surface_coupling_exporter.cpp

@@ -440,10 +440,12 @@ void SurfaceCouplingExporter::compute_eamxx_exports(const double dt, const bool
     const auto qv_i             = ekat::subview(qv, i);
     const auto T_mid_i          = ekat::subview(T_mid, i);
     const auto p_mid_i          = ekat::subview(p_mid, i);
+    const auto p_int_i          = ekat::subview(p_int, i);
     const auto pseudo_density_i = ekat::subview(pseudo_density, i);
     const auto dz_i             = ekat::subview(dz, i);
 
     const auto s_p_mid_i = ekat::scalarize(p_mid_i);
+    const auto s_p_int_i = ekat::scalarize(p_int_i);
     const auto s_T_mid_i = ekat::scalarize(T_mid_i);
     const auto z_int_i = ekat::subview(z_int, i);
     const auto z_mid_i = ekat::subview(z_mid, i);
@@ -487,7 +489,12 @@ void SurfaceCouplingExporter::compute_eamxx_exports(const double dt, const bool
     }
 
     if (export_source(idx_Sa_ptem)==FROM_MODEL) {
-      Sa_ptem(i) = PF::calculate_theta_from_T(s_T_mid_i(num_levs-1), s_p_mid_i(num_levs-1));
+      // WARNING - THE FOLLOWING IS A HACK
+      // To make the flux calculations within the component coupler consistent with EAM we need to
+      // provide theta based on an exner function that evaluates to 1 at the bottom interface.
+      // To accomplish this we calculate a theta that replaces the reference pressure (P0) for exner
+      // with the pressure of the lowest interface level => s_p_int_i(num_levs)
+      Sa_ptem(i) = s_T_mid_i(num_levs-1) / pow( s_p_mid_i(num_levs-1)/s_p_int_i(num_levs), PC::RD*PC::INV_CP);
     }
 
     if (export_source(idx_Sa_pbot)==FROM_MODEL) {

components/eamxx/tests/single-process/surface_coupling/surface_coupling.cpp

@@ -361,10 +361,16 @@ void test_exports(const FieldManager& fm,
     const Real T_int_bot = PF::calculate_surface_air_T(T_mid_i(nlevs-1),z_mid_i(nlevs-1));
 
     Sa_z(i)       = z_mid_i(nlevs-1);
-    Sa_ptem(i)    = PF::calculate_theta_from_T(T_mid_i(nlevs-1), p_mid_i(nlevs-1));
     Sa_dens(i)    = PF::calculate_density(pseudo_density_i(nlevs-1), dz_i(nlevs-1));
     Sa_pslv(i)    = PF::calculate_psl(T_int_bot, p_int_i(nlevs), phis(i));
 
+    // WARNING - THE FOLLOWING IS A HACK
+    // To make the flux calculations within the component coupler consistent with EAM we need to
+    // provide theta based on an exner function that evaluates to 1 at the bottom interface.
+    // To accomplish this we calculate a theta that replaces the reference pressure (P0) for exner
+    // with the pressure of the lowest interface level => p_int_i(nlevs)
+    Sa_ptem(i) = T_mid_i(nlevs-1) / pow( p_mid_i(nlevs-1)/p_int_i(nlevs), PC::RD*PC::INV_CP);
+
     if (not called_directly_after_init) {
       Faxa_rainl(i) = precip_liq_surf_mass(i)/dt*(1000.0/PC::RHO_H2O);
       Faxa_snowl(i) = precip_ice_surf_mass(i)/dt*(1000.0/PC::RHO_H2O);