diff --git a/include_test/Radiation_solver_rt.h b/include_test/Radiation_solver_rt.h
index bb56d27..4d59a97 100644
--- a/include_test/Radiation_solver_rt.h
+++ b/include_test/Radiation_solver_rt.h
@@ -99,6 +99,7 @@ class Radiation_solver_shortwave
         #ifdef __CUDACC__
         void solve_gpu(
                 const bool switch_fluxes,
+                const bool switch_disable_2s,
                 const bool switch_raytracing,
                 const bool switch_independent_column,
                 const bool switch_cloud_optics,
diff --git a/src_cuda_rt/Cloud_optics_rt.cu b/src_cuda_rt/Cloud_optics_rt.cu
index f42e0aa..5c19f60 100644
--- a/src_cuda_rt/Cloud_optics_rt.cu
+++ b/src_cuda_rt/Cloud_optics_rt.cu
@@ -86,6 +86,23 @@ namespace
         }
     }
 
+    __global__
+    void combine_and_store_kernel_single_phase(const int ncol, const int nlay, const Float tmin,
+                  Float* __restrict__ tau,
+                  const Float* __restrict__ l_or_i_tau, const Float* __restrict__ l_or_i_taussa)
+    {
+        const int icol = blockIdx.x*blockDim.x + threadIdx.x;
+        const int ilay = blockIdx.y*blockDim.y + threadIdx.y;
+
+        if ( (icol < ncol) && (ilay < nlay) )
+        {
+            const int idx = icol + ilay*ncol;
+            const Float tau_t = (l_or_i_tau[idx] - l_or_i_taussa[idx]);
+
+            tau[idx] = tau_t;
+        }
+    }
+
     __global__
     void combine_and_store_kernel(const int ncol, const int nlay, const Float tmin,
                   Float* __restrict__ tau, Float* __restrict__ ssa, Float* __restrict__ g,
@@ -108,6 +125,28 @@ namespace
         }
     }
 
+        __global__
+    void combine_and_store_kernel_single_phase(const int ncol, const int nlay, const Float tmin,
+                  Float* __restrict__ tau, Float* __restrict__ ssa, Float* __restrict__ g,
+                  const Float* __restrict__ l_or_i_tau, const Float* __restrict__ l_or_i_taussa, const Float* __restrict__ l_or_i_taussag
+                  )
+    {
+        const int icol = blockIdx.x*blockDim.x + threadIdx.x;
+        const int ilay = blockIdx.y*blockDim.y + threadIdx.y;
+
+        if ( (icol < ncol) && (ilay < nlay) )
+        {
+            const int idx = icol + ilay*ncol;
+            const Float tau_t = l_or_i_tau[idx];
+            const Float taussa = l_or_i_taussa[idx];
+            const Float taussag = l_or_i_taussag[idx];
+
+            tau[idx] = tau_t;
+            ssa[idx] = taussa / max(tau_t, tmin);
+            g[idx]   = taussag/ max(taussa, tmin);
+        }
+    }
+
     __global__
     void set_mask(const int ncol, const int nlay, const Float min_value,
                   Bool* __restrict__ mask, const Float* __restrict__ values)
