Update longwave quadrature, LW testing data

lw_rce.py

@@ -0,0 +1,44 @@
+
+# sigma 
+
+def source(T):
+"""Broadband flux from temperature 
+
+    In:
+    T [K]:  temperature
+
+    Returns:
+    source [W/m2-str]: source of broadband intensity (sigma T^4/pi)
+    """ 
+
+def set_up_rad_equil(Tsurf, tau_array, D = 1./0.6096748751): 
+"""Sets up a temperature profile consistent with grey radiative equillibirum given 
+    a surface temperature and an array of optical depth increments.
+    Assumes the top-of-atmosphere is at index 0
+
+    In:
+    Tsurf [K]: surface temperature
+    tau_array []; optical depth increwments 
+
+    Returns:
+    T [k]: profile of temperature consistent with  grey radiative equillibirum
+    """ 
+
+def olr_in_re(Tsurf, total_tau, OLR, D = 1./0.6096748751):
+"""Checks that the top-of-atmosphere outgoing longwave radiation 
+   is consistent with the optcal depth and surface temperature
+
+    In:
+    Tsurf [K]: surface temperature
+    total_tau []; total_optical_depth 
+    OLR [W/m2]: 
+
+    Returns: True/False 
+    """ 
+    (2._wp * sigma * T**4)/(2 + D total_tau) 
+
+   #
+   # Each check requires surface T and an optical depth distribution 
+   #
+   # Compute fluxes; check OLR and constant net flux 
+   #

pyrte_rrtmgp/kernels/rte.py

@@ -9,22 +9,24 @@
     rte_sw_solver_noscat,
 )
 
-GAUSS_DS = np.array(
+# Weights and angle secants for "Gauss-Jacobi-5" quadrature.
+#   Values from Table 1, R. J. Hogan 2023, doi:10.1002/qj.4598
+GAUSS_DS = 1./np.array(
     [
-        [1.66, 0.0, 0.0, 0.0],  # Diffusivity angle, not Gaussian angle
-        [1.18350343, 2.81649655, 0.0, 0.0],
-        [1.09719858, 1.69338507, 4.70941630, 0.0],
-        [1.06056257, 1.38282560, 2.40148179, 7.15513024],
+        [0.6096748751, np.inf,       np.inf,       np.inf],  
+        [0.2509907356, 0.7908473988, np.inf,       np.inf],
+        [0.1024922169, 0.4417960320, 0.8633751621, np.inf],
+        [0.0454586727, 0.2322334416, 0.5740198775, 0.9030775973]
     ]
 )
 
 
 GAUSS_WTS = np.array(
     [
-        [0.5, 0.0, 0.0, 0.0],
-        [0.3180413817, 0.1819586183, 0.0, 0.0],
-        [0.2009319137, 0.2292411064, 0.0698269799, 0.0],
-        [0.1355069134, 0.2034645680, 0.1298475476, 0.0311809710],
+        [1.,           0.0,          0.0,          0.0], 
+        [0.2300253764, 0.7699746236, 0.0,          0.0],
+        [0.0437820218, 0.3875796738, 0.5686383044, 0.0], 
+        [0.0092068785, 0.1285704278, 0.4323381850, 0.4298845087]
     ]
 )
 

tests/test_python_frontend/test_lw_solver.py

@@ -10,6 +10,7 @@
 
 rte_rrtmgp_dir = download_rrtmgp_data()
 clear_sky_example_files = f"{rte_rrtmgp_dir}/examples/rfmip-clear-sky/inputs"
+clear_sky_ref_files = f"{rte_rrtmgp_dir}/examples/rfmip-clear-sky/reference"
 
 rfmip = xr.load_dataset(
     f"{clear_sky_example_files}/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc"
@@ -18,12 +19,12 @@
 kdist = xr.load_dataset(f"{rte_rrtmgp_dir}/rrtmgp-gas-lw-g256.nc")
 
 rlu = xr.load_dataset(
-    "tests/test_python_frontend/rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"
+    f"{clear_sky_ref_files}/rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"
 )
 ref_flux_up = rlu.isel(expt=0)["rlu"].values
 
 rld = xr.load_dataset(
-    "tests/test_python_frontend/rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"
+    f"{clear_sky_ref_files}/rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"
 )
 ref_flux_down = rld.isel(expt=0)["rld"].values