with --cloud-cam enabled, also compute slanted cloud optical depth an…

…d distance to closest cloudy cell

include_rt/Raytracer_bw.h

@@ -101,9 +101,11 @@ class Raytracer_bw
                 const Vector<int>& grid_cells,
                 const Array_gpu<Float,2>& lwp,
                 const Array_gpu<Float,2>& iwp,
+                const Array_gpu<Float,2>& tau_cloud,
                 const Camera& camera,
-                Array_gpu<Float,2>& lwp_cam,
-                Array_gpu<Float,2>& iwp_cam);
+                Array_gpu<Float,2>& liwp_cam,
+                Array_gpu<Float,2>& tauc_cam,
+                Array_gpu<Float,2>& dist_cam);
 
     private:
 

include_rt_kernels/raytracer_kernels_bw.h

@@ -121,11 +121,13 @@ __global__
 void accumulate_clouds_kernel(
     const Float* __restrict__ lwp, 
     const Float* __restrict__ iwp, 
+    const Float* __restrict__ tau_cloud,
     const Vector<Float> grid_d,
     const Vector<Float> grid_size,
     const Vector<int> grid_cells,
-    Float* __restrict__ lwp_cam, 
-    Float* __restrict__ iwp_cam, 
+    Float* __restrict__ liwp_cam, 
+    Float* __restrict__ tauc_cam, 
+    Float* __restrict__ dist_cam, 
     const Camera camera);
 
 #endif

include_test/Radiation_solver_bw.h

@@ -178,8 +178,9 @@ class Radiation_solver_shortwave
                 const Aerosol_concs_gpu& aerosol_concs,
                 const Camera& camera,
                 Array_gpu<Float,3>& XYZ,
-                Array_gpu<Float,2>& lwp_cam,
-                Array_gpu<Float,2>& iwp_cam);
+                Array_gpu<Float,2>& liwp_cam,
+                Array_gpu<Float,2>& tauc_cam,
+                Array_gpu<Float,2>& dist_cam);
         #endif
 
         #ifdef __CUDACC__
@@ -211,8 +212,9 @@ class Radiation_solver_shortwave
                 const Aerosol_concs_gpu& aerosol_concs,
                 const Camera& camera,
                 Array_gpu<Float,2>& radiance,
-                Array_gpu<Float,2>& lwp_cam,
-                Array_gpu<Float,2>& iwp_cam);
+                Array_gpu<Float,2>& liwp_cam,
+                Array_gpu<Float,2>& tauc_cam,
+                Array_gpu<Float,2>& dist_cam);
 
         int get_n_gpt_gpu() const { return this->kdist_gpu->get_ngpt(); };
         int get_n_bnd_gpu() const { return this->kdist_gpu->get_nband(); };

src_cuda_rt/Raytracer_bw.cu

@@ -665,12 +665,15 @@ void Raytracer_bw::accumulate_clouds(
         const Vector<int>& grid_cells,
         const Array_gpu<Float,2>& lwp,
         const Array_gpu<Float,2>& iwp,
+        const Array_gpu<Float,2>& tau_cloud,
         const Camera& camera,
-        Array_gpu<Float,2>& lwp_cam,
-        Array_gpu<Float,2>& iwp_cam)
+        Array_gpu<Float,2>& liwp_cam,
+        Array_gpu<Float,2>& tauc_cam,
+        Array_gpu<Float,2>& dist_cam)
 {
-    Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, lwp_cam.ptr());
-    Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, iwp_cam.ptr());
+    Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, liwp_cam.ptr());
+    Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, tauc_cam.ptr());
+    Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, dist_cam.ptr());
 
     const int n_pix = camera.nx * camera.ny;
     const int n_block = std::min(n_pix, 512);
@@ -685,11 +688,13 @@ void Raytracer_bw::accumulate_clouds(
     accumulate_clouds_kernel<<<grid, block>>>(
         lwp.ptr(),
         iwp.ptr(),
+        tau_cloud.ptr(),
         grid_d,
         grid_size,
         grid_cells,
-        lwp_cam.ptr(),
-        iwp_cam.ptr(),
+        liwp_cam.ptr(),
+        tauc_cam.ptr(),
+        dist_cam.ptr(),
         camera);
 
 }

src_kernels_cuda_rt/raytracer_kernels_bw.cu

