Finalize Helios CK (#133)

- reads correctly
- requires extra band info
- still needs test to integrate into RT

.github/workflows/main.yml

@@ -351,7 +351,7 @@ jobs:
         working-directory: ${{github.workspace}}/build-netcdf/bin
         # Execute tests defined by the CMake configuration.
         # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail
-        run: python ./run_ktable_simple.py
+        run: python ./run_ktable_earth.py
 
   straka-2d:
     if: github.event.pull_request.draft == false

cmake/hjupiter.cmake

@@ -14,6 +14,6 @@ set(NETCDF OFF)
 set(PNETCDF ON)
 set(MPI ON)
 set(DISORT ON)
-set(PYTHON_BINDINGS OFF)
+set(PYTHON_BINDINGS ON)
 set_if_empty(RSOLVER hllc_transform)
-#set(GLOG ON)
+# set(GLOG ON)

examples/2024-CLi-earth-rt/CMakeLists.txt

@@ -0,0 +1,13 @@
+# ==========================
+# Examples published in XXXX
+# ==========================
+
+string(TOLOWER ${CMAKE_BUILD_TYPE} buildl)
+string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
+
+# 2. Copy input file to run directory
+file(GLOB inputs *.py *.yaml *.inp)
+foreach(input ${inputs})
+  # softlink inp files
+  execute_process(COMMAND ln -sf ${input} ${CMAKE_BINARY_DIR}/bin/${inp})
+endforeach()

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

@@ -0,0 +1,71 @@
+#! /usr/bin/env python3
+import sys, os
+
+sys.path.append("../python")
+sys.path.append(".")
+
+from pyharp import radiation_band
+from utilities import load_configure, find_resource
+from numpy import *
+from rfmlib import *
+
+
+# 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
+
+    # km
+    atm["HGT"] = linspace(0, 20, nlyr, endpoint=False)
+    atm["PRE"] = ps_mbar * exp(-atm["HGT"] / Hscale_km)
+    atm["TEM"] = Ts_K * ones(nlyr)
+    # ppmv
+    atm["CO2"] = 400.0 * 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
+
+    return atm
+
+
+if __name__ == "__main__":
+    hitran_file = find_resource("HITRAN2020.par")
+
+    info = load_configure("example_earth.yaml")
+    opacity = info["opacity-sources"]
+    band_name = str(info["bands"][0])
+
+    band = radiation_band(band_name, info)
+
+    nspec = band.get_num_spec_grids()
+    wmin, wmax = band.get_range()
+    wres = (wmax - wmin) / (nspec - 1)
+
+    num_absorbers = band.get_num_absorbers()
+    absorbers = []
+    for i in range(num_absorbers):
+        absorbers.append(band.get_absorber(i).get_name())
+
+    # create atmosphere dictionary
+    num_layers = 100
+    wav_grid = (wmin, wmax, wres)
+    tem_grid = (5, -20, 20)
+
+    atm = create_rfm_atmosphere(num_layers)
+    driver = create_rfm_driver(wav_grid, tem_grid, absorbers, hitran_file)
+
+    # write rfm atmosphere file to file
+    write_rfm_atm(atm)
+    write_rfm_drv(driver)
+
+    # run rfm and write kcoeff file
+    run_rfm()
+    write_ktable(f"kcoeff-{band_name}.nc", absorbers, atm, wav_grid, tem_grid)

examples/2024-XZhang-cloud-rt/CMakeLists.txt

@@ -8,7 +8,7 @@ endif()
 setup_problem(hjupiter)
 
 # 2. Copy input files to run directory
-file(GLOB inputs *.inp *.yaml)
+file(GLOB inputs *.py *.inp *.yaml)
 foreach(input ${inputs})
-    file(COPY ${input} DESTINATION ${CMAKE_BINARY_DIR}/bin)
+  execute_process(COMMAND ln -sf ${input} ${CMAKE_BINARY_DIR}/bin/${inp})
 endforeach()

examples/2024-XZhang-cloud-rt/hjupiter.yaml

@@ -7,7 +7,18 @@ opacity-sources:
     class: FreedmanSimple2
     parameters: {scale: 1.}
 
-bands: [ir, vis]
+  - name: premix
+    class: HeliosCK
+    data: ck_data_01242024/ck/PM_ck_HELIOSK_cond_11_nOPT_wcia.txt
+    parameters: {band_id: 1}
+
+bands: [ir, vis, ck]
+
+ck:
+  units: cm-1
+  wavenumber-range: [0.2, 0.4]
+  grid-type: cktable
+  opacity: [premix]
 
 ir:
   units: cm-1
@@ -45,3 +56,8 @@ Disort-flags:
   print-intensity: false
   print-transmissivity: false
   print-phase-function: true
+
+thermodynamics:
+  dry:
+    Rd: 3777.
+    gammad: 1.4

examples/2024-XZhang-cloud-rt/hywater.yaml

@@ -56,9 +56,9 @@ Disort-flags:
   print-phase-function: true
 
 thermodynamics:
-  non-condensable:
+  dry:
     Rd: 3777.
-    gammad_ref: 1.4
+    gammad: 1.4
 
 microphysics:
   - water-system

examples/2024-XZhang-cloud-rt/run_disort_hjup.py

@@ -0,0 +1,46 @@
+#! /usr/bin/env python3
+import sys, os
+
+sys.path.append("../python")
+sys.path.append(".")
+
+from pyharp import radiation_band, init_thermo
+from utilities import load_configure
+from numpy import linspace, ones, exp
+from netCDF4 import Dataset
+from pylab import *
+
+
+# create atmosphere dictionary
+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)
+
+    # mbar
+    atm["PRE"] = data["press"][0, :, 0, 0] / 100.0
+
+    # K
+    atm["TEM"] = data["temp"][0, :, 0, 0]
+
+    return atm
+
+
+if __name__ == "__main__":
+    config = load_configure("hjupiter.yaml")
+    init_thermo(config["thermodynamics"])
+    band = radiation_band("ck", config, load_opacity=True)
+
+    ab = band.get_absorber_by_name("premix")
+    print(ab)
+
+    temp = 300.0
+    pres = 10.0e5
+    v1 = 0.0
+    v2 = 0.0
+    v3 = 0.0
+
+    air = [temp, v1, v2, v3, pres]
+    print(ab.get_attenuation(1000.0, 1000.0, air))