@@ -182,11 +221,19 @@ void Cloud_optics_rt::cloud_optics(
         const Array_gpu<Float,2>& reliq, const Array_gpu<Float,2>& reice,
         Optical_props_2str_rt& optical_props)
 {
-    const int ncol = clwp.dim(1);
-    const int nlay = clwp.dim(2);
-
-    Optical_props_2str_rt clouds_liq(ncol, nlay, optical_props);
-    Optical_props_2str_rt clouds_ice(ncol, nlay, optical_props);
+    int ncol = -1;
+    int nlay = -1;
+    if (clwp.ptr() != nullptr)
+    {
+        ncol = clwp.dim(1);
+        nlay = clwp.dim(2);
+        Optical_props_2str_rt clouds_liq(ncol, nlay, optical_props);
+    } else if (ciwp.ptr() != nullptr) 
+    {
+        ncol = ciwp.dim(1);
+        nlay = ciwp.dim(2);
+        Optical_props_2str_rt clouds_ice(ncol, nlay, optical_props);
+    }
 
     // Set the mask.
     constexpr Float mask_min_value = Float(0.);
@@ -199,26 +246,36 @@ void Cloud_optics_rt::cloud_optics(
     dim3 grid_m_gpu(grid_col_m, grid_lay_m);
     dim3 block_m_gpu(block_col_m, block_lay_m);
 
-    Array_gpu<Bool,2> liqmsk({ncol, nlay});
-    set_mask<<<grid_m_gpu, block_m_gpu>>>(
-            ncol, nlay, mask_min_value, liqmsk.ptr(), clwp.ptr());
-
-    Array_gpu<Bool,2> icemsk({ncol, nlay});
-    set_mask<<<grid_m_gpu, block_m_gpu>>>(
-            ncol, nlay, mask_min_value, icemsk.ptr(), ciwp.ptr());
 
     // Temporary arrays for storage.
-    Array_gpu<Float,2> ltau    ({ncol, nlay});
-    Array_gpu<Float,2> ltaussa ({ncol, nlay});
-    Array_gpu<Float,2> ltaussag({ncol, nlay});
-
-    Array_gpu<Float,2> itau    ({ncol, nlay});
-    Array_gpu<Float,2> itaussa ({ncol, nlay});
-    Array_gpu<Float,2> itaussag({ncol, nlay});
-
+    Array_gpu<Bool,2> liqmsk({0, 0});
+    Array_gpu<Float,2> ltau    ({0, 0});
+    Array_gpu<Float,2> ltaussa ({0, 0});
+    Array_gpu<Float,2> ltaussag({0, 0});
+    Array_gpu<Bool,2> icemsk({0, 0});
+    Array_gpu<Float,2> itau    ({0, 0});
+    Array_gpu<Float,2> itaussa ({0, 0});
+    Array_gpu<Float,2> itaussag({0, 0});
+    if (clwp.ptr() != nullptr){
+        liqmsk.set_dims({ncol, nlay});
+        ltau.set_dims({ncol, nlay});
+        ltaussa.set_dims({ncol, nlay});
+        ltaussag.set_dims({ncol, nlay});
+
+        set_mask<<<grid_m_gpu, block_m_gpu>>>(
+                ncol, nlay, mask_min_value, liqmsk.ptr(), clwp.ptr());
+    }
+    if (ciwp.ptr() != nullptr){
+        icemsk.set_dims({ncol, nlay});
+        itau.set_dims({ncol, nlay});
+        itaussa.set_dims({ncol, nlay});
+        itaussag.set_dims({ncol, nlay});
+
+        set_mask<<<grid_m_gpu, block_m_gpu>>>(
+                ncol, nlay, mask_min_value, icemsk.ptr(), ciwp.ptr());
+    }
     const int block_col = 64;
     const int block_lay = 1;
-
     const int grid_col  = ncol/block_col + (ncol%block_col > 0);
     const int grid_lay  = nlay/block_lay + (nlay%block_lay > 0);
 
@@ -226,26 +283,44 @@ void Cloud_optics_rt::cloud_optics(
     dim3 block_gpu(block_col, block_lay);
 
     // Liquid water
-    compute_from_table_kernel<<<grid_gpu, block_gpu>>>(
-            ncol, nlay, ibnd-1, liqmsk.ptr(), clwp.ptr(), reliq.ptr(),
-            this->liq_nsteps, this->liq_step_size, this->radliq_lwr,
-            this->lut_extliq_gpu.ptr(), this->lut_ssaliq_gpu.ptr(),
-            this->lut_asyliq_gpu.ptr(), ltau.ptr(), ltaussa.ptr(), ltaussag.ptr());
+    if (clwp.ptr() != nullptr){
+        compute_from_table_kernel<<<grid_gpu, block_gpu>>>(
+                ncol, nlay, ibnd-1, liqmsk.ptr(), clwp.ptr(), reliq.ptr(),
+                this->liq_nsteps, this->liq_step_size, this->radliq_lwr,
+                this->lut_extliq_gpu.ptr(), this->lut_ssaliq_gpu.ptr(),
+                this->lut_asyliq_gpu.ptr(), ltau.ptr(), ltaussa.ptr(), ltaussag.ptr());
+    }
 
     // Ice.
-    compute_from_table_kernel<<<grid_gpu, block_gpu>>>(
-            ncol, nlay, ibnd-1, icemsk.ptr(), ciwp.ptr(), reice.ptr(),
-            this->ice_nsteps, this->ice_step_size, this->radice_lwr,
-            this->lut_extice_gpu.ptr(), this->lut_ssaice_gpu.ptr(),
-            this->lut_asyice_gpu.ptr(), itau.ptr(), itaussa.ptr(), itaussag.ptr());
-
+    if (ciwp.ptr() != nullptr){
+        compute_from_table_kernel<<<grid_gpu, block_gpu>>>(
+                ncol, nlay, ibnd-1, icemsk.ptr(), ciwp.ptr(), reice.ptr(),
+                this->ice_nsteps, this->ice_step_size, this->radice_lwr,
+                this->lut_extice_gpu.ptr(), this->lut_ssaice_gpu.ptr(),
+                this->lut_asyice_gpu.ptr(), itau.ptr(), itaussa.ptr(), itaussag.ptr());
+    }
     constexpr Float eps = std::numeric_limits<Float>::epsilon();
-
-    combine_and_store_kernel<<<grid_gpu, block_gpu>>>(
+    if ((ciwp.ptr() != nullptr) && (clwp.ptr() != nullptr))
+    {
+        combine_and_store_kernel<<<grid_gpu, block_gpu>>>(
             ncol, nlay, eps,
             optical_props.get_tau().ptr(), optical_props.get_ssa().ptr(), optical_props.get_g().ptr(),
             ltau.ptr(), ltaussa.ptr(), ltaussag.ptr(),
             itau.ptr(), itaussa.ptr(), itaussag.ptr());
+    } else if(ciwp.ptr() == nullptr)
+    {
+        combine_and_store_kernel_single_phase<<<grid_gpu, block_gpu>>>(
+            ncol, nlay, eps,
+            optical_props.get_tau().ptr(), optical_props.get_ssa().ptr(), optical_props.get_g().ptr(),
+            ltau.ptr(), ltaussa.ptr(), ltaussag.ptr());
+    } else if (clwp.ptr() == nullptr)
+    {
+       combine_and_store_kernel_single_phase<<<grid_gpu, block_gpu>>>(
+            ncol, nlay, eps,
+            optical_props.get_tau().ptr(), optical_props.get_ssa().ptr(), optical_props.get_g().ptr(),
+            itau.ptr(), itaussa.ptr(), itaussag.ptr()); 
+    }
+
 }
 
 // 1scl variant of cloud optics.
@@ -255,12 +330,20 @@ void Cloud_optics_rt::cloud_optics(
         const Array_gpu<Float,2>& reliq, const Array_gpu<Float,2>& reice,
         Optical_props_1scl_rt& optical_props)
 {
-    const int ncol = clwp.dim(1);
-    const int nlay = clwp.dim(2);
-
-    Optical_props_1scl_rt clouds_liq(ncol, nlay, optical_props);
-    Optical_props_1scl_rt clouds_ice(ncol, nlay, optical_props);
-
+    int ncol = -1;
+    int nlay = -1;
+    if (clwp.ptr() != nullptr)
+    {
+        ncol = clwp.dim(1);
+        nlay = clwp.dim(2);
+        Optical_props_1scl_rt clouds_liq(ncol, nlay, optical_props);
+    } else if (ciwp.ptr() != nullptr) 
+    {
+        ncol = ciwp.dim(1);
+        nlay = ciwp.dim(2);
+        Optical_props_1scl_rt clouds_ice(ncol, nlay, optical_props);
+    }
+    
     // Set the mask.
     constexpr Float mask_min_value = Float(0.);
     const int block_col_m = 16;
@@ -272,22 +355,36 @@ void Cloud_optics_rt::cloud_optics(
     dim3 grid_m_gpu(grid_col_m, grid_lay_m);
     dim3 block_m_gpu(block_col_m, block_lay_m);
 
-    Array_gpu<Bool,2> liqmsk({ncol, nlay});
-    set_mask<<<grid_m_gpu, block_m_gpu>>>(
-            ncol, nlay, mask_min_value, liqmsk.ptr(), clwp.ptr());
-
-    Array_gpu<Bool,2> icemsk({ncol, nlay});
-    set_mask<<<grid_m_gpu, block_m_gpu>>>(
-            ncol, nlay, mask_min_value, icemsk.ptr(), ciwp.ptr());
-
     // Temporary arrays for storage.
-    Array_gpu<Float,2> ltau    ({ncol, nlay});
-    Array_gpu<Float,2> ltaussa ({ncol, nlay});
-    Array_gpu<Float,2> ltaussag({ncol, nlay});
+    Array_gpu<Bool,2> liqmsk({0, 0});
+    Array_gpu<Float,2> ltau    ({0, 0});
+    Array_gpu<Float,2> ltaussa ({0, 0});
+    Array_gpu<Float,2> ltaussag({0, 0});
+    Array_gpu<Bool,2> icemsk({0, 0});
+    Array_gpu<Float,2> itau    ({0, 0});
+    Array_gpu<Float,2> itaussa ({0, 0});
+    Array_gpu<Float,2> itaussag({0, 0});
+
+    if (clwp.ptr() != nullptr)
+    {
+        liqmsk.set_dims({ncol, nlay});
+        ltau.set_dims({ncol, nlay});
+        ltaussa.set_dims({ncol, nlay});
+        ltaussag.set_dims({ncol, nlay});
+
+        set_mask<<<grid_m_gpu, block_m_gpu>>>(
+            ncol, nlay, mask_min_value, liqmsk.ptr(), clwp.ptr());
+    }
+    if (ciwp.ptr() != nullptr)
+    {
+        icemsk.set_dims({ncol, nlay});
+        itau.set_dims({ncol, nlay});
+        itaussa.set_dims({ncol, nlay});
+        itaussag.set_dims({ncol, nlay});
 
-    Array_gpu<Float,2> itau    ({ncol, nlay});
-    Array_gpu<Float,2> itaussa ({ncol, nlay});
-    Array_gpu<Float,2> itaussag({ncol, nlay});
+        set_mask<<<grid_m_gpu, block_m_gpu>>>(
+                ncol, nlay, mask_min_value, icemsk.ptr(), ciwp.ptr());
+    }
 
     const int block_col = 64;
     const int block_lay = 1;
@@ -299,25 +396,45 @@ void Cloud_optics_rt::cloud_optics(
     dim3 block_gpu(block_col, block_lay);
 
     // Liquid water
-    compute_from_table_kernel<<<grid_gpu, block_gpu>>>(
-            ncol, nlay, ibnd-1, liqmsk.ptr(), clwp.ptr(), reliq.ptr(),
-            this->liq_nsteps, this->liq_step_size, this->radliq_lwr,
-            this->lut_extliq_gpu.ptr(), this->lut_ssaliq_gpu.ptr(),
-            this->lut_asyliq_gpu.ptr(), ltau.ptr(), ltaussa.ptr(), ltaussag.ptr());
+    if (clwp.ptr() != nullptr){
+        compute_from_table_kernel<<<grid_gpu, block_gpu>>>(
+                ncol, nlay, ibnd-1, liqmsk.ptr(), clwp.ptr(), reliq.ptr(),
+                this->liq_nsteps, this->liq_step_size, this->radliq_lwr,
+                this->lut_extliq_gpu.ptr(), this->lut_ssaliq_gpu.ptr(),
+                this->lut_asyliq_gpu.ptr(), ltau.ptr(), ltaussa.ptr(), ltaussag.ptr());
+    }
 
     // Ice.
-    compute_from_table_kernel<<<grid_gpu, block_gpu>>>(
-            ncol, nlay, ibnd-1, icemsk.ptr(), ciwp.ptr(), reice.ptr(),
-            this->ice_nsteps, this->ice_step_size, this->radice_lwr,
-            this->lut_extice_gpu.ptr(), this->lut_ssaice_gpu.ptr(),
-            this->lut_asyice_gpu.ptr(), itau.ptr(), itaussa.ptr(), itaussag.ptr());
+    if (ciwp.ptr() != nullptr){
+        compute_from_table_kernel<<<grid_gpu, block_gpu>>>(
+                ncol, nlay, ibnd-1, icemsk.ptr(), ciwp.ptr(), reice.ptr(),
+                this->ice_nsteps, this->ice_step_size, this->radice_lwr,
+                this->lut_extice_gpu.ptr(), this->lut_ssaice_gpu.ptr(),
+                this->lut_asyice_gpu.ptr(), itau.ptr(), itaussa.ptr(), itaussag.ptr());
+    }
 
     constexpr Float eps = std::numeric_limits<Float>::epsilon();
-
-    combine_and_store_kernel<<<grid_gpu, block_gpu>>>(
+    if ((ciwp.ptr() != nullptr) && (clwp.ptr() != nullptr))
+    {
+        combine_and_store_kernel<<<grid_gpu, block_gpu>>>(
             ncol, nlay, eps,
             optical_props.get_tau().ptr(),
             ltau.ptr(), ltaussa.ptr(),
             itau.ptr(), itaussa.ptr());
+    } else if(ciwp.ptr() == nullptr)
+    {
+        combine_and_store_kernel_single_phase<<<grid_gpu, block_gpu>>>(
+            ncol, nlay, eps,
+            optical_props.get_tau().ptr(),
+            ltau.ptr(), ltaussa.ptr());
+    } else if(clwp.ptr() == nullptr)
+    {
+        combine_and_store_kernel_single_phase<<<grid_gpu, block_gpu>>>(
+            ncol, nlay, eps,
+            optical_props.get_tau().ptr(),
+            itau.ptr(), itaussa.ptr());
+
+    }
+
 }
 
diff --git a/src_test/Radiation_solver_rt.cu b/src_test/Radiation_solver_rt.cu
index 49348a2..883fe8e 100644
--- a/src_test/Radiation_solver_rt.cu
+++ b/src_test/Radiation_solver_rt.cu
@@ -567,6 +567,7 @@ void Radiation_solver_shortwave::load_mie_tables(
 
 void Radiation_solver_shortwave::solve_gpu(
         const bool switch_fluxes,
+        const bool switch_disable_2s,
         const bool switch_raytracing,
         const bool switch_independent_column,
         const bool switch_cloud_optics,
@@ -636,10 +637,13 @@ void Radiation_solver_shortwave::solve_gpu(
 
     if (switch_fluxes)
     {
-        Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_up.ptr());
-        Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn.ptr());
-        Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn_dir.ptr());
-        Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_net.ptr());
+        if (!switch_disable_2s)
+        {
+                Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_up.ptr());
+                Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn.ptr());
+                Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn_dir.ptr());
+                Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_net.ptr());
+        }
         if (switch_raytracing)
         {
             Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_tod_up.ptr());
@@ -836,10 +840,12 @@ void Radiation_solver_shortwave::solve_gpu(
             }
 
             (*fluxes).net_flux();
-
-            Gpt_combine_kernels_cuda_rt::add_from_gpoint(
-                    n_col, n_lev, sw_flux_up.ptr(), sw_flux_dn.ptr(), sw_flux_dn_dir.ptr(), sw_flux_net.ptr(),
-                    (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_dn_dir().ptr(), (*fluxes).get_flux_net().ptr());
+            if (!switch_disable_2s)
+            {
+                Gpt_combine_kernels_cuda_rt::add_from_gpoint(
+                        n_col, n_lev, sw_flux_up.ptr(), sw_flux_dn.ptr(), sw_flux_dn_dir.ptr(), sw_flux_net.ptr(),
+                        (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_dn_dir().ptr(), (*fluxes).get_flux_net().ptr());
+            }
 
             if (switch_raytracing)
             {
diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu
index fc7f204..20d43f7 100644
--- a/src_test/test_rte_rrtmgp_rt.cu
+++ b/src_test/test_rte_rrtmgp_rt.cu
@@ -227,19 +227,26 @@ 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."             }},
+        {"disable-2s"        , { false, "use raytracing onlu for flux computation. must be passed with raytracing" }},
         {"raytracing"        , { true,  "Use raytracing for flux computation. '--raytracing 256': use 256 rays per pixel" }},
         {"independent-column", { false, "run raytracer in independent column mode"}},
-        {"cloud-optics"      , { false, "Enable cloud optics."                      }},
+        {"cloud-optics"      , { false, "Enable cloud optics (both liquid and ice)."}},
+        {"liq-cloud-optics"  , { false, "liquid only cloud optics."                 }},
+        {"ice-cloud-optics"  , { false, "ice only cloud optics."                    }},
         {"cloud-mie"         , { false, "Use Mie tables for cloud scattering in ray tracer"  }},
         {"aerosol-optics"    , { false, "Enable aerosol optics."                    }},
         {"single-gpt"        , { false, "Output optical properties and fluxes for a single g-point. '--single-gpt 100': output 100th g-point" }},
         {"profiling"         , { false, "Perform additional profiling run."         }},
         {"delta-cloud"       , { false, "delta-scaling of cloud optical properties"   }},
-        {"delta-aerosol"     , { false, "delta-scaling of aerosol optical properties"   }} };
+        {"delta-aerosol"     , { false, "delta-scaling of aerosol optical properties"   }} ,
+        {"override-sza"     , { false, "override provided value of sza in input file. IN DEGREES. '--override-sza 50': use a sza of 50 degrees" }},
+        {"override-azi"     , { false, "override provided value of azi in input file. IN DEGREES. '--override-azi 240': use of azi of 240 degrees"   }}};
 
     std::map<std::string, std::pair<int, std::string>> command_line_ints {
         {"raytracing", {32, "Number of rays initialised at TOD per pixel per quadraute."}},
-        {"single-gpt", {1 , "g-point to store optical properties and fluxes of" }} };
+        {"single-gpt", {1 , "g-point to store optical properties and fluxes of" }},
+        {"override-sza", {0, "solar zenith angle (theta) in degrees."}},
+        {"override-azi", {0, "Solar azimuth angle in degrees."}} };
 
     if (parse_command_line_options(command_line_switches, command_line_ints, argc, argv))
         return;
@@ -247,16 +254,20 @@ void solve_radiation(int argc, char** argv)
     const bool switch_shortwave         = command_line_switches.at("shortwave"         ).first;
     const bool switch_longwave          = command_line_switches.at("longwave"          ).first;
     const bool switch_fluxes            = command_line_switches.at("fluxes"            ).first;
+    const bool switch_disable_2s        = command_line_switches.at("disable-2s"        ).first;
     const bool switch_raytracing        = command_line_switches.at("raytracing"        ).first;
     const bool switch_independent_column= command_line_switches.at("independent-column").first;
-    const bool switch_cloud_optics      = command_line_switches.at("cloud-optics"      ).first;
+    bool switch_cloud_optics      = command_line_switches.at("cloud-optics"      ).first;
+    bool switch_liq_cloud_optics  = command_line_switches.at("liq-cloud-optics"  ).first;
+    bool switch_ice_cloud_optics  = command_line_switches.at("ice-cloud-optics"  ).first;
     const bool switch_cloud_mie         = command_line_switches.at("cloud-mie"         ).first;
     const bool switch_aerosol_optics    = command_line_switches.at("aerosol-optics"    ).first;
     const bool switch_single_gpt        = command_line_switches.at("single-gpt"        ).first;
     const bool switch_profiling         = command_line_switches.at("profiling"         ).first;
     const bool switch_delta_cloud       = command_line_switches.at("delta-cloud"       ).first;
     const bool switch_delta_aerosol     = command_line_switches.at("delta-aerosol"     ).first;
-
+    const bool override_sza             = command_line_switches.at("override-sza"    ).first;
+    const bool override_azi             = command_line_switches.at("override-azi"    ).first;
 
     Int photons_per_pixel = Int(command_line_ints.at("raytracing").first);
     if (Float(int(std::log2(Float(photons_per_pixel)))) != std::log2(Float(photons_per_pixel)))
@@ -271,13 +282,40 @@ void solve_radiation(int argc, char** argv)
         throw std::runtime_error(error);
     }
 
+    if (switch_disable_2s && !switch_raytracing) {
+        std::string error = "cannot disable two-stream for flux calculation without turning ray tracing on";
+        throw std::runtime_error(error);
+    }
+
+    if (switch_cloud_optics)
+    {
+        switch_liq_cloud_optics = true;
+        switch_ice_cloud_optics = true;
+    }
+    if (switch_liq_cloud_optics || switch_ice_cloud_optics)
+    {
+        switch_cloud_optics = true;
+    }
+
     // Print the options to the screen.
     print_command_line_options(command_line_switches, command_line_ints);
 
     int single_gpt = command_line_ints.at("single-gpt").first;
+    int sza_deg = Int(command_line_ints.at("override-sza").first);
+    int azi_deg = Int(command_line_ints.at("override-azi").first);
 
     Status::print_message("Using "+ std::to_string(photons_per_pixel) + " rays per pixel");
 
+    if (override_sza) 
+    {
+        Status::print_message("Using SZA of "+ std::to_string(sza_deg) + " degrees");
+    }
+
+    if (override_azi) 
+    {
+        Status::print_message("Using azi of "+ std::to_string(azi_deg) + " degrees");
+    }
+
     ////// READ THE ATMOSPHERIC DATA //////
     Status::print_message("Reading atmospheric input data from NetCDF.");
 
@@ -353,17 +391,22 @@ void solve_radiation(int argc, char** argv)
 
     if (switch_cloud_optics)
     {
-        lwp.set_dims({n_col, n_lay});
-        lwp = std::move(input_nc.get_variable<Float>("lwp", {n_lay, n_col_y, n_col_x}));
+        
+        if(switch_liq_cloud_optics){
+            lwp.set_dims({n_col, n_lay});
+            lwp = std::move(input_nc.get_variable<Float>("lwp", {n_lay, n_col_y, n_col_x}));
 
-        iwp.set_dims({n_col, n_lay});
-        iwp = std::move(input_nc.get_variable<Float>("iwp", {n_lay, n_col_y, n_col_x}));
+            rel.set_dims({n_col, n_lay});
+            rel = std::move(input_nc.get_variable<Float>("rel", {n_lay, n_col_y, n_col_x}));
+        }
 
-        rel.set_dims({n_col, n_lay});
-        rel = std::move(input_nc.get_variable<Float>("rel", {n_lay, n_col_y, n_col_x}));
+        if(switch_ice_cloud_optics){
+            iwp.set_dims({n_col, n_lay});
+            iwp = std::move(input_nc.get_variable<Float>("iwp", {n_lay, n_col_y, n_col_x}));
 
-        rei.set_dims({n_col, n_lay});
-        rei = std::move(input_nc.get_variable<Float>("rei", {n_lay, n_col_y, n_col_x}));
+            rei.set_dims({n_col, n_lay});
+            rei = std::move(input_nc.get_variable<Float>("rei", {n_lay, n_col_y, n_col_x}));
+        }
     }
 
     Array<Float,2> rh;
@@ -639,9 +682,25 @@ void solve_radiation(int argc, char** argv)
             rad_sw.load_mie_tables("mie_lut_broadband.nc");
         }
 
