Updated and simplified LW transport

[RobertPincus]

Introduces two changes to long wave noscattering radiative transfer calculations.

1. Changes first-order Gaussian angular quadrature to "Gauss-Jacobi-5" based on results from Hogan 2023.
2. Removes Planck function evaluated at layer centers, reducing the number of large arrays used in the LW solver by 1.

Results for single-angle calculations are comparable to or slightly more accurate than the original formulation (see below).

Changes kernel API. Reference data for continuous integration is updated and failure thresholds are set to their original, stricter values.

[RobertPincus]

In the first three figures dashed line is mean absolute error and solid line RMS as computed over 100 profiles from the present-day (the RFMIP examples) relative to LBLRTM (which uses Planck functions at both layers and levels). The last plot shows the max error at any level.

Mmm, omitting the level source does increase errors relative to LBLRTM calculations that do include the source.

[RobertPincus]

For the record the script to make the above plots:

#! /usr/bin/env python
#
#
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import os
import urllib.request
import warnings

fluxes = ["rld", "rlu"]

def mae(diff, col_dim):
    #
    # Mean absolute error
    #
    return np.fabs(diff).mean(dim=col_dim)


def rms(diff, col_dim):
    #
    # Root mean square error
    #
    return np.sqrt(np.square(diff).mean(dim=col_dim))

def construct_lbl_esgf_root(var, esgf_node="llnl"):
    #
    # For a given variable name, provide the https URL for the LBLRM
    # line-by-line results
    #
    model = "LBLRTM-12-8"
    prefix = ("http://esgf3.dkrz.de/thredds/fileServer/cmip6/RFMIP/AER/" +
              model + "/rad-irf/r1i1p1f1/Efx/")
    if esgf_node == "llnl":
        prefix = ("http://esgf-data1.llnl.gov/thredds/fileServer/css03_data/"
                  "CMIP6/RFMIP/AER/" + model + "/rad-irf/r1i1p1f1/Efx/")
    return (prefix + var + "/gn/v20190514/" + var +
            "_Efx_" + model + "_rad-irf_r1i1p1f1_gn.nc")


######################
#######################

#
# Comparing reference and test results
#
if __name__ == '__main__':
    warnings.simplefilter("ignore", xr.SerializationWarning)

    lbl_suffix = "_Efx_LBLRTM-12-8_rad-irf_r1i1p1f1_gn.nc"
    gp_suffix  = "_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"

    for v in fluxes:
        try:
            try:
                urllib.request.urlretrieve(construct_lbl_esgf_root(v),
                                           v + lbl_suffix)
            except:
                urllib.request.urlretrieve(
                    construct_lbl_esgf_root(v, esgf_node="dkrz"),
                    v + lbl_suffix)
        except:
            raise Exception("Failed to download {0}".format(v + lbl_suffix))

    lbl = xr.open_mfdataset([                                                                    v + lbl_suffix for v in fluxes], 
        combine="by_coords").sel(expt=0)
    old = xr.open_mfdataset([os.environ["RRTMGP_DATA"] + "/examples/rfmip-clear-sky/reference/"+ v + gp_suffix  for v in fluxes], 
        combine="by_coords").sel(expt=0)
    new = xr.open_mfdataset([                                      "examples/rfmip-clear-sky/" + v + gp_suffix  for v in fluxes], 
        combine="by_coords").sel(expt=0)

    #
    # Mean sbolute and RMS errors in upwelling flux
    #
    fig, axes = plt.subplots(ncols=2, figsize=[1.6 * 6, 6], sharey=True)

    for a, flux in enumerate(fluxes): 
        axes[a].plot(rms((old[flux]-lbl[flux]), col_dim="site"), lbl.level, c = "blue", label="Clough")
        axes[a].plot(rms((new[flux]-lbl[flux]), col_dim="site"), lbl.level, c = "red",  label="G-J-5")

        axes[a].plot(mae((old[flux]-lbl[flux]), col_dim="site"), lbl.level, "--", c = "blue")
        axes[a].plot(mae((new[flux]-lbl[flux]), col_dim="site"), lbl.level, "--", c = "red")
        axes[a].invert_yaxis()

        axes[a].set_title(flux)
        if(a == 1): axes[a].legend(frameon=False)
    
    fig.savefig("Flux-errors.png")

    #
    # Errors in net flux 
    #
    old_err_net = (old.rld - old.rlu) - (lbl.rld - lbl.rlu)
    new_err_net = (new.rld - new.rlu) - (lbl.rld - lbl.rlu)

    fig, axes = plt.subplots(ncols=2, figsize=[1.6 * 6, 6], sharey=True)
    axes[0].plot(rms(old_err_net, col_dim="site"), lbl.level, c = "blue", label="Clough")
    axes[0].plot(rms(new_err_net, col_dim="site"), lbl.level, c = "red", label="G-J-5")

    axes[0].plot(mae(old_err_net, col_dim="site"), lbl.level, "--", c = "blue")
    axes[0].plot(mae(new_err_net, col_dim="site"), lbl.level, "--", c = "red")

    axes[0].legend(frameon=False)
    axes[0].set_title("Net flux")

    #
    # Max error
    #
    axes[1].plot(np.abs(old.rld-lbl.rld).max(dim="site"), lbl.level, c = "blue")
    axes[1].plot(np.abs(new.rld-lbl.rld).max(dim="site"), lbl.level, c = "red")
    axes[1].plot(np.abs(old.rlu-lbl.rlu).max(dim="site"), lbl.level, c = "blue")
    axes[1].plot(np.abs(new.rlu-lbl.rlu).max(dim="site"), lbl.level, c = "red")

    axes[1].set_title("Max. abs. error")
    fig.savefig("Net-and-max-flux-error.png")

[RobertPincus]

#282 introduced the first of these changes. I'll come back to using only level source functions in the future when I'm in better state to compare cost/accuracy tradeoffs.