Close

Service Interruption

Atmospheric and Anisotropic Reflectance Correction Source Code Document



Procedures in this file are for atmospheric and anisotropic reflectance correction of the AVHRR shortwave channel. They are:

call_narr2bb_toa
call_anis_correct
call_anis_correct_blend
call_atmos_correct
narr2bb_toa
anis_correct
anis_correct_blend
atmos_correct

Compilation: See the file README.C for instructions on compiling this file.



#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <math.h>
#include "globals.h"
#include "anisdata.h"

long call_narr2bb_toa(short argc, void *argv[]);
long call_anis_correct(short argc, void *argv[]);
long call_anis_correct_blend(short argc, void *argv[]);
long call_atmos_correct(short argc, void *argv[]);
float narr2bb_toa(float ch1, float ch2, float solzen, char surftype);
float anis_correct(float broadalb, float solzen, float scan, float relaz, char surftype);
float anis_correct_blend(float ch1, float ch1_pure_ocean, float ch1_pure_ice, float broadalb, float solzen, float scan, float relaz, char surftype);
float atmos_correct(float broadalb, float solzen, float wv, float odhaze, char surftype);

/******************** CALL NARROW TO BROADBAND CONVERSION ******************/

long call_narr2bb_toa(short argc, void *argv[])

/*--------------------------------------------------------------------------

This function loops through the input arrays and calls 'surfalb_clear' (below) for each pixel. It is called from an IDL program with the following statement:

ERROR = call_external(objfln,'call_narr2bb_toa',ncols,nrows,rsurf1,rsurf2,$ solzen,surftype,cmask,badpix,broadalbtoa)

Input (as pointers to the specified type):

ncols, nrows : Dimensions of 'tsurf' and the arrays pointed to by 'array' [short integer]. While these arrays are 2D, they are treated here as 1D. If scalar values are passed in, both of these values must be 1. If a 1D array is passed in, then 'nrows' must be 1.

ch1, ch2 : Channels 1 and 2 albedos at TOA, in the range [0,1] and normalized for solar zenith angle and NOT corrected for anisotropy and the atmosphere. [float]

solzen : Solar zenith angle image (degrees) [float]

surftype : Surface type (global variables LAND, OCEAN, ICE, SNOW; see 'globals.pro' for values). LAND means snow-free land (vegetation cover), OCEAN means open water, ICE means sea ice (with or without snow), and SNOW means snow-covered land. However, corrections are currently done for only land and ice surface. For other surface types the channel 1 reflectance is returned. [char]

cmask : Cloud mask array, cloudy=CLOUDY, clear=CLEAR ('cmask.pro' and 'globals.h') [char]

badpix : Valid pixels array, where invalid=TRUE ('globals.h') [char]

broadalbtoa : The albedo array. The values are unimportant, but this array must be defined upon entry. [float]

Output:

broadalbtoa: Broadband TOA albedo in the range [0,1]. For surface types other than snow/ice or vegetation the clear sky albedo is returned. [fltarr or float]. If badpix=TRUE and/or the solar zenith angle is > DARKZEN, MISSING values are returned. For CLOUDY pixels the input values are unchanged. Results outside of the range [0,1] are reset to 0 or 1.

Additional Variables (from globals.h):

TRUE, FALSE : true and false flags to be used for returning error condition

Called by: narr2bb_toa (an IDL procedure)

Calls: narr2bb_toa (the C function below)

----------------------------------------------------------------------------*/
{

/* parameters passed in from IDL */

short ncols, nrows;
char *cmask, *surftype, *badpix;
float *ch1, *ch2, *solzen;
float *broadalbtoa;

/* local parameters */

long i;
long ERROR; /* Error flag, returned to IDL */

/* Check that correct number of parameters was passed */

if (argc != 9) {
fprintf(stdout, "argc = %d\r\n",argc);
fprintf(stdout, "ERROR >>> CALL_NARR2BB_TOA: INCORRECT NUMBER OF ARGUMENTS\r\n");
ERROR = TRUE;
return ERROR;
}

/* Cast passed parameters to local vars */

ncols = *(short *) argv[0];
nrows = *(short *) argv[1];
ch1 = (float *) argv[2];
ch2 = (float *) argv[3];
solzen = (float *) argv[4];
surftype = (char *) argv[5];
cmask = (char *) argv[6];
badpix = (char *) argv[7];
broadalbtoa = (float *) argv[8];

/* Loop through each pixel in the arrays. If it is a "bad" pixel, skip it and assign the fluxes a MISSING value. */

ERROR = FALSE;

for (i=0; i<nrows*ncols; i++) {
if (badpix[i] == (char) FALSE && solzen[i] <= DARKZEN) {
if (cmask[i] == (char) CLEAR)
broadalbtoa[i] = narr2bb_toa(ch1[i],ch2[i],solzen[i],surftype[i]);
}
else
broadalbtoa[i] = MISSING;
}

return ERROR;
}

/***************** CALL ANISOTROPIC REFLECTANCE CORRECTION *****************/

