Oscar/seq drydep set h coeff

src/mam4xx/mo_drydep.hpp

@@ -26,10 +26,8 @@ constexpr Real grav = 9.81;
 constexpr Real karman = 0.4;   // from shr_const_mod.F90
 constexpr Real tmelt = 273.15; // from shr_const_mod.F90 via physconst.F90
 
-// nddvels is used only to determine array sizes, and is not involved in
-// logic below. Because we rely on its constancy for C++ array sizes, we
-// must use its maximum value, which is maxspc above
-constexpr int nddvels = mam4::seq_drydep::maxspc;
+// nddvels is equal to number of species in dry deposition list for gases.
+constexpr int nddvels = mam4::seq_drydep::n_drydep;
 
 KOKKOS_INLINE_FUNCTION
 void calculate_uustar(
@@ -560,7 +558,7 @@ void drydep_xactive(
   // define species-dependent parameters (temperature dependent)
   //-------------------------------------------------------------------------------------
   Real heff[nddvels];
-  mam4::seq_drydep::setHCoeff(sfc_temp, heff);
+  mam4::seq_drydep::set_hcoeff_scalar(sfc_temp, heff);
 
   //-------------------------------------------------------------------------------------
   // 	... set month

src/mam4xx/seq_drydep.hpp

@@ -1,6 +1,7 @@
 #ifndef MAM4XX_SEQ_DRYDEP_HPP
 #define MAM4XX_SEQ_DRYDEP_HPP
 
+#include <haero/math.hpp>
 #include <mam4xx/mam4_types.hpp>
 
 namespace mam4::seq_drydep { // C++ version of E3SM's seq_drydep_mod.F90
@@ -14,6 +15,9 @@ constexpr int NSeas = 5;
 // number of land use types
 constexpr int NLUse = 11;
 
+// number of gas species in dry dep list.
+constexpr int n_drydep = 3;
+
 //=========================================
 // data for E3SM dry deposition of tracers
 //=========================================
@@ -55,13 +59,83 @@ struct Data {
   int so2_ndx;
 };
 
-// This function computes the H coefficients corresponding to the given surface
-// temperature. It must be implemented by the client application as an on-device
-// Kokkos function, and must be defined IN THE TRANSLATION UNIT IN WHICH IT'S
-// CALLED because we don't use relocatable device code. The function can be
-// implemented in the atmospheric host model or in a testing environment, for
-// example.
-KOKKOS_FUNCTION void setHCoeff(Real sfc_temp, Real heff[maxspc]);
+// Define the Species enum for dry deposition species identification
+// BAD CONSTANT
+enum class GasDrydepSpecies { H2O2, H2SO4, SO2, CO2, NH3 };
+/**
+ * Calculates Henry's law coefficients based on surface temperature and other
+ * parameters.
+ *
+ * @param sfc_temp The surface temperature in Kelvin. [input]
+ * @param heff Array where the calculated Henry's law coefficients will be
+ * stored. [output]
+ */
+KOKKOS_INLINE_FUNCTION
+void set_hcoeff_scalar(const Real sfc_temp, Real heff[]) {
+
+  // Define dheff array with size n_species_table*6
+  // NOTE: We are hard-coding the table dheff with only 3 species:
+  // H2O2, H2SO4, SO2.
+  // The original table can be found in the seq_drydep_mod.F90 module.
+  // BAD CONSTANT
+  const GasDrydepSpecies drydep_list[n_drydep] = {
+      GasDrydepSpecies::H2O2, GasDrydepSpecies::H2SO4, GasDrydepSpecies::SO2};
+  // --- data for effective Henry's Law coefficient ---
+  constexpr Real dheff[n_drydep * 6] = {
+      8.70e+04, 7320., 2.2e-12,  -3730., 0.,      0.,   // H2O2
+      1.e+11,   6014., 0.,       0.,     0.,      0.,   // H2SO4
+      1.36e+00, 3100., 1.30e-02, 1960.,  6.6e-08, 1500. // SO2
+  };
+  // NOTE: we are using fortran indexing.
+  constexpr int mapping[n_drydep] = {1, 2, 3};
+  // BAD CONSTANT
+  constexpr Real ph = 1.e-5; // measure of the acidity (dimensionless)
+
+  const Real t0 = 298.0;        // Standard Temperature
+  const Real ph_inv = 1.0 / ph; // Inverse of PH
+
+  const Real wrk = (t0 - sfc_temp) / (t0 * sfc_temp);
+
+  for (int m = 0; m < n_drydep; ++m) {
+    const int l = mapping[m] - 1; // Adjust for 0-based indexing
+    const int id = 6 * l;
+    Real e298 = dheff[id];    // Adjusted for 0-based indexing
+    Real dhr = dheff[id + 1]; // Adjusted for 0-based indexing
+    heff[m] = haero::exp(dhr * wrk) * e298;
+
+    // Calculate coefficients based on the drydep tables
+    if (dheff[id + 2] != 0.0 && dheff[id + 4] == 0.0) {
+      e298 = dheff[id + 2];
+      dhr = dheff[id + 3];
+      Real dk1 = haero::exp(dhr * wrk) * e298;
+      heff[m] =
+          (heff[m] != 0.0) ? heff[m] * (1.0 + dk1 * ph_inv) : dk1 * ph_inv;
+    }
+
+    // For coefficients that are non-zero AND CO2 or NH3 handle things this way
+    if (dheff[id + 4] != 0.0) {
+      GasDrydepSpecies species = drydep_list[m];
+      if (species == GasDrydepSpecies::CO2 ||
+          species == GasDrydepSpecies::NH3 ||
+          species == GasDrydepSpecies::SO2) {
+        e298 = dheff[id + 2];
+        dhr = dheff[id + 3];
+        Real dk1 = haero::exp(dhr * wrk) * e298;
+        e298 = dheff[id + 4];
+        dhr = dheff[id + 5];
+        Real dk2 = haero::exp(dhr * wrk) * e298;
+        if (species == GasDrydepSpecies::CO2 ||
+            species == GasDrydepSpecies::SO2) {
+          heff[m] *= (1.0 + dk1 * ph_inv * (1.0 + dk2 * ph_inv));
+        } else if (species == GasDrydepSpecies::NH3) {
+          heff[m] *= (1.0 + dk1 * ph / dk2);
+        } else {
+          EKAT_KERNEL_ERROR_MSG("ERROR: Bad species encountered.\n");
+        }
+      }
+    }
+  }
+}
 
 } // namespace mam4::seq_drydep
 

src/validation/mo_drydep/drydep_xactive.cpp

@@ -1,5 +1,3 @@
-#include "set_h_coeff.hpp" // <-- implementation of seq_drydep::setHCoeff
-
 #include <mam4xx/mam4.hpp>
 
 #include <skywalker.hpp>