Merge pull request #50 from MennoVeerman/main

bringht weights and secants up to date in LW gpu code

Author: Chiil

src_cuda/Rte_lw.cu

@@ -68,19 +68,24 @@ void Rte_lw_gpu::rte_lw(
         const int n_gauss_angles)
 {
     const int max_gauss_pts = 4;
+
+    // Weights and angle secants for "Gauss-Jacobi-5" quadrature.
+    // Values from Table 1, R. J. Hogan 2023, doi:10.1002/qj.4598
     const Array_gpu<Float,2> gauss_Ds(
-            Array<Float,2>({      1.66,         0.,         0.,         0.,
-             1.18350343, 2.81649655,         0.,         0.,
-             1.09719858, 1.69338507, 4.70941630,         0.,
-             1.06056257, 1.38282560, 2.40148179, 7.15513024},
+            Array<Float,2>(
+            { 1./0.6096748751, 0.            , 0.             , 0.,
+              1./0.2509907356, 1/0.7908473988, 0.             , 0.,
+              1./0.1024922169, 1/0.4417960320, 1./0.8633751621, 0.,
+              1./0.0454586727, 1/0.2322334416, 1./0.5740198775, 1./0.903077597 },
             { max_gauss_pts, max_gauss_pts }));
 
-    const Array<Float,2> gauss_wts(
-            {         0.5,           0.,           0.,           0.,
-             0.3180413817, 0.1819586183,           0.,           0.,
-             0.2009319137, 0.2292411064, 0.0698269799,           0.,
-             0.1355069134, 0.2034645680, 0.1298475476, 0.0311809710},
-            { max_gauss_pts, max_gauss_pts });
+    const Array_gpu<Float,2> gauss_wts(
+            Array<Float,2>(
+            { 1.,           0.,           0.,           0.,
+              0.2300253764, 0.7699746236, 0.,           0.,
+              0.0437820218, 0.3875796738, 0.5686383044, 0.,
+              0.0092068785, 0.1285704278, 0.4323381850, 0.4298845087 },
+            { max_gauss_pts, max_gauss_pts }));
 
     const int ncol = optical_props->get_ncol();
     const int nlay = optical_props->get_nlay();

src_kernels_cuda/rte_solver_kernels.cu

@@ -31,7 +31,7 @@ void lw_secants_array_kernel(
     }
 }
 
- 
+
 __device__
 void lw_transport_noscat_kernel(
         const int icol, const int igpt, const int ncol, const int nlay, const int ngpt, const Bool top_at_1,
@@ -186,9 +186,9 @@ void lw_solver_noscat_step_3_kernel(
         const Float pi = acos(Float(-1.));
 
         const int idx = icol + ilev*ncol + igpt*ncol*(nlay+1);
-        radn_up[idx] *= Float(2.) * pi * weight[0];
-        radn_dn[idx] *= Float(2.) * pi * weight[0];
-        radn_up_jac[idx] *= Float(2.) * pi * weight[0];
+        radn_up[idx] *= pi * weight[0];
+        radn_dn[idx] *= pi * weight[0];
+        radn_up_jac[idx] *= pi * weight[0];
     }
 }
 
@@ -585,7 +585,7 @@ void sw_2stream_function(
         *t_dir = -rt_term2 * ((Float(1.) + k_mu) * (alpha1 + k_gamma4) * t_noscat[0]   -
                                   (Float(1.) - k_mu) * (alpha1 - k_gamma4) * exp_minus2ktau * t_noscat[0] -
                                    Float(2.) * (k_gamma4 + alpha1 * k_mu)  * exp_minusktau);
-        
+
         // fix thanks to peter ukkonen (see https://github.com/earth-system-radiation/rte-rrtmgp/pull/39#issuecomment-1026698541)
         *r_dir = max(tmin<Float>(), min(*r_dir, Float(1.0) - *t_noscat));
         *t_dir = max(tmin<Float>(), min(*t_dir, Float(1.0) - *t_noscat - *r_dir));