@@ -742,11 +742,13 @@ void ray_tracer_kernel_bw(
     void accumulate_clouds_kernel(
             const Float* __restrict__ lwp, 
             const Float* __restrict__ iwp, 
+            const Float* __restrict__ tau_cloud, 
             const Vector<Float> grid_d,
             const Vector<Float> grid_size,
             const Vector<int> grid_cells,
-            Float* __restrict__ lwp_cam, 
-            Float* __restrict__ iwp_cam, 
+            Float* __restrict__ liwp_cam, 
+            Float* __restrict__ tauc_cam, 
+            Float* __restrict__ dist_cam, 
             const Camera camera)
     {
         const int pix = blockDim.x * blockIdx.x + threadIdx.x;
@@ -756,8 +758,10 @@ void ray_tracer_kernel_bw(
 
         if (pix < camera.nx * camera.ny)
         {
-            Float lwp_sum = 0;
-            Float iwp_sum = 0;
+            Float liwp_sum = 0;
+            Float tauc_sum = 0;
+            Float dist = 0;
+            bool reached_cloud = false;
 
             //const int pix = n * pixels_per_thread + ipix;
             const Float i = (Float(pix % camera.nx) + Float(0.5))/ Float(camera.nx);
@@ -807,9 +811,14 @@ void ray_tracer_kernel_bw(
                 const Float sz = abs((direction.z > 0) ? ((k+1) * grid_d.z - position.z)/direction.z : (k*grid_d.z - position.z)/direction.z);
                 const Float s_min = min(sx, min(sy, sz));
 
-                lwp_sum += s_min * lwp[ijk];
-                iwp_sum += s_min * iwp[ijk];
-
+                liwp_sum += s_min * (lwp[ijk] + iwp[ijk]);
+                tauc_sum += s_min * tau_cloud[ijk];
+                if (!reached_cloud)
+                {
+                    dist += s_min;
+                    reached_cloud = tau_cloud[ijk] > 0;
+                }
+
                 position = position + direction * s_min;
 
                 position.x += direction.x >= 0 ? s_eps : -s_eps;
@@ -828,9 +837,10 @@ void ray_tracer_kernel_bw(
 
             }
 
-            // divide out initial layer thicknes, equivalent to first converting lwp (kg/m2) to lwc (kg/m3) 
-            lwp_cam[pix] = lwp_sum / grid_d.z;
-            iwp_cam[pix] = iwp_sum / grid_d.z;
+            // divide out initial layer thicknes, equivalent to first converting lwp (g/m2) to lwc (g/m3) or optical depth to k_ext(1/m)
+            liwp_cam[pix] = liwp_sum / grid_d.z;
+            tauc_cam[pix] = tauc_sum / grid_d.z;
+            dist_cam[pix] = reached_cloud ? dist : Float(-1) ;
 
         }
 

src_test/Radiation_solver_bw.cu

@@ -804,8 +804,9 @@ void Radiation_solver_shortwave::solve_gpu(
         const Aerosol_concs_gpu& aerosol_concs,
         const Camera& camera,
         Array_gpu<Float,3>& XYZ,
-        Array_gpu<Float,2>& lwp_cam,
-        Array_gpu<Float,2>& iwp_cam)
+        Array_gpu<Float,2>& liwp_cam,
+        Array_gpu<Float,2>& tauc_cam,
+        Array_gpu<Float,2>& dist_cam)
 
 {
     const int n_col = p_lay.dim(1);
@@ -1063,14 +1064,24 @@ void Radiation_solver_shortwave::solve_gpu(
 
     if (switch_cloud_cam)
     {
+        cloud_optics_gpu->cloud_optics(
+                11, // 441-615 nm band
+                lwp,
+                iwp,
+                rel,
+                rei,
+                *cloud_optical_props);
+
         raytracer.accumulate_clouds(
                 grid_d,
                 grid_cells,
                 lwp,
                 iwp,
+                dynamic_cast<Optical_props_2str_rt&>(*cloud_optical_props).get_tau(), 
                 camera,
-                lwp_cam,
-                iwp_cam);
+                liwp_cam,
+                tauc_cam,
+                dist_cam);
     }
 }
 
@@ -1102,9 +1113,9 @@ void Radiation_solver_shortwave::solve_gpu_bb(
         const Aerosol_concs_gpu& aerosol_concs,
         const Camera& camera,
         Array_gpu<Float,2>& radiance,
-        Array_gpu<Float,2>& lwp_cam,
-        Array_gpu<Float,2>& iwp_cam)
-
+        Array_gpu<Float,2>& liwp_cam,
+        Array_gpu<Float,2>& tauc_cam,
+        Array_gpu<Float,2>& dist_cam)
 {
     const int n_col = p_lay.dim(1);
     const int n_lay = p_lay.dim(2);
@@ -1309,17 +1320,27 @@ void Radiation_solver_shortwave::solve_gpu_bb(
 
             previous_band = band;
         }
+    }
 
-        if (switch_cloud_cam)
-        {
-            raytracer.accumulate_clouds(
-                    grid_d,
-                    grid_cells,
-                    lwp,
-                    iwp,
-                    camera,
-                    lwp_cam,
-                    iwp_cam);
-        }
+    if (switch_cloud_cam)
+    {
+        cloud_optics_gpu->cloud_optics(
+                11, // 441-615 nm band
+                lwp,
+                iwp,
+                rel,
+                rei,
+                *cloud_optical_props);
+
+        raytracer.accumulate_clouds(
+                grid_d,
+                grid_cells,
+                lwp,
+                iwp,
+                dynamic_cast<Optical_props_2str_rt&>(*cloud_optical_props).get_tau(), 
+                camera,
+                liwp_cam,
+                tauc_cam,
+                dist_cam);
     }
 }

src_test/test_rte_rrtmgp_bw.cu

@@ -238,7 +238,7 @@ void solve_radiation(int argc, char** argv)
         {"shortwave"        , { true,  "Enable computation of shortwave radiation."  }},
         {"longwave"         , { false, "Enable computation of longwave radiation."   }},
         {"fluxes"           , { true,  "Enable computation of fluxes."               }},
-        {"raytracing"       , { false, "Use raytracing for flux computation."        }},
+        {"raytracing"       , { true,  "Use raytracing for flux computation."        }},
         {"cloud-optics"     , { false, "Enable cloud optics."                        }},
         {"cloud-mie"        , { false, "mie cloud droplet scattering."               }},
         {"aerosol-optics"   , { false, "Enable aerosol optics."                      }},
@@ -388,10 +388,7 @@ void solve_radiation(int argc, char** argv)
 
         iwp.set_dims({n_col, n_lay});
         iwp = std::move(input_nc.get_variable<Float>("iwp", {n_lay, n_col_y, n_col_x}));
-    }
-
-    if (switch_cloud_optics)
-    {
+
         rel.set_dims({n_col, n_lay});
         rel = std::move(input_nc.get_variable<Float>("rel", {n_lay, n_col_y, n_col_x}));
 
@@ -697,13 +694,15 @@ void solve_radiation(int argc, char** argv)
                 rad_sw.load_mie_tables("mie_lut_visualisation.nc", switch_broadband);
         }
 
