This is an attempt to document the Landsat 8 (L8) surface reflectance (and top-of-atmosphere reflectance) algorithm delivered to the USGS by Eric Vermote, NASA GSFC. The L8SR algorithm was originally delivered as FORTRAN code, and has since been converted to C (with some speed-ups) by the USGS LSRD staff.

The application uses static input tables in addition to auxiliary files which contain climate information for the date of the scene.

L8SR Inputs/Outputs
===================
Static Inputs:  Each of these (except rationm) are tables of coefficients produced by running the 6S algorithm. The filenames for these specify URBAN_CLEAN, thus the assumption is that 6S was run using the URBAN CLEAN model (??).
anglehdf - angle HDF filename
intrefnm - intrinsic reflectance filename
transmnm - transmission filename
sphereanm - atmosphere spherical albedo filename
cmgdemnm - climate modeling grid DEM filename
rationm - ratio averages filename; "ratio map" used by the aerosol retrieval algorithm

Inputs:
auxnm - auxiliary filename for ozone and water vapor (created using LAADS datasets; usually lag about 2-3 days behind the current date)
xml_infile - Landsat 8 XML file which contains information about the L8 scene and each of the bands for that scene
write_toa - should the TOA values be written for bands 1-7 in addition to the SR values?
process_sr - should SR corrections be applied or just stop at TOA?

Outputs:
TOA bands - top-of-atmosphere values for bands 1-7, 9 (Watts/( m2 * srad * µm)) (TOA values for bands 1-7 will only be written if the write_toa flag was specified, otherwise only SR values are written for those bands.); scale factor is 0.0001 to get to the actual TOA values
brightness temp (BT) bands - brightness temperature values for the thermal bands 10 and 11 (kelvin); scale factor is 0.1 to get to the actual BT values
SR bands - surface reflectance values for bands 1-7 (Watts/( m2 * srad * µm)); scale factor is 0.0001 to get to the actual TOA values. These are only produced if the process_sr flag was specified. SR corrections are not possible for scenes which don't have the TIRS bands.  They are also not possible for scenes with solar zenith angles greater than 76 degrees.)

Overall L8SR Algorithm
======================
Set constants for the overall processing:
    Set MULT_FACTOR = 10000.0  (output multiplier for reflectance bands)
    Set MULT_FACTOR_TH = 10.0  (output multiplier for thermal bands)
    Set MIN_VALID = -2000      (minimum valid scaled reflectance value)
    Set MAX_VALID = 16000      (maximum valid scaled reflectance value)
    Set MIN_VALID_TH = 1500    (minimum valid scaled thermal value)
    Set MAX_VALID_TH = 3500    (maximum valid scaled thermal value)
    Set SCALE_FACTOR = 0.0001  (used to unscale the reflectance bands)
    Set SCALE_FACTOR_TH = 0.1  (used to unscale the thermal bands)

Open and read the information in the XML file
Set xts to the solar zenith angle at the scene center (degrees)
Set xfs to the solar azimuth angle at the scene center (degrees)
Set xmus to the cosine of the solar zenith angle (converted to radians)
Set pixsize to the pixel size of the input reflectance bands
Set nlines to the number of lines in the input reflectance bands
Set nsamps to the number of samples in the input reflectance bands
If this is an OLI-only scene then SR corrections cannot be applied
If the solar zenith angle > 76 degrees then SR corrections cannot be applied
Read the input QA band
If process_sr then make sure the auxiliary files exist
Compute the TOA reflectance and brightness temperatures (BT) (see TOA Algorithm Description)
If write_toa then
    Write the TOA values for bands 1-7 (reflectance bands)
    Append the TOA bands 1-7 to the XML file
Endif
Write the scaled TOA and BT values for bands 9-11 (cirrus and thermal bands)
Append the TOA and BT bands 9-11 to the XML file
If process_sr then
    Apply the surface reflectance corrections (see SR Algorithm Description)