examples/2024-XZhang-cloud-rt/run_ktable_hjup.py

@@ -0,0 +1,63 @@
+#! /usr/bin/env python3
+import sys, os
+
+sys.path.append("../python")
+sys.path.append(".")
+
+from pyharp import radiation_band, subscribe_species
+from utilities import load_configure, find_resource
+from netCDF4 import Dataset
+from numpy import *
+from rfmlib import *
+
+
+# create reference atmosphere dictionary
+# Edit this function to create your own atmosphere
+def create_rfm_atmosphere(nlyr: int) -> dict:
+    atm = {}
+
+    data = Dataset("hjupiter-main.nc", "r")
+    # km
+    atm["HGT"] = (data["x1f"][:-1] + data["x1f"][1:]) / (2.0 * 1.0e3)
+
+    # mbar
+    atm["PRE"] = data["press"][0, :, 0, 0] / 100.0
+
+    # K
+    atm["TEM"] = data["temp"][0, :, 0, 0]
+
+    # ppmv
+    atm["H2"] = 1.0e6
+
+    return atm
+
+
+if __name__ == "__main__":
+    config = load_configure("hjupiter.yaml")
+    band = radiation_band("B1", config)
+
+    nspec = band.get_num_spec_grids()
+    wmin, wmax = band.get_range()
+    wres = (wmax - wmin) / (nspec - 1)
+
+    absorbers = ["H2O", "NH3"]
+
+    # atmospheric properties
+    num_layers = 100
+    wav_grid = (wmin, wmax, wres)
+    tem_grid = (5, -20, 20)
+
+    atm = create_rfm_atmosphere(num_layers)
+    driver = create_rfm_driver(wav_grid, tem_grid, absorbers, hitran_file)
+
+    # write rfm atmosphere file to file
+    write_rfm_atm(atm)
+    write_rfm_drv(driver)
+
+    # run rfm and write kcoeff file
+    run_rfm()
+
+    fname = band.get_absorber_by_name("H2O").get_opacity_file()
+    base_name = os.path.basename(fname)
+    fname, _ = os.path.splitext(base_name)
+    write_ktable(fname, absorbers, atm, wav_grid, tem_grid)

examples/CMakeLists.txt

@@ -32,6 +32,10 @@ if (${TASK} STREQUAL "ktable_jup")
   add_subdirectory(2024-FDing-jupiter-rt)
 endif()
 
+if (${TASK} STREQUAL "ktable_earth")
+  add_subdirectory(2024-CLi-earth-rt)
+endif()
+
 if (${TASK} STREQUAL "exo2" OR
     ${TASK} STREQUAL "exo3")
   add_subdirectory(2023-Chen-exo3)

src/opacity/absorber.cpp

