Possible bug in patchy cloud implementation for reflected light

[Nicholaswogan]

When do_holes == True, then get_reflected_1d is called with input optical properties representing a clear-sky atmosphere. Possible issue is that this call takes in the CLOUDY variants of GCOS2, ftau_cld, and ftau_ray.

This may have been intentional, but maybe not. I don't understand these terms enough to know for sure.

To test if it was a bug I did the following experiment:

• Computed a reflected light spectrum where I never called inputs.clouds.
• Computed a reflected light spectrum where I called inputs.clouds, but with do_holes=True, fhole=1.0 and fthin_cld=0.0

These two calculations do not produce the same result (even though intuition says they should), suggesting it is possibly a bug.

picaso/picaso/justdoit.py

Line 330 in 1930fc8

GCOS2[:,:,ig],ftau_cld[:,:,ig],ftau_ray[:,:,ig],

[Nicholaswogan]

Two other possible issues:

• The variable TAU should be computed after the check for do_holes
• ftau_cld is never updated. Not sure if this matters?

picaso/picaso/optics.py

Line 351 in 7061f8c

TAU = np.zeros((nlayer+1, nwno,ngauss))

diff --git a/picaso/optics.py b/picaso/optics.py
index 02292c2..85570bf 100755
--- a/picaso/optics.py
+++ b/picaso/optics.py
@@ -347,19 +347,20 @@ def compute_opacity(atmosphere, opacityclass, ngauss=1, stream=2, delta_eddingto
     #scattering from clouds still
     W0_no_raman = (TAURAY*0.99999 + TAUCLD*single_scattering_cld) / (TAUGAS + TAURAY + TAUCLD) #TOTAL single scattering 
 
-    #sum up taus starting at the top, going to depth
-    TAU = np.zeros((nlayer+1, nwno,ngauss))
-    for igauss in range(ngauss): TAU[1:,:,igauss]=numba_cumsum(DTAU[:,:,igauss])
-
     # Clearsky case
     if do_holes == True:
         DTAU = TAUGAS + TAURAY + fthin_cld*TAUCLD #fraction of cloud opacity
         COSB = fthin_cld*np.copy(asym_factor_cld) #fraction of cloud asymmetry
+        ftau_cld = (single_scattering_cld * TAUCLD * fthin_cld)/(TAURAY + single_scattering_cld * TAUCLD * fthin_cld)
         ftau_ray = TAURAY/(TAURAY + single_scattering_cld * TAUCLD *fthin_cld)
         GCOS2 = 0.5*ftau_ray # since ftau_ray = 1 without any clouds
         W0 = (TAURAY*raman_factor + fthin_cld*TAUCLD*single_scattering_cld) / DTAU #TOTAL single scattering
         W0_no_raman = (TAURAY*0.99999 + TAUCLD*single_scattering_cld* fthin_cld) / DTAU #TOTAL single scattering
 
+    #sum up taus starting at the top, going to depth
+    TAU = np.zeros((nlayer+1, nwno,ngauss))
+    for igauss in range(ngauss): TAU[1:,:,igauss]=numba_cumsum(DTAU[:,:,igauss])
+
     if plot_opacity:
         opt_figure.line(1e4/opacityclass.wno, DTAU[plot_layer,:,0], legend_label='TOTAL', line_width=4, color=colors[0],
             muted_color=colors[c], muted_alpha=0.2)

[natashabatalha]

Hmm.. it's also not clear to me asymmetry should be adjust for a thinning factor.

COSB = fthin_cld*np.copy(asym_factor_cld) #fraction of cloud asymmetry

Ill think about this and readjust the code