ldndc::EcHy Class Reference

Watercycle model EcosystemHydrology - EcHy. More...

Inherits ldndc::MBE_LegacyModel.

Public Member Functions
lerr_t	solve ()

Private Member Functions
lerr_t	EcHyIrrigation ()
	Irrigation.
lerr_t	EcHyFlood ()
	sets hydrologic conditions during flooding events, e.g., More...
lerr_t	EcHyPercolation (size_t, size_t)
	Calculates water percolation within the soil profile.
lerr_t	EcHyBypassFlow (double, double &)
	Calculates water percolation within the soil profile.
lerr_t	EcHyEvapotranspiration ()
	Calculates evapotranspiration within the soil profile.
lerr_t	EcHySubsl2 (double const &, double const &, int const &, double const &, double &)
	... More...
lerr_t	EcHyBalanceCheck (double &)
	... More...
void	EcHyreset ()
	...
lerr_t	EcHyGroundwater ()
	... More...
lerr_t	EcHyIntegration ()
	... More...
lerr_t	EcHyStepInit ()
	...
lerr_t	EcHyStepExit ()
	...
double	EcHyGetInterceptionCapacity ()
	...
lerr_t	EcHyCalculateLeafWaterDistribution (lvector_t< double > &, lvector_t< double > &)
	...
lerr_t	EcHyPotentialEvapotranspiration ()
	Calculates potential evapotranspiration. Specific concept can be given as model option.
lerr_t	EcHySnowIce ()
	Calls SnowDNDC for the calculation of snowpack and soil ice formation.
double	EcHyGetMinimumWater (size_t)
	...
double	EcHyGetWiltingPoint (size_t)
	...
double	EcHyGetAvailableWaterTranspiration (size_t)
	...
double	EcHyGetAvailableWaterEvaporation (size_t)
	...
double	EcHyGetRootLimitation (size_t _sl, double)
	...
double	EcHyGetWaterLimitationTranspiration (double, size_t)
	...
double	EcHySoilWaterChange (size_t)

Private Attributes
double	kst_bottom
double	gw_depth_static
double	accumulated_potentialtranspiration_old
double	accumulated_irrigation_old
lvector_t< double >	trwl_sl
lvector_t< double >	wlfl_sl
lvector_t< double >	cr_fill_sl
lvector_t< double >	gw_fill_sl
lvector_t< double >	bypass_fill_sl
double	ev_surfacewater
lvector_t< double >	evsws_sl
lvector_t< double >	kst_sl
lvector_t< double >	wl_sl
lvector_t< double >	wlfc_sl
lvector_t< double >	wlwp_sl
double	ev_leaf
double	gw_fill_surface
double	cr_fill_groundwater
double	runoff
double	throughfall_snow
double	throughfall_water
double	canopy_snowmelt
double	canopy_water
double	thornthwaite_heat_index
double	daily_potential_evapotranspiration
double	daily_potential_leaf_evaporation
double	daily_potential_soil_evaporation
double	daily_potential_transpiration
WaterCycleSnowDNDC::IceContentStateIn	m_icecontent_in
cbm::string_t	evapotranspiration_method
cbm::string_t	runoff_method

Detailed Description

Watercycle model EcosystemHydrology - EcHy.

Member Function Documentation

EcHyBalanceCheck()

lerr_t ldndc::EcHy::EcHyBalanceCheck ( double &  _balance )

private

...

Checks balance between all incoming and outgoing water fluxes.

References canopy_snowmelt, canopy_water, cr_fill_groundwater, ev_leaf, ev_surfacewater, evsws_sl, gw_fill_sl, gw_fill_surface, runoff, throughfall_snow, throughfall_water, trwl_sl, wl_sl, and wlfl_sl.

Referenced by solve().

{
    double  balance( canopy_water + surface_water + canopy_snow + surface_snow);
    for ( size_t sl = 0;  sl < soillayers_in->soil_layer_cnt();  ++sl)
    {
        balance += wl_sl[sl] + wc_.ice_sl[sl] * sc_.h_sl[sl];
    }

    if ( _balance > 0.0)
    {
        balance +=  wlfl_sl[soillayers_in->soil_layer_cnt()]
                    + ev_leaf + ev_surfacewater + evsws_sl.sum() + trwl_sl.sum() + runoff
                    - gw_fill_surface - cr_fill_groundwater -  gw_fill_sl.sum()
                    - irrigation - throughfall_water - throughfall_snow;

        double const balance_delta( std::abs( _balance - balance));
        if ( cbm::flt_greater( balance_delta, 1.0e-4))
        {
            KLOGWARN( "Water leakage in: ", name(),
                      " Difference: ", balance - _balance);
            KLOGWARN("throughfall_water ",throughfall_water);
            KLOGWARN("throughfall_snow ",throughfall_snow);

            KLOGWARN("canopy_snowmelt ",canopy_snowmelt);
            KLOGWARN("surface_snowmelt ",surface_snowmelt);
            return LDNDC_ERR_FAIL;
        }
    }
    else
    {
        _balance = balance;
    }
    
    return  LDNDC_ERR_OK;
}