-        Array_gpu<Float,2> lwp_cam;
-        Array_gpu<Float,2> iwp_cam;
+        Array_gpu<Float,2> liwp_cam;
+        Array_gpu<Float,2> tauc_cam;
+        Array_gpu<Float,2> dist_cam;
 
         if (switch_cloud_cam)
         {
-            lwp_cam.set_dims({camera.nx, camera.ny});
-            iwp_cam.set_dims({camera.nx, camera.ny});
+            liwp_cam.set_dims({camera.nx, camera.ny});
+            tauc_cam.set_dims({camera.nx, camera.ny});
+            dist_cam.set_dims({camera.nx, camera.ny});
         }
 
         // Solve the radiation.
@@ -767,8 +766,9 @@ void solve_radiation(int argc, char** argv)
                     aerosol_concs,
                     camera,
                     radiance,
-                    lwp_cam,
-                    iwp_cam);
+                    liwp_cam,
+                    tauc_cam,
+                    dist_cam);
 
             cudaEventRecord(stop, 0);
             cudaEventSynchronize(stop);
@@ -840,8 +840,9 @@ void solve_radiation(int argc, char** argv)
                     aerosol_concs,
                     camera,
                     XYZ,
-                    lwp_cam,
-                    iwp_cam);
+                    liwp_cam,
+                    tauc_cam,
+                    dist_cam);
 
             cudaEventRecord(stop, 0);
             cudaEventSynchronize(stop);
@@ -921,21 +922,23 @@ void solve_radiation(int argc, char** argv)
 
         if (switch_cloud_cam)
         {
-            Array<Float,2> lwp_cam_cpu(lwp_cam);
-            Array<Float,2> iwp_cam_cpu(iwp_cam);
+            Array<Float,2> liwp_cam_cpu(liwp_cam);
+            Array<Float,2> tauc_cam_cpu(tauc_cam);
+            Array<Float,2> dist_cam_cpu(dist_cam);
 
-            auto nc_var_liq = output_nc.add_variable<Float>("lwp_cam", {"y","x"});
-            nc_var_liq.insert(lwp_cam_cpu.v(), {0, 0});
-            nc_var_liq.add_attribute("long_name", "accumulated liquid water path");
+            auto nc_var_liwp = output_nc.add_variable<float>("liq_ice_wp_cam", {"y","x"});
+            nc_var_liwp.insert(liwp_cam_cpu.v(), {0, 0});
+            nc_var_liwp.add_attribute("long_name", "accumulated liquid+ice water path");
 
-            auto nc_var_ice = output_nc.add_variable<Float>("iwp_cam", {"y","x"});
-            nc_var_ice.insert(iwp_cam_cpu.v(), {0, 0});
-            nc_var_ice.add_attribute("long_name", "accumulated ice water path");
-
-
+            auto nc_var_tauc = output_nc.add_variable<float>("tau_cld_cam", {"y","x"});
+            nc_var_tauc.insert(tauc_cam_cpu.v(), {0, 0});
+            nc_var_tauc.add_attribute("long_name", "accumulated cloud optical depth (441-615nm band)");
+
+            auto nc_var_dist = output_nc.add_variable<float>("dist_cld_cam", {"y","x"});
+            nc_var_dist.insert(dist_cam_cpu.v(), {0, 0});
+            nc_var_dist.add_attribute("long_name", "distance to first cloudy cell");
         }
 
-
         auto nc_mu0 = output_nc.add_variable<Float>("sza");
         nc_mu0.insert(acos(mu0({1}))/M_PI * Float(180.), {0});