Endif
Write the scaled SR values for bands 1-7 (reflectance bands)
Append the SR bands 1-7 to the XML file
Write the cloud mask band
Append the cloud mask band to the XML file
   ==> Cloud mask bit values representation (LSB to MSB):
       bit 0 - cirrus cloud bit
               ==> 0 - not cirrus cloud
               ==> 1 - cirrus cloud
       bit 1 - cloud bit
               ==> 0 - not cloud
               ==> 1 - cloud
       bit 2 - adjacent cloud bit
               ==> 0 - not adjacentcloud
               ==> 1 - adjacent cloud
       bit 3 - cloud shadow bit
               ==> 0 - not cloud shadow
               ==> 1 - cloud shadow
       bit 4 - first aerosol bit
       bit 5 - second aerosol bit
               ==> 00 - climatology (corrections based on aerosol climatology)
                   01 - low aerosols
                   10 - average aerosols
                   11 - high aerosols (not recommended for use, leads to
                        lower quality surface reflectance product)
       bit 6 - not used
       bit 7 - water bit (internal)
               ==> 0 - not water
               ==> 1 - water
    

TOA Algorithm Description
=========================
Inputs:
    Scene metadata
    QA band
    xmus  (cosine of the solar zenith angle at scene center)

Outputs:
    sband (output scaled TOA reflectance and brightness temp values for bands 1-11, except band 8)

Set constants for the TOA and BT computations:
    Set refl_mult = 2.0E-05  (reflectance multiplier for bands 1-9)
    Set refl_add = -0.1      (reflectance additive for bands 1-9)
    Set xcals = 3.3420E-04   (radiance multiplier for bands 10 and 11)
    Set xcalo = 0.10000      (radiance additive for bands 10 and 11)
    Set k1b10 = 774.89       (temperature constant for band 10)
    Set k1b11 = 480.89       (temperature constant for band 11)
    Set k2b10 = 1321.08      (temperature constant for band 10)
    Set k2b11 = 1201.14      (temperature constant for band 11)

Loop through bands 1 through 11
    If band 8 (pan band) then don't process
    If band 1 through 9
        Read the input DN values for this band
        Loop through each pixel in the band
            If the QA for this pixel isn't FILL
                /* Compute the TOA reflectance (rotoa) based on the scene
                   center sun angle.  Scale the value for output. */
                Set rotoa = (DN_value[pixel] * refl_mult) + refl_add
                Set rotoa = rotoa * MULT_FACTOR / xmus

                /* Save the scaled TOA reflectance value (sband), but make sure
                   it falls within the defined valid range. */
                If (rotoa < MIN_VALID)
                    Set sband[band][pixel] = MIN_VALID
                Else if (rotoa > MAX_VALID)
                    Set sband[band][pixel] = MAX_VALID
                Else
                    Set sband[band][pixel] = (int) rotoa
                Endif
            Else
                Set sband[band][pixel] = FILL
            Endif    
        Endloop
    Else if band 10
        Read the input DN values for this band
        Loop through each pixel in the band
            If the QA for this pixel isn't FILL
                /* Compute the TOA spectral radiance */
                Set tmpf = xcals * DN_value[pixel] + xcalo

                /* Compute the at-satellite brightness temp (K) and scale for
                   output */
                Set tmpf = k2b10 / log (k1b10 / tmpf + 1.0)
                Set tmpf = tmpf * MULT_FACTOR_TH

                /* Save the scaled BT value (sband), but make sure the
                   brightness temp falls within the specified range */
                If (tmpf < MIN_VALID_TH)
                    Set sband[band][pixel] = MIN_VALID_TH
                Else if (tmpf > MAX_VALID_TH)
                    Set sband[band][pixel] = MAX_VALID_TH
                Else
                    Set sband[band][pixel] = (int) (tmpf + 0.5)
                Endif
            Else
                Set sband[band][pixel] = FILL
            Endif    
        Endloop
    Else if band 11
        Read the input DN values for this band
        Loop through each pixel in the band
            If the QA for this pixel isn't FILL
                /* Compute the TOA spectral radiance */
                Set tmpf = xcals * DN_value[pixel] + xcalo

                /* Compute the at-satellite brightness temp (K) and scale for
                   output */
                Set tmpf = k2b11 / log (k1b11 / tmpf + 1.0)
                Set tmpf = tmpf * MULT_FACTOR_TH

                /* Save the scaled BT value (sband), but make sure the
                   brightness temp falls within the specified range */
                If (tmpf < MIN_VALID_TH)
                    Set sband[band][pixel] = MIN_VALID_TH
                Else if (tmpf > MAX_VALID_TH)
                    Set sband[band][pixel] = MAX_VALID_TH
                Else
                    Set sband[band][pixel] = (int) (tmpf + 0.5)
                Endif
            Else
                Set sband[band][pixel] = FILL
            Endif    
        Endloop
    Endif