Here is the caller graph for this function:

EcHyFlood()

lerr_t ldndc::EcHy::EcHyFlood ( )

private

sets hydrologic conditions during flooding events, e.g.,

Flooding.

• surface water table
• bund height

References kst_sl, and wl_sl.

Referenced by solve().

{
    lerr_t  rc = m_eventflood.solve();
    if ( rc != LDNDC_ERR_OK){ return  rc; }

    bund_height = m_eventflood.get_bund_height();

    /* todo */
    //irrigation_amount = m_eventflood.get_irrigation_amount();

    /* todo */
    //double const irrigation_height = m_eventflood.get_irrigation_height();
    double const max_percolation = m_eventflood.get_maximum_percolation();
    double const water_table_flooding = m_eventflood.get_water_table();


    irrigation_switch = NONE;
    minimum_watertable_height = 0.0;
    if ( cbm::is_valid( water_table_flooding))
    {
        if ( cbm::flt_greater_zero( water_table_flooding))
        {
            irrigation_switch = CONSTANT_POSITIVE_WATER_TABLE;
            minimum_watertable_height = water_table_flooding;
            bund_height = water_table_flooding;
        }
        else
        {
            irrigation_switch = CONSTANT_NEGATIVE_WATER_TABLE;
            minimum_watertable_height = water_table_flooding;
            bund_height = 0.0;
        }
    }

    bool cracked( false);
    for (size_t sl = 0; sl < soillayers_in->soil_layer_cnt()-1; sl++)
    {
        if ( cbm::flt_greater_equal( sc_.depth_sl[sl], 0.3))
        {
            if ( cbm::is_valid( max_percolation))
            {
                if ((kst_eq_sl[sl] > max_percolation / lclock()->time_resolution()))
                {
                    
                    kst_sl[sl] = cbm::bound_min( 0.0,
                                                kst_eq_sl[sl] -
                                                (kst_eq_sl[sl] - max_percolation / lclock()->time_resolution()) *
                                                (sc_.depth_sl[sl] - 0.2) / (sc_.depth_sl[soillayers_in->soil_layer_cnt()-1] - 0.2));
                }
                kst_sl[sl] = max_percolation / lclock()->time_resolution();
            }
            else if ( wl_sl[sl] < -0.8 * pore_space( sl))
            {
                cracked = true;
            }

            if ( cracked)
            {
                double const time_rate( 0.1 / lclock()->time_resolution());
                kst_sl[sl] = cbm::bound_min( 0.0,
                                             kst_sl[sl] - (kst_sl[sl] - kst_eq_sl[sl]) * time_rate);
            }
            break;
        }
    }

    return  LDNDC_ERR_OK;
}

Here is the caller graph for this function:

EcHyGroundwater()

lerr_t ldndc::EcHy::EcHyGroundwater ( )

private

...

Negative groundwater table represents water on the soil surface

Negative groundwater table represents water on the soil surface

References cr_fill_groundwater, cr_fill_sl, EcHySoilWaterChange(), gw_depth_static, gw_fill_sl, gw_fill_surface, kst_sl, and wl_sl.

Referenced by solve().

