fix indexing error to close LW balance

include_test/radiation_solver_rt.h

@@ -27,9 +27,9 @@
 #include "aerosol_optics_rt.h"
 #include "rte_lw_rt.h"
 #include "rte_sw_rt.h"
-#include "raytracer.h"
+#include "raytracer_sw.h"
 #include "raytracer_lw.h"
-#include "raytracer_kernels.h"
+#include "raytracer_kernels_sw.h"
 #include "source_functions_rt.h"
 #include <curand_kernel.h>
 

les_cloudfield/README.md

@@ -0,0 +1,11 @@
+# Example LES cloudfield
+Example cloudfield from a 2560 x 2560 x 4000 m3 (dx=dy=dz=20m) large-eddy simulations based on the Rain In Cumulus over the Ocean (RICO) campaign (Van Zanten, 2011).
+Input and reference output files are stored on Zenodo: https://doi.org/10.5281/zenodo.18757088
+
+How to run:
+
+1. `Compile code (following the basic instructions) with -DUSECUDA`
+2. `ln -s {BUILD_DIRECTORY}/{EXECUTABLE_NAME} test_rte_rrtmgp_*` (link the executables)
+3. `./make_links.sh`                                             (link the coefficients)
+4. `download input and reference output files from https://doi.org/10.5281/zenodo.18757088`
+5. `run "./test_rte_rrtmgp_rt test" (forward) or "test_rte_rrtmgp_bw test" (backward)`

les_cloudfield/make_links.sh

@@ -0,0 +1,9 @@
+#! /bin/sh
+ln -sf ../rrtmgp-data/rrtmgp-clouds-sw.nc cloud_coefficients_sw.nc
+ln -sf ../rrtmgp-data/rrtmgp-clouds-lw.nc cloud_coefficients_lw.nc
+ln -sf ../rrtmgp-data/rrtmgp-gas-sw-g112.nc coefficients_sw.nc
+ln -sf ../rrtmgp-data/rrtmgp-gas-lw-g128.nc coefficients_lw.nc
+ln -sf ../data/aerosol_optics_sw.nc
+ln -sf ../data/aerosol_optics_lw.nc
+
+echo "Don't forget to link the raytracering executables 'test_rte_rrtmgp_rt' and 'test_rte_rrtmgp_bw'"

les_cloudfield/test.ini

@@ -0,0 +1,77 @@
+[switches]
+##### General switches
+shortwave = true
+longwave = true
+fluxes = true
+
+cloud-optics = true
+delta-cloud = false
+
+aerosol-optics = true
+delta-aerosol = false
+
+##### test_rte_rrtmgp only (CPU & GPU)
+output-optical = false
+output-bnd-fluxes = false
+
+##### test_rte_rrtmgp only (GPU)
+timings = false
+
+##### Forward raytracer only
+
+# switch on ray tracer. Specifiy ray counts under [ints]: samples per pixel for sw, for lw provide n to use 2^n samples
+sw-raytracing = true
+lw-raytracing = true
+
+# Also perform two-stream computations for shortwave.
+sw-two-stream = false
+
+# Independent column mode
+independent-column = false
+# Tilted independent column mode
+tica = false
+
+# Enable scattering for longwave. If enabled, 1D solver will also switch from no-scat solution to two-stream formulation
+lw-scattering = true
+
+# Only perform ray tracing when the minimum gasous mean free path is larger than (min-mfp-grid-ratio * grid spacing) and use a 1D solutions otherwise.
+# Specify threshold ratio under [floats]
+min-mfp-grid-ratio = true
+
+# liq and ice flags are automatically switch on when cloud-optics is set to true
+# to run either liquid or ice cloud optics
+liq-cloud-optics = false
+ice-cloud-optics = false
+
+##### Common raytracer (forward & backward) switches
+cloud-mie = false
+
+# Save optical properties and fluxes for a single g-point. Specify under [ints]
+single-gpt = false
+
+# Enable profiling
+profiling = false
+
+##### Backward raytracer only
+raytracing = true
+
+# Read variable "surface_type" to determine scattering type and spectral albedo for each surface grid cell.
+# Set either to 0 (water; spectrally constant albedo; specular reflection) or a number between 1 (soil-like spectral albedo; lambertian) and 2 (grass-line spectral albedo; lambertian).
+lu-albedo = true
+
+# broadband radiance computation (sum over all g-points)
+broadband = true
+
+# image computation: run only g-points in visible spectrum and convert to XYZ tristimulus values
+image = true
+
+[ints]
+single-gpt = 1
+
+sw-raytracing = 256
+lw-raytracing = 22
+
+raytracing = 128
+
+[floats]
+min-mfp-grid-ratio = 1.0