@@ -41,7 +41,7 @@ void Absorber::LoadOpacity(int bid) {
 
 std::string Absorber::ToString() const {
   std::stringstream ss;
-  ss << "Absorber: " << GetName();
+  ss << "Absorber: " << GetName() << std::endl;
   ss << "Opacity file: " << opacity_filename_;
 
   return ss.str();

src/opacity/absorber_ck.cpp

@@ -61,13 +61,3 @@ Real AbsorberCK::GetAttenuation(Real g1, Real g2, AirParcel const& var) const {
   Real dens = var.q[IPR] / (Constants::kBoltz * var.q[IDN]);
   return exp(val) * dens;  // ln(m*2/kmol) -> 1/m
 }
-
-void HeliosCK::ModifySpectralGrid(std::vector<SpectralBin>& spec) const {
-  spec.resize(weights_.size());
-
-  for (size_t i = 0; i < weights_.size(); ++i) {
-    spec[i].wav1 = axis_[len_[0] + len_[1] + i];
-    spec[i].wav2 = axis_[len_[0] + len_[1] + i];
-    spec[i].wght = weights_[i];
-  }
-}

src/opacity/absorber_ck.hpp

@@ -32,7 +32,7 @@ class AbsorberCK : public Absorber {
 
 class HeliosCK : public AbsorberCK {
  public:
-  HeliosCK(std::string name) : AbsorberCK(name) {}
+  HeliosCK(std::string name);
   virtual ~HeliosCK() {}
 
   void LoadCoefficient(std::string fname, int bid) override;

src/opacity/helios_ck_premix.cpp → src/opacity/helios_ck.cpp

@@ -8,29 +8,53 @@
 #include <air_parcel.hpp>
 #include <constants.hpp>
 
+// application
+#include <application/exceptions.hpp>
+
 // climath
 #include <climath/interpolation.h>
 
+// harp
+#include <harp/spectral_grid.hpp>
+
 // opacity
 #include "absorber_ck.hpp"
 
+HeliosCK::HeliosCK(std::string name) : AbsorberCK(name) {}
+
+void HeliosCK::ModifySpectralGrid(std::vector<SpectralBin>& spec) const {
+  spec.resize(weights_.size());
+
+  for (size_t i = 0; i < weights_.size(); ++i) {
+    spec[i].wav1 = axis_[len_[0] + len_[1] + i];
+    spec[i].wav2 = axis_[len_[0] + len_[1] + i];
+    spec[i].wght = weights_[i];
+  }
+}
+
 void HeliosCK::LoadCoefficient(std::string fname, int bid) {
   std::ifstream file(fname);
   if (!file.is_open()) {
     throw std::runtime_error("Failed to open file: " + fname);
   }
 
+  // over ride band id if it is set
+  if (HasPar("band_id")) {
+    bid = static_cast<int>(GetPar<Real>("band_id"));
+  }
+
   // skip the first line
   std::string junk_line;
   std::getline(file, junk_line);
 
-  size_t num_bands;
+  int num_bands;
 
   // temperature, pressure, band, g-points
   file >> len_[0] >> len_[1] >> num_bands >> len_[2];
 
-  if (bid >= num_bands) {
-    throw std::runtime_error("Band index out of range: " + std::to_string(bid));
+  if (bid >= num_bands || bid < 0) {
+    throw RuntimeError("HeliosCK::LoadCoefficient",
+                       "Band index out of range: " + std::to_string(bid));
   }
 
   axis_.resize(len_[0] + len_[1] + len_[2]);
@@ -62,18 +86,22 @@ void HeliosCK::LoadCoefficient(std::string fname, int bid) {
 
   // skip unimportant wavelengths
   Real dummy;
-  for (int i = 0; i < num_bands - bid; ++i) {
+  for (int i = 1; i < num_bands - bid; ++i) {
     file >> dummy;
   }
 
   // g-points and weights
   weights_.resize(len_[2]);
   for (int g = 0; g < len_[2]; ++g) {
     Real gpoint;
-    file >> gpoint >> weights_[g];
+    file >> gpoint;
     axis_[len_[0] + len_[1] + g] = wmin + (wmax - wmin) * gpoint;
   }
 
+  for (int g = 0; g < len_[2]; ++g) {
+    file >> weights_[g];
+  }
+
   int n = 0;
   for (int i = 0; i < len_[0]; ++i)
     for (int j = 0; j < len_[1]; ++j)
@@ -84,7 +112,7 @@ void HeliosCK::LoadCoefficient(std::string fname, int bid) {
             kcoeff_[n] = log(std::max(kcoeff_[n], 1.0e-99));
           }
         } else {
-          for (int g = 0; g < len_[2]; ++g, ++n) {
+          for (int g = 0; g < len_[2]; ++g) {
             file >> dummy;
           }
         }