{
    cbm::invalidate( kst_bottom_gw_vertical);
    if ( cbm::flt_greater( sc_.depth_sl[soillayers_in->soil_layer_cnt()-1], gw_depth_static))
    {
        if ( cbm::flt_greater( -gw_depth_static, surface_water))
        {
            gw_fill_surface = -gw_depth_static - surface_water;
        }

        //groundwater access
        {
            size_t sl( 0);
            double const layer_midpoint_depth( 0.5 * sc_.h_sl[sl]);
            if ( cbm::flt_greater_equal( layer_midpoint_depth, gw_depth_static))
            {
                gw_fill_sl[sl] = cbm::bound_min( 0.0,
                                                 0.999 * pore_space( sl) - wc_.ice_sl[sl] * sc_.h_sl[sl]
                                                 - (wl_sl[sl] + EcHySoilWaterChange( sl)));
                wc_.accumulated_groundwater_access += gw_fill_sl[sl];
            }
            else
            {
                double const wl_new( wl_sl[sl] + EcHySoilWaterChange( sl));
                double const wc_new( cbm::bound( wc_.wc_res_sl[sl], wl_new / sc_.h_sl[sl], wc_.wc_sat_sl[sl]));
                kst_bottom_gw_vertical = cbm::bound_min( 0.0,
                                                    ldndc::hydraulic_conductivity(
                                                                              wc_new,
                                                                              wc_.vgm_sl[sl],
                                                                              wc_.mvg_tau_sl[sl],
                                                                              wc_.wc_sat_sl[sl],
                                                                              wc_.wc_res_sl[sl],
                                                                              kst_sl[sl]));
                gw_fill_sl[sl] = 0.0;
            }
        }

        for (size_t sl = 1; sl < soillayers_in->soil_layer_cnt()-1; sl++)
        {
            double const layer_midpoint_depth( sc_.depth_sl[sl] - 0.5 * sc_.h_sl[sl]);
            if ( cbm::flt_greater_equal( layer_midpoint_depth, gw_depth_static))
            {
                gw_fill_sl[sl] = cbm::bound_min( 0.0,
                                                 0.999 * pore_space( sl) - wc_.ice_sl[sl] * sc_.h_sl[sl]
                                                 - (wl_sl[sl] + EcHySoilWaterChange( sl)));
                wc_.accumulated_groundwater_access += gw_fill_sl[sl];
            }
            else
            {
                double const wl_new( wl_sl[sl] + EcHySoilWaterChange( sl));
                double const wc_new( cbm::bound( wc_.wc_res_sl[sl], wl_new / sc_.h_sl[sl], wc_.wc_sat_sl[sl]));
                kst_bottom_gw_vertical = cbm::bound_min( 0.0,
                                                    ldndc::hydraulic_conductivity(
                                                                              wc_new,
                                                                              wc_.vgm_sl[sl],
                                                                              wc_.mvg_tau_sl[sl],
                                                                              wc_.wc_sat_sl[sl],
                                                                              wc_.wc_res_sl[sl],
                                                                              kst_sl[sl]));
                gw_fill_sl[sl] = 0.0;
            }
        }

        {
            size_t sl( soillayers_in->soil_layer_cnt()-1);
            double const layer_midpoint_depth( sc_.depth_sl[sl] - 0.5 * sc_.h_sl[sl]);
            if ( cbm::flt_greater_equal( layer_midpoint_depth, gw_depth_static))
            {
                gw_fill_sl[sl] = cbm::bound_min( 0.0,
                                                 0.999 * pore_space( sl) - wc_.ice_sl[sl] * sc_.h_sl[sl]
                                                 - (wl_sl[sl] + EcHySoilWaterChange( sl)));
                wc_.accumulated_groundwater_access += gw_fill_sl[sl];
            }
            else
            {
                double const wl_new( wl_sl[sl] + EcHySoilWaterChange( sl));
                double const wc_new( cbm::bound( wc_.wc_res_sl[sl], wl_new / sc_.h_sl[sl], wc_.wc_sat_sl[sl]));
                kst_bottom_gw_vertical = cbm::bound_min( 0.0,
                                                    ldndc::hydraulic_conductivity(
                                                                              wc_new,
                                                                              wc_.vgm_sl[sl],
                                                                              wc_.mvg_tau_sl[sl],
                                                                              wc_.wc_sat_sl[sl],
                                                                              wc_.wc_res_sl[sl],
                                                                              kst_sl[sl]));
                gw_fill_sl[sl] = 0.0;
            }
        }
    }

    //capillary rise
    if ( m_param->WCDNDC_HAVE_CAPILLARY_ACTION())
    {
        size_t const STEPS( 10);
        for (size_t steps = 0; steps < STEPS; steps++)
        {
            for (size_t sl = 0; sl < soillayers_in->soil_layer_cnt()-1; sl++)
            {
                double const layer_midpoint_depth( sc_.depth_sl[sl] - 0.5 * sc_.h_sl[sl]);
                if ( cbm::flt_greater_equal( layer_midpoint_depth, gw_depth_static))
                {
                    break;
                }
                
                /* Water contents with until now calculated water balance */
                double const wl_new_0( wl_sl[sl] + EcHySoilWaterChange( sl));
                double const wl_new_1( wl_sl[sl+1] + EcHySoilWaterChange( sl+1));
                
                double const wc_vg_0( cbm::bound( wc_.wc_res_sl[sl], wl_new_0 / sc_.h_sl[sl], wc_.wc_sat_sl[sl]));
                double const wc_vg_1( cbm::bound( wc_.wc_res_sl[sl+1], wl_new_1 / sc_.h_sl[sl+1], wc_.wc_sat_sl[sl+1]));
                
                /* Unsaturated hydraulic conductivity with until now calculated water balance */
                double const kust_0( cbm::bound_min( 0.0,
                                                ldndc::hydraulic_conductivity(
                                                                          wc_vg_0, wc_.vgm_sl[sl], wc_.mvg_tau_sl[sl],
                                                                          wc_.wc_sat_sl[sl], wc_.wc_res_sl[sl],
                                                                          1.0/double(STEPS)*kst_sl[sl])));
                double const kust_1( cbm::bound_min( 0.0,
                                                ldndc::hydraulic_conductivity(
                                                                          wc_vg_1, wc_.vgm_sl[sl+1], wc_.mvg_tau_sl[sl],
                                                                          wc_.wc_sat_sl[sl+1], wc_.wc_res_sl[sl+1],
                                                                          1.0/double(STEPS)*kst_sl[sl+1])));

                double const capillary_pressure_0( ldndc::capillary_pressure(
                                                                         wl_new_0/sc_.h_sl[sl], wc_.vga_sl[sl],
                                                                         wc_.vgn_sl[sl], wc_.vgm_sl[sl],
                                                                         wc_.wc_sat_sl[sl], wc_.wc_res_sl[sl]));
                double const capillary_pressure_1( ldndc::capillary_pressure(
                                                                         wl_new_1/sc_.h_sl[sl+1], wc_.vga_sl[sl+1],
                                                                         wc_.vgn_sl[sl+1], wc_.vgm_sl[sl+1],
                                                                         wc_.wc_sat_sl[sl], wc_.wc_res_sl[sl+1]));

                double const kust( cbm::flt_greater( capillary_pressure_0, capillary_pressure_1) ? kust_1 : kust_0 );
                double const flow( (capillary_pressure_1 - capillary_pressure_0) / (0.5 * (sc_.h_sl[sl] + sc_.h_sl[sl+1])) * kust);

                /* downward flow*/
                if ( cbm::flt_greater_zero( flow))
                {
                    double const left( wl_sl[sl] + EcHySoilWaterChange( sl) - (wc_.wc_res_sl[sl] * sc_.h_sl[sl]));
                    double const space( 0.999 * pore_space( sl+1) - wc_.ice_sl[sl+1] * sc_.h_sl[sl+1]
                                        - (wl_sl[sl+1] + EcHySoilWaterChange( sl+1)));

                    if ( cbm::flt_greater_zero( left) &&
                         cbm::flt_greater_zero( space))
                    {
                        double const bound_flow( cbm::bound_max( flow, std::min( left, space)) / static_cast<double>(STEPS));
                        cr_fill_sl[sl] -= bound_flow;
                        cr_fill_sl[sl+1] += bound_flow;
                    }
                }
                /* upward flow */
                else if ( cbm::flt_less( flow, 0.0))
                {
                    double const left( wl_sl[sl+1] + EcHySoilWaterChange( sl+1) - (wc_.wc_res_sl[sl+1] * sc_.h_sl[sl+1]));
                    double const space( 0.999 * pore_space( sl) - wc_.ice_sl[sl] * sc_.h_sl[sl]
                                        - (wl_sl[sl] + EcHySoilWaterChange( sl)));

                    if ( cbm::flt_greater_zero( left) &&
                         cbm::flt_greater_zero( space))
                    {
                        double const bound_flow( cbm::bound_max( -flow, std::min( left, space)) / static_cast<double>(STEPS));
                        cr_fill_sl[sl] += bound_flow;
                        cr_fill_sl[sl+1] -= bound_flow;
                    }
                }
            }

            size_t sl( soillayers_in->soil_layer_cnt()-1);
            double const layer_midpoint_depth( sc_.depth_sl[sl] - 0.5 * sc_.h_sl[sl]);
            if ( cbm::flt_greater( gw_depth_static, layer_midpoint_depth))
            {
                /* Water contents with until now calculated water balance */
                double const wl_new( wl_sl[sl] + EcHySoilWaterChange( sl));
                
                double const wc_vg( cbm::bound( wc_.wc_res_sl[sl], wl_new / sc_.h_sl[sl], wc_.wc_sat_sl[sl]));
                
                /* Unsaturated hydraulic conductivity with until now calculated water balance */
                double const kust( cbm::bound_min( 0.0,
                                              ldndc::hydraulic_conductivity(
                                                                        wc_vg, wc_.vgm_sl[sl], wc_.mvg_tau_sl[sl],
                                                                        wc_.wc_sat_sl[sl], wc_.wc_res_sl[sl],
                                                                        1.0/double(STEPS)*kst_sl[sl])));

                double const capillary_pressure( ldndc::capillary_pressure(
                                                                       wl_new/sc_.h_sl[sl], wc_.vga_sl[sl],
                                                                       wc_.vgn_sl[sl], wc_.vgm_sl[sl],
                                                                       wc_.wc_sat_sl[sl], wc_.wc_res_sl[sl]));

                /* upward flow */
                double const flow( -capillary_pressure / (gw_depth_static - layer_midpoint_depth) * kust);
                if ( cbm::flt_less( flow, 0.0))
                {
                    double const space( 0.999 * pore_space( sl) - wc_.ice_sl[sl] * sc_.h_sl[sl]
                                       - (wl_sl[sl] + EcHySoilWaterChange( sl)));
                    if ( cbm::flt_greater_zero( space))
                    {
                        double const bound_flow( cbm::bound_max( -flow, 1.0 / static_cast<double>(STEPS) * space));
                        cr_fill_sl[sl] += bound_flow;
                        cr_fill_groundwater += bound_flow;
                    }
                }
            }
        }
    }

    return LDNDC_ERR_OK;
}

