/******************************************************************************
* @section DESCRIPTION
*
* This routine computes all surface fluxes, and solves the snow accumulation
* and ablation algorithm. Solutions are for the current snow band and
* vegetation type.
******************************************************************************/

#include <vic_run.h>

/******************************************************************************
* @brief        This routine computes all surface fluxes
******************************************************************************/
int
surface_fluxes(bool                 overstory,
               double               BareAlbedo,
               double               ice0,
               double               moist0,
               double               surf_atten,
               double              *Melt,
               double              *Le,
               double              *aero_resist,
               double              *displacement,
               double              *gauge_correction,
               double              *out_prec,
               double              *out_rain,
               double              *out_snow,
               double              *ref_height,
               double              *roughness,
               double              *snow_inflow,
               double              *wind,
               double              *root,
               size_t               Nlayers,
               size_t               Nveg,
               unsigned short       band,
               double               dp,
               unsigned short       iveg,
               unsigned short       veg_class,
               force_data_struct   *force,
               dmy_struct          *dmy,
               energy_bal_struct   *energy,
               global_param_struct *gp,
               cell_data_struct    *cell,
               snow_data_struct    *snow,
               soil_con_struct     *soil_con,
               veg_var_struct      *veg_var,
               double               lag_one,
               double               sigma_slope,
               double               fetch,
               double              *CanopLayerBnd)
{
    extern veg_lib_struct   *vic_run_veg_lib;
    extern option_struct     options;
    extern parameters_struct param;

    int                      ErrorFlag;
    int                      INCLUDE_SNOW = false;
    int                      UNSTABLE_CNT;
    int                      UNSTABLE_SNOW = false;
    int                      N_steps;
    int                      UnderStory;
    size_t                   hidx; // index of initial element of atmos array
    size_t                   step_inc; // number of atmos array elements to skip per surface fluxes step
    size_t                   endhidx; // index of final element of atmos array
    double                   step_dt; // time length of surface fluxes step (in seconds)
    size_t                   lidx;
    int                      over_iter;
    int                      under_iter;
    int                      q;
    double                   Ls;
    double                   LongUnderIn; // inmoing LW to ground surface
    double                   LongUnderOut; // outgoing LW from ground surface
    double                   NetLongSnow; // net LW over snowpack
    double                   NetShortSnow; // net SW over understory
    double                   NetShortGrnd; // net SW over snow-free surface
    double                   OldTSurf; // previous snow surface temperature
    double                   ShortUnderIn; // incoming SW to understory
    double                   Tair; // air temperature
    double                   Tcanopy; // canopy air temperature
    double                   Tgrnd; // soil surface temperature
    double                   Tsurf; // ground surface temperature
    double                   VPDcanopy; // vapor pressure deficit in canopy/atmos
    double                   VPcanopy; // vapor pressure in canopy/atmos
    double                   coverage; // mid-step snow cover fraction
    double                   delta_coverage; // change in snow cover fraction
    double                   delta_snow_heat; // change in snowpack heat storage
    double                   last_Tcanopy;
    double                   last_snow_coverage; // previous snow covered area
    double                   last_snow_flux;
    double                   last_tol_under; // previous surface iteration tol
    double                   last_tol_over; // previous overstory iteration tol
    double                   ppt; // precipitation/melt reaching soil surface
    double                   rainfall; // rainfall
    double                   snowfall; // snowfall
    double                   snow_flux; // heat flux through snowpack
    double                   snow_grnd_flux; // ground heat flux into snowpack
    double                   new_snow_albedo; // new snow albedo
    double                   tol_under;
    double                   tol_over;
    double                  *aero_resist_used;
    double                  *inflow;
    layer_data_struct       *layer;

    // Step-specific quantities
    double                   step_Wdew;
    double                   step_melt;
    double                   step_melt_energy; /* energy used to reduce snow coverage */
    double                   step_out_prec;
    double                   step_out_rain;
    double                   step_out_snow;
    double                   step_ppt;
    double                   step_prec;

    // Quantities that need to be summed or averaged over multiple snow steps
    // energy structure
    double            store_AlbedoOver;
    double            store_AlbedoUnder;
    double            store_AtmosLatent;
    double            store_AtmosLatentSub;
    double            store_AtmosSensible;
    double            store_LongOverIn;
    double            store_LongUnderIn;
    double            store_LongUnderOut;
    double            store_NetLongAtmos;
    double            store_NetLongOver;
    double            store_NetLongUnder;
    double            store_NetShortAtmos;
    double            store_NetShortGrnd;
    double            store_NetShortOver;
    double            store_NetShortUnder;
    double            store_ShortOverIn;
    double            store_ShortUnderIn;
    double            store_advected_sensible;
    double            store_advection;
    double            store_canopy_advection;
    double            store_canopy_latent;
    double            store_canopy_latent_sub;
    double            store_canopy_sensible;
    double            store_canopy_refreeze;
    double            store_deltaCC;
    double            store_deltaH;
    double            store_fusion;
    double            store_grnd_flux;
    double            store_latent;
    double            store_latent_sub;
    double            store_melt_energy;
    double            store_refreeze_energy;
    double            store_sensible;
    double            store_snow_flux;
    // snow structure
    double            store_canopy_vapor_flux;
    double            store_melt;
    double            store_vapor_flux;
    double            store_blowing_flux;
    double            store_surface_flux;
    // veg_var structure
    double            store_canopyevap;
    double            store_throughfall;
    // cell structure
    double            store_layerevap[MAX_LAYERS];
    double            store_ppt;
    double            store_aero_cond_used[2];
    double            store_pot_evap;

    // Structures holding values for current snow step
    energy_bal_struct snow_energy;    // energy fluxes at snowpack surface
    energy_bal_struct soil_energy;    // energy fluxes at soil surface
    veg_var_struct    snow_veg_var;    // veg fluxes/storages in presence of snow
    veg_var_struct    soil_veg_var;    // veg fluxes/storages in soil energy balance
    snow_data_struct  step_snow;
    layer_data_struct step_layer[MAX_LAYERS];

    // Structures holding values for current iteration
    energy_bal_struct iter_snow_energy;    // energy fluxes at snowpack surface
    energy_bal_struct iter_soil_energy;    // energy fluxes at soil surface
    veg_var_struct    iter_snow_veg_var;    // veg fluxes/storages in presence of snow
    veg_var_struct    iter_soil_veg_var;    // veg fluxes/storages in soil energy balance
    snow_data_struct  iter_snow;
    layer_data_struct iter_layer[MAX_LAYERS];
    double            iter_aero_resist[3];
    double            iter_aero_resist_veg[3];
    double            iter_aero_resist_used[2];
    double            iter_pot_evap;

    // handle bisection of understory solution
    double            store_tol_under;
    double            A_tol_under;

    // handle bisection of overstory solution
    double            store_tol_over;

    // Carbon cycling
    double            dryFrac;
    double           *LAIlayer = NULL;
    double           *faPAR = NULL;
    size_t            cidx;
    double            store_gc;
    double           *store_gsLayer = NULL;
    double            store_Ci;
    double            store_GPP;
    double            store_Rdark;
    double            store_Rphoto;
    double            store_Rmaint;
    double            store_Rgrowth;
    double            store_Raut;
    double            store_NPP;

    if (options.CARBON) {
        store_gsLayer = calloc(options.Ncanopy, sizeof(*store_gsLayer));
        check_alloc_status(store_gsLayer, "Memory allocation error.");
    }

    /***********************************************************************
       Set temporary variables for convenience
    ***********************************************************************/
    aero_resist_used = cell->aero_resist;
    inflow = &(cell->inflow);
    layer = cell->layer;

    /***********************************************************************
       Set temporary variables - preserves original values until iterations
       are completed
    ***********************************************************************/

    energy->advection = 0;
    energy->deltaCC = 0;
    if (snow->swq > 0) {
        snow_flux = energy->snow_flux;
    }
    else {
        snow_flux = -(energy->grnd_flux + energy->deltaH + energy->fusion);
    }
    energy->refreeze_energy = 0;
    coverage = snow->coverage;
    snow_energy = (*energy);
    soil_energy = (*energy);
    iter_soil_energy = (*energy);
    snow_veg_var = (*veg_var);
    soil_veg_var = (*veg_var);
    step_snow = (*snow);
    for (lidx = 0; lidx < Nlayers; lidx++) {
        step_layer[lidx] = layer[lidx];
    }
    for (lidx = 0; lidx < Nlayers; lidx++) {
        step_layer[lidx].evap = 0;
    }
    soil_veg_var.canopyevap = 0;
    snow_veg_var.canopyevap = 0;
    soil_veg_var.throughfall = 0;
    snow_veg_var.throughfall = 0;

    /********************************
       Set-up sub-time step controls
       (May eventually want to set this up so that it is also true
       if frozen soils are present)
    ********************************/

    if (snow->swq > 0 || snow->snow_canopy > 0 || force->snowflag[NR]) {
        hidx = 0;
        step_inc = 1;
        endhidx = hidx + NF;
        step_dt = gp->snow_dt;
    }
    else {
        hidx = NR;
        step_inc = 1;
        endhidx = hidx + step_inc;
        step_dt = gp->dt;
    }

    /*******************************************
       Initialize sub-model time step variables
    *******************************************/

    // energy structure
    store_AlbedoOver = 0;
    store_AlbedoUnder = 0;
    store_AtmosLatent = 0;
    store_AtmosLatentSub = 0;
    store_AtmosSensible = 0;
    store_LongOverIn = 0;
    store_LongUnderIn = 0;
    store_LongUnderOut = 0;
    store_NetLongAtmos = 0;
    store_NetLongOver = 0;
    store_NetLongUnder = 0;
    store_NetShortAtmos = 0;
    store_NetShortGrnd = 0;
    store_NetShortOver = 0;
    store_NetShortUnder = 0;
    store_ShortOverIn = 0;
    store_ShortUnderIn = 0;
    store_advected_sensible = 0;
    store_advection = 0;
    store_canopy_advection = 0;
    store_canopy_latent = 0;
    store_canopy_latent_sub = 0;
    store_canopy_sensible = 0;
    store_canopy_refreeze = 0;
    store_deltaCC = 0;
    store_deltaH = 0;
    store_fusion = 0;
    store_grnd_flux = 0;
    store_latent = 0;
    store_latent_sub = 0;
    store_melt_energy = 0;
    store_refreeze_energy = 0;
    store_sensible = 0;
    store_snow_flux = 0;
    // snow structure
    last_snow_coverage = snow->coverage;
    store_canopy_vapor_flux = 0;
    store_melt = 0;
    store_vapor_flux = 0;
    store_surface_flux = 0;
    store_blowing_flux = 0;
    // veg_var and cell structures
    store_throughfall = 0.;
    store_canopyevap = 0.;
    for (lidx = 0; lidx < options.Nlayer; lidx++) {
        store_layerevap[lidx] = 0.;
    }
    step_Wdew = veg_var->Wdew;
    // misc
    store_ppt = 0;
    store_aero_cond_used[0] = 0;
    store_aero_cond_used[1] = 0;
    (*snow_inflow) = 0;
    store_pot_evap = 0;
    N_steps = 0;

    // Carbon cycling
    if (options.CARBON) {
        store_gc = 0;
        for (cidx = 0; cidx < options.Ncanopy; cidx++) {
            store_gsLayer[cidx] = 0;
        }
        store_Ci = 0;
        store_GPP = 0;
        store_Rdark = 0;
        store_Rphoto = 0;
        store_Rmaint = 0;
        store_Rgrowth = 0;
        store_Raut = 0;
        store_NPP = 0;
    }

    /*************************
       Compute surface fluxes
    *************************/

    do
    {
        /** Solve energy balance for all sub-model time steps **/


        /* set air temperature and precipitation for this snow band */
        Tair = force->air_temp[hidx] + soil_con->Tfactor[band];
        step_prec = force->prec[hidx] * soil_con->Pfactor[band];

        // initialize ground surface temperaure
        Tgrnd = energy->T[0];

        // initialize canopy terms
        Tcanopy = Tair;
        VPcanopy = force->vp[hidx];
        VPDcanopy = force->vpd[hidx];

        over_iter = 0;
        tol_over = 999;

        last_Tcanopy = 999;
        last_snow_flux = 999;

        // compute LAI and absorbed PAR per canopy layer
        if (options.CARBON && iveg < Nveg) {
            LAIlayer = calloc(options.Ncanopy, sizeof(*LAIlayer));
            check_alloc_status(LAIlayer, "Memory allocation error.");
            faPAR = calloc(options.Ncanopy, sizeof(*faPAR));
            check_alloc_status(faPAR, "Memory allocation error.");

            /* Compute absorbed PAR per ground area per canopy layer (W/m2)
               normalized to PAR = 1 W, i.e. the canopy albedo in the PAR
               range (alb_total ~ 0.45*alb_par + 0.55*alb_other) */
            faparl(CanopLayerBnd,
                   veg_var->LAI,
                   soil_con->AlbedoPar,
                   force->coszen[hidx],
                   force->fdir[hidx],
                   LAIlayer,
                   faPAR);

            /* Convert to absolute (unnormalized) absorbed PAR per leaf area per canopy layer
               (umol(photons)/m2 leaf area / s); dividing by Epar converts PAR from W to umol(photons)/s */
            veg_var->aPAR = 0;
            for (cidx = 0; cidx < options.Ncanopy; cidx++) {
                if (LAIlayer[cidx] > 1e-10) {
                    veg_var->aPARLayer[cidx] =
                        (force->par[hidx] /
                         param.PHOTO_EPAR) * faPAR[cidx] / LAIlayer[cidx];
                    veg_var->aPAR += force->par[hidx] * faPAR[cidx] /
                                     LAIlayer[cidx];
                }
                else {
                    veg_var->aPARLayer[cidx] = force->par[hidx] /
                                               param.PHOTO_EPAR *
                                               faPAR[cidx] / 1e-10;
                    veg_var->aPAR += force->par[hidx] * faPAR[cidx] / 1e-10;
                }
            }
            free((char*) LAIlayer);
            free((char*) faPAR);
        }

        // Compute mass flux of blowing snow
        if (!overstory && options.BLOWING && step_snow.swq > 0.) {
            Ls = calc_latent_heat_of_sublimation(step_snow.surf_temp);
            step_snow.blowing_flux = CalcBlowingSnow(step_dt, Tair,
                                                     step_snow.last_snow,
                                                     step_snow.surf_water,
                                                     wind[2], Ls,
                                                     force->density[hidx],
                                                     force->vp[hidx],
                                                     roughness[2],
                                                     ref_height[2],
                                                     step_snow.depth,
                                                     lag_one, sigma_slope,
                                                     step_snow.surf_temp, iveg,
                                                     Nveg, fetch,
                                                     displacement[1],
                                                     roughness[1],
                                                     &step_snow.transport);
            if ((int) step_snow.blowing_flux == ERROR) {
                return (ERROR);
            }
            step_snow.blowing_flux *= step_dt / CONST_RHOFW; /* m/time step */
        }
        else {
            step_snow.blowing_flux = 0.0;
        }

        do
        {
            /** Iterate for overstory solution **/

            over_iter++;
            last_tol_over = tol_over;

            under_iter = 0;
            tol_under = 999;
            UnderStory = 999;

            UNSTABLE_CNT = 0;

            // bisect understory
            A_tol_under = 999;
            store_tol_under = 999;

            store_tol_over = 999;

            do
            {
                /** Iterate for understory solution - itererates to find snow flux **/

                under_iter++;
                last_tol_under = tol_under;

                if (last_Tcanopy != 999) {
                    Tcanopy = (last_Tcanopy + Tcanopy) / 2.;
                }
                last_Tcanopy = Tcanopy;

                // update understory energy balance terms for iteration
                if (last_snow_flux != 999) {
                    if ((fabs(store_tol_under) > fabs(A_tol_under) &&
                         A_tol_under != 999 &&
                         fabs(store_tol_under - A_tol_under) > 1.) ||
                        tol_under < 0) { // stepped the correct way
                        UNSTABLE_CNT++;
                        if (UNSTABLE_CNT > 3 || tol_under < 0) {
                            UNSTABLE_SNOW = true;
                        }
                    }
                    else if (!INCLUDE_SNOW) { // stepped the wrong way
                        snow_flux =
                            (last_snow_flux + iter_soil_energy.snow_flux) / 2.;
                    }
                }
                last_snow_flux = snow_flux;
                A_tol_under = store_tol_under;

                snow_grnd_flux = -snow_flux;

                // Initialize structures for new iteration
                iter_snow_energy = snow_energy;
                iter_soil_energy = soil_energy;
                iter_snow_veg_var = snow_veg_var;
                iter_soil_veg_var = soil_veg_var;
                iter_snow = step_snow;
                for (lidx = 0; lidx < Nlayers; lidx++) {
                    iter_layer[lidx] = step_layer[lidx];
                }
                iter_snow_veg_var.Wdew = step_Wdew;
                iter_soil_veg_var.Wdew = step_Wdew;
                iter_snow_veg_var.canopyevap = 0;
                iter_soil_veg_var.canopyevap = 0;
                for (lidx = 0; lidx < Nlayers; lidx++) {
                    iter_layer[lidx].evap = 0;
                }
                for (q = 0; q < 3; q++) {
                    iter_aero_resist[q] = aero_resist[q];
                    iter_aero_resist_veg[q] = aero_resist[q];
                }
                iter_aero_resist_used[0] = aero_resist_used[0];
                iter_aero_resist_used[1] = aero_resist_used[1];
                iter_snow.canopy_vapor_flux = 0;
                iter_snow.vapor_flux = 0;
                iter_snow.surface_flux = 0;
                /* iter_snow.blowing_flux has already been reset to step_snow.blowing_flux */
                LongUnderOut = iter_soil_energy.LongUnderOut;
                dryFrac = -1;

                // set new snow albedo
                if (options.MAX_SNOW_ALBEDO) {
                    // use maximum snow albedo from parameter file
                    if (iveg != Nveg) {
                        new_snow_albedo =
                            vic_run_veg_lib[veg_class].max_snow_albedo;
                    }
                    else {
                        new_snow_albedo = param.SNOW_NEW_SNOW_ALB;
                    }
                }
                else {
                    new_snow_albedo = param.SNOW_NEW_SNOW_ALB;
                }


                /** Solve snow accumulation, ablation and interception **/
                step_melt = solve_snow(overstory, BareAlbedo, LongUnderOut,
                                       param.SNOW_MIN_RAIN_TEMP,
                                       param.SNOW_MAX_SNOW_TEMP,
                                       new_snow_albedo,
                                       Tcanopy, Tgrnd, Tair,
                                       step_prec, snow_grnd_flux,
                                       &energy->AlbedoUnder, Le,
                                       &LongUnderIn, &NetLongSnow,
                                       &NetShortGrnd,
                                       &NetShortSnow, &ShortUnderIn, &OldTSurf,
                                       iter_aero_resist, iter_aero_resist_used,
                                       &coverage, &delta_coverage,
                                       &delta_snow_heat, displacement,
                                       gauge_correction, &step_melt_energy,
                                       &step_out_prec, &step_out_rain,
                                       &step_out_snow,
                                       &step_ppt, &rainfall, ref_height,
                                       roughness, snow_inflow, &snowfall,
                                       &surf_atten,
                                       wind, root, UNSTABLE_SNOW,
                                       Nveg, iveg, band, step_dt, hidx,
                                       veg_class,
                                       &UnderStory, CanopLayerBnd, &dryFrac,
                                       dmy, force, &(iter_snow_energy),
                                       iter_layer, &(iter_snow),
                                       soil_con,
                                       &(iter_snow_veg_var));

                if (step_melt == ERROR) {
                    return (ERROR);
                }

                /* Check that the snow surface temperature was estimated, if not
                   prepare to include thin snowpack in the estimation of the
                   snow-free surface energy balance */
                if ((iter_snow.surf_temp == 999 || UNSTABLE_SNOW) &&
                    iter_snow.swq > 0) {
                    INCLUDE_SNOW = UnderStory + 1;
                    iter_soil_energy.advection = iter_snow_energy.advection;
                    iter_snow.surf_temp = step_snow.surf_temp;
                    step_melt_energy = 0;
                }
                else {
                    INCLUDE_SNOW = false;
                }

                if (iter_snow.snow) {
                    iter_aero_resist_veg[0] = iter_aero_resist_used[0];
                    iter_aero_resist_veg[1] = iter_aero_resist_used[1];
                }

                /**************************************************
                   Solve Energy Balance Components at Soil Surface
                **************************************************/

                Tsurf = calc_surf_energy_bal((*Le), LongUnderIn, NetLongSnow,
                                             NetShortGrnd, NetShortSnow,
                                             OldTSurf,
                                             ShortUnderIn, iter_snow.albedo,
                                             iter_snow_energy.latent,
                                             iter_snow_energy.latent_sub,
                                             iter_snow_energy.sensible,
                                             Tcanopy, VPDcanopy,
                                             VPcanopy,
                                             delta_coverage, dp,
                                             ice0, step_melt_energy, moist0,
                                             iter_snow.coverage,
                                             (step_snow.depth + iter_snow.depth) / 2.,
                                             BareAlbedo, surf_atten,
                                             iter_aero_resist, iter_aero_resist_veg, iter_aero_resist_used,
                                             displacement, &step_melt, &step_ppt,
                                             rainfall, ref_height, roughness,
                                             snowfall, wind, root, INCLUDE_SNOW,
                                             UnderStory, options.Nnode, Nveg,
                                             step_dt, hidx, iveg,
                                             (int) overstory, veg_class,
                                             CanopLayerBnd, &dryFrac, force,
                                             dmy, &iter_soil_energy,
                                             iter_layer,
                                             &(iter_snow), soil_con,
                                             &iter_soil_veg_var);

                if ((int) Tsurf == ERROR) {
                    // Return error flag to skip rest of grid cell
                    return (ERROR);
                }

                if (INCLUDE_SNOW) {
                    /* store melt from thin snowpack */
                    step_ppt += step_melt;
                }

                /*****************************************
                   Compute energy balance with atmosphere
                *****************************************/
                if (iter_snow.snow && overstory) {
                    // do this if overstory is active and energy balance is closed
                    Tcanopy = calc_atmos_energy_bal(
                        iter_snow_energy.canopy_sensible,
                        iter_soil_energy.sensible,
                        iter_snow_energy.canopy_latent,
                        iter_soil_energy.latent,
                        iter_snow_energy.canopy_latent_sub,
                        iter_soil_energy.latent_sub,
                        iter_snow_energy.NetLongOver,
                        iter_soil_energy.NetLongUnder,
                        iter_snow_energy.NetShortOver,
                        iter_soil_energy.NetShortUnder,
                        iter_aero_resist_veg[1], Tair,
                        force->density[hidx],
                        &iter_soil_energy.AtmosError,
                        &iter_soil_energy.AtmosLatent,
                        &iter_soil_energy.AtmosLatentSub,
                        &iter_soil_energy.NetLongAtmos,
                        &iter_soil_energy.NetShortAtmos,
                        &iter_soil_energy.AtmosSensible,
                        &iter_soil_energy.Tcanopy_fbflag,
                        &iter_soil_energy.Tcanopy_fbcount);

                    /* iterate to find Tcanopy which will solve the atmospheric energy
                       balance.  Since I do not know vp in the canopy, use the
                       sum of latent heats from the ground and foliage, and iterate
                       on the temperature used for the sensible heat flux from the
                       canopy air to the mixing level */
                    if ((int) Tcanopy == ERROR) {
                        // Return error flag to skip rest of grid cell
                        return (ERROR);
                    }
                }
                else {
                    // else put surface fluxes into atmospheric flux storage so that
                    // the model will continue to function
                    iter_soil_energy.AtmosLatent = iter_soil_energy.latent;
                    iter_soil_energy.AtmosLatentSub =
                        iter_soil_energy.latent_sub;
                    iter_soil_energy.AtmosSensible = iter_soil_energy.sensible;
                    iter_soil_energy.NetLongAtmos =
                        iter_soil_energy.NetLongUnder;
                    iter_soil_energy.NetShortAtmos =
                        iter_soil_energy.NetShortUnder;
                }
                iter_soil_energy.Tcanopy = Tcanopy;
                iter_snow_energy.Tcanopy = Tcanopy;

                /*****************************************
                   Compute iteration tolerance statistics
                *****************************************/

                // compute understory tolerance
                if (INCLUDE_SNOW ||
                    (iter_snow.swq == 0 && delta_coverage == 0)) {
                    store_tol_under = 0;
                    tol_under = 0;
                }
                else {
                    store_tol_under = snow_flux - iter_soil_energy.snow_flux;
                    tol_under = fabs(store_tol_under);
                }
                if (fabs(tol_under - last_tol_under) < param.TOL_GRND &&
                    tol_under >
                    1.) {
                    tol_under = -999;
                }

                // compute overstory tolerance
                if (overstory && iter_snow.snow) {
                    store_tol_over = Tcanopy - last_Tcanopy;
                    tol_over = fabs(store_tol_over);
                }
                else {
                    store_tol_over = 0;
                    tol_over = 0;
                }
            }
            while ((fabs(tol_under - last_tol_under) > param.TOL_GRND) &&
                   (tol_under != 0) &&
                   (under_iter < param.MAX_ITER_GRND_CANOPY));
        }
        while ((fabs(tol_over - last_tol_over) > param.TOL_OVER &&
                overstory) && (tol_over != 0) &&
               (over_iter < param.MAX_ITER_GRND_CANOPY));

        /**************************************
           Compute GPP, Raut, and NPP
        **************************************/
        if (options.CARBON) {
            if (iveg < Nveg && !step_snow.snow && dryFrac > 0) {
                canopy_assimilation(vic_run_veg_lib[veg_class].Ctype,
                                    vic_run_veg_lib[veg_class].MaxCarboxRate,
                                    vic_run_veg_lib[veg_class].MaxETransport,
                                    vic_run_veg_lib[veg_class].CO2Specificity,
                                    iter_soil_veg_var.NscaleFactor,
                                    Tair,
                                    force->shortwave[hidx],
                                    iter_soil_veg_var.aPARLayer,
                                    soil_con->elevation,
                                    force->Catm[hidx],
                                    CanopLayerBnd,
                                    veg_var->LAI,
                                    "rs",
                                    iter_soil_veg_var.rsLayer,
                                    &(iter_soil_veg_var.rc),
                                    &(iter_soil_veg_var.Ci),
                                    &(iter_soil_veg_var.GPP),
                                    &(iter_soil_veg_var.Rdark),
                                    &(iter_soil_veg_var.Rphoto),
                                    &(iter_soil_veg_var.Rmaint),
                                    &(iter_soil_veg_var.Rgrowth),
                                    &(iter_soil_veg_var.Raut),
                                    &(iter_soil_veg_var.NPP));
                /* Adjust by fraction of canopy that was dry and account for any other inhibition`*/
                dryFrac *= iter_soil_veg_var.NPPfactor;
                iter_soil_veg_var.GPP *= dryFrac;
                iter_soil_veg_var.Rdark *= dryFrac;
                iter_soil_veg_var.Rphoto *= dryFrac;
                iter_soil_veg_var.Rmaint *= dryFrac;
                iter_soil_veg_var.Rgrowth *= dryFrac;
                iter_soil_veg_var.Raut *= dryFrac;
                iter_soil_veg_var.NPP *= dryFrac;
                /* Adjust by veg cover fraction */
                iter_soil_veg_var.GPP *= iter_soil_veg_var.fcanopy;
                iter_soil_veg_var.Rdark *= iter_soil_veg_var.fcanopy;
                iter_soil_veg_var.Rphoto *= iter_soil_veg_var.fcanopy;
                iter_soil_veg_var.Rmaint *= iter_soil_veg_var.fcanopy;
                iter_soil_veg_var.Rgrowth *= iter_soil_veg_var.fcanopy;
                iter_soil_veg_var.Raut *= iter_soil_veg_var.fcanopy;
                iter_soil_veg_var.NPP *= iter_soil_veg_var.fcanopy;
            }
            else {
                iter_soil_veg_var.rc = param.HUGE_RESIST;
                for (cidx = 0; cidx < options.Ncanopy; cidx++) {
                    iter_soil_veg_var.rsLayer[cidx] = param.HUGE_RESIST;
                }
                iter_soil_veg_var.Ci = 0;
                iter_soil_veg_var.GPP = 0;
                iter_soil_veg_var.Rdark = 0;
                iter_soil_veg_var.Rphoto = 0;
                iter_soil_veg_var.Rmaint = 0;
                iter_soil_veg_var.Rgrowth = 0;
                iter_soil_veg_var.Raut = 0;
                iter_soil_veg_var.NPP = 0;
            }
        }

        /**************************************
           Compute Potential Evap
        **************************************/

        compute_pot_evap(gp->model_steps_per_day,
                         vic_run_veg_lib[veg_class].rmin,
                         iter_soil_veg_var.albedo, force->shortwave[hidx],
                         iter_soil_energy.NetLongAtmos,
                         vic_run_veg_lib[veg_class].RGL, Tair, VPDcanopy,
                         iter_soil_veg_var.LAI, soil_con->elevation,
                         iter_aero_resist_veg,
                         vic_run_veg_lib[veg_class].overstory,
                         vic_run_veg_lib[veg_class].rarc,
                         iter_soil_veg_var.fcanopy, iter_aero_resist_used[0],
                         &iter_pot_evap);

        /**************************************
           Store sub-model time step variables
        **************************************/

        snow_energy = iter_snow_energy;
        soil_energy = iter_soil_energy;
        snow_veg_var = iter_snow_veg_var;
        soil_veg_var = iter_soil_veg_var;
        step_snow = iter_snow;
        for (lidx = 0; lidx < options.Nlayer; lidx++) {
            step_layer[lidx] = iter_layer[lidx];
        }

        if (iveg != Nveg) {
            if (step_snow.snow) {
                store_throughfall += snow_veg_var.throughfall;
                store_canopyevap += snow_veg_var.canopyevap;
                soil_veg_var.Wdew = snow_veg_var.Wdew;
            }
            else {
                store_throughfall += soil_veg_var.throughfall;
                store_canopyevap += soil_veg_var.canopyevap;
                snow_veg_var.Wdew = soil_veg_var.Wdew;
            }
            step_Wdew = soil_veg_var.Wdew;
            if (options.CARBON) {
                store_gc += 1 / soil_veg_var.rc;
                for (cidx = 0; cidx < options.Ncanopy; cidx++) {
                    store_gsLayer[cidx] += 1 / soil_veg_var.rsLayer[cidx];
                }
                store_Ci += soil_veg_var.Ci;
                store_GPP += soil_veg_var.GPP;
                store_Rdark += soil_veg_var.Rdark;
                store_Rphoto += soil_veg_var.Rphoto;
                store_Rmaint += soil_veg_var.Rmaint;
                store_Rgrowth += soil_veg_var.Rgrowth;
                store_Raut += soil_veg_var.Raut;
                store_NPP += soil_veg_var.NPP;
            }
        }
        for (lidx = 0; lidx < options.Nlayer; lidx++) {
            store_layerevap[lidx] += step_layer[lidx].evap;
        }
        store_ppt += step_ppt;
        if (iter_aero_resist_used[0] > 0) {
            store_aero_cond_used[0] += 1 / iter_aero_resist_used[0];
        }
        else {
            store_aero_cond_used[0] += param.HUGE_RESIST;
        }
        if (iter_aero_resist_used[1] > 0) {
            store_aero_cond_used[1] += 1 / iter_aero_resist_used[1];
        }
        else {
            store_aero_cond_used[1] += param.HUGE_RESIST;
        }

        if (iveg != Nveg) {
            store_canopy_vapor_flux += step_snow.canopy_vapor_flux;
        }
        store_melt += step_melt;
        store_vapor_flux += step_snow.vapor_flux;
        store_surface_flux += step_snow.surface_flux;
        store_blowing_flux += step_snow.blowing_flux;

        out_prec[0] += step_out_prec;
        out_rain[0] += step_out_rain;
        out_snow[0] += step_out_snow;

        if (INCLUDE_SNOW) {
            /* copy needed flux terms to the snowpack */
            snow_energy.advected_sensible = soil_energy.advected_sensible;
            snow_energy.advection = soil_energy.advection;
            snow_energy.deltaCC = soil_energy.deltaCC;
            snow_energy.latent = soil_energy.latent;
            snow_energy.latent_sub = soil_energy.latent_sub;
            snow_energy.refreeze_energy = soil_energy.refreeze_energy;
            snow_energy.sensible = soil_energy.sensible;
            snow_energy.snow_flux = soil_energy.snow_flux;
        }

        store_AlbedoOver += snow_energy.AlbedoOver;
        store_AlbedoUnder += soil_energy.AlbedoUnder;
        store_AtmosLatent += soil_energy.AtmosLatent;
        store_AtmosLatentSub += soil_energy.AtmosLatentSub;
        store_AtmosSensible += soil_energy.AtmosSensible;
        store_LongOverIn += snow_energy.LongOverIn;
        store_LongUnderIn += LongUnderIn;
        store_LongUnderOut += soil_energy.LongUnderOut;
        store_NetLongAtmos += soil_energy.NetLongAtmos;
        store_NetLongOver += snow_energy.NetLongOver;
        store_NetLongUnder += soil_energy.NetLongUnder;
        store_NetShortAtmos += soil_energy.NetShortAtmos;
        store_NetShortGrnd += NetShortGrnd;
        store_NetShortOver += snow_energy.NetShortOver;
        store_NetShortUnder += soil_energy.NetShortUnder;
        store_ShortOverIn += snow_energy.ShortOverIn;
        store_ShortUnderIn += soil_energy.ShortUnderIn;
        store_canopy_advection += snow_energy.canopy_advection;
        store_canopy_latent += snow_energy.canopy_latent;
        store_canopy_latent_sub += snow_energy.canopy_latent_sub;
        store_canopy_sensible += snow_energy.canopy_sensible;
        store_canopy_refreeze += snow_energy.canopy_refreeze;
        store_deltaH += soil_energy.deltaH;
        store_fusion += soil_energy.fusion;
        store_grnd_flux += soil_energy.grnd_flux;
        store_latent += soil_energy.latent;
        store_latent_sub += soil_energy.latent_sub;
        store_melt_energy += step_melt_energy;
        store_sensible += soil_energy.sensible;
        if (step_snow.swq == 0 && INCLUDE_SNOW) {
            if (last_snow_coverage == 0 && step_prec > 0) {
                last_snow_coverage = 1;
            }
            store_advected_sensible += snow_energy.advected_sensible *
                                       last_snow_coverage;
            store_advection += snow_energy.advection * last_snow_coverage;
            store_deltaCC += snow_energy.deltaCC * last_snow_coverage;
            store_snow_flux += soil_energy.snow_flux * last_snow_coverage;
            store_refreeze_energy += snow_energy.refreeze_energy *
                                     last_snow_coverage;
        }
        else if (step_snow.snow || INCLUDE_SNOW) {
            store_advected_sensible += snow_energy.advected_sensible *
                                       (step_snow.coverage + delta_coverage);
            store_advection += snow_energy.advection *
                               (step_snow.coverage + delta_coverage);
            store_deltaCC += snow_energy.deltaCC *
                             (step_snow.coverage + delta_coverage);
            store_snow_flux += soil_energy.snow_flux *
                               (step_snow.coverage + delta_coverage);
            store_refreeze_energy += snow_energy.refreeze_energy *
                                     (step_snow.coverage + delta_coverage);
        }
        store_pot_evap += iter_pot_evap;

        /* increment time step */
        N_steps++;
        hidx += step_inc;
    }
    while (hidx < endhidx);

    /************************************************
       Store snow variables for sub-model time steps
    ************************************************/

    (*snow) = step_snow;
    snow->vapor_flux = store_vapor_flux;
    snow->blowing_flux = store_blowing_flux;
    snow->surface_flux = store_surface_flux;
    snow->canopy_vapor_flux = store_canopy_vapor_flux;
    (*Melt) = store_melt;
    snow->melt = store_melt;
    ppt = store_ppt;

    /******************************************************
       Store energy flux averages for sub-model time steps
    ******************************************************/

    (*energy) = soil_energy;
    energy->AlbedoOver = store_AlbedoOver / (double) N_steps;
    energy->AlbedoUnder = store_AlbedoUnder / (double) N_steps;
    energy->AtmosLatent = store_AtmosLatent / (double) N_steps;
    energy->AtmosLatentSub = store_AtmosLatentSub / (double) N_steps;
    energy->AtmosSensible = store_AtmosSensible / (double) N_steps;
    energy->LongOverIn = store_LongOverIn / (double) N_steps;
    energy->LongUnderIn = store_LongUnderIn / (double) N_steps;
    energy->LongUnderOut = store_LongUnderOut / (double) N_steps;
    energy->NetLongAtmos = store_NetLongAtmos / (double) N_steps;
    energy->NetLongOver = store_NetLongOver / (double) N_steps;
    energy->NetLongUnder = store_NetLongUnder / (double) N_steps;
    energy->NetShortAtmos = store_NetShortAtmos / (double) N_steps;
    energy->NetShortGrnd = store_NetShortGrnd / (double) N_steps;
    energy->NetShortOver = store_NetShortOver / (double) N_steps;
    energy->NetShortUnder = store_NetShortUnder / (double) N_steps;
    energy->ShortOverIn = store_ShortOverIn / (double) N_steps;
    energy->ShortUnderIn = store_ShortUnderIn / (double) N_steps;
    energy->advected_sensible = store_advected_sensible / (double) N_steps;
    energy->canopy_advection = store_canopy_advection / (double) N_steps;
    energy->canopy_latent = store_canopy_latent / (double) N_steps;
    energy->canopy_latent_sub = store_canopy_latent_sub / (double) N_steps;
    energy->canopy_refreeze = store_canopy_refreeze / (double) N_steps;
    energy->canopy_sensible = store_canopy_sensible / (double) N_steps;
    energy->deltaH = store_deltaH / (double) N_steps;
    energy->fusion = store_fusion / (double) N_steps;
    energy->grnd_flux = store_grnd_flux / (double) N_steps;
    energy->latent = store_latent / (double) N_steps;
    energy->latent_sub = store_latent_sub / (double) N_steps;
    energy->melt_energy = store_melt_energy / (double) N_steps;
    energy->sensible = store_sensible / (double) N_steps;
    if (snow->snow || INCLUDE_SNOW) {
        energy->advection = store_advection / (double) N_steps;
        energy->deltaCC = store_deltaCC / (double) N_steps;
        energy->refreeze_energy = store_refreeze_energy / (double) N_steps;
        energy->snow_flux = store_snow_flux / (double) N_steps;
    }
    energy->Tfoliage = snow_energy.Tfoliage;
    energy->Tfoliage_fbflag = snow_energy.Tfoliage_fbflag;
    energy->Tfoliage_fbcount = snow_energy.Tfoliage_fbcount;

    /**********************************************************
       Store vegetation variable sums for sub-model time steps
    **********************************************************/

    if (iveg != Nveg) {
        veg_var->throughfall = store_throughfall;
        veg_var->canopyevap = store_canopyevap;
        if (snow->snow) {
            veg_var->Wdew = snow_veg_var.Wdew;
        }
        else {
            veg_var->Wdew = soil_veg_var.Wdew;
        }
    }

    /**********************************************************
       Store soil layer variables for sub-model time steps
    **********************************************************/

    for (lidx = 0; lidx < Nlayers; lidx++) {
        layer[lidx] = step_layer[lidx];
        layer[lidx].evap = store_layerevap[lidx];
    }
    if (store_aero_cond_used[0] > 0 && store_aero_cond_used[0] <
        param.HUGE_RESIST) {
        aero_resist_used[0] = 1 / (store_aero_cond_used[0] / (double) N_steps);
    }
    else if (store_aero_cond_used[0] >= param.HUGE_RESIST) {
        aero_resist_used[0] = 0;
    }
    else {
        aero_resist_used[0] = param.HUGE_RESIST;
    }
    if (store_aero_cond_used[1] > 0 && store_aero_cond_used[1] <
        param.HUGE_RESIST) {
        aero_resist_used[1] = 1 / (store_aero_cond_used[1] / (double) N_steps);
    }
    else if (store_aero_cond_used[1] >= param.HUGE_RESIST) {
        aero_resist_used[1] = 0;
    }
    else {
        aero_resist_used[1] = param.HUGE_RESIST;
    }
    cell->pot_evap = store_pot_evap;

    /**********************************************************
       Store carbon cycle variable sums for sub-model time steps
    **********************************************************/

    if (options.CARBON && iveg != Nveg) {
        veg_var->rc = 1 / store_gc / (double) N_steps;
        for (cidx = 0; cidx < options.Ncanopy; cidx++) {
            veg_var->rsLayer[cidx] = 1 / store_gsLayer[cidx] / (double) N_steps;
        }
        veg_var->Ci = store_Ci / (double) N_steps;
        veg_var->GPP = store_GPP / (double) N_steps;
        veg_var->Rdark = store_Rdark / (double) N_steps;
        veg_var->Rphoto = store_Rphoto / (double) N_steps;
        veg_var->Rmaint = store_Rmaint / (double) N_steps;
        veg_var->Rgrowth = store_Rgrowth / (double) N_steps;
        veg_var->Raut = store_Raut / (double) N_steps;
        veg_var->NPP = store_NPP / (double) N_steps;

        free((char *) (store_gsLayer));

        soil_carbon_balance(soil_con, energy, cell, veg_var);

        // Update running total annual NPP
        if (veg_var->NPP > 0) {
            veg_var->AnnualNPP += veg_var->NPP * CONST_MWC / MOLE_PER_KMOLE *
                                  gp->dt;
        }
    }

    /********************************************************
       Compute Runoff, Baseflow, and Soil Moisture Transport
    ********************************************************/

    (*inflow) = ppt;

    ErrorFlag = runoff(cell, energy, soil_con, ppt, soil_con->frost_fract,
                       options.Nnode);

    return(ErrorFlag);
}