Endloop


SR Algorithm Description
========================
Inputs:
    Scene metadata
    QA band
    sband  (TOA and BT values for bands 1-11 (except band 8))
    xmus   (cosine of the solar zenith angle at scene center)
    xts    (solar zenith angle at the scene center)
    xfs    (solar azimuth angle at the scene center)

Outputs:
    sband  (output scaled surface reflectance values for bands 1-7)

Set constants for the SR computations:
    Set cfac = 6.0               (cloud factor)
    Set aot550nm[22] =           (AOT look-up table)
        {0.01, 0.05, 0.10, 0.15, 0.20, 0.30, 0.40, 0.60, 0.80, 1.00, 1.20,
         1.40, 1.60, 1.80, 2.00, 2.30, 2.60, 3.00, 3.50, 4.00, 4.50, 5.00}
    Set float tpres[7] =         (surface pressure table)
        {1050.0, 1013.0, 900.0, 800.0, 700.0, 600.0, 500.0}
    Set tauray[NSR_BANDS] =      (molecular optical thickness coefficients --
        produced by running 6S)
        {0.23638, 0.16933, 0.09070, 0.04827, 0.01563, 0.00129, 0.00037,
         0.07984}
    Set oztransa[NSR_BANDS] = (ozone transmission coeff)
        {-0.00255649, -0.0177861, -0.0969872, -0.0611428, 0.0001, 0.0001,
          0.0001, -0.0834061}
    Set wvtransa[NSR_BANDS] = (water vapor transmission coeff)
        {2.29849e-27, 2.29849e-27, 0.00194772, 0.00404159, 0.000729136,
         0.00067324, 0.0177533, 0.00279738}
    Set wvtransb[NSR_BANDS] = (water vapor transmission coeff)
        {0.999742, 0.999742, 0.775024, 0.774482, 0.893085, 0.939669, 0.65094,
         0.759952}
    Set ogtransa1[NSR_BANDS] = (other gases transmission coeff)
        {4.91586e-20, 4.91586e-20, 4.91586e-20, 1.04801e-05, 1.35216e-05,
         0.0205425, 0.0256526, 0.000214329}
    Set ogtransb0[NSR_BANDS] = (other gases transmission coeff)
        {0.000197019, 0.000197019, 0.000197019, 0.640215, -0.195998, 0.326577,
         0.243961, 0.396322}
    Set ogtransb1[NSR_BANDS] = (other gases transmission coeff)
        {9.57011e-16, 9.57011e-16, 9.57011e-16, -0.348785, 0.275239, 0.0117192,
         0.0616101, 0.04728}

Setup input space to geolocation space mapping coefficients
Initialize the look up tables and atmospheric correction variables (see Init SR Vars Algorithm Description)

