Closed RobertPincus closed 3 months ago
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")
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).
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 andvariant
is the new code (weights, angles).