Update pyharp (#134)

- change unit to SI
- change get_tau to get_dtau

examples/2024-CLi-earth-rt/run_ktable_earth.py

@@ -13,25 +13,24 @@
 # create atmosphere dictionary
 def create_rfm_atmosphere(nlyr: int) -> dict:
     atm = {}
-    ps_mbar = 1.0e3
-    Ts_K = 300.0
-    grav_ms2 = 9.8
-    mu_kgmol = 29.0e-3
-    Rgas_JKmol = 8.314
-    Rgas_Jkg = Rgas_JKmol / mu_kgmol
-    Hscale_km = Rgas_Jkg * Ts_K / grav_ms2 * 1.0e-3
+    Ps = 1.0e5
+    Ts = 300.0
+    grav = 9.8
+    mu = 29.0e-3
+    Rgas = 8.314 / mu
+    Hscale = Rgas * Ts / grav
 
     # km
     atm["HGT"] = linspace(0, 20, nlyr, endpoint=False)
-    atm["PRE"] = ps_mbar * exp(-atm["HGT"] / Hscale_km)
-    atm["TEM"] = Ts_K * ones(nlyr)
+    atm["PRE"] = Ps * exp(-atm["HGT"] / Hscale)
+    atm["TEM"] = Ts * ones(nlyr)
     # ppmv
-    atm["CO2"] = 400.0 * ones(nlyr)
+    atm["CO2"] = 4.0e-4 * ones(nlyr)
     # ppmv
-    atm["H2O"] = 4000.0 * exp(-atm["HGT"] / (Hscale_km / 5.0))
-    atm["O2"] = (1e6 - atm["CO2"] - atm["H2O"]) * 0.21
-    atm["N2"] = (1e6 - atm["CO2"] - atm["H2O"]) * 0.78
-    atm["Ar"] = (1e6 - atm["CO2"] - atm["H2O"]) * 0.01
+    atm["H2O"] = 4.0e-3 * exp(-atm["HGT"] / (Hscale / 5.0))
+    atm["O2"] = (1.0 - atm["CO2"] - atm["H2O"]) * 0.21
+    atm["N2"] = (1.0 - atm["CO2"] - atm["H2O"]) * 0.78
+    atm["Ar"] = (1.0 - atm["CO2"] - atm["H2O"]) * 0.01
 
     return atm
 

examples/2024-FDing-jupiter-rt/run_disort_jup.py

@@ -16,20 +16,20 @@ def create_atmosphere(nlyr: int) -> dict:
     atm = {}
 
     data = Dataset("jupiter1d-main.nc", "r")
-    # km
-    atm["HGT"] = (data["x1f"][:-1] + data["x1f"][1:]) / (2.0 * 1.0e3)
+    # m
+    atm["HGT"] = (data["x1f"][:-1] + data["x1f"][1:]) / 2.0
 
-    # mbar
-    atm["PRE"] = data["press"][0, :, 0, 0] / 100.0
+    # pa
+    atm["PRE"] = data["press"][0, :, 0, 0]
 
     # K
     atm["TEM"] = data["temp"][0, :, 0, 0]
 
-    # ppmv
-    atm["H2O"] = data["vapor1"][0, :, 0, 0] / 18.0 * 2.2 * 1.0e6
+    # mole fraction
+    atm["H2O"] = data["vapor1"][0, :, 0, 0] / 18.0 * 2.2
 
-    # ppmv
-    atm["NH3"] = data["vapor2"][0, :, 0, 0] / 17.0 * 2.2 * 1.0e6
+    # mole fraction
+    atm["NH3"] = data["vapor2"][0, :, 0, 0] / 17.0 * 2.2
 
     return atm
 
@@ -66,8 +66,9 @@ def create_atmosphere(nlyr: int) -> dict:
     num_layers = 100
     band.resize(num_layers)
 
-    # atm = create_atmosphere(num_layers)
-    # band.set_spectral_properties(atm)
+    atm = create_atmosphere(num_layers)
+    band.set_spectral_properties(atm)
 
-    tau = band.get_tau()
-    print(tau.shape)
+    dtau = band.get_dtau()
+    print(dtau.shape)
+    print(dtau)

examples/2024-FDing-jupiter-rt/run_ktable_jup.py

@@ -17,23 +17,23 @@ def create_rfm_atmosphere(nlyr: int) -> dict:
     atm = {}
 
     data = Dataset("jupiter1d-main.nc", "r")
-    # km
-    atm["HGT"] = (data["x1f"][:-1] + data["x1f"][1:]) / (2.0 * 1.0e3)
+    # m
+    atm["HGT"] = (data["x1f"][:-1] + data["x1f"][1:]) / 2.0
 
-    # mbar
-    atm["PRE"] = data["press"][0, :, 0, 0] / 100.0
+    # pa
+    atm["PRE"] = data["press"][0, :, 0, 0]
 
     # K
     atm["TEM"] = data["temp"][0, :, 0, 0]
 
-    # ppmv
-    atm["H2O"] = data["vapor1"][0, :, 0, 0] / 18.0 * 2.2 * 1.0e6
+    # mole fraction
+    atm["H2O"] = data["vapor1"][0, :, 0, 0] / 18.0 * 2.2
 
-    # ppmv
-    atm["NH3"] = data["vapor2"][0, :, 0, 0] / 17.0 * 2.2 * 1.0e6
+    # mole fraction
+    atm["NH3"] = data["vapor2"][0, :, 0, 0] / 17.0 * 2.2
 
-    # ppmv
-    atm["H2"] = 1.0e6 - atm["H2O"] - atm["NH3"]
+    # mole fraction
+    atm["H2"] = 1.0 - atm["H2O"] - atm["NH3"]
 
     return atm
 

python/pyharp.cpp

@@ -128,12 +128,14 @@ PYBIND11_MODULE(pyharp, m) {
              return ndarray;
            })
 