/* Apply first level of atmospheric corrections */
Loop through reflectance bands 1 through 7
    Call atmcorlamb2 for the current band to implement the Lambertian atmospheric correction (passing in the solar angles, observation angles, the SR constant variables, the look up tables, ozone, water vapor, and TOA reflectance).  lambertian SR (roslamb), other gaseous transmittance (tgo), atmospheric reflectance (roatm), total atmospheric transmission (ttatmg), atmosphere spherical albedo (satm), and molecular reflectance (xrorayp) are all returned.

    /* Save the band-related parameters for later use */
    Set btgo[band] = tgo
    Set broatm[band] = roatm
    Set bttatmg[band] = ttatmg
    Set bsatm[band] = satm

    /* Perform atmospheric corrections for bands 1-7 */
    Loop through each pixel in the band
        If the QA for this pixel isn't FILL
            /* Save the scaled TOA values for later use before adding
               atmospheric corrections */
            If (ib == DN_BAND1)
                Set aerob1[pixel] = sband[band][pixel]
            Else if (ib == DN_BAND2)
                Set aerob2[pixel] = sband[band][pixel]
            Else if (ib == DN_BAND4)
                Set aerob4[pixel] = sband[band][pixel]
            Else if (ib == DN_BAND5)
                Set aerob5[pixel] = sband[band][pixel]
            Else if (ib == DN_BAND7)
                Set aerob7[pixel] = sband[band][pixel]
            Endif

            /* Apply the atmospheric corrections, and store the scaled value
               for further corrections */
            Set rotoa = sband[band][pixel] * SCALE_FACTOR
            Set roslamb = rotoa / tgo
            Set roslamb = roslamb - roatm
            Set roslamb = roslamb / ttatmg
            Set roslamb = roslamb / (1.0 + satm * roslamb)
            Set sband[band][pixel] = (int) (roslamb * MULT_FACTOR)
        Endif
    Endloop
Endloop

Loop through reflectance bands 1 through 7
    Set erelc[band] = -1.0    (band ratio variable)
    Set troatm[band] = 0.0    (atmospheric reflectance table)
Endloop