rcemip/make_links.sh

@@ -3,5 +3,5 @@ ln -sf ../rrtmgp-data/rrtmgp-clouds-sw.nc cloud_coefficients_sw.nc
 ln -sf ../rrtmgp-data/rrtmgp-clouds-lw.nc cloud_coefficients_lw.nc
 ln -sf ../rrtmgp-data/rrtmgp-gas-sw-g224.nc coefficients_sw.nc
 ln -sf ../rrtmgp-data/rrtmgp-gas-lw-g256.nc coefficients_lw.nc
-ln -sf ../data/aerosol_optics.nc
+ln -sf ../data/aerosol_optics_sw.nc
 echo "Don't forget to link your executable as 'test_rte_rrtmgp'"

src_cuda_rt/raytracer_lw.cu

@@ -356,17 +356,4 @@ void Raytracer_lw::trace_rays(
 }
 
 Raytracer_lw::Raytracer_lw()
-{
-//    curandDirectionVectors32_t* qrng_vectors;
-//    curandGetDirectionVectors32(
-//                &qrng_vectors,
-//                CURAND_SCRAMBLED_DIRECTION_VECTORS_32_JOEKUO6);
-//    unsigned int* qrng_constants;
-//    curandGetScrambleConstants32(&qrng_constants);
-//
-//    this->qrng_vectors_gpu = allocate_gpu<curandDirectionVectors32_t>(2);
-//    this->qrng_constants_gpu = allocate_gpu<unsigned int>(2);
-//
-//    copy_to_gpu(qrng_vectors_gpu, qrng_vectors, 2);
-//    copy_to_gpu(qrng_constants_gpu, qrng_constants, 2);
-}
+{}

src_cuda_rt/raytracer.cu → src_cuda_rt/raytracer_sw.cu

@@ -1,12 +1,12 @@
 #include <curand_kernel.h>
 
-#include "raytracer.h"
+#include "raytracer_sw.h"
 #include "array.h"
 #include "optical_props_rt.h"
 #include <cub/cub.cuh>
 #include "raytracer_definitions.h"
 #include "raytracer_functions.h"
-#include "raytracer_kernels.h"
+#include "raytracer_kernels_sw.h"
 
 #include "gas_optics_rrtmgp_kernels_cuda_rt.h"
 

src_kernels_cuda_rt/raytracer_kernels_lw.cu

@@ -11,68 +11,6 @@ namespace
 
     constexpr Float w_thres = 0.5;
 
