Skip to article frontmatterSkip to article content

All-sky example (based on RCE)

This notebook demonstrates how to use the PyRTE-RRTMGP package to solve an idealized problem including clouds and clear skies. PyRTE-RRTMGP is a Python implementation of the Radiative Transfer for Energetics (RTE).

Overview

PyRTE-RRTMGP provides a flexible and efficient framework for computing radiative fluxes in planetary atmospheres. This example shows an end-to-end problem with both clear skies and clouds.

  1. Loading data for cloud and gas optics
  2. Computing gas and cloud optical properties and combining them to produce an all-sky problem
  3. Solving the radiative transfer equation to obtain upward and downward fluxes
  4. Validating results against reference solutions generated with the original RTE fortran code

The package leverages xarray to represent data.

Key Components

This example demonstrates the workflow for both longwave and shortwave radiative transfer calculations.

See the documentation for more information.

Preliminaries

Installing dependencies

import numpy as np

Importing pyRTE components

from pyrte_rrtmgp.rrtmgp_data_files import (
    CloudOpticsFiles,
    GasOpticsFiles,
)
from pyrte_rrtmgp.examples import (
    compute_RCE_clouds,
    compute_RCE_profiles,
    load_example_file,
    ALLSKY_EXAMPLES,
)
from pyrte_rrtmgp import rte
from pyrte_rrtmgp.rrtmgp import GasOptics, CloudOptics

Setting up the problem

The routine compute_RCE_profiles() packaged with pyRTE_RRTMGP computes temperature, pressure, and humidity profiles following a moist adibat. The concentrations of other gases are also needed. Clouds are distributed in 2/3 of the columns

def make_profiles(ncol=24, nlay=72):
    # Create atmospheric profiles and gas concentrations
    atmosphere = compute_RCE_profiles(300, ncol, nlay)

    # Add other gas values
    gas_values = {
        "co2": 348e-6,
        "ch4": 1650e-9,
        "n2o": 306e-9,
        "n2": 0.7808,
        "o2": 0.2095,
        "co": 0.0,
    }

    for gas_name, value in gas_values.items():
        atmosphere[gas_name] = value

    return atmosphere

Longwave calculations

In this example datasets are saved to intermediate variables at each step

Initialize the cloud and gas optics data

cloud_optics_lw = CloudOptics(
    cloud_optics_file=CloudOpticsFiles.LW_BND
)
gas_optics_lw = GasOptics(
    gas_optics_file=GasOpticsFiles.LW_G256
)

cloud_optics_lw, gas_optics_lw

Atmospheric profiles - clear sky, then clouds

atmosphere = make_profiles()
#
# Temporary workaround - compute_RCE_clouds() needs to know the particle size;
#   that's set as the mid-point of the valid range from cloud_optics
#
cloud_props = compute_RCE_clouds(
    cloud_optics_lw, 
    atmosphere["pres_layer"], 
    atmosphere["temp_layer"]
)

atmosphere = atmosphere.merge(cloud_props)
atmosphere

Clear-sky (gases) optical properties; surface boundary conditions

optical_props = gas_optics_lw.compute(
    atmosphere, 
    add_to_input=False,
)
optical_props["surface_emissivity"] = 0.98
optical_props

Calculate cloud properties; create all-sky optical properties

First compute the optical properties of the clouds

clouds_optical_props = cloud_optics_lw.compute(
    atmosphere,
    problem_type = rte.OpticsTypes.ABSORPTION,
)
# The optical properties of the clouds alone
clouds_optical_props

Then add the optical properties of the clouds to the clear sky to get the total

# add_to() changes the value of `optical_props`
clouds_optical_props.rte.add_to(optical_props)
optical_props

Compute broadband fluxes

fluxes = optical_props.rte.solve(add_to_input=False)
fluxes

Compare to reference results

# Load reference data and verify results
ref_data = load_example_file(ALLSKY_EXAMPLES.REF_LW_NO_AEROSOL)
assert np.isclose(
    fluxes["lw_flux_up"],
    ref_data["lw_flux_up"].T,
    atol=1e-7,
).all()

assert np.isclose(
    fluxes["lw_flux_down"],
    ref_data["lw_flux_dn"].T,
    atol=1e-7,
).all()

print("All-sky longwave calculations validated successfully")

Shortwave calculations

In this example steps are combined where possible

Initialize optics data

cloud_optics_sw = CloudOptics(
    cloud_optics_file=CloudOpticsFiles.SW_BND
)
gas_optics_sw = GasOptics(
    gas_optics_file=GasOpticsFiles.SW_G224
)

cloud_optics_sw, gas_optics_sw

Atmospheric profiles - clear sky, then clouds

atmosphere = make_profiles()
#
# Temporary workaround - compute_RCE_clouds() needs to know the particle size;
#    that's set as the mid-point of the valid range from cloud_optics
#
atmosphere = atmosphere.merge(
    compute_RCE_clouds(
        cloud_optics_sw, atmosphere["pres_layer"], atmosphere["temp_layer"]
    )
)
atmosphere

Compute gas and cloud optics and combine in one step

# compute_cloud_optics() returns two-stream properties by default?
optical_props = gas_optics_sw.compute(
    atmosphere, 
    add_to_input=False,
)
# add_to() changes the values in optical_props
cloud_optics_sw.compute(atmosphere).rte.add_to(
    optical_props, 
    delta_scale=True,
)
#
# Add SW-specific surface and angle properties
#
optical_props["surface_albedo_direct"] = 0.06
optical_props["surface_albedo_diffuse"] = 0.06
# Could also specific a single "surface_albedo"
optical_props["mu0"] = 0.86
optical_props

Compute fluxes

fluxes = optical_props.rte.solve(add_to_input=False)
fluxes

Compare to reference results

The fluxes computed here have a column dimension while the reference fluxes have a site dimension. But np.close() is happy to compare the arrays since they are the same size.

ref_data = load_example_file(ALLSKY_EXAMPLES.REF_SW_NO_AEROSOL)
assert np.isclose(
    fluxes["sw_flux_up"],
    ref_data["sw_flux_up"].T,
    atol=1e-7,
).all()
assert np.isclose(
    fluxes["sw_flux_down"],
    ref_data["sw_flux_dn"].T,
    atol=1e-7,
).all()
assert np.isclose(
    fluxes["sw_flux_dir"],
    ref_data["sw_flux_dir"].T,
    atol=1e-7,
).all()

print("All-sky shortwave calculations validated successfully")
pyRTE user documentation
Clear sky example (RFMIP)