long call_anis_correct(short argc, void *argv[])

/*--------------------------------------------------------------------------

This function loops through the input arrays and calls 'anis_correct' (below) for each pixel. It is called from an IDL program with the following statement:

ERROR = call_external(objfln,'call_anis_correct',ncols,nrows,broadalb,$ solzen,scan,relaz,surftype,cmask,badpix)

Input (as pointers to the specified type):

ncols, nrows : Dimensions of 'tsurf' and the arrays pointed to by 'array' [short integer]. While these arrays are 2D, they are treated here as 1D. If scalar values are passed in, both of these values must be 1. If a 1D array is passed in, then 'nrows' must be 1.

broadalb : Broadband albedo at TOA, in the range [0,1]. [float]

solzen : Solar zenith angle image (degrees) [float]

scan : Satellite scan angle (degrees) [float]

relaz : Satellite-sun relative azimuth angle (degrees), where 180 is looking into the sun and 0 is looking away (0..180) [float]

surftype : Surface type (global variables LAND, OCEAN, ICE, SNOW; see 'globals.pro' for values). LAND means snow-free land (vegetation cover), OCEAN means open water, ICE means sea ice (with or without snow), and SNOW means snow-covered land. However, corrections are currently done for only land and ice surface. Any other surface type uses the ice procedure. [char]

cmask : Cloud mask array, cloudy=CLOUDY, clear=CLEAR ('globals.h') [char]

badpix : Valid pixels array, where invalid=TRUE ('globals.h') [char]

Output:

broadalb : The broadband albedo corrected for anisotropy, in the range [0,1]. If badpix=TRUE and/or the solar zenith angle is > DARKZEN, MISSING values are returned. For CLOUDY pixels the input values are unchanged. Results outside of the range [0,1] are reset to 0 or 1. [fltarr or float]

Additional Variables (from globals.h):

TRUE, FALSE: true and false flags to be used for returning error condition

Called by: anis_correct (an IDL procedure)

Calls: anis_correct (the C function below)

----------------------------------------------------------------------------*/
{

/* parameters passed in from IDL */

short ncols, nrows;
char *cmask, *surftype, *badpix;
float *broadalb, *solzen, *scanang, *relaz;

/* local parameters */

long i;
long ERROR; /* Error flag, returned to IDL */
short j, k, l, m, index;

/* Check that correct number of parameters was passed */

if (argc != 9) {
fprintf(stdout, "argc = %d\r\n",argc);
fprintf(stdout, "ERROR >>> CALL_ANIS_CORRECT: INCORRECT NUMBER OF ARGUMENTS\r\n");
ERROR = TRUE;
return ERROR;
}

/* Cast passed parameters to local vars. Copy the data table to fixed size
arrays. */

ncols = *(short *) argv[0];
nrows = *(short *) argv[1];
broadalb = (float *) argv[2];
solzen = (float *) argv[3];
scanang = (float *) argv[4];
relaz = (float *) argv[5];
surftype = (char *) argv[6];
cmask = (char *) argv[7];
badpix = (char *) argv[8];

/* Loop through each pixel in the arrays. If it is a "bad" pixel, skip it and assign the fluxes a MISSING value. */

ERROR = FALSE;

for (i=0; i<nrows*ncols; i++) {
if (badpix[i] == (char) FALSE && solzen[i] <= DARKZEN) {
if (cmask[i] == (char) CLEAR)
broadalb[i] = anis_correct(broadalb[i],solzen[i],scanang[i],
relaz[i],surftype[i]);
}
else {
broadalb[i] = MISSING;
}
}

return ERROR;
}

/******CALL ANISOTROPIC REFLECTANCE CORRECTION WITH BLENDED OCEAN/ICE *******/

long call_anis_correct_blend(short argc, void *argv[])

