From 9de04e87095f462c5251484194de4eb8adc37d72 Mon Sep 17 00:00:00 2001 From: Cheng Li <69489965+chengcli@users.noreply.github.com> Date: Mon, 18 Mar 2024 23:35:34 -0400 Subject: [PATCH] Finalize Helios CK (#133) - reads correctly - requires extra band info - still needs test to integrate into RT --- .github/workflows/main.yml | 2 +- cmake/hjupiter.cmake | 4 +- examples/2024-CLi-earth-rt/CMakeLists.txt | 13 ++ .../2024-CLi-earth-rt/earth_flux.yaml | 0 .../2024-CLi-earth-rt/earth_rt.yaml | 0 .../2024-CLi-earth-rt}/example_earth.yaml | 0 .../2024-CLi-earth-rt/run_disort_earth.py | 0 .../run_disort_earth_flux.py | 0 .../2024-CLi-earth-rt/run_ktable_earth.py | 71 +++++++ examples/2024-XZhang-cloud-rt/CMakeLists.txt | 4 +- examples/2024-XZhang-cloud-rt/hjupiter.yaml | 18 +- examples/2024-XZhang-cloud-rt/hywater.yaml | 4 +- .../2024-XZhang-cloud-rt/run_disort_hjup.py | 46 ++++ .../2024-XZhang-cloud-rt/run_ktable_hjup.py | 63 ++++++ examples/CMakeLists.txt | 4 + src/opacity/absorber.cpp | 2 +- src/opacity/absorber_ck.cpp | 10 - src/opacity/absorber_ck.hpp | 2 +- .../{helios_ck_premix.cpp => helios_ck.cpp} | 40 +++- src/snap/thermodynamics/thermodynamics.cpp | 4 +- tools/run_ktable_simple.py | 199 ------------------ 21 files changed, 259 insertions(+), 227 deletions(-) create mode 100644 examples/2024-CLi-earth-rt/CMakeLists.txt rename tools/example_earth_flux.yaml => examples/2024-CLi-earth-rt/earth_flux.yaml (100%) rename tools/example_earth_rt.yaml => examples/2024-CLi-earth-rt/earth_rt.yaml (100%) rename {tools => examples/2024-CLi-earth-rt}/example_earth.yaml (100%) rename tools/run_disort_rt.py => examples/2024-CLi-earth-rt/run_disort_earth.py (100%) rename tools/run_disort_flux.py => examples/2024-CLi-earth-rt/run_disort_earth_flux.py (100%) create mode 100755 examples/2024-CLi-earth-rt/run_ktable_earth.py create mode 100755 examples/2024-XZhang-cloud-rt/run_disort_hjup.py create mode 100755 examples/2024-XZhang-cloud-rt/run_ktable_hjup.py rename src/opacity/{helios_ck_premix.cpp => helios_ck.cpp} (69%) delete mode 100755 tools/run_ktable_simple.py diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 5618ec42..b8c76197 100644 --- a/.github/workflows/main.yml +++ b/.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 diff --git a/cmake/hjupiter.cmake b/cmake/hjupiter.cmake index f29106b7..f63363be 100644 --- a/cmake/hjupiter.cmake +++ b/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) diff --git a/examples/2024-CLi-earth-rt/CMakeLists.txt b/examples/2024-CLi-earth-rt/CMakeLists.txt new file mode 100644 index 00000000..af12c377 --- /dev/null +++ b/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() diff --git a/tools/example_earth_flux.yaml b/examples/2024-CLi-earth-rt/earth_flux.yaml similarity index 100% rename from tools/example_earth_flux.yaml rename to examples/2024-CLi-earth-rt/earth_flux.yaml diff --git a/tools/example_earth_rt.yaml b/examples/2024-CLi-earth-rt/earth_rt.yaml similarity index 100% rename from tools/example_earth_rt.yaml rename to examples/2024-CLi-earth-rt/earth_rt.yaml diff --git a/tools/example_earth.yaml b/examples/2024-CLi-earth-rt/example_earth.yaml similarity index 100% rename from tools/example_earth.yaml rename to examples/2024-CLi-earth-rt/example_earth.yaml diff --git a/tools/run_disort_rt.py b/examples/2024-CLi-earth-rt/run_disort_earth.py similarity index 100% rename from tools/run_disort_rt.py rename to examples/2024-CLi-earth-rt/run_disort_earth.py diff --git a/tools/run_disort_flux.py b/examples/2024-CLi-earth-rt/run_disort_earth_flux.py similarity index 100% rename from tools/run_disort_flux.py rename to examples/2024-CLi-earth-rt/run_disort_earth_flux.py diff --git a/examples/2024-CLi-earth-rt/run_ktable_earth.py b/examples/2024-CLi-earth-rt/run_ktable_earth.py new file mode 100755 index 00000000..324d28ce --- /dev/null +++ b/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) diff --git a/examples/2024-XZhang-cloud-rt/CMakeLists.txt b/examples/2024-XZhang-cloud-rt/CMakeLists.txt index 0a3509e0..a4f57ed9 100644 --- a/examples/2024-XZhang-cloud-rt/CMakeLists.txt +++ b/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() diff --git a/examples/2024-XZhang-cloud-rt/hjupiter.yaml b/examples/2024-XZhang-cloud-rt/hjupiter.yaml index ec2a6cd8..616f8219 100644 --- a/examples/2024-XZhang-cloud-rt/hjupiter.yaml +++ b/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 diff --git a/examples/2024-XZhang-cloud-rt/hywater.yaml b/examples/2024-XZhang-cloud-rt/hywater.yaml index cca45592..d2f31cd3 100644 --- a/examples/2024-XZhang-cloud-rt/hywater.yaml +++ b/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 diff --git a/examples/2024-XZhang-cloud-rt/run_disort_hjup.py b/examples/2024-XZhang-cloud-rt/run_disort_hjup.py new file mode 100755 index 00000000..1311c9d1 --- /dev/null +++ b/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)) diff --git a/examples/2024-XZhang-cloud-rt/run_ktable_hjup.py b/examples/2024-XZhang-cloud-rt/run_ktable_hjup.py new file mode 100755 index 00000000..db14382e --- /dev/null +++ b/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) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index d54d747f..9394182c 100644 --- a/examples/CMakeLists.txt +++ b/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) diff --git a/src/opacity/absorber.cpp b/src/opacity/absorber.cpp index 7ac48443..ab1e3f48 100644 --- a/src/opacity/absorber.cpp +++ b/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(); diff --git a/src/opacity/absorber_ck.cpp b/src/opacity/absorber_ck.cpp index 59191a90..1ee2b106 100644 --- a/src/opacity/absorber_ck.cpp +++ b/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& 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]; - } -} diff --git a/src/opacity/absorber_ck.hpp b/src/opacity/absorber_ck.hpp index 8ec7d189..4ae3488a 100644 --- a/src/opacity/absorber_ck.hpp +++ b/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; diff --git a/src/opacity/helios_ck_premix.cpp b/src/opacity/helios_ck.cpp similarity index 69% rename from src/opacity/helios_ck_premix.cpp rename to src/opacity/helios_ck.cpp index 193ae820..9a4676fd 100644 --- a/src/opacity/helios_ck_premix.cpp +++ b/src/opacity/helios_ck.cpp @@ -8,29 +8,53 @@ #include #include +// application +#include + // climath #include +// harp +#include + // opacity #include "absorber_ck.hpp" +HeliosCK::HeliosCK(std::string name) : AbsorberCK(name) {} + +void HeliosCK::ModifySpectralGrid(std::vector& 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(GetPar("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,7 +86,7 @@ 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; } @@ -70,10 +94,14 @@ void HeliosCK::LoadCoefficient(std::string fname, int bid) { 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; } } diff --git a/src/snap/thermodynamics/thermodynamics.cpp b/src/snap/thermodynamics/thermodynamics.cpp index 5c4555e4..c4bd155f 100644 --- a/src/snap/thermodynamics/thermodynamics.cpp +++ b/src/snap/thermodynamics/thermodynamics.cpp @@ -55,8 +55,8 @@ Thermodynamics* Thermodynamics::fromYAMLInput(YAML::Node const& node) { mythermo_ = new Thermodynamics(); // non-condensable - mythermo_->Rd_ = node["non-condensable"]["Rd"].as(); - mythermo_->gammad_ref_ = node["non-condensable"]["gammad_ref"].as(); + mythermo_->Rd_ = node["dry"]["Rd"].as(); + mythermo_->gammad_ref_ = node["dry"]["gammad"].as(); // saturation adjustment if (node["saturation_adjustment"]) { diff --git a/tools/run_ktable_simple.py b/tools/run_ktable_simple.py deleted file mode 100755 index 0efb630e..00000000 --- a/tools/run_ktable_simple.py +++ /dev/null @@ -1,199 +0,0 @@ -#! /usr/bin/env python3 -import os, subprocess -from canoe import * -from numpy import * - -from pyharp import radiation_band -from utilities import load_file -from collections import OrderedDict - -hitran_file = f"{cmake_source_dir}/data/HITRAN2020.par" - - -def check_file_exist(filename: str) -> bool: - if not os.path.isfile(filename): - print("File %s does not exist!" % filename) - sys.exit(1) - return True - - -# create rfm atmosphere file -def create_rfm_atm(atm: dict) -> None: - print("# Creating rfm.atm ...") - num_layers = atm["HGT"].shape[0] - 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]) - file.write("\n*PRE [mb]\n") - for i in range(num_layers): - file.write("%.8g " % atm["PRE"][i]) - 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]) - file.write("\n*END") - print("# rfm.atm written.") - - -# create rfm driver file -def create_rfm_drv(driver: dict) -> None: - print("# Creating rfm.drv ...") - with open("rfm.drv", "w") as file: - for sec in driver: - if driver[sec] != None: - file.write(sec + "\n") - file.write(" " * 4 + driver[sec] + "\n") - print("# rfm.drv written.") - - -def run_rfm() -> None: - process = subprocess.Popen( - ["./rfm.release"], stdout=subprocess.PIPE, stderr=subprocess.STDOUT - ) - - for line in iter(process.stdout.readline, b""): - # decode the byte string and end='' to avoid double newlines - print(line.decode(), end="") - - process.communicate() - - -def create_netcdf_input( - bane_name: str, - absorbers: list, - atm: dict, - wmin: float, - wmax: float, - wres: float, - tnum: int, - tmin: float, - tmax: float, -) -> str: - print(f"# Creating kcoeff.inp-{band_name} ...") - fname = "kcoeff.inp-%s" % band_name - with open("kcoeff.inp-%s" % band_name, "w") as file: - file.write("# Molecular absorber\n") - file.write("%d\n" % len(absorbers)) - file.write(" ".join(absorbers) + "\n") - file.write("# Molecule data files\n") - for ab in absorbers: - file.write("%-40s\n" % ("./tab_" + ab.lower() + ".txt",)) - file.write("# Wavenumber range\n") - file.write( - "%-14.6g%-14.6g%-14.6g\n" % (wmin, wmax, int((wmax - wmin) / wres) + 1) - ) - file.write("# Relative temperature range\n") - file.write("%-14.6g%-14.6g%-14.6g\n" % (tmin, tmax, tnum)) - file.write("# Number of vertical levels\n") - file.write("%d\n" % len(atm["TEM"])) - file.write("# Temperature\n") - for i in range(len(atm["TEM"])): - file.write("%-14.6g" % atm["TEM"][-(i + 1)]) - if (i + 1) % 10 == 0: - file.write("\n") - if (i + 1) % 10 != 0: - file.write("\n") - file.write("# Pressure\n") - for i in range(len(atm["PRE"])): - file.write("%-14.6g" % atm["PRE"][-(i + 1)]) - if (i + 1) % 10 == 0: - file.write("\n") - if (i + 1) % 10 != 0: - file.write("\n") - print(f"# kcoeff.inp-{band_name} written.") - return fname - - -def run_kcoeff(inpfile: str, ncfile: str) -> None: - process = subprocess.Popen( - ["./kcoeff.release", "-i", inpfile, "-o", ncfile], - stdout=subprocess.PIPE, - stderr=subprocess.STDOUT, - ) - - for line in iter(process.stdout.readline, b""): - # decode the byte string and end='' to avoid double newlines - print(line.decode(), end="") - - process.communicate() - - -# create atmosphere dictionary -def create_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__": - info = load_file("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 - atm = create_atmosphere(num_layers) - - # create rfm atmosphere file - create_rfm_atm(atm) - - # create spectral grid - tem_grid = (5, -20, 20) - wav_grid = (wmin, wmax, wres) - - # create rfm driver file - driver = OrderedDict( - [ - ("*HDR", "Header for rfm"), - ("*FLG", "TAB CTM"), - ("*SPC", "%.4f %.4f %.4f" % wav_grid), - ("*GAS", " ".join(absorbers)), - ("*ATM", "rfm.atm"), - ("*DIM", "PLV \n %d %.4f %.4f" % tem_grid), - ("*TAB", "tab_*.txt"), - ("*HIT", hitran_file), - ("*END", ""), - ] - ) - check_file_exist(hitran_file) - create_rfm_drv(driver) - run_rfm() - inp = create_netcdf_input(band_name, absorbers, atm, *wav_grid, *tem_grid) - run_kcoeff(inp, f"kcoeff-{band_name}.nc")