/* Interpolate the auxiliary data for each pixel location in the scene */
Loop through all the lines in the scene
    Loop through all the samples in the scene
        Get the lat/long (in degrees) for the current pixel
        Determine the x/y location of the current pixel in the CMG-related lookup tables
        Interpolate the water vapor (twvi) and ozone (tozi) based on the weighted sum of the current and window water vapor and ozone CMG pixels.  For the ozone values, if any of the window values are 0 then set them to 120.
        Set twvi[pixel] = twvi[pixel] / 100.0
        Set tozi[pixel] = tozi[pixel] / 400.0
        Interpolate the atmospheric pressure (tp) based on the weighted sum of the current pixel and window CMG pixels. The pressure for any CMG pixel in the window is computed as follows:
            If the CMGDEM value for this CMG pixel is not fill
                Pressure[cmg_pixel] = 1013.0 * exp (-dem[cmg_pixel] / 8500.0)
            Else
                Pressure[cmg_pixel] = 1013.0 
            Endif

        /* Use the CMGDEM to flag potential water pixels */
        If the CMGDEM for the current CMG pixel or any of the surrounding pixels are fill then
            Set cloud_qa[pixel] = 128 (turns on water bit)
            Set tresi[pixel] = -1.0   (residual for the current pixel;
                                       tresi < 0 flags water pixels and pixels
                                       with high residuals)
        Endif
        
        /* Invert the aerosols, but filter the cirrus pixels */
        If (sband[BAND9][pixel] > 100.0 / tp[pixel] / 1013.0) then
            Add 1 to cloud_qa[pixel]  (sets the cirrus bit)
        Else
            /* Determine the band ratios */
            If ratiob1[cmg_pixel] == 0
                /* Average the valid ratio around the location */
                Set erelc[BAND1] = 0.4817
                Set erelc[BAND2] = erelc[DN_BAND1] / 0.844239
                Set erelc[BAND4] = 1.0
                Set erelc[BAND7] = 1.79
            Else
                /* Use the NDWI to calculate the band ratio */
                xndwi = (sband[BAND5][pixel] - (sband[BAND7][pixel] * 0.5)) /
                        (sband[BAND5][pixel] + (sband[BAND7][pixel] * 0.5))

                Set th1 = (andwi[cmg_pixel] + 2.0 * sndwi[cmg_pixel]) * 0.001
                Set th2 = (andwi[cmg_pixel] - 2.0 * sndwi[cmg_pixel]) * 0.001
                If (xndwi > th1) then Set xndwi = th1
                If (xndwi < th2) then Set xndwi = th2

                Set erelc[BAND1] = (xndwi * slpratiob1[cmg_pixel] +
                    intratiob1[cmg_pixel]) * 0.001
                Set erelc[BAND2] = (xndwi * slpratiob2[cmg_pixel] +
                    intratiob2[cmg_pixel]) * 0.001
                Set erelc[BAND4] = 1.0
                Set erelc[BAND7] = (xndwi * slpratiob7[cmg_pixel] +
                    intratiob7[cmg_pixel]) * 0.001
            Endif

            /* Retrieve the TOA reflectance values for the current pixel */
            Set troatm[BAND1] = aerob1[pixel] * SCALE_FACTOR
            Set troatm[BAND2] = aerob2[pixel] * SCALE_FACTOR
            Set troatm[BAND4] = aerob4[pixel] * SCALE_FACTOR
            Set troatm[BAND7] = aerob7[pixel] * SCALE_FACTOR

            If cloud_qa[pixel] has water bit set from the CMGDEM then
                Compute the NDVI from sband BAND5 and BAND4 to validate if this really is water
                If NDVI < 0.1 then
                    /* This is water */
                    Set taero[pixel] = 0.0     (aerosol value for current pixel)
                    Set tresi[pixel] = -0.01
                    Continue to next pixel
                Else
                    /* Remove the water designation */
                    Subtract the 128 water bit from cloud_qa[pixel]
                Endif
            Endif

            Call subaeroret with bands 4 and 1 to retrieve the aerosol information (passing in the band ratios, solar angles, observation angles, the SR constant variables, the look up tables, ozone, water vapor, and surface pressure).  AOT reflectance (raot) and model residual (residual) are both returned.
            Set corf = raot / xmus   (aerosol impact; higher values represent high aerosol)

            /* Test the quality of the aerosol inversion */
            If residual < (0.015 + 0.005 * corf) then
                /* Test the AOT reflectance with BAND5 and BAND4 */
                Call atmcorlamb2 for BAND5 to implement the Lambertian atmospheric correction (passing in the calculated raot as raot550nm, solar angles, observation angles, the SR constant variables, the look up tables, ozone, water vapor, and TOA reflectance for BAND5).  lambertian SR (roslamb), other gaseous transmittance (tgo), atmospheric reflectance (roatm), total atmospheric transmission (ttatmg), atmosphere spherical albedo (satm), and molecular reflectance (xrorayp) are all returned.
                Set ros5 = roslamb  (surface reflectance for band 5)

                Call atmcorlamb2 for BAND4 to implement the Lambertian atmospheric correction (passing in the calculated raot as raot550nm, solar angles, observation angles, the SR constant variables, the look up tables, ozone, water vapor, and TOA reflectance for BAND4).  lambertian SR (roslamb), other gaseous transmittance (tgo), atmospheric reflectance (roatm), total atmospheric transmission (ttatmg), spherical atmosphere albedo (satm), and molecular reflectance (xrorayp) are all returned.
                Set ros4 = roslamb  (surface reflectance for band 4)

                If (ros5 > 0.1) and ((ros5 - ros4) / (ros5 + ros4) > 0) then
                    Set taero[pixel] = raot
                    Set tresi[pixel] = residual
                Else
                    Set taero[pixel] = 0.0
                    Set tresi[pixel] = -0.01
                Endif 
            Else
                Set taero[pixel] = 0.0
                Set tresi[pixel] = -0.01
            Endif
        Endif
    Endloop
Endloop

/* Refine the cloud mask using the average temperature of the clear pixels versus the average temperature of all the pixels */
Set nbval = 0
Set nbclear = 0
Set mclear = 0.0
Set mall = 0.0
Loop through each pixel in the QA band
    If the QA for this pixel isn't FILL
        /* Keep track of the number of total (non-fill) pixels in addition
           to the sum of the unscaled thermal values */
        Increment nbval by 1
        Set mall = mall + sband[BAND10][pixel] * SCALE_FACTOR_TH

        /* Check for clear pixels */
        If (cloud_qa[pixel] is not cirrus) and (sband[BAND5][pixel] > 300) then
            Set anom = sband[BAND2][pixel] - sband[BAND4][pixel] * 0.5
            If (anom < 300) then
                /* Keep track of the number of clear pixels in addition to
                   the sum of the unscaled thermal values */
                Increment nbclear by 1
                Set mclear = mclear + sband[BAND10][pixel] * SCALE_FACTOR_TH
            Endif
        Endif
    Endif