/*--------------------------------------------------------------------------

This function loops through the input arrays and calls 'anis_correct_blend' (below) for each pixel. It is called from an IDL program with the following statement:

ERROR = call_external(objfln,'call_anis_correct_blend',ncols,nrows,$ ch1, ch1_pure_ocean, ch1_pure_ice,$broadalb,solzen,scan,relaz,surftype,cmask,badpix)

Input (as pointers to the specified type):

ncols, nrows : Dimensions of 'tsurf' and the arrays pointed to by 'array' [short integer]. While these arrays are 2D, they are treated here as 1D. If scalar values are passed in, both of these values must be 1. If a 1D array is passed in, then 'nrows' must be 1.

ch1 : Channels 1 albedos at TOA, in the range [0,1] and normalized for solar zenith angle and NOT corrected for anisotropy and the atmosphere. [float]

ch1_pure_ocean: The ch1 value which will be considered pure ocean. [scalar float]

ch1_pure_ice : The ch1 value which will be considered pure ocean. [scalar float]

broadalb : Broadband albedo at TOA, in the range [0,1]. [float]

solzen : Solar zenith angle image (degrees) [float]

scan : Satellite scan angle (degrees) [float]

relaz : Satellite-sun relative azimuth angle (degrees), where 180 is looking into the sun and 0 is looking away (0..180) [float]

surftype : Surface type (global variables LAND, OCEAN, ICE, SNOW; see 'globals.pro' for values). LAND means snow-free land (vegetation cover), OCEAN means open water, ICE means sea ice (with or without snow), and SNOW means snow-covered land. However, corrections are currently done for only land and ice surface. Any other surface type uses the ice procedure. [char]

cmask : Cloud mask array, cloudy=CLOUDY, clear=CLEAR ('globals.h') [char]

badpix : Valid pixels array, where invalid=TRUE ('globals.h') [char]

Output:

broadalb : The broadband albedo corrected for anisotropy with blending of albedos for OCEAN and ICE pixels, in the range [0,1]. If badpix=TRUE and/or the solar zenith angle is > DARKZEN, MISSING values are returned. For CLOUDY pixels the input values are unchanged. Results outside of the range [0,1] are reset to 0 or 1. [fltarr or float]

Additional Variables (from globals.h):

TRUE, FALSE : true and false flags to be used for returning error condition

Called by: anis_correct_blend (an IDL procedure)

Calls: anis_correct_blend (the C function below)

----------------------------------------------------------------------------*/


{

/* parameters passed in from IDL */

short ncols, nrows;
char *cmask, *surftype, *badpix;
float *ch1, ch1_pure_ocean, ch1_pure_ice;
float *broadalb, *solzen, *scanang, *relaz;

/* local parameters */

long i;
long ERROR; /* Error flag, returned to IDL */
short j, k, l, m, index;

/* Check that correct number of parameters was passed */

if (argc != 12) {
fprintf(stdout, "argc = %d\r\n",argc);
fprintf(stdout, "ERROR >>> CALL_ANIS_CORRECT_BLEND: INCORRECT NUMBER OF ARGUMENTS\r\n");
ERROR = TRUE;
return ERROR;
}

/* Cast passed parameters to local vars. Copy the data table to fixed size
arrays. */

ncols = *(short *) argv[0];
nrows = *(short *) argv[1];
ch1 = (float *) argv[2];
ch1_pure_ocean = *(float *) argv[3];
ch1_pure_ice = *(float *) argv[4];
broadalb = (float *) argv[5];
solzen = (float *) argv[6];
scanang = (float *) argv[7];
relaz = (float *) argv[8];
surftype = (char *) argv[9];
cmask = (char *) argv[10];
badpix = (char *) argv[11];

/* Loop through each pixel in the arrays. If it is a "bad" pixel, skip
it and assign the fluxes a MISSING value. */

ERROR = FALSE;

for (i=0; i<nrows*ncols; i++) {
if (badpix[i] == (char) FALSE && solzen[i] <= DARKZEN) {
if (cmask[i] == (char) CLEAR)
broadalb[i] = anis_correct_blend(ch1[i],ch1_pure_ocean,ch1_pure_ice,
broadalb[i],solzen[i],scanang[i],
relaz[i],surftype[i]);
}
else {
broadalb[i] = MISSING;
}
}

return ERROR;
}

/************************ CALL ATMOSPHERIC CORRECTION **********************/

long call_atmos_correct(short argc, void *argv[])

/*--------------------------------------------------------------------------

This function loops through the input arrays and calls 'atmos_correct' (below) for each pixel. It is called from an IDL program with the following statement:

ERROR = call_external(objfln,'call_atmos_correct',ncols,nrows,broadalb,$ solzen,wv,odhaze,surftype,cmask,badpix)

Input (as pointers to the specified type):

ncols, nrows : Dimensions of 'tsurf' and the arrays pointed to by 'array' [short integer]. While these arrays are 2D, they are treated here as 1D. If scalar values are passed in, both of these values must be 1. If a 1D array is passed in, then 'nrows' must be 1.

broadalb : Broadband TOA albedo, in the range [0,1]. [float]

solzen : Solar zenith angle image (degrees) [float]

wv: Water vapor amount image (0.5 - 5.0 g/cm^2 = precipitable water in cm) [float]

odhaze : Aerosol optical depth (0.05 - 0.5) [float]

surftype : Surface type (global variables LAND, OCEAN, ICE, SNOW; see 'globals.pro' for values). LAND means snow-free land (vegetation cover), OCEAN means open water, ICE means sea ice (with or without snow), and SNOW means snow-covered land. However, corrections are currently done for only land and ice surface. Any other surface type uses the ice procedure. [char]

cmask : Cloud mask array, cloudy=CLOUDY, clear=CLEAR ('cmask.pro' and 'globals.h') [char]

badpix : Valid pixels array, where invalid=TRUE ('globals.h') [char]

Output:

broadalb : "Apparent" surface albedo corrected for the atmosphere in the range [0,1], which are directional, hemispherical (i.e., are not scan angle or azimuthally-dependent but do vary with the atmospheric conditions, hence "apparent"). If badpix=TRUE or the solar zenith angle > DARKZEN, MISSING values are returned. If the pixel is CLOUDY then the values that these variables have on input are unchanged. Results outside of the range [0,1] are reset to 0 or 1. [fltarr or float]

Additional Variables (from globals.h):

TRUE, FALSE : true and false flags to be used for returning error condition

Called by: atmos_correct (an IDL procedure)

Calls: atmos_correct (the C function below)

----------------------------------------------------------------------------*/