Here is the call graph for this function:

Here is the caller graph for this function:

EcHyIntegration()

lerr_t ldndc::EcHy::EcHyIntegration ( )

private

...

Integrates all water state variables.

References bypass_fill_sl, canopy_snowmelt, canopy_water, EcHySoilWaterChange(), ev_surfacewater, gw_fill_surface, runoff, throughfall_snow, throughfall_water, wl_sl, and wlfl_sl.

Referenced by solve().

{
    canopy_snow = cbm::bound_min( 0.0,
                                  canopy_snow - canopy_snowmelt);
    surface_snow = cbm::bound_min( 0.0,
                                   surface_snow
                                   + throughfall_snow
                                   - surface_snowmelt);

    canopy_water = cbm::bound_min( 0.0,
                                   canopy_water
                                   + canopy_snowmelt);
    surface_water = cbm::bound_min( 0.0,
                                    surface_water
                                    + throughfall_water + surface_snowmelt + gw_fill_surface + irrigation
                                    - ev_surfacewater - wlfl_sl[0] - runoff - bypass_fill_sl.sum());

    for (size_t sl = 0; sl < soillayers_in->soil_layer_cnt(); sl++)
    {
        wl_sl[sl] = cbm::bound_min( 0.0, wl_sl[sl] + EcHySoilWaterChange( sl));
        wc_.mskpa_sl[sl] = ldndc::hydrology::capillary_pressure( wl_sl[sl]/sc_.h_sl[sl], wc_.vga_sl[sl],
                                                                 wc_.vgn_sl[sl], wc_.vgm_sl[sl],
                                                                 wc_.wc_sat_sl[sl], wc_.wc_res_sl[sl]);
    }

    return  LDNDC_ERR_OK;
}

