From 2996019d8c16d42821e0f939cf5d4da1788fa021 Mon Sep 17 00:00:00 2001
From: Sriharsha Kandala <sriharsha.kvs@gmail.com>
Date: Sat, 15 Jun 2024 00:01:02 +0530
Subject: [PATCH] Update longwave secants and weights (#498)

---
 src/optics/AngularDiscretizations.jl | 21 ++++++++++-----------
 src/rte/longwave1scalar.jl           |  2 +-
 test/runtests.jl                     |  1 +
 3 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/optics/AngularDiscretizations.jl b/src/optics/AngularDiscretizations.jl
index 8e3e53350..666d752e0 100644
--- a/src/optics/AngularDiscretizations.jl
+++ b/src/optics/AngularDiscretizations.jl
@@ -8,9 +8,8 @@ export AngularDiscretization
 """
     AngularDiscretization{FT,FTA1D}
 
-Weights and angle secants for first order (k=1) Gaussian quadrature.
-Values from Table 2, Clough et al, 1992, doi:10.1029/92JD01419
-after Abramowitz & Stegun 1972, page 921
+Weights and angle secants for "Gauss-Jacobi-5" quadrature.
+Values from Table 1, R. J. Hogan 2023, doi:10.1002/qj.4598
 
 # Fields
 $(DocStringExtensions.FIELDS)
@@ -29,17 +28,17 @@ function AngularDiscretization(::Type{FT}, ::Type{DA}, n_gauss_angles::Int) wher
     max_gauss_pts = 4
     @assert 1 ≤ n_gauss_angles ≤ max_gauss_pts
     if n_gauss_angles == 1
-        gauss_Ds = DA{FT}([1.660000000000])
-        gauss_wts = DA{FT}([0.500000000000])
+        gauss_Ds = DA{FT}(FT(1) ./ [0.6096748751])
+        gauss_wts = DA{FT}([1])
     elseif n_gauss_angles == 2
-        gauss_Ds = DA{FT}(FT[1.183503430000 2.816496550000])
-        gauss_wts = DA{FT}([0.318041381700 0.181958618300])
+        gauss_Ds = DA{FT}(FT(1) ./ [0.2509907356, 0.7908473988])
+        gauss_wts = DA{FT}([0.2300253764, 0.7699746236])
     elseif n_gauss_angles == 3
-        gauss_Ds = DA{FT}([1.097198580000 1.693385070000 4.709416300000])
-        gauss_wts = DA{FT}([0.200931913700 0.229241106400 0.069826979900])
+        gauss_Ds = DA{FT}(FT(1) ./ [0.1024922169, 0.4417960320, 0.8633751621])
+        gauss_wts = DA{FT}([0.0437820218, 0.3875796738, 0.5686383044])
     else
-        gauss_Ds = DA{FT}([1.060562570000 1.382825600000 2.401481790000 7.155130240000])
-        gauss_wts = DA{FT}([0.135506913400 0.203464568000 0.129847547600 0.031180971000])
+        gauss_Ds = DA{FT}(FT(1) ./ [0.0454586727, 0.2322334416, 0.5740198775, 0.9030775973])
+        gauss_wts = DA{FT}([0.0092068785, 0.1285704278, 0.4323381850, 0.4298845087])
     end
 
     return AngularDiscretization{eltype(gauss_Ds), typeof(gauss_Ds)}(
diff --git a/src/rte/longwave1scalar.jl b/src/rte/longwave1scalar.jl
index b0858841e..e31708f3d 100644
--- a/src/rte/longwave1scalar.jl
+++ b/src/rte/longwave1scalar.jl
@@ -133,7 +133,7 @@ Transport for no-scattering longwave problem.
     FT = eltype(τ)
     τ_thresh = 100 * eps(FT) # or abs(eps(FT))?
 
-    intensity_to_flux = FT(2) * FT(π) * w_μ[1]
+    intensity_to_flux = FT(π) * w_μ[1]
     flux_to_intensity = FT(1) / intensity_to_flux
 
     # Transport is for intensity
diff --git a/test/runtests.jl b/test/runtests.jl
index 870b9d013..b0fd53a15 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,7 @@ printstyled("==========================\n\n", color = color1)
     context = ClimaComms.context()
     include("clear_sky_utils.jl")
     for FT in (Float32, Float64)
+        clear_sky(context, OneScalar, TwoStream, NoScatLWRTE, TwoStreamSWRTE, VmrGM, FT)
         clear_sky(context, TwoStream, TwoStream, TwoStreamLWRTE, TwoStreamSWRTE, VmrGM, FT)
     end
 end