{

/* parameters passed in from IDL */

short ncols, nrows;
char *surftype, *cmask, *badpix;
float *solzen, *wv, odhaze;
float *broadalb;

/* local parameters */

long i;
long ERROR; /* Error flag, returned to IDL */

/* Check that correct number of parameters was passed */

if (argc != 9) {
fprintf(stdout, "argc = %d\r\n",argc);
fprintf(stdout, "ERROR >>> CALL_ATMOS_CORRECT: INCORRECT NUMBER OF ARGUMENTS\r\n");
ERROR = TRUE;
return ERROR;
}

/* Cast passed parameters to local vars */

ncols = *(short *) argv[0];
nrows = *(short *) argv[1];
broadalb = (float *) argv[2];
solzen = (float *) argv[3];
wv = (float *) argv[4];
odhaze = *(float *) argv[5];
surftype = (char *) argv[6];
cmask = (char *) argv[7];
badpix = (char *) argv[8];

/* Loop through each pixel in the arrays. If it is a "bad" pixel, skip
it and assign the fluxes a MISSING value. */

ERROR = FALSE;

for (i=0; i<nrows*ncols; i++) {
if (badpix[i] == (char) FALSE && solzen[i] <= DARKZEN) {
if (cmask[i] == (char) CLEAR)
broadalb[i] = atmos_correct(broadalb[i],solzen[i],wv[i],odhaze,
surftype[i]);
}
else {
broadalb[i] = MISSING;
}
}

return ERROR;
}

/*********************** NARROW TO BROADBAND CONVERSION ********************/

float narr2bb_toa(float ch1, float ch2, float solzen, char surftype)

/*--------------------------------------------------------------------------

Converts the channels 1 and 2 reflectance/albedo values to a broadband albedo for clear conditions over snow/ice, vegetation (generic), and open water. This procedure operates on scalar values only.

Input:

ch1, ch2 : Channels 1 and 2 albedos at TOA, in the range [0,1], normalized for solar zenith angle, NOT corrected for the atmospheric attenuation, and corrected for anisotropic reflectance. [float]

solzen : The solar zenith angle (degrees) [float]

surftype : One byte (char) value representing the surface type (global variables LAND, OCEAN, ICE, SNOW* see 'globals.pro' for values). LAND means snow-free land (vegetation cover), OCEAN means open water, ICE means sea ice (with or without snow), and SNOW means snow-covered land. However, corrections are currently done for only land and ice surface. For other surface types the channel 1 reflectance is returned. [char]

Output as value of function:

Broadband TOA albedo in the range [0,1]. For surface types other than snow/ice or vegetation the clear sky albedo is returned. [fltarr or float] If badpix=TRUE and/or the solar zenith angle is > DARKZEN, MISSING values are returned. For CLOUDY pixels the input values are unchanged. Results outside of the range [0,1] are reset to 0 or 1.

Algorithm Summary:

The clear sky broadband surface albedo is computed as a function of the channels 1 and 2 reflectances. The regression relationship has the form:

broadalbtoa = a + b*rsurf1 + c*rsurf2

Two sets of coefficients were developed: one for snow-free land (generic vegetation) and one for snow/ice. Streamer was used to simulate AVHRR channels 1 and 2 albedos and a TOA broadband albedo for a variety of surface albedo, atmospheric, and viewing angle relationships. A linear regression procedures was then applied to the Streamer data. For open water we assume that the broadband albedo is the same as the channel 1 albedo (to a certain extent after Brieglieb et al. JCAM 1986, where the albedo is the same for all wavelengths).

----------------------------------------------------------------------------*/


{
float broadalbtoa;

/* If it's dark then the input reflectances were previously set to zero (in 'correct12'). Do the same with the broadband albedo.
Otherwise, compute it. */

if (ch1 == MISSING || ch2 == MISSING || solzen >= DARKZEN)
broadalbtoa = MISSING;

else if (surftype == (char) LAND)
broadalbtoa = 0.0404522 + 0.545025*ch1 + 0.299113*ch2;

else if ((surftype == (char) SNOW) || (surftype == (char) ICE))
broadalbtoa = 0.0215773 + 0.277479*ch1 + 0.506755*ch2;

else
broadalbtoa = ch1; /* Open water */

if ((broadalbtoa < 0.0) && (broadalbtoa != MISSING)) broadalbtoa = 0.0;

return broadalbtoa;
}

/************************ ANISOTROPIC CORRECTION *****************************/

float anis_correct(float broadalb, float solzen, float scan, float relaz, char surftype)