Here is the call graph for this function:

Here is the caller graph for this function:

EcHySoilWaterChange()

double ldndc::EcHy::EcHySoilWaterChange ( size_t  _sl )

private

Integrates all water state variables.

References cr_fill_sl, evsws_sl, gw_fill_sl, trwl_sl, and wlfl_sl.

Referenced by EcHyGetAvailableWaterEvaporation(), EcHyGetAvailableWaterTranspiration(), EcHyGroundwater(), EcHyIntegration(), and EcHyIrrigation().

{
    return wlfl_sl[_sl] - wlfl_sl[_sl+1] + cr_fill_sl[_sl] - trwl_sl[_sl] - evsws_sl[_sl] + gw_fill_sl[_sl];
}

Here is the caller graph for this function:

EcHySubsl2()

lerr_t ldndc::EcHy::EcHySubsl2 ( double const &  _cp, double const &  _delta_z, int const &  _sl, double const &  _saturation, double &  _flow  )

private

...

This routine calculates the rate of capillary flow or percolation between groundwater table and root zone. The stationary flow is found by integration of dZL = K.d(MH)/(K + FLW), where Z= height above groundwater, MH= matric head, K= conductivity and FLW= chosen flow. In an iteration loop the correct flow is found. The integration goes at most over four intervals: [0,45], [45,170], [170,330] and [330,MH-rootzone] (last one on logarithmic scale).

Author

C. Rappoldt M. Wopereis

Date

January 1986, revised June 1990

Chapter 15 in documentation WOFOST Version 4.1 (1988)

References kst_sl.