Endloop
If (nbclear > 0) then
    /* Compute the average clear temp */
    Set mclear = mclear / nbclear
Else
    Set mclear = 275.0
Endif
If (nbval > 0) then
    /* Compute the average/mean temperature of the clear pixels */
    Set mall = mall / nbval    (NOTE: not used except to print for info sake)
Endif
Loop through each pixel in the reflectance bands
    If (tresi[pixel] < 0.0) then
        If (((sband[BAND2][pixel] - sband[BAND4][pixel] * 0.5) > 500) and
            ((sband[BAND10][pixel] * SCALE_FACTOR_TH) < (mclear - 2.0))) then
            Set the cloud bit for cloud_qa[pixel]
        Endif
    Endif
Endloop

/* Set up the adjacent to something (cloud or cirrus) bit */
Loop through each pixel in the reflectance bands
    If cloud_qa[pixel] is cloud or cirrus then
        Loop through a 5x5 window around the current pixel
            If cloud_qa[window_pixel] is not already cloud or cirrus then
                Set the adjacent to cloud bit for cloud_qa[window_pixel]
            Endif
        Endloop
    Endif
Endloop

/* Compute the cloud shadow */
Set facl = cos(xfs) * tan(xts) / pixsize  (factor for lines)
Set fack = sin(xfs) * tan(xts) / pixsize  (factor for samples)
Loop through each line in the QA band
    Loop through each sample in the QA band
        If cloud_qa[pixel] is cloud or cirrus then
            Set tcloud = sband[BAND10][pixel] * SCALE_FACTOR_TH
            Set cldh = (mclear - tcloud) * 1000.0 / cfac
            If (cldh < 0.0) then Set cldh = 0.0
            Set cldhmin = cldh - 1000.0
            Set cldhmax = cldh + 1000.0
            Set mband5 = 9999
            Set mband5k = -9999
            Set mband5l = -9999
            If (cldhmin < 0) then Set cldhmin = 0.0
            Loop from (icldh = cldhmin * 0.1) until (icldh <= cldhmax * 0.1) by steps of 1 for icldh
                Set cldh = icldh * 10.0
                Set k = i + facl * cldh    (lines)
                Set l = j - fack * cldh    (samps)

                /* Make sure the line and sample is valid */
                If (k < 0 or k >= nlines or l < 0 or l >= nsamps) then
                    Skip to the next icldh

                Set win_pix = k * nsamps + l
                If ((sband[BAND6][win_pix] < 800) and
                   ((sband[BAND3][win_pix] - sband[BAND4][win_pix]) < 100)) then
                    If cloud_qa[win_pix] is cloud or cirrus or cloud shadow then
                        Skip to the next icldh
                    Else
                        /* Store value of BAND6 as well as l and k values */
                        If (sband[BAND6][win_pix] < mband5) then
                            Set mband5 = sband[BAND6][win_pix]
                            Set mband5k = k
                            Set mband5l = l
                        Endif
                    Endif
                Endif
            Endloop

            If (mband5 < 9999) then
                Set cloud shadow bit for cloud_qa[mband5k*nsamps + mband5l]
            Endif
        Endif
    Endloop
Endloop

/* Expand the cloud shadow using the residual */
Loop through each pixel in the QA band
    If cloud_qa[pixel] is cloud shadow then
        Loop through a 6x6 window around the current pixel
            If cloud_qa[window_pixel] is cloud or cloud shadow then
                Skip to the next pixel in the window
            Else
                If cloud_qa[window_pixel] has the temporary bit set then
                    Skip to the next pixel in the window
                Else 
                    If tresi[window_pixel] < 0 then
                        Set the temporary bit for cloud_qa[window_pixel]
                    Endif
                Endif
            Endif
        Endloop
    Endif