/*--------------------------------------------------------------------------

Corrects the broadband TOA albedo data for anisotropic reflectance. This procedure operates on scalar values only.

Input:

broadalb : TOA broadband albedo, in the range [0,1]. [float]

solzen : Solar zenith angle image (degrees) [float]

scan : Satellite scan angle (degrees) [float]

relaz : Satellite-sun relative azimuth angle (degrees), where 180 is looking into the sun and 0 is looking away (0..180) [float]

surftype : One byte (char) value representing the surface type (global variables LAND, OCEAN, ICE, SNOW see 'globals.h' for values). LAND means snow-free land (vegetation cover), OCEAN means open water, ICE means sea ice (with or without snow), and SNOW means snow-covered land. However, corrections are currently done for only snow-free land, ocean, and snow/ice (the same) surfaces. Any other surface type uses the snow procedure.

Static data in anisdata.h:

rmatrx : Anisotropic reflectance values (not quite factors) based on ERBE and GOES, dimension: (3,10,7,8) for eight relative azimuth bins, seven viewing zenith bins, ten solar zenith bins, and three scene types. [float]

albmn : Normalizing factors for 'albmn', dimension: (3,10), for ten solar zenith bins and three scene types. [float]

Output as value of function:

The broadband albedo corrected for anisotropy, in the range [0,1]. If badpix=TRUE and/or the solar zenith angle is > DARKZEN, MISSING values are returned. For CLOUDY pixels the input values are unchanged. Results outside of the range [0,1] are reset to 0 or 1. [fltarr or float]

Algorithm Summary:

The anisotropic reflectance factors are based on ERBE and GOES data for clear ice (clear sky over ice), clear land, and clear ocean as given in NASA Reference Publication 1184, July 1988, "Angular Radiation Models for Earth-Atmosphere System, Volume I-Shortwave Radiation", by Suttles et al. The anisotropic reflectance factors are obtained via interpolation in the table, and are divided into the TOA reflectance to produce a TOA albedo. Data were provided by L. Smith (NASA/Langley).

----------------------------------------------------------------------------*/
{
short ISCENE, ITEMP, IVZEN, ISZEN, IRAZ;
float VZEN, CSUN, RAZ, FVZEN, CVZEN, FSZEN, CSZEN, FRAZ, CRAZ, TEMP, F, bbalb;

float rearth = 6357.0; /* Radius of the earth */
float rehsat = 7207.0; /* Rearth + Hsat(850 km) */

short MVZEN[] = {0,0,0,0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,5,5,5};
float GVZEN[] = {0.,21.,33.,45.,57.,69.,90};
short MSZEN[] = {0,0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,8};
float GSZEN[] = {0.,.15,.25,.35,.45,.55,.65,.75,.85,1.00};
short MRAZ[] = {0,0,0,0,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,5,5,5,5,5,6,6,6,6};
float GRAZ[] = {0.,20.,45.,75.,105.,135.,160.,180.};

/* If it's dark or if the input value is missing, return MISSING. */

if (solzen >= DARKZEN || broadalb == MISSING) return MISSING;

/* Convert values to what the interpolation code expects: VZEN is viewing
zenith angle (0-90 deg), CSUN is the solar zenith angle cosine (0-1),
RAZ is the relative azimuth angle (0-180 deg, 0 is looking into the sun
as in Streamer, but opposite of CASPR input), ISCENE is the scene type
(0=clear ocean, 1=clear land, 2=clear snow/ice). */

if (surftype == (char) OCEAN)
ISCENE = 0;
else if (surftype == (char) LAND)
ISCENE = 1;
else
ISCENE = 2;
RAZ = 180.0 - relaz;
VZEN = asin((sin(scan*DEGRAD)*rehsat)/rearth)*RADDEG;
CSUN = cos(solzen*DEGRAD);

/* Perform tri-linear interpolation in the table. First determine the indices
and fractions. Then linearly interpolate. */

ITEMP = VZEN/3.000001;
IVZEN = MVZEN[ITEMP];
FVZEN = (VZEN - GVZEN[IVZEN]) / (GVZEN[IVZEN+1] - GVZEN[IVZEN]);
CVZEN = 1. - FVZEN;

TEMP = 1. - CSUN;
ITEMP = TEMP/.05000001;
ISZEN = MSZEN[ITEMP];
FSZEN = (TEMP - GSZEN[ISZEN]) / (GSZEN[ISZEN+1] - GSZEN[ISZEN]);
CSZEN = 1. - FSZEN;

ITEMP = RAZ/5.000001;
IRAZ = MRAZ[ITEMP];
FRAZ = (RAZ - GRAZ[IRAZ]) / (GRAZ[IRAZ+1] - GRAZ[IRAZ]);
CRAZ = 1. - FRAZ;

F = (CVZEN*CRAZ*rmatrx[ISCENE][ISZEN][IVZEN][IRAZ]
+ FVZEN*CRAZ*rmatrx[ISCENE][ISZEN][IVZEN+1][IRAZ]
+ CVZEN*FRAZ*rmatrx[ISCENE][ISZEN][IVZEN][IRAZ+1]
+ FVZEN*FRAZ*rmatrx[ISCENE][ISZEN][IVZEN+1][IRAZ+1] )
* CSZEN / albmn[ISCENE][ISZEN]
+(CVZEN*CRAZ*rmatrx[ISCENE][ISZEN+1][IVZEN][IRAZ]
+ FVZEN*CRAZ*rmatrx[ISCENE][ISZEN+1][IVZEN+1][IRAZ]
+ CVZEN*FRAZ*rmatrx[ISCENE][ISZEN+1][IVZEN][IRAZ+1]
+ FVZEN*FRAZ*rmatrx[ISCENE][ISZEN+1][IVZEN+1][IRAZ+1] )
* FSZEN / albmn[ISCENE][ISZEN+1];

/* Convert the albedo. */

bbalb = broadalb / F;

if (bbalb < 0.) bbalb = 0.;

return bbalb;
}