{
    //calculation of pf from matrix head
    double pf( std::log10( _cp));

    //in case of small matric head (high water contents)
    if ( cbm::flt_greater( 1.0, pf))
    {
        _flow = 0.0;
        return LDNDC_ERR_OK;
    }

    double elog10( 2.302585);
    double logst4( 2.518514);

    double start[4] = {0.0, 45.0, 170.0, 330.0};

    //number and width of integration intervals
    int iint( 0);
    double del[4] = { 0.0, 0.0, 0.0, 0.0};
    for (int i1 = 0; i1 < 4; i1++)
    {
        if (i1 < 3)
        {
            del[i1] = std::min( start[i1+1], _cp) - start[i1];
        }
        else
        {
            del[i1] = pf - logst4;
        }
        if (del[i1] <= 0.0)
        {
            break;
        }
        iint += 1;
    }

    //preparation of three-point gaussian integration
    double hulp[12] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    double conduc[12] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    for (int i1 = 0; i1 < iint; i1++)
    {
        for (int i2 = 0; i2 < 3; i2++)
        {
            double pgau[3] = {0.1127016654, 0.5, 0.8872983346};
            double wgau[3] = {0.2777778, 0.4444444, 0.2777778};

            int i3( 3 * i1 + i2);

            //the three points in the full-width intervals are standard
            double pfstan[9] = {0.705143, 1.352183, 1.601282, 1.771497, 2.031409, 2.192880, 2.274233, 2.397940, 2.494110};
            double pfgau( pfstan[i3]);

            //the three points in the last interval are calculated
            if ( i1 == (iint-1))
            {
                if ( iint < 4)
                {
                    pfgau = std::log10( start[iint-1] + pgau[i2] * del[iint-1]);
                }
                else
                {
                    pfgau = logst4 + pgau[i2] * del[iint-1];
                }
            }

            double const wcl_loc( ldndc::hydrology::water_content( std::exp(elog10 * pfgau), wc_.vga_sl[_sl], wc_.vgn_sl[_sl], wc_.vgm_sl[_sl], _saturation, 0.1));
            conduc[i3] = ldndc::hydraulic_conductivity( wcl_loc, wc_.vgm_sl[_sl], wc_.mvg_tau_sl[_sl],
                                                        _saturation, wc_.wc_wp_sl[_sl], kst_sl[_sl] * cbm::CM_IN_M);
            hulp[i3] = del[i1] * wgau[i2] * conduc[i3];
            if (i3 > 8)
            {
                hulp[i3] = hulp[i3] * elog10 * std::exp( elog10 * pfgau);
            }
        }
    }


    //setting upper and lower limit
    double const wcl_loc( ldndc::hydrology::water_content( _cp, wc_.vga_sl[_sl], wc_.vgn_sl[_sl], wc_.vgm_sl[_sl], _saturation, 0.1));
    double const kms( ldndc::hydraulic_conductivity( wcl_loc, wc_.vgm_sl[_sl], wc_.mvg_tau_sl[_sl],
                                                     _saturation, wc_.wc_wp_sl[_sl], kst_sl[_sl] * cbm::CM_IN_M));
    double fu( (_cp <= _delta_z) ? 0.0 : 1.27);
    double fl( (_cp >= _delta_z) ? 0.0 : -kms);


    if ( !cbm::flt_equal( _cp, _delta_z))
    {
        //iteration loop
        int imax( 3 * iint);
        for ( int i1 = 0; i1 < 14; i1++)
        {
            double flw( (fu + fl) / 2.0);
            double df( (fu - fl) / 2.0);
            if ( (df < 0.01) &&
                ((df / std::abs(flw)) < 0.1))
            {
                break;
            }

            double z( 0.0);
            for ( int i2 = 0; i2 < imax; i2++)
            {
                z  = z + hulp[i2] / (conduc[i2] + flw);
            }
            if (z >= _delta_z)
            {
                fl = flw;
            }
            if (z <= _delta_z)
            {
                fu = flw;
            }
        }
    }

    //flow output in mm
    _flow = 10.0 * (fu + fl) / 2.0;

    return LDNDC_ERR_OK;
}

solve()

lerr_t ldndc::EcHy::solve

(

)

Kicks off computation for one time step.

References EcHyBalanceCheck(), EcHyBypassFlow(), EcHyEvapotranspiration(), EcHyFlood(), EcHyGroundwater(), EcHyIntegration(), EcHyIrrigation(), EcHyPercolation(), EcHyreset(), EcHySnowIce(), EcHyStepExit(), EcHyStepInit(), kst_bottom, kst_sl, runoff, runoff_method, throughfall_water, and wlfl_sl.