-    struct Quasi_random_number_generator_2d
-    {
-        __device__ Quasi_random_number_generator_2d(
-                curandDirectionVectors32_t* vectors, unsigned int* constants, unsigned int offset)
-        {
-            curand_init(vectors[0], constants[0], offset, &state_x);
-            curand_init(vectors[1], constants[1], offset, &state_y);
-        }
-
-        __device__ void xy(unsigned int* x, unsigned int* y,
-                           const Vector<int>& grid_cells,
-                           const Int qrng_grid_x, const Int qrng_grid_y,
-                           Int& photons_shot)
-        {
-            *x = curand(&state_x);
-            *y = curand(&state_y);
-
-            while (true)
-            {
-                const int i = *x / static_cast<unsigned int>((1ULL << 32) / qrng_grid_x);
-                const int j = *y / static_cast<unsigned int>((1ULL << 32) / qrng_grid_y);
-
-                ++photons_shot;
-                if (i < grid_cells.x && j < grid_cells.y)
-                {
-                    return;
-                }
-                else
-                {
-                    *x = curand(&state_x);
-                    *y = curand(&state_y);
-                }
-            }
-        }
-
-        curandStateScrambledSobol32_t state_x;
-        curandStateScrambledSobol32_t state_y;
-    };
-
-    __device__
-
-    inline int find_source_index(const Float* weights, int n, const Float r)
-    {
-        int left = 0;
-        int right = n;
-
-        while (left < right)
-        {
-            int mid = left + (right - left) / 2;
-            if (weights[mid] <= r)
-            {
-                left = mid + 1;
-            }
-            else
-            {
-                right = mid;
-            }
-        }
-
-        return min(left, n-1);
-    }
-
     __device__
     inline void write_emission(
             Photon& photon,
@@ -115,6 +53,8 @@ namespace
             Float& total_absorbed_weight,
             int& src_type)
     {
+        __syncwarp();
+
         ++photons_shot;
 
         if (photons_shot < photons_to_shoot)
@@ -131,12 +71,13 @@ namespace
 
             const int ij = i + j * grid_cells.x;
 
+            photon.position.x = (i + rng()) * grid_d.x;
+            photon.position.y = (j + rng()) * grid_d.y;
+
             if (k == 0) // surface
             {
                 src_type = 1;
 
-                photon.position.x = (i + rng()) * grid_d.x;
-                photon.position.y = (j + rng()) * grid_d.y;
                 photon.position.z = Float(0.);
 
                 mu = sqrt(rng());
@@ -147,8 +88,6 @@ namespace
             else if (k == grid_cells.z+1) // top-of-domain
             {
                 src_type = 2;
-                photon.position.x = (i + rng()) * grid_d.x;
-                photon.position.y = (j + rng()) * grid_d.y;
                 photon.position.z = grid_size.z - Float_epsilon;
 
                 mu = Float(-1.)*sqrt(rng());
@@ -162,8 +101,6 @@ namespace
                 src_type = 0;
 
                 const int km = k - 1;
-                photon.position.x = (i + rng()) * grid_d.x;
-                photon.position.y = (j + rng()) * grid_d.y;
                 photon.position.z = (km + rng()) * grid_d.z;
 
                 mu = rng()*Float(2.) - Float(1.);
@@ -330,10 +267,11 @@ void ray_tracer_lw_kernel(
 
                 const Float mu_surface = sqrt(rng());
                 const Float azimuth_surface = Float(2.*M_PI)*rng();
+                const Float sin_theta = sqrt(Float(1.) - mu_surface*mu_surface + Float_epsilon);
 
-                photon.direction.x = mu_surface*sin(azimuth_surface);
-                photon.direction.y = mu_surface*cos(azimuth_surface);
-                photon.direction.z = sqrt(Float(1.) - mu_surface*mu_surface + Float_epsilon);
+                photon.direction.x = sin_theta*sin(azimuth_surface);
+                photon.direction.y = sin_theta*cos(azimuth_surface);
+                photon.direction.z = mu_surface;
             }
 
             // TOD exit

src_kernels_cuda_rt/raytracer_kernels.cu → src_kernels_cuda_rt/raytracer_kernels_sw.cu

@@ -63,6 +63,7 @@ namespace
             const Vector<Float> sun_direction,
             Float& weight)
     {
+        __syncwarp();
         unsigned int random_number_x;
         unsigned int random_number_y;
         qrng.xy(&random_number_x, &random_number_y, grid_cells, qrng_grid_x, qrng_grid_y, photons_shot);

src_test/radiation_solver.cpp

@@ -335,26 +335,26 @@ namespace
         Netcdf_file coef_nc(coef_file, Netcdf_mode::Read);
 
         // Read look-up table coefficient dimensions
-        int n_band     = coef_nc.get_dimension_size("band_sw");
+        int n_band     = coef_nc.get_dimension_size("band");
         int n_hum      = coef_nc.get_dimension_size("relative_humidity");
         int n_philic = coef_nc.get_dimension_size("hydrophilic");
         int n_phobic = coef_nc.get_dimension_size("hydrophobic");
 
         Array<Float,2> band_lims_wvn({2, n_band});
 
         Array<Float,2> mext_phobic(
-                coef_nc.get_variable<Float>("mass_ext_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("mass_ext_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
         Array<Float,2> ssa_phobic(
-                coef_nc.get_variable<Float>("ssa_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("ssa_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
         Array<Float,2> g_phobic(
-                coef_nc.get_variable<Float>("asymmetry_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("asymmetry_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
 
         Array<Float,3> mext_philic(
-                coef_nc.get_variable<Float>("mass_ext_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("mass_ext_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
         Array<Float,3> ssa_philic(
-                coef_nc.get_variable<Float>("ssa_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("ssa_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
         Array<Float,3> g_philic(
-                coef_nc.get_variable<Float>("asymmetry_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("asymmetry_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
 
         Array<Float,1> rh_upper(
                 coef_nc.get_variable<Float>("relative_humidity2", {n_hum}), {n_hum});

src_test/radiation_solver.cu

@@ -370,26 +370,26 @@ namespace
         Netcdf_file coef_nc(coef_file, Netcdf_mode::Read);
 
         // Read look-up table coefficient dimensions
-        int n_band     = coef_nc.get_dimension_size("band_sw");
+        int n_band     = coef_nc.get_dimension_size("band");
         int n_hum      = coef_nc.get_dimension_size("relative_humidity");
         int n_philic = coef_nc.get_dimension_size("hydrophilic");
         int n_phobic = coef_nc.get_dimension_size("hydrophobic");
 
         Array<Float,2> band_lims_wvn({2, n_band});
 
         Array<Float,2> mext_phobic(
-                coef_nc.get_variable<Float>("mass_ext_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("mass_ext_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
         Array<Float,2> ssa_phobic(
-                coef_nc.get_variable<Float>("ssa_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("ssa_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
         Array<Float,2> g_phobic(
-                coef_nc.get_variable<Float>("asymmetry_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("asymmetry_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
 
         Array<Float,3> mext_philic(
-                coef_nc.get_variable<Float>("mass_ext_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("mass_ext_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
         Array<Float,3> ssa_philic(
-                coef_nc.get_variable<Float>("ssa_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("ssa_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
         Array<Float,3> g_philic(
-                coef_nc.get_variable<Float>("asymmetry_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("asymmetry_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
 
         Array<Float,1> rh_upper(
                 coef_nc.get_variable<Float>("relative_humidity2", {n_hum}), {n_hum});

src_test/radiation_solver_bw.cu

@@ -445,26 +445,26 @@ namespace
         Netcdf_file coef_nc(coef_file, Netcdf_mode::Read);
 
         // Read look-up table coefficient dimensions
-        int n_band     = coef_nc.get_dimension_size("band_sw");
+        int n_band     = coef_nc.get_dimension_size("band");
         int n_hum      = coef_nc.get_dimension_size("relative_humidity");
         int n_philic = coef_nc.get_dimension_size("hydrophilic");
         int n_phobic = coef_nc.get_dimension_size("hydrophobic");
 
         Array<Float,2> band_lims_wvn({2, n_band});
 
         Array<Float,2> mext_phobic(
-                coef_nc.get_variable<Float>("mass_ext_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("mass_ext_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
         Array<Float,2> ssa_phobic(
-                coef_nc.get_variable<Float>("ssa_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("ssa_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
         Array<Float,2> g_phobic(
-                coef_nc.get_variable<Float>("asymmetry_sw_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
+                coef_nc.get_variable<Float>("asymmetry_hydrophobic", {n_phobic, n_band}), {n_band, n_phobic});
 
         Array<Float,3> mext_philic(
-                coef_nc.get_variable<Float>("mass_ext_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("mass_ext_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
         Array<Float,3> ssa_philic(
-                coef_nc.get_variable<Float>("ssa_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("ssa_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
         Array<Float,3> g_philic(
-                coef_nc.get_variable<Float>("asymmetry_sw_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
+                coef_nc.get_variable<Float>("asymmetry_hydrophilic", {n_philic, n_hum, n_band}), {n_band, n_hum, n_philic});
 
         Array<Float,1> rh_upper(
                 coef_nc.get_variable<Float>("relative_humidity2", {n_hum}), {n_hum});

src_test/radiation_solver_rt.cu

@@ -99,7 +99,7 @@ namespace
         const int grid_col = ncol/block_col + (ncol%block_col > 0);
 
         convert_1d_to_rt_flx_kernels<<<grid_col, block_col>>>(ncol, 0, flux_dn.ptr(), flux_up.ptr(), flux_sfc_dn.ptr(), flux_sfc_up.ptr());
-        convert_1d_to_rt_flx_kernels<<<grid_col, block_col>>>(ncol, nz+1, flux_dn.ptr(), flux_up.ptr(), flux_tod_dn.ptr(), flux_tod_up.ptr());
+        convert_1d_to_rt_flx_kernels<<<grid_col, block_col>>>(ncol, nz, flux_dn.ptr(), flux_up.ptr(), flux_tod_dn.ptr(), flux_tod_up.ptr());
 
         const dim3 block_2d(block_col, 1, 1);
         const dim3 grid_2d(grid_col, nz, 1);