/********** ANISOTROPIC CORRECTION WITH OCEAN/ICE BLENDING *****************/

float anis_correct_blend(float ch1, float ch1_pure_ocean, float ch1_pure_ice,
float broadalb, float solzen, float scan,
float relaz, char surftype)

/*--------------------------------------------------------------------------

Corrects the broadband TOA albedo data for anisotropic reflectance using a weighted sum blend of ocean and ice albedos for ocean and ice surface types. This procedure operates on scalar values only.

Input:

ch1: Channel 1 albedo at TOA , in the range [0,1], normalized for solar zenith angle, NOT corrected for the atmospheric attenuation, and corrected for anisotropic reflectance. [float]

ch1_pure_ocean: The ch1 value which will be considered pure ocean. [float]

ch1_pure_ice : The ch1 value which will be considered pure ocean. [float]


broadalb,
solzen,
scan,
relaz,
surftype : See anis_correct.

Algorithm Summary:
If surface type is LAND, or ch1 <= ch1_pure_ocean,
or ch1 >= ch1_pure_ice, anis_correct is called
with the given broadalb, solzen, scan, relaz, and surftype
parmameters and the return value from anis_correct is
returned without modification. Otherwise, anis_correct is
called twice with the given broadalb, solzen, scan, and relaz
parameters, but with the surftype parameter forced to OCEAN
on the first call and ICE on the second call. A weighted
sum average of the returned OCEAN and ICE albedos is then
calculated and returned.

The weighted sum is computed as follows:

range = ch1_pure_ice - ch1_pure_ocean
ocean_weight = (ch1_pure_ice - ch1) / range
ice_weight = 1.0 - ocean_weight
blended_albedo = ocean_albedo * ocean_weight +
ice_albedo * ice_weight
------------------------------------------------------------------------*/
{
float range;
float ocean_weight;
float ice_weight;
float ocean_albedo;
float ice_albedo;
float blended_albedo;

/* If it's dark or if the input value is missing, return MISSING. */

if (solzen >= DARKZEN || broadalb == MISSING) return MISSING;

if (surftype == (char)LAND || surftype == (char)SNOW ||
ch1 <= ch1_pure_ocean || ch1 >= ch1_pure_ice)
return(anis_correct(broadalb, solzen, scan, relaz, surftype));

range = ch1_pure_ice - ch1_pure_ocean;
ocean_weight = (ch1_pure_ice - ch1) / range;
ice_weight = 1.0 - ocean_weight;

ocean_albedo = anis_correct(broadalb, solzen, scan, relaz, (char)OCEAN);
ice_albedo = anis_correct(broadalb, solzen, scan, relaz, (char)ICE);
blended_albedo = ocean_albedo * ocean_weight + ice_albedo * ice_weight;

if (blended_albedo < 0.0)
blended_albedo = 0.0;

return blended_albedo;
}

/**************************** ATMOSPHERIC CORRECTION ************************/