{
    /* reset of all fluxes */
    EcHyreset();

    /* flooding events */
    lerr_t rc = EcHyFlood();
    if ( rc != LDNDC_ERR_OK)
    {
        KLOGERROR("Handling flooding event failed in: ", name());
        return  rc;
    }

    /* tasks:
     *  - update internal state (e.g., throughfall_water, ...)
     *  - calculate evapotranspiration on demand
     */
    rc = EcHyStepInit();
    if ( rc != LDNDC_ERR_OK)
    { return  rc; }
    
    /* begin water balance */
    double balance( -1.0);
    EcHyBalanceCheck( balance);

    /* snow and ice calculated */
    EcHySnowIce();

    /* 
     * apply irrigation triggered by
     *  - irrigation events
     *  - flooding events
     */
    EcHyIrrigation();

    double const potential_infiltration( surface_water + irrigation + throughfall_water + surface_snowmelt);
    double bypass_flow( m_param->CRACK_FRACTION() * potential_infiltration);

    /* bypass flow */
    EcHyBypassFlow( m_param->CRACK_DEPTH(), bypass_flow);
    double const infiltration( potential_infiltration - bypass_flow);

    /* downward water flow */
    wlfl_sl[0] = cbm::bound_max( infiltration, kst_sl[0]);
    for ( size_t sl = 0; sl < soillayers_in->soil_layer_cnt();  ++sl)
    {
        EcHyPercolation( 1, sl);
    }

    if ( cbm::is_valid( kst_bottom))
    {
        wlfl_sl[soillayers_in->soil_layer_cnt()] = cbm::bound_max( wlfl_sl[soillayers_in->soil_layer_cnt()],
                                                                   kst_bottom);
    }
    if ( cbm::is_valid( kst_bottom_gw_lateral))
    {
        wlfl_sl[soillayers_in->soil_layer_cnt()] = cbm::bound_max( wlfl_sl[soillayers_in->soil_layer_cnt()],
                                                                   kst_bottom_gw_lateral);
    }
    if ( cbm::is_valid( kst_bottom_gw_vertical))
    {
        wlfl_sl[soillayers_in->soil_layer_cnt()] = cbm::bound_max( wlfl_sl[soillayers_in->soil_layer_cnt()],
                                                                   kst_bottom_gw_vertical);
    }

    EcHyEvapotranspiration();

    /* upward water flow */
    for ( int sl = soillayers_in->soil_layer_cnt()-1;  sl >= 0;  --sl)
    {
        EcHyPercolation( 2, sl);
    }

    /* runoff */
    if ( runoff_method == "curvenumber")
    {
        double const surface_water_in_mm( surface_water * cbm::MM_IN_M);
        double const curvenumber( 70.0);

        //retention parameter [mm]
        double const retention_factor( 25.4 * (1000.0 / curvenumber - 10));

        //runoff [mm]
        double const q_surf( cbm::sqr( surface_water_in_mm) /
                             ((surface_water_in_mm + retention_factor) * lclock()->time_resolution()));

        runoff = cbm::bound( 0.0,
                             q_surf * cbm::M_IN_MM,
                             0.99 * (surface_water + irrigation + throughfall_water + surface_snowmelt
                                     - wlfl_sl[0] - bypass_flow));
    }
    else
    {
        double const f_runoff( cbm::bound_max( m_param->FRUNOFF() / lclock()->time_resolution(), 1.0));
        runoff = cbm::bound_min( 0.0, f_runoff * (surface_water + irrigation + throughfall_water + surface_snowmelt
                                                  - wlfl_sl[0] - bypass_flow - bund_height));
    }

    /* groundwater interaction:
     *  - set soil layers water-saturated
     *  - capillary rise
     */
    EcHyGroundwater();

    /* integration of state variables */
    EcHyIntegration();

    /* update external state */
    EcHyStepExit();

    /* perform water balance */
    rc = EcHyBalanceCheck( balance);
    if ( rc){ return rc; }

    rc = EcHy_check_for_negative_value("exit");
    if ( rc){ return rc; }

    return  LDNDC_ERR_OK;
}

Here is the call graph for this function:

Member Data Documentation

accumulated_irrigation_old

double ldndc::EcHy::accumulated_irrigation_old

private

Accumulated irrigation of last time step [m]

Referenced by EcHyIrrigation().

accumulated_potentialtranspiration_old

double ldndc::EcHy::accumulated_potentialtranspiration_old

private

Accumulated potential transpiration of last time step [m]

Referenced by EcHyPotentialEvapotranspiration().

bypass_fill_sl

lvector_t< double > ldndc::EcHy::bypass_fill_sl

private

Soillayer water input due to bypass flow

Referenced by EcHyIntegration(), EcHyPotentialEvapotranspiration(), and EcHyreset().

canopy_snowmelt

double ldndc::EcHy::canopy_snowmelt

private

Snow melting

Referenced by EcHyBalanceCheck(), EcHyIntegration(), and EcHySnowIce().

canopy_water

double ldndc::EcHy::canopy_water

private

Amount of canopy intercepted water

Referenced by EcHyBalanceCheck(), EcHyCalculateLeafWaterDistribution(), EcHyIntegration(), EcHyPotentialEvapotranspiration(), EcHySnowIce(), and EcHyStepInit().

cr_fill_groundwater

double ldndc::EcHy::cr_fill_groundwater

private

Water addition to last soil layer from groundwater water by cappillary rise (only used for water balance calculation)

Referenced by EcHyBalanceCheck(), EcHyGroundwater(), and EcHyreset().

cr_fill_sl

lvector_t< double > ldndc::EcHy::cr_fill_sl

private

Soillayer water input due to capillary rise

Referenced by EcHyGroundwater(), EcHyreset(), and EcHySoilWaterChange().

daily_potential_evapotranspiration

double ldndc::EcHy::daily_potential_evapotranspiration

private

...

Referenced by EcHyPotentialEvapotranspiration(), and EcHyStepExit().

daily_potential_leaf_evaporation

double ldndc::EcHy::daily_potential_leaf_evaporation

private

...

Referenced by EcHyPotentialEvapotranspiration().

daily_potential_soil_evaporation

double ldndc::EcHy::daily_potential_soil_evaporation

private

...

Referenced by EcHyPotentialEvapotranspiration().

daily_potential_transpiration

double ldndc::EcHy::daily_potential_transpiration

private

...

Referenced by EcHyPotentialEvapotranspiration().

ev_leaf

double ldndc::EcHy::ev_leaf

private

