28 add more usage examples (#78)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

examples/all_sky_example.ipynb

@@ -0,0 +1,197 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "\n",
+    "from pyrte_rrtmgp.rrtmgp_gas_optics import GasOpticsFiles, load_gas_optics\n",
+    "from pyrte_rrtmgp.rrtmgp_data import download_rrtmgp_data\n",
+    "from pyrte_rrtmgp.rte_solver import RTESolver\n",
+    "from pyrte_rrtmgp.all_skys_funcs import (\n",
+    "    compute_clouds,\n",
+    "    compute_cloud_optics,\n",
+    "    combine_optical_props,\n",
+    ")\n",
+    "from pyrte_rrtmgp.utils import compute_profiles, create_gas_dataset\n",
+    "\n",
+    "# Set up paths and dimensions\n",
+    "rte_rrtmgp_dir = download_rrtmgp_data()\n",
+    "rfmip_dir = os.path.join(rte_rrtmgp_dir, \"examples\", \"all-sky\")\n",
+    "ref_dir = os.path.join(rfmip_dir, \"reference\")\n",
+    "lw_clouds = os.path.join(rte_rrtmgp_dir, \"rrtmgp-clouds-lw-bnd.nc\")\n",
+    "\n",
+    "ncol = 24\n",
+    "nlay = 72\n",
+    "\n",
+    "# Create atmospheric profiles and gas concentrations\n",
+    "profiles = compute_profiles(300, ncol, nlay)\n",
+    "gas_values = {\n",
+    "    \"co2\": 348e-6,\n",
+    "    \"ch4\": 1650e-9,\n",
+    "    \"n2o\": 306e-9,\n",
+    "    \"n2\": 0.7808,\n",
+    "    \"o2\": 0.2095,\n",
+    "    \"co\": 0.0,\n",
+    "}\n",
+    "gases = create_gas_dataset(gas_values, dims={\"site\": ncol, \"layer\": nlay})\n",
+    "\n",
+    "# Set up atmosphere dataset\n",
+    "atmosphere = xr.merge([profiles, gases])\n",
+    "top_at_1 = (\n",
+    "    atmosphere[\"pres_layer\"].values[0, 0] < atmosphere[\"pres_layer\"].values[0, -1]\n",
+    ")\n",
+    "t_sfc = profiles[\"temp_level\"][:, nlay if top_at_1 else 0]\n",
+    "atmosphere[\"surface_temperature\"] = xr.DataArray(t_sfc, dims=[\"site\"])\n",
+    "\n",
+    "# Calculate gas optical properties\n",
+    "gas_optics_lw = load_gas_optics(gas_optics_file=GasOpticsFiles.LW_G256)\n",
+    "clear_sky_optical_props = gas_optics_lw.gas_optics.compute(\n",
+    "    atmosphere, problem_type=\"absorption\", add_to_input=False\n",
+    ")\n",
+    "clear_sky_optical_props[\"surface_emissivity\"] = 0.98\n",
+    "\n",
+    "# Calculate cloud properties and optical properties\n",
+    "cloud_optics = xr.load_dataset(lw_clouds)\n",
+    "cloud_properties = compute_clouds(\n",
+    "    cloud_optics, ncol, nlay, profiles[\"pres_layer\"], profiles[\"temp_layer\"]\n",
+    ")\n",
+    "clouds_optical_props = compute_cloud_optics(cloud_properties, cloud_optics)\n",
+    "\n",
+    "# Combine optical properties and solve RTE\n",
+    "combined_optical_props = combine_optical_props(clouds_optical_props, clear_sky_optical_props)\n",
+    "solver = RTESolver()\n",
+    "fluxes = solver.solve(combined_optical_props, add_to_input=False)\n",
+    "\n",
+    "# Load reference data and verify results\n",
+    "ref_data = xr.load_dataset(\n",
+    "    os.path.join(ref_dir, \"rrtmgp-allsky-lw-no-aerosols.nc\"),\n",
+    "    decode_cf=False,\n",
+    ")\n",
+    "assert np.isclose(fluxes[\"lw_flux_up\"].values, ref_data[\"lw_flux_up\"].values.T, atol=1e-7).all()\n",
+    "assert np.isclose(fluxes[\"lw_flux_down\"].values, ref_data[\"lw_flux_dn\"].values.T, atol=1e-7).all()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "\n",
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "from pyrte_rrtmgp.rrtmgp_gas_optics import GasOpticsFiles, load_gas_optics\n",
+    "from pyrte_rrtmgp.rrtmgp_data import download_rrtmgp_data\n",
+    "from pyrte_rrtmgp.rte_solver import RTESolver\n",
+    "from pyrte_rrtmgp.all_skys_funcs import (\n",
+    "    compute_clouds,\n",
+    "    compute_cloud_optics,\n",
+    "    combine_optical_props,\n",
+    "    delta_scale_optical_props,\n",
+    ")\n",
+    "from pyrte_rrtmgp.utils import compute_profiles, create_gas_dataset\n",
+    "\n",
+    "# Set up paths and dimensions\n",
+    "rte_rrtmgp_dir = download_rrtmgp_data()\n",
+    "rfmip_dir = os.path.join(rte_rrtmgp_dir, \"examples\", \"all-sky\")\n",
+    "ref_dir = os.path.join(rfmip_dir, \"reference\")\n",
+    "sw_clouds = os.path.join(rte_rrtmgp_dir, \"rrtmgp-clouds-sw-bnd.nc\")\n",
+    "\n",
+    "ncol = 24\n",
+    "nlay = 72\n",
+    "\n",
+    "# Create atmospheric profiles and gas concentrations\n",
+    "profiles = compute_profiles(300, ncol, nlay)\n",
+    "gas_values = {\n",
+    "    \"co2\": 348e-6,\n",
+    "    \"ch4\": 1650e-9,\n",
+    "    \"n2o\": 306e-9,\n",
+    "    \"n2\": 0.7808,\n",
+    "    \"o2\": 0.2095,\n",
+    "    \"co\": 0.0,\n",
+    "}\n",
+    "gases = create_gas_dataset(gas_values, dims={\"site\": ncol, \"layer\": nlay})\n",
+    "\n",
+    "# Set up atmosphere dataset\n",
+    "atmosphere = xr.merge([profiles, gases])\n",
+    "top_at_1 = (\n",
+    "    atmosphere[\"pres_layer\"].values[0, 0] < atmosphere[\"pres_layer\"].values[0, -1]\n",
+    ")\n",
+    "\n",
+    "# Add SW-specific surface and angle properties\n",
+    "gas_optics_sw = load_gas_optics(gas_optics_file=GasOpticsFiles.SW_G224)\n",
+    "ngpt = gas_optics_sw.sizes[\"gpt\"]\n",
+    "\n",
+    "atmosphere[\"surface_albedo_dir\"] = xr.DataArray(\n",
+    "    0.06,\n",
+    "    dims=[\"site\", \"gpt\"],\n",
+    "    coords={\"site\": np.arange(ncol), \"gpt\": np.arange(ngpt)},\n",
+    ")\n",
+    "atmosphere[\"surface_albedo_dif\"] = xr.DataArray(\n",
+    "    0.06,\n",
+    "    dims=[\"site\", \"gpt\"],\n",
+    "    coords={\"site\": np.arange(ncol), \"gpt\": np.arange(ngpt)},\n",
+    ")\n",
+    "atmosphere[\"mu0\"] = xr.DataArray(\n",
+    "    0.86,\n",
+    "    dims=[\"site\", \"layer\"],\n",
+    "    coords={\"site\": np.arange(ncol), \"layer\": np.arange(nlay)},\n",
+    ")\n",
+    "\n",
+    "# Calculate gas optical properties\n",
+    "clear_sky_optical_props = gas_optics_sw.gas_optics.compute(\n",
+    "    atmosphere, problem_type=\"two-stream\", add_to_input=False\n",
+    ")\n",
+    "\n",
+    "# Calculate cloud properties and optical properties\n",
+    "cloud_optics = xr.load_dataset(sw_clouds)\n",
+    "cloud_properties = compute_clouds(\n",
+    "    cloud_optics, ncol, nlay, profiles[\"pres_layer\"], profiles[\"temp_layer\"]\n",
+    ")\n",
+    "clouds_optical_props = compute_cloud_optics(cloud_properties, cloud_optics, lw=False)\n",
+    "clouds_optical_props = delta_scale_optical_props(clouds_optical_props)\n",
+    "\n",
+    "# Combine optical properties and solve RTE\n",
+    "combined_optical_props = combine_optical_props(clouds_optical_props, clear_sky_optical_props)\n",
+    "solver = RTESolver()\n",
+    "fluxes = solver.solve(combined_optical_props, add_to_input=False)\n",
+    "\n",
+    "# Load reference data and verify results\n",
+    "ref_data = xr.load_dataset(\n",
+    "    os.path.join(ref_dir, \"rrtmgp-allsky-sw-no-aerosols.nc\"),\n",
+    "    decode_cf=False,\n",
+    ")\n",
+    "assert np.isclose(fluxes[\"sw_flux_up\"].values, ref_data[\"sw_flux_up\"].values.T, atol=1e-7).all()\n",
+    "assert np.isclose(fluxes[\"sw_flux_down\"].values, ref_data[\"sw_flux_dn\"].values.T, atol=1e-7).all()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "pyrte",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

examples/lw_example.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -40,10 +40,10 @@
     "rld = xr.load_dataset(rld_reference, decode_cf=False)\n",
     "\n",
     "assert np.isclose(\n",
-    "    fluxes[\"lw_flux_up_broadband\"], rlu[\"rlu\"], atol=ERROR_TOLERANCE\n",
+    "    fluxes[\"lw_flux_up\"], rlu[\"rlu\"], atol=ERROR_TOLERANCE\n",
     ").all()\n",
     "assert np.isclose(\n",
-    "    fluxes[\"lw_flux_down_broadband\"], rld[\"rld\"], atol=ERROR_TOLERANCE\n",
+    "    fluxes[\"lw_flux_down\"], rld[\"rld\"], atol=ERROR_TOLERANCE\n",
     ").all()"
    ]
   }