float atmos_correct(float broadalb, float solzen, float wv, float odhaze,
char surftype)
/*------------------------------------------------------------------------
Converts a TOA broadband albedo to a surface broadband albedo. The
apparent (varies with atmospheric conditions) albedo will be calculated.
The apparent albedo is what would be measured with a pair of up/down
radiometers given arbitrary atmospheric conditions. The inherent albedo,
which is not calculated here (but could be), is an albedo without an
atmosphere; i.e., it doesn't depend upon gases and aerosols.

This procedure operates on scalar values only.

Input:
broadalb : TOA broadband albedo (corrected for anisotropy), in the range
[0,1].
solzen : Solar zenith angle, in degrees.
wv: Water vapor amount (0.5 - 5.0 g/cm^2
= precipitable water in cm)
odhaze: Aerosol optical depth (0.05 - 0.5)
surftype : One byte (char) value representing the surface type (global
variables LAND, OCEAN, ICE, SNOW see 'globals.h' for
values). LAND means snow-free land (vegetation cover),
OCEAN means open water, ICE means sea ice (with or without
snow), and SNOW means snow-covered land.
However, corrections are currently done for only snow-free
land, ocean, and snow/ice (the same) surfaces. Any other
surface type uses the snow procedure.

Output as value of function:
"Apparent" surface albedo corrected for the atmosphere in the
range [0,1], which are directional, hemispherical (i.e., are
not scan angle or azimuthally-dependent but do vary with the
atmospheric conditions, hence "apparent"). If badpix=TRUE or
the solar zenith angle > DARKZEN, MISSING values are returned.
If the pixel is CLOUDY then the values that these variables
have on input are unchanged. Results outside of the range
[0,1] are reset to 0 or 1. [fltarr or float]

Algorithm Summary:

The correction for atmospheric attenuation in in the broadband albedo
was developed with data generated by Streamer for two surface types
(snow and generic vegetation combined) and a variety of atmospheric
conditions (global). It is a function of solar zenith angle,
water vapor, and aerosol optical depth. Aerosol optical depth
and water vapor amounts are accounted for explicitly through linear
interpolation of coefficients (which probably introduces some error);
ozone amount is prescribed (0.24 cm). The form of the correction is:

surf_alb = (toa_alb - a) / b

where surf_alb is the atmospherically corrected broadband surface albedo,
toa_alb is the TOA broadband albedo, and a and b are determined through
regression using the Streamer data.

Unfrozen ocean uses a simpler approach, where the surface albedo is
related to the TOA albedo, zenith angle, aerosol optical depth,
and precipitable water by

surf_alb = a + b*toa_alb+ c*cos(zen) + d*pw + e*aertau

where a, b, c, d, and e are determined empirically using modeled
albedos. Ocean surface albedos less than 0.04 are set to that value.

The retrieved surface albedo is an "apparent" albedo; i.e., that which
would be measured with radiometers and which vary with changes in
atmospheric conditions, especially for bright surfaces.
------------------------------------------------------------------------*/
{
int nsols = 20;
float bbalb;
float sqrtpw, broadalb2, alpha1, alpha2, beta1, beta2, alpha, beta;

/* These coefficients are for 'nsols' solar zenith angle cosines from 1 down
to 0.05 in increments of 0.05. The first dimension corresponds to the
solar zenith angles in the order shown below. The second dimension holds
the 'a' and 'b' coefficients, in that order.

All calculations are for an ozone amount of 325 Dobson units.

The solar Zenith angles are as follows.
Cosines:
1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55,
0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05
Degrees:
0.00, 18.19, 25.84, 31.79, 36.87, 41.41, 45.57, 49.46, 53.13, 56.63,
60.00, 63.26, 66.42, 69.51, 72.54, 75.52, 78.46, 81.37, 84.26, 87.13 */

float odmin,odmax;
float wvmin,wvmax;
float odrange,oddiff,odscale;
float wvrange,wvdiff,wvscale;
int k;
float a1,b1,a2,b2,awv1,bwv1,awv2,bwv2;
float a,b;
float coszen;
int zen_index;

/*------ COEFFICIENTS FOR APPARENT ALBEDO OF SNOW/ICE AND VEGETATION -------*/

/* First for a water vapor amount of 0.5 g/cm^2 and aerosol optical
depth 0.05 (aer1_wv1), then wv=0.5, aer=0.5 (aer2_wv1), then
wv=5.0, aer=0.05 (aer1_wv2), then wv=5.0, aer=0.5 (aer2_wv2). */

float aer1_wv1[20][2] =
{ { 0.0369468, 0.7716573}, { 0.0393092, 0.7681513},
{ 0.0420182, 0.7638956}, { 0.0449436, 0.7596443},
{ 0.0480998, 0.7549872}, { 0.0515801, 0.7498252},
{ 0.0556449, 0.7435693}, { 0.0598915, 0.7376858},
{ 0.0648356, 0.7302867}, { 0.0702911, 0.7221886},
{ 0.0770464, 0.7122877}, { 0.0843299, 0.7015929},
{ 0.0933098, 0.6879621}, { 0.1038388, 0.6723217},
{ 0.1166860, 0.6532161}, { 0.1327160, 0.6291631},
{ 0.1532862, 0.5983471}, { 0.1807689, 0.5560154},
{ 0.2196831, 0.4951470}, { 0.2790609, 0.3983846} };

float aer2_wv1[20][2] =
{ { 0.0440627, 0.4562957}, { 0.0491487, 0.4497921},
{ 0.0544951, 0.4427959}, { 0.0607027, 0.4344482},
{ 0.0671587, 0.4260309}, { 0.0741401, 0.4170585},
{ 0.0819679, 0.4068565}, { 0.0903022, 0.3963189},
{ 0.0995805, 0.3844509}, { 0.1100535, 0.3704767},
{ 0.1216115, 0.3557419}, { 0.1344762, 0.3394325},
{ 0.1491109, 0.3203476}, { 0.1654393, 0.2998796},
{ 0.1841938, 0.2767967}, { 0.2054262, 0.2506457},
{ 0.2296357, 0.2223978}, { 0.2571991, 0.1917805},
{ 0.2878296, 0.1623356}, { 0.3209661, 0.1382332} };

float aer1_wv2[20][2] =
{ { 0.0301825, 0.7455254}, { 0.0328166, 0.7420131},
{ 0.0358166, 0.7377133}, { 0.0390927, 0.7331624},
{ 0.0426720, 0.7282085}, { 0.0464936, 0.7231386},
{ 0.0507808, 0.7172813}, { 0.0556268, 0.7107663},
{ 0.0611204, 0.7033693}, { 0.0675476, 0.6945344},
{ 0.0745468, 0.6851625}, { 0.0830090, 0.6737533},
{ 0.0926029, 0.6612120}, { 0.1043418, 0.6455948},
{ 0.1183378, 0.6274127}, { 0.1360318, 0.6041222},
{ 0.1584471, 0.5754020}, { 0.1884491, 0.5367222},
{ 0.2302195, 0.4839584}, { 0.2928393, 0.4020659} };

float aer2_wv2[20][2]=
{ { 0.0146276, 0.6222795}, { 0.0217345, 0.6142944},
{ 0.0296301, 0.6050117}, { 0.0379796, 0.5953020},
{ 0.0471900, 0.5848151}, { 0.0569098, 0.5736121},
{ 0.0678248, 0.5610381}, { 0.0795343, 0.5476506},
{ 0.0924229, 0.5327814}, { 0.1066507, 0.5160899},
{ 0.1222040, 0.4982739}, { 0.1395290, 0.4782666},
{ 0.1589199, 0.4553649}, { 0.1801602, 0.4307515},
{ 0.2038677, 0.4030677}, { 0.2301649, 0.3719538},
{ 0.2588158, 0.3378740}, { 0.2892141, 0.3007274},
{ 0.3194277, 0.2621159}, { 0.3452438, 0.2208999} };

/*-------------------------------------------------------------------------*/

/* If it's dark or if the input value is missing or if water vapor is
missing, return MISSING. */

if (solzen >= DARKZEN || broadalb == MISSING || wv == MISSING) return MISSING;

coszen = cos(solzen*DEGRAD);
if (coszen < 0.0501) coszen = .0501; /* Limit of solar angles */

/* SNOW, ICE, or LAND */

if (surftype != (char) OCEAN) {

zen_index = (int)(((1.0 - coszen)/0.05)+0.5);

/* The min and max aerosol optical depth and water vapor amounts
for the coefficients specified above. */

odmin=0.05;
odmax=0.5;
wvmin=0.5;
wvmax=5.0;

/* Determine the proper scale factor for the specified aerosol optical
depth and the water vapor amount. */

odrange = odmax - odmin;
oddiff = odhaze - odmin;
odscale = oddiff/odrange;
wvrange = wvmax - wvmin;
wvdiff = wv - wvmin;
wvscale = wvdiff/wvrange;

/* Get the estimated coefficients for the specified haze optical depth by
linearly interpolating between those for the optical depths in the
data statements above. This is done for each of the two water
vapor amounts. The final coefficients are then obtained by
interpolating (again linearly) between those computed for the desired
optical depth at the two water vapor extremes. */

a1 = aer1_wv1[zen_index][0];
b1 = aer1_wv1[zen_index][1];
a2 = aer2_wv1[zen_index][0];
b2 = aer2_wv1[zen_index][1];
awv1 = a1 + odscale*(a2 - a1);
bwv1 = b1 + odscale*(b2 - b1);

a1 = aer1_wv2[zen_index][0];
b1 = aer1_wv2[zen_index][1];
a2 = aer2_wv2[zen_index][0];
b2 = aer2_wv2[zen_index][1];
awv2 = a1 + odscale*(a2 - a1);
bwv2 = b1 + odscale*(b2 - b1);

a = awv1 + wvscale*(awv2 - awv1);
b = bwv1 + wvscale*(bwv2 - bwv1);

bbalb = (broadalb - a) / b;
}

/* OCEAN */

else {
bbalb = -0.112236 + broadalb*0.948389 + coszen*0.108496 + wv*0.00242575 -
odhaze*0.125026;
}

/*---------- ALTERNATE METHOD BY LI & GARAND, ALL SURFACE TYPES ------------*/
/*---------- (hasn't been tested in this function) Uncomment to use.

coszen = cos(solzen*DEGRAD);
sqrtpw = sqrt(wv)
broadalb2 = broadalb * 100.

alpha1 = -0.96882 + 0.718*sqrtpw;
alpha2 = -4.11460 - 0.76347*sqrtpw;
beta1 = 1.16711 + 0.05963*sqrtpw;
beta2 = 0.07514 + 0.04105*sqrtpw;
alpha = alpha1 + alpha2/coszen;
beta = beta1 + beta2/coszen;
bbalb = alpha + beta*broadalb2;
bbalb = bbalb/100.0;
-------------------------------------------------------------------------*/

if (bbalb < 0.04) bbalb = 0.04;

return bbalb;
}