Evaporation from plant surface

Referenced by EcHyBalanceCheck(), EcHyPotentialEvapotranspiration(), and EcHyStepExit().

ev_surfacewater

double ldndc::EcHy::ev_surfacewater

private

Evaporation from surface water

Referenced by EcHyBalanceCheck(), EcHyIntegration(), EcHyPotentialEvapotranspiration(), and EcHyreset().

evapotranspiration_method

cbm::string_t ldndc::EcHy::evapotranspiration_method

private

...

Referenced by EcHyPotentialEvapotranspiration().

evsws_sl

lvector_t< double > ldndc::EcHy::evsws_sl

private

Evaporation from soil layer

Referenced by EcHyBalanceCheck(), EcHyPotentialEvapotranspiration(), EcHyreset(), EcHySoilWaterChange(), and EcHyStepExit().

gw_depth_static

double ldndc::EcHy::gw_depth_static

private

Depth of groundwater table [m]

Referenced by EcHyGroundwater(), and EcHyStepInit().

gw_fill_sl

lvector_t< double > ldndc::EcHy::gw_fill_sl

private

Soillayer water input due to groundwater flow

Referenced by EcHyBalanceCheck(), EcHyGroundwater(), EcHyreset(), and EcHySoilWaterChange().

gw_fill_surface

double ldndc::EcHy::gw_fill_surface

private

Water addition to surface water from groundwater water boundary condition

Referenced by EcHyBalanceCheck(), EcHyGroundwater(), EcHyIntegration(), and EcHyreset().

kst_bottom

double ldndc::EcHy::kst_bottom

private

Saturated hydraulic conductivity below last soil layer [cm:min-1]

Referenced by EcHyStepInit(), and solve().

kst_sl

lvector_t< double > ldndc::EcHy::kst_sl

private

Saturated hydraulic conductivity

Referenced by EcHyFlood(), EcHyGroundwater(), EcHySubsl2(), and solve().

m_icecontent_in

WaterCycleSnowDNDC::IceContentStateIn ldndc::EcHy::m_icecontent_in

private

...

Referenced by EcHySnowIce().

runoff

double ldndc::EcHy::runoff

private

Lateral surface runoff

Referenced by EcHyBalanceCheck(), EcHyIntegration(), EcHyreset(), EcHyStepExit(), and solve().

runoff_method

cbm::string_t ldndc::EcHy::runoff_method

private

...

Referenced by solve().

thornthwaite_heat_index

double ldndc::EcHy::thornthwaite_heat_index

private

...

Referenced by EcHyPotentialEvapotranspiration().

throughfall_snow

double ldndc::EcHy::throughfall_snow

private

Snowfall

Referenced by EcHyBalanceCheck(), EcHyIntegration(), EcHyreset(), and EcHySnowIce().

throughfall_water

double ldndc::EcHy::throughfall_water

private

Throughfall to soil

Referenced by EcHyBalanceCheck(), EcHyIntegration(), EcHyPotentialEvapotranspiration(), EcHySnowIce(), EcHyStepExit(), EcHyStepInit(), and solve().

trwl_sl

lvector_t< double > ldndc::EcHy::trwl_sl

private

Water withdrawl by transpiration [m d−1]

Referenced by EcHyBalanceCheck(), EcHyPotentialEvapotranspiration(), EcHyreset(), EcHySoilWaterChange(), and EcHyStepExit().

wl_sl

lvector_t< double > ldndc::EcHy::wl_sl

private

Amount of water in soil layer

Referenced by EcHyBalanceCheck(), EcHyFlood(), EcHyGetAvailableWaterEvaporation(), EcHyGetAvailableWaterTranspiration(), EcHyGroundwater(), EcHyIntegration(), EcHyIrrigation(), EcHyPotentialEvapotranspiration(), EcHySnowIce(), EcHyStepExit(), and EcHyStepInit().

wlfc_sl

lvector_t< double > ldndc::EcHy::wlfc_sl

private

Amount of water in soil layer at field capacity

Referenced by EcHyGetWaterLimitationTranspiration(), and EcHyIrrigation().

wlfl_sl

lvector_t< double > ldndc::EcHy::wlfl_sl

private

Water flux at boundaries of soil layer [m d-1] the flux of water between soil layer is tracked by wlfl (mm d−1). in total, there are nl+1 flow rates, where nl is the number of soil layers, wlfl[1] is the flow rate between the ponded water layer and the soil surface, wlfl[2] is the flow rate between soil layer 1 and soil layer 2, etc.

Referenced by EcHyBalanceCheck(), EcHyIntegration(), EcHyPotentialEvapotranspiration(), EcHyreset(), EcHySoilWaterChange(), EcHyStepExit(), and solve().

wlwp_sl

lvector_t< double > ldndc::EcHy::wlwp_sl

private

Amount of water in soil layer at wilting point

Referenced by EcHyGetWiltingPoint().