+        Array<Float,1> mu0({n_col});
+        Array<Float,1> azi({n_col});
+
+        if (override_sza) {
+            Float mu0_in = cosf(sza_deg * 3.14159f / 180.0f);
+            for (int icol=1; icol<=n_col; ++icol)
+                mu0({icol}) = mu0_in;
+        } else {
+            mu0 = input_nc.get_variable<Float>("mu0", {n_col_y, n_col_x});
+        }
+
+        if (override_azi) {
+            Float azi_in = azi_deg * 3.14159f / 180.0f;
+            for (int icol=1; icol<=n_col; ++icol)
+                azi({icol}) = azi_in;
+        } else {
+            azi = input_nc.get_variable<Float>("azi", {n_col_y, n_col_x});
+        }
 
-        Array<Float,1> mu0(input_nc.get_variable<Float>("mu0", {n_col_y, n_col_x}), {n_col});
-        Array<Float,1> azi(input_nc.get_variable<Float>("azi", {n_col_y, n_col_x}), {n_col});
         Array<Float,2> sfc_alb_dir(input_nc.get_variable<Float>("sfc_alb_dir", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col});
         Array<Float,2> sfc_alb_dif(input_nc.get_variable<Float>("sfc_alb_dif", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col});
 
@@ -702,10 +761,13 @@ void solve_radiation(int argc, char** argv)
 
         if (switch_fluxes)
         {
-            sw_flux_up    .set_dims({n_col, n_lev});
-            sw_flux_dn    .set_dims({n_col, n_lev});
-            sw_flux_dn_dir.set_dims({n_col, n_lev});
-            sw_flux_net   .set_dims({n_col, n_lev});
+            if(!switch_disable_2s)
+            {
+                sw_flux_up    .set_dims({n_col, n_lev});
+                sw_flux_dn    .set_dims({n_col, n_lev});
+                sw_flux_dn_dir.set_dims({n_col, n_lev});
+                sw_flux_net   .set_dims({n_col, n_lev});
+            }
 
             if (switch_raytracing)
             {
@@ -766,6 +828,7 @@ void solve_radiation(int argc, char** argv)
 
             rad_sw.solve_gpu(
                     switch_fluxes,
+                    switch_disable_2s,
                     switch_raytracing,
                     switch_independent_column,
                     switch_cloud_optics,
@@ -903,27 +966,31 @@ void solve_radiation(int argc, char** argv)
 
         if (switch_fluxes)
         {
-            auto nc_sw_flux_up     = output_nc.add_variable<Float>("sw_flux_up"    , {"lev", "y", "x"});
-            auto nc_sw_flux_dn     = output_nc.add_variable<Float>("sw_flux_dn"    , {"lev", "y", "x"});
-            auto nc_sw_flux_dn_dir = output_nc.add_variable<Float>("sw_flux_dn_dir", {"lev", "y", "x"});
-            auto nc_sw_flux_net    = output_nc.add_variable<Float>("sw_flux_net"   , {"lev", "y", "x"});
+            if (!switch_disable_2s)
+            {
+                auto nc_sw_flux_up     = output_nc.add_variable<Float>("sw_flux_up"    , {"lev", "y", "x"});
+                auto nc_sw_flux_dn     = output_nc.add_variable<Float>("sw_flux_dn"    , {"lev", "y", "x"});
+                auto nc_sw_flux_dn_dir = output_nc.add_variable<Float>("sw_flux_dn_dir", {"lev", "y", "x"});
+                auto nc_sw_flux_net    = output_nc.add_variable<Float>("sw_flux_net"   , {"lev", "y", "x"});
+
+                nc_sw_flux_up    .insert(sw_flux_up_cpu    .v(), {0, 0, 0});
+                nc_sw_flux_dn    .insert(sw_flux_dn_cpu    .v(), {0, 0, 0});
+                nc_sw_flux_dn_dir.insert(sw_flux_dn_dir_cpu.v(), {0, 0, 0});
+                nc_sw_flux_net   .insert(sw_flux_net_cpu   .v(), {0, 0, 0});
 
-            nc_sw_flux_up    .insert(sw_flux_up_cpu    .v(), {0, 0, 0});
-            nc_sw_flux_dn    .insert(sw_flux_dn_cpu    .v(), {0, 0, 0});
-            nc_sw_flux_dn_dir.insert(sw_flux_dn_dir_cpu.v(), {0, 0, 0});
-            nc_sw_flux_net   .insert(sw_flux_net_cpu   .v(), {0, 0, 0});
+                nc_sw_flux_up.add_attribute("long_name","Upwelling shortwave fluxes (TwoStream solver)");
+                nc_sw_flux_up.add_attribute("units","W m-2");
 
-            nc_sw_flux_up.add_attribute("long_name","Upwelling shortwave fluxes (TwoStream solver)");
-            nc_sw_flux_up.add_attribute("units","W m-2");
+                nc_sw_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes (TwoStream solver)");
+                nc_sw_flux_dn.add_attribute("units","W m-2");
 
-            nc_sw_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes (TwoStream solver)");
-            nc_sw_flux_dn.add_attribute("units","W m-2");
+                nc_sw_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes (TwoStream solver)");
+                nc_sw_flux_dn_dir.add_attribute("units","W m-2");
 
-            nc_sw_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes (TwoStream solver)");
-            nc_sw_flux_dn_dir.add_attribute("units","W m-2");
+                nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)");
+                nc_sw_flux_net.add_attribute("units","W m-2");
 
-            nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)");
-            nc_sw_flux_net.add_attribute("units","W m-2");
+            }        
 
             if (switch_raytracing)
             {
@@ -964,28 +1031,31 @@ void solve_radiation(int argc, char** argv)
 
 
             if (switch_single_gpt)
-            {
-                auto nc_sw_gpt_flux_up     = output_nc.add_variable<Float>("sw_gpt_flux_up"    , {"lev", "y", "x"});
-                auto nc_sw_gpt_flux_dn     = output_nc.add_variable<Float>("sw_gpt_flux_dn"    , {"lev", "y", "x"});
-                auto nc_sw_gpt_flux_dn_dir = output_nc.add_variable<Float>("sw_gpt_flux_dn_dir", {"lev", "y", "x"});
-                auto nc_sw_gpt_flux_net    = output_nc.add_variable<Float>("sw_gpt_flux_net"   , {"lev", "y", "x"});
-
-                nc_sw_gpt_flux_up    .insert(sw_gpt_flux_up_cpu    .v(), {0, 0, 0});
-                nc_sw_gpt_flux_dn    .insert(sw_gpt_flux_dn_cpu    .v(), {0, 0, 0});
-                nc_sw_gpt_flux_dn_dir.insert(sw_gpt_flux_dn_dir_cpu.v(), {0, 0, 0});
-                nc_sw_gpt_flux_net   .insert(sw_gpt_flux_net_cpu   .v(), {0, 0, 0});
-
-                nc_sw_gpt_flux_up.add_attribute("long_name","Upwelling shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)");
-                nc_sw_gpt_flux_up.add_attribute("units","W m-2");
-
-                nc_sw_gpt_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)");
-                nc_sw_gpt_flux_dn.add_attribute("units","W m-2");
-
-                nc_sw_gpt_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)");
-                nc_sw_gpt_flux_dn_dir.add_attribute("units","W m-2");
-
-                nc_sw_gpt_flux_net.add_attribute("long_name","Net shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)");
-                nc_sw_gpt_flux_net.add_attribute("units","W m-2");
+            { 
+                if (!switch_disable_2s)
+                {
+                    auto nc_sw_gpt_flux_up     = output_nc.add_variable<Float>("sw_gpt_flux_up"    , {"lev", "y", "x"});
+                    auto nc_sw_gpt_flux_dn     = output_nc.add_variable<Float>("sw_gpt_flux_dn"    , {"lev", "y", "x"});
+                    auto nc_sw_gpt_flux_dn_dir = output_nc.add_variable<Float>("sw_gpt_flux_dn_dir", {"lev", "y", "x"});
+                    auto nc_sw_gpt_flux_net    = output_nc.add_variable<Float>("sw_gpt_flux_net"   , {"lev", "y", "x"});
+
+                    nc_sw_gpt_flux_up    .insert(sw_gpt_flux_up_cpu    .v(), {0, 0, 0});
+                    nc_sw_gpt_flux_dn    .insert(sw_gpt_flux_dn_cpu    .v(), {0, 0, 0});
+                    nc_sw_gpt_flux_dn_dir.insert(sw_gpt_flux_dn_dir_cpu.v(), {0, 0, 0});
+                    nc_sw_gpt_flux_net   .insert(sw_gpt_flux_net_cpu   .v(), {0, 0, 0});
+
+                    nc_sw_gpt_flux_up.add_attribute("long_name","Upwelling shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)");
+                    nc_sw_gpt_flux_up.add_attribute("units","W m-2");
+
+                    nc_sw_gpt_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)");
+                    nc_sw_gpt_flux_dn.add_attribute("units","W m-2");
+
+                    nc_sw_gpt_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)");
+                    nc_sw_gpt_flux_dn_dir.add_attribute("units","W m-2");
+
+                    nc_sw_gpt_flux_net.add_attribute("long_name","Net shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)");
+                    nc_sw_gpt_flux_net.add_attribute("units","W m-2");
+                }
             }
         }
     }