examples/sw_example.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [

pybind_interface.cpp

@@ -1960,7 +1960,7 @@ PYBIND11_MODULE(pyrte_rrtmgp, m) {
     [](
         int ncol,
         int nlay,
-        int nbnd,
+        int ngpt,
         py::array_t<Bool>  mask,
         py::array_t<Float> lwp,
         py::array_t<Float> re,
@@ -1974,19 +1974,19 @@ PYBIND11_MODULE(pyrte_rrtmgp, m) {
         py::array_t<Float> taussa,
         py::array_t<Float> taussag
     ) {
-        if (ncol <= 0 || nlay <= 0 || nbnd <= 0 || nsteps <= 0) {
-            throw std::runtime_error("ncol, nlay, nbnd and nsteps must be positive integers");
+        if (ncol <= 0 || nlay <= 0 || ngpt <= 0 || nsteps <= 0) {
+            throw std::runtime_error("ncol, nlay, ngpt and nsteps must be positive integers");
         }
 
         if (mask.size()      != ncol * nlay)        throw std::runtime_error("Invalid size for input array 'mask'");
         if (lwp.size()       != ncol * nlay)        throw std::runtime_error("Invalid size for input array 'lwp'");
         if (re.size()        != ncol * nlay)        throw std::runtime_error("Invalid size for input array 're'");
-        if (tau_table.size() != nsteps * nbnd)      throw std::runtime_error("Invalid size for input array 'tau_table'");
-        if (ssa_table.size() != nsteps * nbnd)      throw std::runtime_error("Invalid size for input array 'ssa_table'");
-        if (asy_table.size() != nsteps * nbnd)      throw std::runtime_error("Invalid size for input array 'asy_table'");
-        if (tau.size()       != ncol * nlay * nbnd) throw std::runtime_error("Invalid size for input array 'tau'");
-        if (taussa.size()    != ncol * nlay * nbnd) throw std::runtime_error("Invalid size for input array 'taussa'");
-        if (taussag.size()   != ncol * nlay * nbnd) throw std::runtime_error("Invalid size for input array 'taussag'");
+        if (tau_table.size() != nsteps * ngpt)      throw std::runtime_error("Invalid size for input array 'tau_table'");
+        if (ssa_table.size() != nsteps * ngpt)      throw std::runtime_error("Invalid size for input array 'ssa_table'");
+        if (asy_table.size() != nsteps * ngpt)      throw std::runtime_error("Invalid size for input array 'asy_table'");
+        if (tau.size()       != ncol * nlay * ngpt) throw std::runtime_error("Invalid size for input array 'tau'");
+        if (taussa.size()    != ncol * nlay * ngpt) throw std::runtime_error("Invalid size for input array 'taussa'");
+        if (taussag.size()   != ncol * nlay * ngpt) throw std::runtime_error("Invalid size for input array 'taussag'");
 
         py::buffer_info buf_mask      = mask.request();
         py::buffer_info buf_lwp       = lwp.request();
@@ -2001,7 +2001,7 @@ PYBIND11_MODULE(pyrte_rrtmgp, m) {
         fortran::rrtmgp_compute_cld_from_table(
             ncol,
             nlay,
-            nbnd,
+            ngpt,
             reinterpret_cast<Bool*>(buf_mask.ptr),
             reinterpret_cast<Float*>(buf_lwp.ptr),
             reinterpret_cast<Float*>(buf_re.ptr),

pyproject.toml

@@ -82,7 +82,7 @@ extend_skip_glob = ["!pyrte_rrtmgp/tests/test_python_frontend/*"]
 
 [tool.flake8]
 max-line-length = 88
-extend-ignore = ["E203", "W503", "D100", "D101", "D102", "D103", "D104", "D107", "E501"]
+extend-ignore = ["E203", "W503", "F401", "D100", "D101", "D102", "D103", "D104", "D107", "E501"]
 exclude = [
     ".git",
     "__pycache__",