Endloop

/* Update the cloud shadow */
Loop through each pixel in the QA band
    If the temporary bit was set for cloud_qa[pixel] then
        Remove the temporary bit and set the cloud shadow bit for cloud_qa[pixel]
    Endif
Endloop

/* Handle the aerosol interpolation. Does not use water, cloud, or cirrus
   pixels. */
Initialize hole = true
Initialize step = 10
Loop while hole is true and step < 1000
    Set hole = false
    Loop through the lines in the reflectance band, by increments of step
        Loop through the samples in the reflectance band, by increments of step
            Initialize the nbaot = 0 (number of AOT)
            Initialize the aaot = 0.0 (average of AOT)
            Initialize the sresi = 0.0 (1/residual)
            Loop through the line, sample to line+step-1, sample+step-1 window one pixel at a time
                If tresi[window_pixel] > 0 and cloud_qa is clear then
                    Increment the nbaot count by one
                    Add taero[window_pixel] / tresi[window_pixel] to aaot
                    Add 1 / tresi[window_pixel] to sresi
                Endif
            Endloop

            /* Check to see if any clear pixels with positive residuals were
               found in the window */
            If nbaot != 0 then
                Set aaot = aaot / sresi
                Loop through the same line, sample to line+step-1, sample+step-1 window as above, one pixel at a time
                    If tresi[window_pixel] < 0 and cloud_qa[window_pixel] is not cirrus, not cloud, and not water then
                        Set taero[window_pixel] = aaot
                        Set tresi[window_pixel] = 1.0
                    Endif
                Endloop
            Else
                Set hole = true
            Endif
        Endloop
    Endloop

    Increase step = step * 2
Endloop

/* Perform the second level of atmospheric correction for the aerosols. This
   correction is not applied to water, cirrus, or cloud pixels. */
Loop through reflectance bands 1 through 7
    Loop through each pixel in the band
        If the QA for this pixel isn't FILL then
            If tresi[pixel] > 0 and cloud_qa[pixel] is not cirrus and is not cloud then
                Set rsurf = sband[band][pixel] * SCALE_FACTOR
                Set rotoa = (rsurf * bttatmg[band] / (1.0 - bsatm[band] * rsurf)
                    + broatm[band]) * btgo[band]
                Set raot550nm = taero[pixel]
                Set pres = tp[pixel]
                Set uwv = twvi[pixel]
                Set uoz = tozi[pixel]
                Call atmcorlamb2 for the current band to implement the Lambertian atmospheric correction (passing in the solar angles, observation angles, the SR constant variables, the look up tables, ozone, water vapor, and TOA reflectance).  lambertian SR (roslamb), other gaseous transmittance (tgo), atmospheric reflectance (roatm), total atmospheric transmission (ttatmg), spherical atmosphere albedo (satm), and molecular reflectance (xrorayp) are all returned.

                /* Handle the aerosol computation in the cloud mask if this
                   is the coastal aerosol band */
                If band is BAND1
                    If roslamb < -0.005 then
                        Set taero[pixel] = 0.05
                        Set roat550nm = 0.05
                        Call atmcorlamb2 for the current band to implement the Lambertian atmospheric correction (passing in the solar angles, observation angles, the SR constant variables, the look up tables, ozone, water vapor, and TOA reflectance).  lambertian SR (roslamb), other gaseous transmittance (tgo), atmospheric reflectance (roatm), total atmospheric transmission (ttatmg), spherical atmosphere albedo (satm), and molecular reflectance (xrorayp) are all returned.
                    Else
                        If absolute value (rsurf - roslamb) <= 0.015 then
                            Set first aerosol bit for cloud_qa[pixel] (designates low aerosols)
                        Else
                            If absolute value (rsurf - roslamb) <= 0.03 then
                                Set second aerosol bit for cloud_qa[pixel] (designates average aerosols)
                            Else
                                Set both aerosol bits for cloud_qa[pixel] (designates high aerosols)
                            Endif
                        Endif
                    Endif
                Endif

                /* Save the scaled surface reflectance value (sband), but make
                   sure it falls within the defined valid range */
                Set roslamb = roslamb * MULT_FACTOR
                If roslamb < MIN_VALID
                    Set sband[band][pixel] = MIN_VALID
                Else if roslamb > MAX_VALID
                    Set sband[band][pixel] = MAX_VALID
                Else
                    Set sband[band][pixel] = (int) roslamb
                Endif
            Endif
        Endif
    Endloop
