Change LW quadrature angles

[RobertPincus]

This updates the weight and secants for LW quadrature to the "Gauss-Jacobi-5" described in R. J. Hogan 2023, doi:10.1002/qj.4598

The data repository has been updated to match so tolerances for comparisons against examples have also been updated. The names for the rfmip-clear-sky files have not changed but perhaps they should since the data repository is no longer consistent with the RTE+RRTMGP contribution to CMIP6.

The comparison with LBLRTM, which uses the older integration angles and weights, is not uniformly better. In these figures dashed lines are mean absolute error and solid lines RMS; "baseline" is the current version of develop for both code and data and variant is the new code (weights, angles).

[RobertPincus]

Figures were made with this script:

#! /usr/bin/env python
#
# General purpose comparison script -- compare all variables in a set of files, 
#   write output if differences exceed some threshold, 
#   error exit if differences exceed a different threshold 
#
# Currently thresholds are specified as absolute differences 
#    worth revising but could change development practice 
# Thresholds come from environement variables when set? 
# 
#
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 absolute 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="Baseline")
        axes[a].plot(rms((new[flux]-lbl[flux]), col_dim="site"), lbl.level, c = "red",  label="Variant")

        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].set_title(flux)
    
    axes[1].legend(frameon=False)
    axes[1].invert_yaxis()

    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="Baseline")
    axes[0].plot(rms(new_err_net, col_dim="site"), lbl.level, c = "red", label="Variant")

    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")
    axes[1].invert_yaxis()

    fig.savefig("Net-and-max-flux-error.png")

[Chiil]

Edits look good. The results are not obviously better judging from the figure, but the tests from Hogan's paper are convincing. Maybe add some comments in the code on why the multiplications in the weights go from 2._wp * pi to pi.

[RobertPincus]

To address @Chiil's comment: the weights as provided by Hogan are normalized to 1; the weights from Clough et al. 1922 in use prior were normalized to 1/2 (flux from isotropic intensity).