atmocor2.c

Go to the documentation of this file.

#include "l12_proto.h"
#include "atrem_corl1.h"
 
/* --------------------------------------------------------------------------------------- */
/* atmocor2() - atmospheric correction, converts Lt -> nLw                                 */
/*                                                                                         */
/* C-version, September 2004, B. Franz                                                     */
 
/* --------------------------------------------------------------------------------------- */
 
int atmocor2(l2str *l2rec, aestr *aerec, int32_t ip) {
    static int firstCall = 1;
    static int want_nirRrs = 0;
    static int want_nirLw = 0;
    static int want_mumm = 0;
    static int want_ramp = 1;
    static int32_t aer_iter_max = 1;
    static int32_t aer_iter_min = 1;
 
    static float pi = PI;
    static float radeg = RADEG;
    static float badval = BAD_FLT;
    static float badchl = BAD_FLT;
    static float p0 = STDPR;
    static float df = 0.33;
    static float cbot = 0.7;
    static float ctop = 1.3;
    static float seed_chl = 0.0;
    static float seed_green = 0.0;
    static float seed_red = 0.0;
    static float nir_chg = 0.02;
    static float glint_min = GLINT_MIN;
    static float cslp;
    static float cint;
 
    static int32_t green;
    static int32_t red;
    static int32_t nir_s;
    static int32_t nir_l;
    static int32_t swir_s;
    static int32_t swir_l;
    static int32_t aer_s;
    static int32_t aer_l;
    static int32_t daer;
    static int32_t nwvis;
    static int32_t nwave;
    static float *wave;
 
    l1str *l1rec = l2rec->l1rec;
    filehandle *l1file = l1rec->l1file;
 
    int32_t sensorID = l1file->sensorID;
    int32_t brdf_opt = input->brdf_opt;
    int32_t aer_opt = input->aer_opt;
    int32_t glint_opt = input->glint_opt;
    int32_t cirrus_opt = input->cirrus_opt;
 
    float *Fo = l1rec->Fo;
    float *Fobar = l1file->Fobar;
    float fsol = l1rec->fsol;
    float solz = l1rec->solz [ip];
    float senz = l1rec->senz [ip];
    float delphi = l1rec->delphi[ip];
    float ws = l1rec->ws [ip];
    float gc = l1rec->glint_coef[ip];
 
    int32_t *aermodmin = &l2rec->aermodmin[ip];
    int32_t *aermodmax = &l2rec->aermodmax[ip];
    float *aermodrat = &l2rec->aerratio [ip];
    int32_t *aermodmin2 = &l2rec->aermodmin2[ip];
    int32_t *aermodmax2 = &l2rec->aermodmax2[ip];
    float *aermodrat2 = &l2rec->aerratio2 [ip];
    float *eps = &l2rec->eps [ip];
 
    float *TLg = &l1rec->TLg [ip * l1file->nbands];
    float *Lw = &l2rec->Lw [ip * l1file->nbands];
    float *nLw = &l2rec->nLw [ip * l1file->nbands];
    float *La = &l2rec->La [ip * l1file->nbands];
    float *taua = &l2rec->taua [ip * l1file->nbands];
    float *t_sol = &l1rec->t_sol[ip * l1file->nbands];
    float *t_sen = &l1rec->t_sen[ip * l1file->nbands];
    float *brdf = &l2rec->brdf [ip * l1file->nbands];
    float *Rrs = &l2rec->Rrs [ip * l1file->nbands];
 
    float *nLw_unc = &l2rec->nLw_unc[ip * l1file->nbands];
    float *Rrs_unc = &l2rec->Rrs_unc[ip * l1file->nbands];
 
    float *tg_sol = &l1rec->tg_sol[ip * l1file->nbands];
 
    float *taur;
    float *tLw;
    float *rhown_nir;
    float *tLw_nir, *last_tLw_nir;
    int32_t ib, ipb;
    int32_t status;
    int32_t iter, iter_max, iter_min, last_iter, iter_reset;
    float mu, mu0;
    int32_t nneg;
    float airmass;
    float *Ltemp;
    float *Lunc;
    float chl;
    int want_glintcorr;
    float refl_nir = 0, last_refl_nir = 0;
    float tindx;
    float Ka = 0.8; /* 1.375 um channel transmittance for water vapor (<=1)  Gao et al. 1998 JGR */
 
    if (firstCall == 1) {
        firstCall = 0;
        nwave = l1file->nbands;
        if ((wave = (float *) calloc(nwave, sizeof (float))) == NULL) {
            printf("-E- : Error allocating memory to wave\n");
            exit(FATAL_ERROR);
        }
        for (ib = 0; ib < nwave; ib++) {
            wave[ib] = l1file->fwave[ib];
        }
        if (sensorID == CZCS) {
            want_ramp = 0;
            seed_chl = 0.01;
            seed_green = 0.003;
            seed_red = 0.00125;
        }
        nwvis = rdsensorinfo(l1file->sensorID, l1_input->evalmask, "NbandsVIS", NULL);
        nir_s = bindex_get(input->aer_wave_short);
        nir_l = bindex_get(input->aer_wave_long);
        if (nir_s < 0 || nir_l < 0) {
            printf("Aerosol selection bands %d and %d not available for this sensor\n",
                    input->aer_wave_short, input->aer_wave_long);
            exit(1);
        }
        if (nir_l < nir_s) {
            printf("Invalid aerosol selection bands: long (%d) must be greater than short (%d).\n",
                    input->aer_wave_long, input->aer_wave_short);
            exit(1);
        }
        if (wave[nir_s] < 600) {
            printf("Aerosol selection band(s) must be greater than 600nm");
            exit(1);
        }
 
        aer_s = nir_s;
        aer_l = nir_l;
        daer = MAX(nir_l - nir_s, 1);
        cslp = 1. / (ctop - cbot);
        cint = -cslp * cbot;
        printf("Aerosol selection bands %d and %d\n", l1file->iwave[aer_s], l1file->iwave[aer_l]);
 
        switch (aer_opt) {
        case AERWANGNIR:
        case AERRHNIR:
        case FIXANGSTROMNIR:
        case FIXMODPAIRNIR:
            want_nirLw = 1;
            aer_iter_min = 1;
            aer_iter_max = input->aer_iter_max;
            printf("NIR correction enabled.\n");
            break;
        case AERMUMM:
        case AERRHMUMM:
            want_mumm = 1;
            aer_iter_min = 3;
            aer_iter_max = input->aer_iter_max;
            printf("MUMM correction enabled.\n");
            break;
        case AERWANGSWIR:
        case AERRHSWIR:
            want_nirLw = 1;
            aer_iter_min = 1;
            aer_iter_max = input->aer_iter_max;
            swir_s = bindex_get(input->aer_swir_short);
            swir_l = bindex_get(input->aer_swir_long);
            if (swir_s < 0 || swir_l < 0) {
                printf("Aerosol selection bands %d and %d not available for this sensor\n",
                        input->aer_swir_short, input->aer_swir_long);
                exit(1);
            }
            printf("NIR/SWIR switching correction enabled.\n");
            break;
        case AERRHMSEPS:
            want_nirLw = 1;
            aer_iter_min = 1;
            aer_iter_max = input->aer_iter_max;
            printf("NIR correction enabled --> for multi-scattering epsilon.\n");
            break;
        case AERRHSM:
            want_nirLw = 1; //This needs to be turned on, but a new SWIR water model is needed for it to work
            aer_iter_min = 1;
            aer_iter_max = input->aer_iter_max;
            printf("NIR correction enabled --> for spectral matching.\n");
            break;
        default:
            if (input->aer_rrs_short >= 0.0 && input->aer_rrs_long >= 0.0) {
                want_nirRrs = 1;
                aer_iter_min = 3;
                aer_iter_max = input->aer_iter_max;
                printf("NIR correction via input Rrs enabled.\n");
            }
            break;
        }
        if (input->aer_iter_max < 1)
            want_nirLw = 0;
        if ( want_nirLw || want_nirRrs) {
            if ((red = bindex_get(670)) < 0) {
                if ((red = bindex_get(680)) < 0)
                    if ((red = bindex_get(620)) < 0)
                        if ((red = bindex_get(765)) < 0)
                            if ((red = bindex_get(655)) < 0)
                                if ((red = bindex_get(664)) < 0) /* added for MSI */ {
                                    printf("%s line %d: can't find red band\n", __FILE__, __LINE__);
                                    exit(1);
                                }
            }
            if ((green = bindex_get(550)) < 0) {
                if ((green = bindex_get(555)) < 0)
                    if ((green = bindex_get(560)) < 0)
                        if ((green = bindex_get(565)) < 0) {
                            printf("%s line %d: can't find green band\n", __FILE__, __LINE__);
                            exit(1);
                        }
            }
        }
 
    }
 
    if ((taur = (float *) calloc(nwave, sizeof (float))) == NULL) {
        printf("-E- : Error allocating memory to taur\n");
        exit(FATAL_ERROR);
    }
    if ((tLw = (float *) calloc(nwave, sizeof (float))) == NULL) {
        printf("-E- : Error allocating memory to tLw\n");
        exit(FATAL_ERROR);
    }
    if ((rhown_nir = (float *) calloc(nwave, sizeof (float))) == NULL) {
        printf("-E- : Error allocating memory to rhown_nir\n");
        exit(FATAL_ERROR);
    }
    if ((tLw_nir = (float *) calloc(nwave, sizeof (float))) == NULL) {
        printf("-E- : Error allocating memory to tLw_nir\n");
        exit(FATAL_ERROR);
    }
    if ((last_tLw_nir = (float *) calloc(nwave, sizeof (float))) == NULL) {
        printf("-E- : Error allocating memory to last_tLw_nir\n");
        exit(FATAL_ERROR);
    }
    if ((Ltemp = (float *) calloc(nwave, sizeof (float))) == NULL) {
        printf("-E- : Error allocating memory to Ltemp\n");
        exit(FATAL_ERROR);
    }
    if ((Lunc = (float *) calloc(nwave, sizeof (float))) == NULL) {
        printf("-E- : Error allocating memory to Lunc\n");
        exit(FATAL_ERROR);
    }
 
    /* Initialize output values. If any input radiances are negative, just return. */
 
    status = 0;
    iter = 0;
    nneg = 0;
 
    for (ib = 0; ib < nwave; ib++) {
 
        ipb = ip * nwave + ib;
 
        //t_sol [ib]  = 1.0;    leave them as rayleigh only
        //t_sen [ib]  = 1.0;
        La [ib] = badval;
        tLw [ib] = badval;
        Lw [ib] = badval;
        nLw [ib] = badval;
        taua [ib] = badval;
        if (glint_opt != 2) TLg [ib] = 0.0;
        brdf [ib] = 1.0;
        Rrs [ib] = badval;
 
        l2rec->eps[ip] = badval;
 
        if (l2rec->l1rec->Lt[ipb] <= 0.0)
            if (wave[ib] < 1000) nneg++; /* don't test SWIR */
 }
 
    /* If any expected channels are negative */
    if (nneg > 0) {
        free(taur);
        free(tLw);
        free(rhown_nir);
        free(tLw_nir);
        free(last_tLw_nir);
        free(Ltemp);
        free(Lunc);
 
        status = 1;
        return (status);
    }
 
    mu0 = cos(solz / radeg);
    mu = cos(senz / radeg);
    airmass = 1.0 / mu0 + 1.0 / mu;
 
    /* Remove pre-computed atmospheric effects */
    /* --------------------------------------- */
    for (ib = 0; ib < nwave; ib++) {
 
        ipb = ip * nwave + ib;
 
        /* Pressure correct the Rayleigh optical thickness */
        taur[ib] = l1rec->pr[ip] / p0 * l1file->Tau_r[ib];
 
        /* Copy TOA radiance to temp var, eliminate missing bands */
        Ltemp[ib] = l2rec->l1rec->Lt[ipb];
        Lunc [ib] = l2rec->l1rec->Lt_unc[ipb];
 
        /* Correct for ozone absorption.  We correct for inbound and outbound here, */
        /* then we put the inbound back when computing Lw.                          */
        Ltemp[ib] = Ltemp[ib] / l1rec->tg_sol[ipb] / l1rec->tg_sen[ipb];
        Lunc [ib] = Lunc [ib] / l1rec->tg_sol[ipb] / l1rec->tg_sen[ipb];
 
        /* Do Cirrus correction - subtract off cirrus reflectance from Ltemp */
        /* Ka is the 1.375 um transmittance of water vapor above cirrus clouds */
        /* Add cirrus_opt to input */
        if (cirrus_opt) Ltemp[ib] -= l1rec->rho_cirrus[ip] / Ka * Fo[ib] * mu0 / pi;
 
        /*  Apply polarization correction */
        Ltemp[ib] /= l1rec->polcor[ipb];
        Lunc [ib] /= l1rec->polcor[ipb];
 
        /* Remove whitecap radiance */
        Ltemp[ib] -= l1rec->tLf[ipb];
 
        /* Subtract the Rayleigh contribution for this geometry. */
        Ltemp[ib] -= l1rec->Lr[ipb];
 
        /* If selected, correct for Ding and Gordon O2 absorption for the aerosol component (replace O2 losses) */
        if (input->oxaband_opt == 1) {
            Ltemp[ib] /= l1rec->t_o2[ipb];
            Lunc [ib] /= l1rec->t_o2[ipb];
        }
    }
 
 
    /* Compute BRDF correction, if not dependent on radiance */
    if (brdf_opt < FOQMOREL && brdf_opt > NOBRDF)
        ocbrdf(l2rec, ip, brdf_opt, wave, nwvis, solz, senz, delphi, ws, -1.0, NULL, NULL, brdf);
 
    /* Initialize iteration loop */
    chl = seed_chl;
    iter = 0;
    last_iter = 0;
    iter_max = aer_iter_max;
    iter_min = aer_iter_min;
    iter_reset = 0;
    last_refl_nir = 100.;
    want_glintcorr = 0;
 
    /*  new glint_opt usage - a 2 will use the simple TLg from atmocor1 */
    if (glint_opt == 1 && l1rec->glint_coef[ip] > glint_min) {
        iter_max = MAX(2, iter_max);
        iter_min = MAX(2, iter_min);
        want_glintcorr = 1;
    }
    if (glint_opt == 2) {
        want_glintcorr = 1;
    }
 
    if (aer_opt == AERWANGSWIR || aer_opt == AERRHSWIR) {
        tindx_shi(l2rec, ip, &tindx);
        if (tindx >= 1.05) {
            iter_max = 1;
            aer_s = swir_s;
            aer_l = swir_l;
            want_nirLw = 0;
        } else {
            aer_s = nir_s;
            aer_l = nir_l;
            want_nirLw = 1;
        }
        daer = MAX(aer_l - aer_s, 1);
    }
NIRSWIR:
 
    if (want_nirLw || want_nirRrs) {
        for (ib = 0; ib < nwave; ib++) {
            last_tLw_nir[ib] = 0.0;
            tLw_nir[ib] = 0.0;
            Rrs[ib] = 0.0;
        }
        Rrs[green] = seed_green;
        Rrs[red ] = seed_red;
    }
 
 
    /* -------------------------------------------------------------------- */
    /* Begin interations for aerosol with corrections for non-zero nLw(NIR) */
    /* -------------------------------------------------------------------- */
 
    while (!last_iter) {
 
        iter++;
        status = 0;
 
        /* Initialize tLw as surface + aerosol radiance */
        for (ib = 0; ib < nwave; ib++) {
            tLw[ib] = Ltemp[ib];
        }
 
        /*  Compute and subtract glint radiance */
        if (want_glintcorr) {
 
            if (glint_opt == 1)
                glint_rad(iter, nwave, aer_s, aer_l, gc, airmass, mu0, Fo, taur, taua, tLw, TLg);
 
 
            for (ib = 0; ib < nwave; ib++) {
                tLw[ib] -= TLg[ib];
            }
        }
 
        /* Adjust for non-zero NIR water-leaving radiances using input Rrs */
        if (want_nirRrs) {
 
            rhown_nir[aer_s] = pi * input->aer_rrs_short;
            rhown_nir[aer_l] = pi * input->aer_rrs_long;
 
            for (ib = aer_s; ib <= aer_l; ib += daer) {
 
                /* Convert NIR reflectance to TOA W-L radiance */
                tLw_nir[ib] = rhown_nir[ib] / pi * Fo[ib] * mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
 
                /* Avoid over-subtraction */
                if (tLw_nir[ib] > tLw[ib] && tLw[ib] > 0.0)
                    tLw_nir[ib] = tLw[ib];
 
                /* Remove estimated NIR water-leaving radiance */
                tLw[ib] -= tLw_nir[ib];
            }
        }
 
 
        /* Adjust for non-zero NIR water-leaving radiances using MUMM */
        if (want_mumm) {
 
            get_rhown_mumm(l2rec, ip, aer_s, aer_l, rhown_nir);
 
            for (ib = aer_s; ib <= aer_l; ib += daer) {
 
                /* Convert NIR reflectance to TOA W-L radiance */
                tLw_nir[ib] = rhown_nir[ib] / pi * Fo[ib] * mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
 
                /* Remove estimated NIR water-leaving radiance */
                tLw[ib] -= tLw_nir[ib];
            }
        }
 
 
        /* Adjust for non-zero NIR water-leaving radiances using IOP model */
        if (want_nirLw) {
 
            ipb = ip*nwave;
            get_rhown_eval(input->fqfile, Rrs, wave, aer_s, aer_l, nwave, &l1rec->sw_a_avg[ipb], &l1rec->sw_bb_avg[ipb], chl, solz, senz, delphi, rhown_nir);
 
            for (ib = aer_s; ib <= aer_l; ib += daer) {
 
                /* Convert NIR reflectance to TOA W-L radiance */
                tLw_nir[ib] = rhown_nir[ib] / pi * Fo[ib] * mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
 
                /* Iteration damping */
                tLw_nir[ib] = (1.0 - df) * tLw_nir[ib] + df * last_tLw_nir[ib];
 
                /* Ramp-up ?*/
                if (want_ramp) {
                    if (chl > 0.0 && chl <= cbot)
                        tLw_nir[ib] = 0.0;
                    else if ((chl > cbot) && (chl < ctop))
                        tLw_nir[ib] *= (cslp * chl + cint);
                }
 
                /* Remove estimated NIR water-leaving radiance */
                tLw[ib] -= tLw_nir[ib];
            }
        }
 
        /*  Compute the aerosol contribution */
        if (status == 0) {
            if (aer_opt != AERNULL)
                status = aerosol(l2rec, aer_opt, aerec, ip, wave, nwave, aer_s, aer_l, Fo, tLw,
                    La, t_sol, t_sen, eps, taua, aermodmin, aermodmax, aermodrat,
                    aermodmin2, aermodmax2, aermodrat2);
            else {
                for (ib = 0; ib < nwave; ib++) {
 
                    ipb = ip * nwave + ib;
 
                    t_sol [ib] = 1.0;
                    t_sen [ib] = 1.0;
                    La [ib] = 0.0;
                    taua [ib] = 0.0;
                    *eps = 1.0;
                }
            }
        }
 
 
        /* Subtract aerosols and normalize */
        if (status == 0) {
 
            for (ib = 0; ib < nwave; ib++) {
 
                /* subtract aerosol and normalize */
                tLw[ib] = tLw[ib] - La[ib];
                Lw [ib] = tLw[ib] / t_sen[ib] * tg_sol[ib];
                nLw[ib] = Lw [ib] / t_sol[ib] / tg_sol[ib] / mu0 / fsol * brdf[ib];
            }
 
            /* Compute new estimated chlorophyll */
            if (want_nirLw) {
                refl_nir = Rrs[red];
                for (ib = aer_s; ib <= aer_l; ib += daer)
                    last_tLw_nir[ib] = tLw_nir[ib];
                for (ib = 0; ib < nwvis; ib++) {
                    Rrs[ib] = nLw[ib] / Fobar[ib];
                }
                chl = get_default_chl(l2rec, Rrs);
 
                // if we passed atmospheric correction but the spectral distribution of
                // Rrs is bogus (chl failed), assume this is a turbid-water case and
                // reseed iteration as if all 670 reflectance is from water.
 
                if (chl == badchl && iter_reset == 0 && iter < iter_max) {
                    chl = 10.0;
                    Rrs[red] = 1.0 * (Ltemp[red] - TLg[red]) / t_sol[red] / tg_sol[red] / mu0 / fsol / Fobar[red];
                    iter_reset = 1;
                }
 
                // if we already tried a reset, and still no convergence, force one last
                // pass with an assumption that all red radiance is water component, and
                // force iteration to end.  this will be flagged as atmospheric correction
                // failure, but a qualitatively useful retrieval may still result.
 
                if (chl == badchl && iter_reset == 1 && iter < iter_max) {
                    chl = 10.0;
                    Rrs[red] = 1.0 * (Ltemp[red] - TLg[red]) / t_sol[red] / tg_sol[red] / mu0 / fsol / Fobar[red];
                    iter = iter_max; // so iter will trigger maxiter flag and ATMFAIL
                    iter_reset = 2;
                }
            }
 
        } else {
 
            /* Aerosol determination failure */
            for (ib = 0; ib < nwave; ib++) {
                La [ib] = badval;
                tLw[ib] = badval;
                Lw[ib] = badval;
                nLw[ib] = badval;
                Rrs[ib] = badval;
            }
        }
 
        /* If NIR/SWIR switching, do secondary test for turbidity and reset if needed */
        if (iter == 1 && (aer_opt == AERWANGSWIR || aer_opt == AERRHSWIR)) {
            if (tindx >= 1.05 && status == 0) {
                for (ib = 0; ib < nwvis; ib++) {
                    Rrs[ib] = nLw[ib] / Fobar[ib];
                }
                chl = get_default_chl(l2rec, Rrs);
                //printf("Checking turbidity %d %f %f %f\n",ip,tindx,chl,nLw[nir_l]);
                if (chl > 0 && (chl <= 1.0 || nLw[nir_l] < 0.08)) {
                    iter_max = aer_iter_max;
                    aer_s = nir_s;
                    aer_l = nir_l;
                    daer = MAX(aer_l - aer_s, 1);
                    want_nirLw = 1;
                    //printf("Reverting to NIR %d %f %f %f\n",ip,tindx,chl,nLw[nir_l]);
                    goto NIRSWIR;
                } else
                    l1rec->flags[ip] |= TURBIDW;
            }
        }
 
 
        /* Shall we continue iterating */
        if (status != 0) {
            last_iter = 1;
        } else if (iter < iter_min) {
            last_iter = 0;
        } else if (want_nirLw && (fabs(refl_nir - last_refl_nir) < fabs(nir_chg * refl_nir) || refl_nir < 0.0)) {
            last_iter = 1;
        } else if (want_mumm || want_nirRrs) {
            last_iter = 1;
        } else if (iter > iter_max) {
            last_iter = 1;
        }
 
        last_refl_nir = refl_nir;
 
    } /* end of iteration loop */
 
    l2rec->num_iter[ip] = iter;
 
 
    /* If the atmospheric correction failed, we don't need to do more. */
    if (status != 0) {
        free(taur);
        free(tLw);
        free(rhown_nir);
        free(tLw_nir);
        free(last_tLw_nir);
        free(Ltemp);
        free(Lunc);
        return (status);
 
    }
 
    /* If we used a NIR Lw correction, we record the tLw as it was predicted. */
    if (want_nirLw || want_nirRrs || want_mumm) {
        for (ib = aer_s; ib <= aer_l; ib += daer) {
            tLw[ib] = tLw_nir[ib];
            Lw [ib] = tLw[ib] / t_sen[ib] * tg_sol[ib];
            nLw[ib] = Lw [ib] / t_sol[ib] / tg_sol[ib] / mu0 / fsol * brdf[ib];
            Rrs[ib] = nLw[ib] / Fobar[ib];
        }
    }
 
    /* Convert water-leaving radiances from band-averaged to band-centered.  */
    /* Switch mean solar irradiance from band-averaged to band centered also.*/
    if (l1_input->outband_opt >= 2) {
        nlw_outband(l1_input->evalmask, sensorID, wave, nwave, Lw, nLw, &l2rec->outband_correction[ip*nwave]);
        Fobar = l1file->Fonom;
    }
 
    /* Compute f/Q correction and apply to nLw */
    if (brdf_opt >= FOQMOREL) {
        ocbrdf(l2rec, ip, brdf_opt, wave, nwvis, solz, senz, delphi, ws, -1.0, nLw, Fobar, brdf);
        for (ib = 0; ib < nwvis; ib++) {
            nLw[ib] = nLw[ib] * brdf[ib];
        }
    }
 
    /* Compute final Rrs */
    for (ib = 0; ib < nwave; ib++) {
        if (ib != aer_s && ib != aer_l) {
            Rrs[ib] = nLw[ib] / Fobar[ib];
            l2rec->Rrs[ip * nwave + ib] = Rrs[ib];
 
            nLw_unc[ib] = Lunc[ib] / t_sen[ib] / t_sol[ib] / mu0 / fsol * brdf[ib];
                    Rrs_unc[ib] = nLw_unc[ib] / Fobar[ib];
            }
        }
 
    /* Compute final chl from final nLw (needed for flagging) */
    l2rec->chl[ip] = get_default_chl(l2rec, Rrs);
 
    /*Determine Raman scattering contribution to Rrs*/
    run_raman_cor(l2rec, ip);
 
    free(taur);
    free(tLw);
    free(rhown_nir);
    free(tLw_nir);
    free(last_tLw_nir);
    free(Ltemp);
    free(Lunc);
 
    return (status);
}