-      .def("get_tau",
+      .def("get_dtau",
            [](RadiationBand &band) {
-             py::array_t<Real> ndarray(
+             py::capsule tau_capsule(band._GetTauPtr(), [](void *) {});
+             py::array_t<double> tau(
                  {band.GetNumSpecGrids(), band.GetNumLayers()},
-                 band._GetTauPtr());
-             return ndarray;
+                 band._GetTauPtr(), tau_capsule);
+
+             return tau;
            })
 
       .def("set_spectral_properties",
@@ -148,9 +150,9 @@ PYBIND11_MODULE(pyharp, m) {
 
              for (int i = 0; i < nlayer; ++i) {
                ac[i].SetType(AirParcel::Type::MoleFrac);
-               x1v[i] = atm["HGT"][i] * 1e3;  // km -> m
+               x1v[i] = atm["HGT"][i];
                ac[i].w[IDN] = atm["TEM"][i];
-               ac[i].w[IPR] = atm["PRE"][i] * 100.;  // mbar -> Pa
+               ac[i].w[IPR] = atm["PRE"][i];
              }
 
              for (int i = 1; i < nlayer; ++i) {
@@ -169,26 +171,22 @@ PYBIND11_MODULE(pyharp, m) {
 
                if (pindex->HasVapor(pair.first))
                  for (int i = 0; i < nlayer; ++i) {
-                   ac[i].w[pindex->GetVaporId(pair.first)] =
-                       pair.second[i] * 1e-6;
+                   ac[i].w[pindex->GetVaporId(pair.first)] = pair.second[i];
                  }
 
                if (pindex->HasCloud(pair.first))
                  for (int i = 0; i < nlayer; ++i) {
-                   ac[i].c[pindex->GetCloudId(pair.first)] =
-                       pair.second[i] * 1e-6;
+                   ac[i].c[pindex->GetCloudId(pair.first)] = pair.second[i];
                  }
 
                if (pindex->HasChemistry(pair.first))
                  for (int i = 0; i < nlayer; ++i) {
-                   ac[i].q[pindex->GetChemistryId(pair.first)] =
-                       pair.second[i] * 1e-6;
+                   ac[i].q[pindex->GetChemistryId(pair.first)] = pair.second[i];
                  }
 
                if (pindex->HasTracer(pair.first))
                  for (int i = 0; i < nlayer; ++i) {
-                   ac[i].x[pindex->GetTracerId(pair.first)] =
-                       pair.second[i] * 1e-6;
+                   ac[i].x[pindex->GetTracerId(pair.first)] = pair.second[i];
                  }
              }
 

tools/rfmlib.py

@@ -34,20 +34,20 @@ def write_rfm_atm(atm: dict) -> None:
     with open("rfm.atm", "w") as file:
         file.write("%d\n" % num_layers)
         file.write("*HGT [km]\n")
-        for i in range(num_layers):
-            file.write("%.8g " % atm["HGT"][i])
+        for i in range(num_layers):  # m -> km
+            file.write("%.8g " % (atm["HGT"][i] / 1.0e3,))
         file.write("\n*PRE [mb]\n")
-        for i in range(num_layers):
-            file.write("%.8g " % atm["PRE"][i])
+        for i in range(num_layers):  # pa -> mb
+            file.write("%.8g " % (atm["PRE"][i] / 100.0,))
         file.write("\n*TEM [K]\n")
         for i in range(num_layers):
             file.write("%.8g " % atm["TEM"][i])
         for name, val in atm.items():
             if name in ["HGT", "PRE", "TEM"]:
                 continue
             file.write("\n*" + name + " [ppmv]\n")
-            for j in range(num_layers):
-                file.write("%.8g " % val[j])
+            for j in range(num_layers):  # mol/mol -> ppmv
+                file.write("%.8g " % (val[j] * 1.0e6,))
         file.write("\n*END")
     print("# rfm.atm written.")