Endloop


Init SR Vars Algorithm Description
==================================
Inputs:
    Scene metadata
    Input auxiliary filenames

Outputs:
    xtv,         (observation zenith angle (deg))
    xmuv,        (cosine of observation zenith angle)
    xfi,         (azimuthal difference between sun and observation (deg))
    cosxfi,      (cosine of azimuthal difference)
    raot550nm,   (nearest value of AOT)
    pres,        (surface pressure)
    uoz,         (total column ozone)
    uwv,         (total column water vapor (precipital water vapor))
    xtsstep,     (solar zenith step value)
    xtsmin,      (minimum solar zenith value)
    xtvstep,     (observation step value)
    xtvmin,      (minimum observation value)
    tsmax,       (maximum scattering angle table)
    tsmin,       (minimum scattering angle table)
    tts,         (sun angle table)
    ttv,         (view angle table)
    indts,       (index for the sun angle table)
    rolutt,      (intrinsic reflectance table)
    transt,      (transmission table)
    sphalbt,     (atmosphere spherical albedo table)
    normext,     (aerosol extinction coefficient at the current wavelength (normalized at 550nm))
    nbfic,       (communitive number of azimuth angles)
    nbfi,        (number of azimuth angles)
    dem,         (CMG DEM data array)
    andwi,       (avg NDWI)
    sndwi,       (standard NDWI)
    ratiob1,     (mean band1 ratio)
    ratiob2,     (mean band2 ratio)
    ratiob7,     (mean band7 ratio)
    intratiob1,  (integer band1 ratio)
    intratiob2,  (integer band2 ratio)
    intratiob7,  (integer band7 ratio)
    slpratiob1,  (slope band1 ratio)
    slpratiob2,  (slope band2 ratio)
    slpratiob7,  (slope band7 ratio)
    wv,          (water vapor values)
    oz           (ozone values)

    Set xtv = 0.0  (NOTE: this is static and remains 0.0)
    Set xmuv = cosine (xtv)
    Set xfi = 0.0
    Set cosxfi = cosine (xfi)
    Set xtsmin = 0
    Set xtsstep = 4.0
    Set xtvmin = 2.84090
    Set xtvstep = 6.52107 - xtvmin

    Read the {clean case v2.0} look-up tables for tsmax, tsmin, ttv, tts, nbfi, nbfic, indts, rolutt, transt, sphalbt, and normext
    Read the auxiliary data files for dem, andwi, sndwi, ratiob1, ratiob2, ratiob7, intratiob1, intratiob2, intratiob7, slpratiob1, slprationb2, slpratiob7, wv, and oz

    Set raot550nm = 0.12
    Set pres = 1013.0         (atmospheric pressure at sea level)
    Set uoz = 0.30
    Set uwv = 0.5

    /* Use the scene center lat/long to determine the line/sample in the CMG-related lookup tables */
    Get the lat/long for the scene center line/sample
    Set ycmg = (89.975 - center_lat) / 0.05
    Set xcmg = (179.975 + center_lon) / 0.05
    Set lcmg = (int) (ycmg + 0.5)
    Set scmg = (int) (xcmg + 0.5)

    If wv[lcmg][scmg] != 0 then
        Set uwv = wv[lcmg][scmg] / 200.0
    Else
        Set uwv = 0.5
    Endif

    If oz[lcmg][scmg] != 0 then
        Set uoz = oz[lcmg][scmg] / 400.0
    Else
        Set uoz = 0.3
    Endif

    If dem[lcmg][scmg] is not FILL
        Set pres = 1013.0 * exp (-dem[lcmg][scmg] / 8500)
    Else
        Set pres = 1013.0
    Endif
    Set raot550nm = 0.05

