ldndc::WatercycleDNDC Class Reference

Watercycle model WatercycleDNDC. More...

#include <models/watercycle/dndc/watercycle-dndc.h>

Public Member Functions
lerr_t	solve ()
	Kicks off computation for one time step.

Private Member Functions
lerr_t	event_irrigate ()
	considers water input due to irrigation and/or liquid manure evetns
lerr_t	event_flood ()
	sets hydrologic conditions during flooding events, e.g.,
double	balance_check (double *)
	checks each time step for correct water balance
void	CalcRaining (unsigned int)
	apportion daily rainfall to subdaily rainfall
double	rainfall_intensity ()
void	CalcIrrigation (unsigned int)
	apportion daily irrigation to subdaily irrigation
lerr_t	CalcIceContent ()
	Call of ice-content related processes.
lerr_t	CalcPotEvapoTranspiration (size_t)
	Potential evaporation.
void	CalcInterception ()
	Calculates water quantity that is retained on leaf and stem bark surfaces.
double	get_leaf_water ()
	Collects leaf water from all canopy layers.
double	get_impedance_factor (size_t)
	Reduces water fluxes out of a layer if ice lenses have formed.
lerr_t	CalcTranspiration ()
lerr_t	CalcTranspirationCouvreur (double const &hr_potTransinm)
void	CalcSnowEvaporation ()
void	CalcSurfaceWaterEvaporation ()
void	CalcSurfaceFlux ()
void	WatercycleDNDC_preferential_flow ()
void	CalcSoilEvapoPercolation ()
void	SoilWaterEvaporation (size_t)
double	WaterFlowCapillaryRise (size_t, double const &)
	Similar to WaterFlowBelowFieldCapacity(...) but restricts water flow by water availability from soil layer below.
double	WaterFlowAboveFieldCapacity (size_t, double const &)
double	get_fsl (size_t)
	Returns factor accounting for soil layer depth. Originally, a constant soil layer depth of 0.02 m has been assumed. Due to flexible layer depth, an additional factor has been introduced that accounts for a higher resistance of larger layers.
double	substep_day_fraction ()
	Length of a substep as fraction of a day.
double	get_clay_limitation (size_t _sl, double _water_avail, double _water_ref)
	Restrains water availability due to clay content.
double	hours_per_time_step ()

Private Attributes
double	thornthwaite_heat_index
	heat index for potential evaporation using approach after Thornthwaite
double	rainfall_intensity_cache

Detailed Description

Watercycle model WatercycleDNDC.

Member Function Documentation

balance_check()

double ldndc::WatercycleDNDC::balance_check ( double * _balance )

private

checks each time step for correct water balance

Balance check of water.

{
    double  balance = get_leaf_water() + wc_.surface_ice + wc_.surface_water;
    for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
    {
        balance += (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl];
    }
    balance += (   - m_wc.accumulated_irrigation_event_executed
                   - wc_.accumulated_irrigation_automatic
                   - wc_.accumulated_irrigation_ggcmi
                   - wc_.accumulated_irrigation_reservoir_withdrawal
                   - wc_.accumulated_precipitation
                   - wc_.accumulated_groundwater_access
                   + wc_.accumulated_percolation
                   + wc_.accumulated_runoff
                   + wc_.accumulated_wateruptake_sl.sum()
                   + wc_.accumulated_interceptionevaporation
                   + wc_.accumulated_soilevaporation
                   + wc_.accumulated_surfacewaterevaporation);
 
    if ( _balance)
    {
        double const  balance_delta( std::abs( *_balance - balance));
        if ( balance_delta > cbm::DIFFMAX)
        {
            KLOGWARN( "Water leakage:  difference is ", balance - *_balance);
        }
    }
 
    for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
    {
        if ( cbm::flt_less( wc_.wc_sl[sl], 0.0))
        {
            KLOGWARN( "Negative water content of: ", wc_.wc_sl[sl], " in layer: ", sl, ". Water content set to zero.");
            wc_.wc_sl[sl] = 0.0;
        }
        if ( cbm::flt_less( wc_.ice_sl[sl], 0.0))
        {
            KLOGWARN( "Negative ice content of: ", wc_.ice_sl[sl], " in layer: ", sl, ". Ice content set to zero.");
            wc_.ice_sl[sl] = 0.0;
        }
    }
 
    return  balance;
}

References balance_check(), and get_leaf_water().

Referenced by balance_check(), and solve().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcInterception()

void ldndc::WatercycleDNDC::CalcInterception ( )

private

Calculates water quantity that is retained on leaf and stem bark surfaces.

• Rainfall and interception are assumed to happen at the end of the timestep and evaporation is realised at once at the start of the timestep.

• The model assumes a homogeneous horizontal distribution of LAI although a non-linear dependency to crown cover (relative foliage clumping) has been reported (e.g. Teklehaimanot et al. 1991).

• The default water-interception for foliage almost equal to the largest species specific value

References get_impedance_factor(), and get_leaf_water().

Referenced by solve().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcIrrigation()

void ldndc::WatercycleDNDC::CalcIrrigation ( unsigned int _substep )

private

apportion daily irrigation to subdaily irrigation

Automatic irrigation to prevent water stress

GGCMI style irrigation

• fill up the water content of each soil layer to field capacity at the start of each day
• water supplied just injected into soil layer
• no interception, no run off, no surface-water
• only irrigate if a crop is present

Irrigation due to flooding

irrigation triggered by flooding management

Irrigation due to explicit irrigation event

{
    //reset irrigation
    m_wc.irrigation = 0.0;

    if (cbm::flt_in_range_lu(0.0, m_cfg.automaticirrigation, 1.0))
    {
        if (m_veg->size() > 0) //check crop is present
        {
            if (cbm::flt_greater(mc_.nd_airtemperature,10.0)) // only add water if temperature > 10C)
            {
                for (size_t l = 0; l < m_soillayers->soil_layer_cnt(); ++l)
                {
                    if (cbm::flt_greater_zero(m_veg->mfrt_sl(l))) //only count water need for soil layers that contain roots
                    {
                        // trigger irrigation events as soon as drought stress appears
                        double const wc_target(wc_.wc_wp_sl[l] + m_cfg.automaticirrigation * (wc_.wc_fc_sl[l] - wc_.wc_wp_sl[l]));
                        if (cbm::flt_less(wc_.wc_sl[l], wc_target) &&
                            cbm::flt_equal_zero(wc_.ice_sl[l])) //only add water if no ice is present
                        {
                            // irrigation assumed until target water content
                            double const irrigate_layer((wc_target - wc_.wc_sl[l]) * sc_.h_sl[l]);
                            m_wc.irrigation += irrigate_layer;
                            wc_.accumulated_irrigation_automatic += irrigate_layer;
                        }
                    }
                }
 
            } 
            else 
            {
            }
 
        }
    }
 
    if ( cbm::flt_in_range_lu( 0.0, m_eventflood.get_saturation_level(), 1.0) &&
         cbm::is_valid( m_eventflood.get_irrigation_height()) &&
         cbm::is_valid( m_eventflood.get_soil_depth()))
    {
        double wc_sum( 0.0);
        double wc_fc_sum( 0.0);
        double wc_wp_sum( 0.0);
        for ( size_t l = 0;  l < m_soillayers->soil_layer_cnt();  ++l)
        {
            wc_sum += wc_.wc_sl[l] * sc_.h_sl[l];
            wc_fc_sum += wc_.wc_fc_sl[l] * sc_.h_sl[l];
            wc_wp_sum += wc_.wc_wp_sl[l] * sc_.h_sl[l];
 
            if ( cbm::flt_greater_equal( sc_.depth_sl[l], m_eventflood.get_soil_depth()))
            {
                break;
            }
        }
 
        // trigger irrigation events as soon as drought stress appears
        double const wc_target( wc_wp_sum + m_eventflood.get_saturation_level() * (wc_fc_sum -  wc_wp_sum));
        if ( cbm::flt_less( wc_sum, wc_target))
        {
            m_wc.irrigation += m_eventflood.get_irrigation_height();
            wc_.accumulated_irrigation_automatic += m_eventflood.get_irrigation_height();
        }
    }

    if ( m_cfg.ggcmi_irrigation)
    {
        if( (lclock()->subday() == 1) && (_substep == 0))
        {
            if( m_veg->size() > 0 ) //check crop is present
            {
                double water_supply( 0.0);
                for ( size_t l = 0;  l < m_soillayers->soil_layer_cnt();  ++l)
                {
                    if ( cbm::flt_greater_zero( m_veg->mfrt_sl( l))) //only irrigate soil layers that contain roots
                    {
                        if( cbm::flt_less( wc_.wc_sl[l], wc_.wc_fc_sl[l]) &&
                            cbm::flt_equal_zero( wc_.ice_sl[l])) //only add water if no ice is present
                        {
                            water_supply += (wc_.wc_fc_sl[l] - wc_.wc_sl[l]) * sc_.h_sl[l]; //amount of water added (m)
                            wc_.wc_sl[l] = wc_.wc_fc_sl[l];
                        }
                    }
                }
                wc_.accumulated_irrigation_ggcmi += water_supply;
            }
        }
    }

    if ( cbm::is_valid( m_eventflood.get_water_table()))
    {
        /* Dynamic flooding */
        if ( cbm::is_valid( m_eventflood.get_irrigation_height()))
        {
            bool irrigate( false);
 
            /* Irrigate after surface water table drop */
            if ( cbm::flt_greater_equal_zero( m_eventflood.get_water_table()))
            {
                if ( cbm::flt_greater( m_eventflood.get_water_table(), wc_.surface_water))
                {
                    irrigate = true;
                }
            }
            else
            {
                if ( !cbm::flt_greater_zero( wc_.surface_water))
                {
                    double water_tot_sum( 0.0);
                    double water_fc_sum( 0.0);
                    for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
                    {
                        water_tot_sum += wc_.wc_sl[sl] * sc_.h_sl[sl];
                        water_fc_sum += wc_.wc_fc_sl[sl] * sc_.h_sl[sl];
                        if ( cbm::flt_greater_equal( sc_.depth_sl[sl], -m_eventflood.get_water_table()))
                        {
                            break;
                        }
                    }
                    // Note: The tipping bucket approach may block percolation into deeper soil layers
                    // if the upper layers have higher water content. This can lead to water reduction
                    // from evapotranspiration in deeper layers without replenishment from above.
                    // Workaround: AWD is triggered when the integrated topsoil water drops below
                    // 90% of the integrated topsoil field capacity.
                    if ( cbm::flt_less( water_tot_sum, 0.9 * water_fc_sum))
                    {
                        irrigate = true;
                    }
                }
            }
            // trigger AWD irrigation event
            if ( irrigate)
            {
                double water_supply( cbm::bound_min( 0.0, m_eventflood.get_irrigation_height() - wc_.surface_water));
                for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
                {
                    // set soil fully water saturated
                    if ( cbm::flt_greater( sc_.poro_sl[sl], wc_.wc_sl[sl]))
                    {
                        water_supply += ((sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]);
                    }
                }
 
                //remove rainfall (less irrigation water is supplied if it if raining)
                water_supply = cbm::bound_min( 0.0, water_supply - m_mc.precipitation_ts);
 
                if ( m_eventflood.have_unlimited_water())
                {
                    wc_.accumulated_irrigation_automatic += water_supply;
                }
                else
                {
                    if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
                    {
                        wc_.irrigation_reservoir -= water_supply;
                    }
                    else
                    {
                        water_supply = wc_.irrigation_reservoir;
                        wc_.irrigation_reservoir = 0.0;
                    }
                    wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
                }
                m_wc.irrigation += water_supply;
            }
        }
        /* Static flooding */
        else
        {
            if ( !cbm::flt_greater_zero( m_eventflood.get_water_table()))
            {
                double water_supply( 0.0);
                for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
                {
                    if ( cbm::flt_greater( sc_.depth_sl[sl], -m_eventflood.get_water_table()) &&
                         cbm::flt_greater( sc_.poro_sl[sl], wc_.wc_sl[sl]))
                    {
                        double const scale( cbm::flt_less( sc_.depth_sl[sl] - sc_.h_sl[sl], -m_eventflood.get_water_table()) ?
                                            (-m_eventflood.get_water_table() - (sc_.depth_sl[sl] - sc_.h_sl[sl])) / sc_.h_sl[sl] :
                                            1.0);
 
                        water_supply += scale * (sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl];
                    }
                }
 
                //remove rainfall (less irrigation water is supplied if it is raining)
                //and limit by saturated water flow into the first soil layer during the substep
                water_supply = cbm::bound( 0.0, water_supply - m_mc.precipitation_ts,
                                           WaterFlowSaturated( 0, get_impedance_factor( 0)));
 
                /* Surface water will be actively drained */
                if ( cbm::flt_greater_equal( wc_.surface_water, water_supply))
                {
                    wc_.accumulated_runoff += wc_.surface_water - water_supply;
                    wc_.surface_water = water_supply;
                    water_supply = 0.0;
                }
                else if ( cbm::flt_greater_zero( wc_.surface_water))
                {
                    water_supply -=  wc_.surface_water;
                }
 
                if ( m_eventflood.have_unlimited_water())
                {
                    wc_.accumulated_irrigation_automatic += water_supply;
                }
                else
                {
                    if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
                    {
                        wc_.irrigation_reservoir -= water_supply;
                    }
                    else
                    {
                        water_supply = wc_.irrigation_reservoir;
                        wc_.irrigation_reservoir = 0.0;
                    }
                    wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
                }
 
                m_wc.irrigation += water_supply;
            }
            else if ( cbm::flt_greater( m_eventflood.get_water_table(), wc_.surface_water))
            {
                //adjust surface and soil water in case of flooding event
                //water required to fill up surface water
                double water_supply( m_eventflood.get_water_table() - wc_.surface_water);
 
                //add water required to fill top x layers up to saturation
                unsigned long no_layers(3);
                if( no_layers > m_soillayers->soil_layer_cnt())
                {
                    no_layers = m_soillayers->soil_layer_cnt();
                }
 
                for ( size_t sl = 0;  sl < no_layers;  ++sl)
                {
                    water_supply += ((sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]);
                }
 
                //remove rainfall (less irrigation water is supplied if it if raining)
                water_supply = cbm::bound_min( 0.0, water_supply - m_mc.precipitation_ts);
 
                if ( m_eventflood.have_unlimited_water())
                {
                    wc_.accumulated_irrigation_automatic += water_supply;
                }
                else
                {
                    if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
                    {
                        wc_.irrigation_reservoir -= water_supply;
                    }
                    else
                    {
                        water_supply = wc_.irrigation_reservoir;
                        wc_.irrigation_reservoir = 0.0;
                    }
                    wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
                }
                m_wc.irrigation += water_supply;
            }
        }
    }

    double const irrigation_scheduled( cbm::bound_min( 0.0, wc_.accumulated_irrigation - m_wc.accumulated_irrigation_event_executed));
 
    // note: amounts are given per substep (substep_hours hours)
    if ( cbm::flt_greater( irrigation_scheduled, 24.0 * rainfall_intensity()))
    {
        m_wc.irrigation += irrigation_scheduled * substep_hours / 24.0;
        m_wc.accumulated_irrigation_event_executed += irrigation_scheduled * substep_hours / 24.0;
    }
    else if ( cbm::flt_greater( irrigation_scheduled, substep_hours * rainfall_intensity()))
    {
        m_wc.irrigation += substep_hours * rainfall_intensity();
        m_wc.accumulated_irrigation_event_executed += substep_hours * rainfall_intensity();
    }
    else
    {
        m_wc.irrigation += irrigation_scheduled;
        m_wc.accumulated_irrigation_event_executed += irrigation_scheduled;
    }
}

References CalcIrrigation(), get_impedance_factor(), and rainfall_intensity().

Referenced by CalcIrrigation(), and solve().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcPotEvapoTranspiration()

lerr_t ldndc::WatercycleDNDC::CalcPotEvapoTranspiration ( size_t _substep )

private

Potential evaporation.

there are two methods of calculating potential evaporation (choose via macro POTENTIAL_EVAPORATION above)

• [0] Procedure according to Thornthwaite 1948 with modifications by Camargo et al. 1999 and Peireira and Pruitt 2004 (implemented from Pereira and Pruitt 2004) to better account for dry environments and daylength.
  NOTE:
  • monthly temperatures are assumed to be equal to daily average temperature of the middle of the month, calculated from continous equations equal to those used in the weather generator.
  • the model uses always daily average temperature as input even if it is run in sub-daily mode
  • the original PnET-N-DNDC code which contains modification that could not be found in the literature has been outcommented.

• [1] uses the Priestley Taylor equation (Priestly and Taylor 1972) to calculate daily and hourly potential evapotranspiration (implemented from internet) vps calculation according to Allen et al. 1998, cit. in Cai et al. 2007

{
    // Daily evaportation demand is derived from the Thornthwaite, the Penman, or the Priestley-Taylor method (standard is 'thornthwaite).
    day_potevapo = 0.0;
    if ( m_cfg.potentialevapotranspiration_method == "thornthwaite")
    {
        if ( lclock()->is_position( TMODE_PRE_YEARLY))
        {
            thornthwaite_heat_index = ldndc::thornthwaite_heat_index( m_setup->latitude(),
                                                                      lclock()->days_in_year(),
                                                                      m_climate->annual_temperature_average(),
                                                                      m_climate->annual_temperature_amplitude());
        }
 
        double const daylength(ldndc::meteo::daylength(m_setup->latitude(), lclock()->yearday()));
        double const avg_dayl(12.0);
 
        day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
                            * ldndc::thornthwaite(
                                           mc_.nd_airtemperature,
                                           daylength, avg_dayl, lclock()->days_in_month(),
                                           thornthwaite_heat_index);
    }
    else if ( (m_cfg.potentialevapotranspiration_method == "penman") ||
              (m_cfg.potentialevapotranspiration_method == "priestleytaylor"))
    {
        cbm::sclock_t const &  clk( lclock_ref());
 
        /* vapour pressure [10^3:Pa] */
        double vps( 0.0);
        ldndc::meteo::vps( mc_.nd_airtemperature, &vps);
        double vp( vps - m_climate->vpd_day( clk));
 
        if( m_cfg.potentialevapotranspiration_method == "priestleytaylor")
        {
            day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
                                * ldndc::priestleytaylor( clk.yearday(),
                                                          clk.days_in_year(),
                                                          mc_.albedo,
                                                          m_setup->latitude(),
                                                          mc_.nd_airtemperature,
                                                          mc_.nd_shortwaveradiation_in,
                                                          vp,
                                                          m_sipar->PT_ALPHA());
        }
        else if( m_cfg.potentialevapotranspiration_method == "penman")
        {
            /* leaf area index */
            double const lai( m_veg->lai());
            ldndc::surface_type surface( cbm::flt_greater_zero( lai) ? short_grass : bare_soil);
 
            day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
                                * ldndc::penman( surface,
                                                 clk.yearday(),
                                                 clk.days_in_year(),
                                                 mc_.albedo,
                                                 m_setup->latitude(),
                                                 mc_.nd_shortwaveradiation_in * cbm::SEC_IN_DAY,
                                                 mc_.nd_airtemperature,
                                                 mc_.nd_windspeed,
                                                 vp);
        }
    }
    else
    {
        KLOGERROR("[BUG] huh : it seems that this potentialevapotranspiration method is not registered",
                  "tell the maintainers refering to this error message");
        return  LDNDC_ERR_FAIL;
    }
 
    // The daily evaporation demand is distributed throughout the day (standard is 'previouspattern)
    // TODO: design a routine that distributes evaporation demand according to temperature or humidity
    if ( m_cfg.evapotranspiration_method == "uniform")
    {
        // hourly evaporation demand is calculated from daily demand equally throughout the day
        hr_pot_evapotrans = day_potevapo * substep_hours / cbm::HR_IN_DAY;
    }
    else if ( m_cfg.evapotranspiration_method == "previouspattern")
    {
        // evaporation demand of the time step is derived from the weighted expression of previous values
        // and equally distributed over the substeps of the time step
        hr_pot_evapotrans = day_potevapo * pot_evapotranspiration_ts[lclock()->subday()-1] / substeps_per_timestep;
    }
    else if (m_cfg.evapotranspiration_method == "daytimeonly")
    {
        // hourly evaporation demand is distributed throughout the sunlight hours only (similar as radiation)
        // hour of the day at the end of the current substep
        double const hour( double(lclock()->subday() - 1) * substep_hours * substeps_per_timestep + double(_substep + 1) * substep_hours);
        double const day( lclock()->yearday());
        double const hsun( 0.833 * cbm::PI / 180.0);
        double const lat( m_setup->latitude() * cbm::PI / 180.0);
        double const dec( -23.44 * (cbm::PI / 180.0) * std::cos(2.0 * cbm::PI * (day - 172.0) / lclock()->days_in_year()));
        double const help( (std::sin(hsun) - std::sin(lat) * std::sin(dec)) / (std::cos(lat) * std::cos(dec)));
        double dayl(0.0);
        if (help < -0.999) dayl = 0.0;
        else if (help > 0.999) dayl = 24.0;
        else dayl = 24.0 - 24.0 / cbm::PI * std::acos((std::sin(hsun) - std::sin(lat) * std::sin(dec)) / (std::cos(lat) * std::cos(dec)));
        if (dayl < 1.0) dayl = 0.0;
        double const rise(12.0 - dayl * 0.5);
 
        double factor(0.0);
        if (dayl > 0.0) factor = (-std::cos(2.0 * cbm::PI * (hour - rise) / (cbm::HR_IN_DAY - 2.0 * rise)) + 1.0) / (cbm::HR_IN_DAY - 2.0 * rise);
 
        double fday(0.0);
        if (hour > rise && hour < (24.0 - rise)) fday = factor;
 
        // fday is given per hour
        hr_pot_evapotrans = day_potevapo * fday * substep_hours;
    }
    else
    {
        KLOGERROR("[BUG] huh : it seems that this evapotranspiration method option is not registered ",
            "tell the maintainers refering to this error message");
        return  LDNDC_ERR_FAIL;
    }
 
    // the evaporation limit for all timesteps is the daily demand
    wc_.accumulated_potentialevaporation += hr_pot_evapotrans;
 
    return  LDNDC_ERR_OK;
}

References CalcPotEvapoTranspiration(), ldndc::thornthwaite(), ldndc::thornthwaite_heat_index(), and thornthwaite_heat_index.

Referenced by CalcPotEvapoTranspiration(), and solve().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcSnowEvaporation()

void ldndc::WatercycleDNDC::CalcSnowEvaporation ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    //Average limit of snowEvaporation is 0.00013m day-1 (Granger and Male 1978)
    double hr_potESnow( cbm::bound(0.0, hr_pot_evapotrans, 0.00013 * substep_hours/24.0));
 
    if ( cbm::flt_greater_zero( wc_.surface_ice) &&
         cbm::flt_greater_zero( hr_potESnow))
    {
        //If there is a snowpack and there is a potential to evaporate
        if ( wc_.surface_ice > hr_potESnow)
        {
            //Partial evaporation from the snow
            //note: snow evaporation is booked as surface water evaporation (no separate state variable)
            hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - hr_potESnow);
            wc_.accumulated_surfacewaterevaporation += hr_potESnow;
            wc_.surface_ice -= hr_potESnow;
        }
        else
        {
            //Total evaporation of the snowpack
            hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - wc_.surface_ice);;
            wc_.accumulated_surfacewaterevaporation += wc_.surface_ice;
            wc_.surface_ice = 0.0;
        }
    }
}

References CalcSnowEvaporation().

Referenced by CalcSnowEvaporation(), CalcTranspiration(), and solve().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcSoilEvapoPercolation()

void ldndc::WatercycleDNDC::CalcSoilEvapoPercolation ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    // reset all water fluxes
    for (size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
    {
        m_wc.percolation_sl[sl] = 0.0;
    }
 
    // fill groundwater layer with water
    for (int sl = m_soillayers->soil_layer_cnt() - 1; sl >= 0; --sl)
    {
        if (cbm::flt_greater(sc_.depth_sl[sl], ground_water_table))
        {
            // ice volume is excepted from groundwater
            double const ice_vol(wc_.ice_sl[sl] / cbm::DICE);
            if (cbm::flt_greater(sc_.poro_sl[sl], ice_vol))
            {
                // new water in the layer is generated from groundwater 
                double const gw_access(cbm::bound_min(0.0, (sc_.poro_sl[sl] - ice_vol - wc_.wc_sl[sl]) * sc_.h_sl[sl] ));
                wc_.wc_sl[sl] += gw_access / sc_.h_sl[sl];
                wc_.accumulated_groundwater_access += gw_access;
            }
        }
        else
        {
            break;
        }
    }
 
    for (size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
    {
        // the first soil (litter) layer is filled with surface water up to its capacity
        if (sl == 0)
        {
            double delta_wc(0.0);
            double const wc_old(wc_.wc_sl[sl]);
            double const wc_cap( (sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]); // uptake capacity in first soil layer [mm]
            
            if (cbm::flt_greater(wc_.surface_water, wc_cap))
            {
                wc_.wc_sl[sl] = sc_.poro_sl[sl];
                delta_wc = cbm::bound_max(wc_.surface_water, (cbm::bound_min(0.0, wc_.wc_sl[sl] - wc_old) * sc_.h_sl[sl]));
                wc_.surface_water -= delta_wc;
            }
            else
            {
                wc_.wc_sl[sl] = (wc_.wc_sl[sl] * sc_.h_sl[sl] + wc_.surface_water) / sc_.h_sl[sl];
                delta_wc = cbm::bound_max(wc_.surface_water, (cbm::bound_min(0.0, wc_.wc_sl[sl] - wc_old) * sc_.h_sl[sl]));
                wc_.surface_water = 0.0;
            }
 
            wc_.accumulated_infiltration += delta_wc;
        }
        // other soil layers get updated by percolation rates from the layer above
        else
        {
            wc_.wc_sl[sl] += m_wc.percolation_sl[sl - 1] / sc_.h_sl[sl];
        }
 
        // percolation for all layers below groundwater 
        if (cbm::flt_greater(sc_.depth_sl[sl], ground_water_table))
        {
 
            // in case groundwater reaches the soil surface, the first outflux is estimated by saturated water flow
            if (sl == 0)
            {
                m_wc.percolation_sl[sl] = WaterFlowSaturated(sl, get_impedance_factor(sl));
            }
            // other soil layers are assumed to percolate by the same rate as the last non-groundwater layer
            else
            {
                m_wc.percolation_sl[sl] = m_wc.percolation_sl[sl - 1];
            }
 
        }
        // in all other cases percolation 
        else if (sl + 1 == m_soillayers->soil_layer_cnt())
        {
            m_wc.percolation_sl[sl] = 0.0;
            if (cbm::flt_greater_zero(wc_.sks_sl[sl]))
            {
                if (cbm::flt_greater(wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
                    cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_fc_sl[sl])) //security, should be always true if first condition is met
                {
                    //maximum allowed flow until field capacity
                    double const max_flow((wc_.wc_sl[sl] - wc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
                    m_wc.percolation_sl[sl] = cbm::bound_max(WaterFlowSaturated(sl, get_impedance_factor(sl)), max_flow);
                }
 
                //water flux above field capacity
                else if (cbm::flt_greater_equal(wc_.wc_sl[sl] + wc_.ice_sl[sl], wc_.wc_fc_sl[sl]))
                {
                    m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity(sl, get_impedance_factor(sl));
                }
 
                //water flux below field capacity and above wilting point
                else if ((sl > 0) && cbm::flt_greater_zero(m_wc.percolation_sl[sl - 1]) && 
                         cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_wp_sl[sl]))
                {
                    m_wc.percolation_sl[sl] = m_sipar->FPERCOL() * m_wc.percolation_sl[sl - 1];
                }
 
                //bound water flux by defined maximal percolation rate (e.g., during flooding events)
                m_wc.percolation_sl[sl] = cbm::bound_max(m_wc.percolation_sl[sl], max_percolation);
            }
        }
 
        //infiltrating water in current time step is exceeding available pore space of last time step
        //second condition: security, should be always true if first condition is met
        else if (cbm::flt_greater(wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
            cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_fc_sl[sl]))
        {
            //maximum allowed flow until field capacity
            double const max_flow((wc_.wc_sl[sl] - wc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
            m_wc.percolation_sl[sl] = cbm::bound_max(WaterFlowSaturated(sl, get_impedance_factor(sl)),
                max_flow);
        }
        else if (cbm::flt_greater_equal(wc_.wc_sl[sl] + wc_.ice_sl[sl], wc_.wc_fc_sl[sl]))
        {
            m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity(sl, get_impedance_factor(sl));
        }
        else if (cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_wp_sl[sl]))
        {
            m_wc.percolation_sl[sl] = WaterFlowCapillaryRise(sl, get_impedance_factor(sl));
        }
        else
        {
            m_wc.percolation_sl[sl] = 0.0;
        }
 
 
        // water content must be greater or equal wilting point
        if (cbm::flt_less((wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl] - m_wc.percolation_sl[sl], wc_.wc_wp_sl[sl] * sc_.h_sl[sl]) &&
            cbm::flt_greater_zero(m_wc.percolation_sl[sl]))
        {
            m_wc.percolation_sl[sl] = cbm::bound(0.0,
                (wc_.wc_sl[sl] + wc_.ice_sl[sl] - wc_.wc_wp_sl[sl]) * sc_.h_sl[sl],
                m_wc.percolation_sl[sl]);
        }
 
        // Update of the current soil layer water content
        wc_.wc_sl[sl] -= m_wc.percolation_sl[sl] / sc_.h_sl[sl];
 
        SoilWaterEvaporation(sl);
    }
 
    /* Correct water content and percolation if water content is above porosity
     * If the layer below cannot take up all of the water assigned for percolation,
     * so if the air filled pore space is not sufficient,
     * the water remains in the upper layer.
     */
    for (size_t sl = m_soillayers->soil_layer_cnt() - 1; sl > 0; sl--)
    {
        double const ice_vol(wc_.ice_sl[sl] / cbm::DICE);
        if (cbm::flt_greater(wc_.wc_sl[sl] + ice_vol, sc_.poro_sl[sl]))
        {
            double delta_wc(wc_.wc_sl[sl] + ice_vol - sc_.poro_sl[sl]);
            if (cbm::flt_greater_equal(ice_vol, sc_.poro_sl[sl]))
            {
                delta_wc = wc_.wc_sl[sl];
                wc_.wc_sl[sl] = 0.0;
            }
            else
            {
                wc_.wc_sl[sl] = sc_.poro_sl[sl] - ice_vol;
            }
            m_wc.percolation_sl[sl - 1] -= delta_wc * sc_.h_sl[sl];
            wc_.wc_sl[sl - 1] += delta_wc * sc_.h_sl[sl] / sc_.h_sl[sl - 1];
        }
    }
 
    double const ice_vol(wc_.ice_sl[0] / cbm::DICE);
    if (cbm::flt_greater(wc_.wc_sl[0] + ice_vol, sc_.poro_sl[0]))
    {
 
        double delta_wc(wc_.wc_sl[0] + ice_vol - sc_.poro_sl[0]);
        if (cbm::flt_greater_equal(ice_vol, sc_.poro_sl[0]))
        {
            delta_wc = wc_.wc_sl[0];
            wc_.wc_sl[0] = 0.0;
        }
        else
        {
            wc_.wc_sl[0] = sc_.poro_sl[0] - ice_vol;
        }
        wc_.accumulated_infiltration -= delta_wc * sc_.h_sl[0];
        wc_.surface_water += delta_wc * sc_.h_sl[0];
    }
 
    if (cbm::flt_greater_zero(wc_.surface_water))
    {
        CalcSurfaceFlux();
    }
 
    wc_.accumulated_waterflux_sl += m_wc.percolation_sl;
    wc_.accumulated_percolation += m_wc.percolation_sl[ m_soillayers->soil_layer_cnt()-1];
}

References CalcSoilEvapoPercolation(), CalcSurfaceFlux(), get_impedance_factor(), SoilWaterEvaporation(), WaterFlowAboveFieldCapacity(), and WaterFlowCapillaryRise().

Referenced by CalcSoilEvapoPercolation(), CalcTranspiration(), and solve().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcSurfaceFlux()

void ldndc::WatercycleDNDC::CalcSurfaceFlux ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    if ( cbm::flt_greater_zero( m_eventflood.get_bund_height()))
    {
        if ( cbm::flt_greater( wc_.surface_water, m_eventflood.get_bund_height()))
        {
            wc_.accumulated_runoff += wc_.surface_water - m_eventflood.get_bund_height();
            wc_.surface_water = m_eventflood.get_bund_height();
        }
    }
    else
    {
        /* runoff fraction per time step equals RO fraction per hour * length of a time step */
        double const runoff_fraction(m_sipar->FRUNOFF() * substep_day_fraction());
 
        // - not all surface water is removed within one timestep
        if ( cbm::flt_less( runoff_fraction, 1.0))
        {
            double const runoff( runoff_fraction * wc_.surface_water);
            wc_.surface_water -= runoff;
            wc_.accumulated_runoff += runoff;
        }
        // - all surface water is immediatly removed
        else
        {
            wc_.accumulated_runoff += wc_.surface_water;
            wc_.surface_water = 0.0;
        }
    }
}

References CalcSurfaceFlux(), and substep_day_fraction().

Referenced by CalcSoilEvapoPercolation(), CalcSurfaceFlux(), and CalcTranspiration().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    if ( cbm::flt_greater_zero( wc_.surface_water) &&
         cbm::flt_greater_zero( hr_pot_evapotrans))
    {
        if ( hr_pot_evapotrans > wc_.surface_water)
        {
            wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
            hr_pot_evapotrans -= wc_.surface_water;
            wc_.surface_water = 0.0;
        }
        else
        {
            wc_.surface_water -= hr_pot_evapotrans;
            wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
            hr_pot_evapotrans = 0.0;
        }
    }
}

References CalcSurfaceWaterEvaporation().

Referenced by CalcSurfaceWaterEvaporation(), CalcTranspiration(), and solve().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcTranspiration()

lerr_t ldndc::WatercycleDNDC::CalcTranspiration ( )

private

Parameters

[in]	None
[out]	None

Returns

None

References CalcSnowEvaporation(), CalcSoilEvapoPercolation(), CalcSurfaceFlux(), CalcSurfaceWaterEvaporation(), CalcTranspirationCouvreur(), get_clay_limitation(), get_fsl(), hours_per_time_step(), SoilWaterEvaporation(), substep_day_fraction(), WatercycleDNDC_preferential_flow(), WaterFlowAboveFieldCapacity(), and WaterFlowCapillaryRise().

Referenced by solve().

Here is the call graph for this function:

Here is the caller graph for this function:

CalcTranspirationCouvreur()

lerr_t ldndc::WatercycleDNDC::CalcTranspirationCouvreur ( double const & hr_potTransinm )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    // hourly water uptake
    double hr_uptake( 0.0);
 
    // Calculate the root water uptake and hence the transpiration for every plant individually
    for( PlantIterator vt = m_veg->begin(); vt != m_veg->end(); ++vt)
    {
        MoBiLE_Plant *p = *vt;
 
        // split up the potential transpiration equally to the various plants
        double const transpiration_possible_vt( _transpiration_possible * p->mFrt / m_veg->dw_frt() * cbm::CM_IN_M);  // in cm
 
        // If some soil layers are frozen and contain only ice, no fluid water is present and the water potential diverges.
        // The uptake from these layers is prevented by considering only roots in layers with more water than ice.
        // In this case the root distribution needs to be rescaled (otherwise the non-contributing layers would effectively contribute a potential = 0).
        double unfrozen_root_fraction( 1.0);
 
        /* TODO: include ice_water_ratio_critical as proper parameter */
        double const ice_water_ratio_critical( 0.1);
        double eff_soil_water_potential( 0.0);
        for( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
        {
            // fine root conditions to prevent transpiration without uptake organs, wateruptake from frozen ground
            if( cbm::flt_greater_zero( p->fFrt_sl[sl] * p->mFrt))
            {
                if( cbm::flt_greater_zero( wc_.wc_sl[sl]) &&
                    cbm::flt_less( wc_.ice_sl[sl] / wc_.wc_sl[sl], ice_water_ratio_critical))
                {
                    // overall water potential in this layer
                    eff_soil_water_potential += Soillayerwaterhead( sl) * p->fFrt_sl[sl];  // negative
                }
                else
                {
                    unfrozen_root_fraction -= p->fFrt_sl[sl];
                }
            }
            else
            {
                // roots are not wireless ;-)
                break;
            }
        }
        // normalize the root distribution in the potential, if some soil is frozen
        if( cbm::flt_greater_zero( unfrozen_root_fraction))
        {
            eff_soil_water_potential /= unfrozen_root_fraction;
        }
 
        // hydraulic conductivities of plant system
        double const k_rs( vt->KRSINIT() * p->mFrt / vt->FRTMASSINIT() * substep_hours);       // cm(water)/cm(from pressure)/timestep
        double const k_comp( vt->KCOMPINIT() * p->mFrt / vt->FRTMASSINIT() * substep_hours);   // cm(water)/cm(from pressure)/timestep
 
        // water potential in the leafs, depending on the potential transpiration
        double const leaf_water_potential( cbm::bound_min( vt->HLEAFCRIT(),
                                                           eff_soil_water_potential - (transpiration_possible_vt / k_rs)));    // cm
 
        if( cbm::flt_greater_zero( leaf_water_potential))
        {
            KLOGERROR( "Failed to determine transpiration: Non-negative leaf water potential!");
            return  LDNDC_ERR_FAIL;
        }
        if( !cbm::flt_greater_zero( eff_soil_water_potential - leaf_water_potential))
        {
            KLOGERROR( "Failed to determine transpiration!");
            return  LDNDC_ERR_FAIL;
        }
 
        double const effuptake( k_rs * (eff_soil_water_potential - leaf_water_potential));   // cm/timestep
        double uptake_tot( 0.0);
        double uptakecomp_tot( 0.0);
        //double uptakecompabs_tot( 0.0);
 
        for( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
        {
            // fine root condition to prevent transpiration without uptake organs
            if( cbm::flt_greater_zero( p->fFrt_sl[sl] * p->mFrt))
            {
                if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], ice_water_ratio_critical))
                {
                    double const soil_water_potential( Soillayerwaterhead( sl));
                    double const compuptake( k_comp * (soil_water_potential - eff_soil_water_potential));
                    double uptake_layer( (effuptake + compuptake) * p->fFrt_sl[sl]);    // absolute value in cm
 
                    wc_.accumulated_wateruptake_sl[sl] += uptake_layer * cbm::M_IN_CM;
                    uptake_tot += uptake_layer;                                         // absolute value in cm
                    uptakecomp_tot += compuptake * p->fFrt_sl[sl];                      // absolute value in cm
                    //uptakecompabs_tot += fabs(compuptake * p->fFrt_sl[sl]);           // absolute value in cm
 
                    wc_.wc_sl[sl] -= uptake_layer * cbm::M_IN_CM /sc_.h_sl[sl];         // fraction of water in this layer, in m
 
                    if( !cbm::flt_greater_zero( wc_.wc_sl[sl]))
                    {
                        KLOGERROR( "In soil layer ", sl, " the water content is negative: ", wc_.wc_sl[sl]);
                        return  LDNDC_ERR_FAIL;
                    }
                }
            }
            else
            {
                // roots are not wireless ;-)
                break;
            }
        }
 
        hr_uptake += uptake_tot * cbm::M_IN_CM;                 // absolute value in m for all species
 
        if( cbm::flt_greater_zero( fabs(uptakecomp_tot))){
            KLOGWARN( "Compensatory water uptake ", uptakecomp_tot, " > 0.");
        }
        if( cbm::flt_greater( uptake_tot, transpiration_possible_vt, 0.0000000001))
        {
            KLOGERROR( "Total water uptake ", uptake_tot, " larger than potential transpiration ", transpiration_possible_vt);
        }
        if(cbm::flt_equal_eps(unfrozen_root_fraction, 1.0, 0.0000000001) && !cbm::flt_equal_eps( uptake_tot, transpiration_possible_vt, 0.0000000001) && !cbm::flt_equal_eps( leaf_water_potential, vt->HLEAFCRIT(), 0.0000000001))
        {
            KLOGERROR( "T ", uptake_tot, "cm < Tpot ", transpiration_possible_vt, "cm even though leaf_water_potential ", leaf_water_potential, " > HLEAFCRIT", vt->HLEAFCRIT());
        }
    }
 
    // total evaporation demand reduced to apply on other evaporating sources (ground, snow, surface water)
    hr_pot_evapotrans = cbm::bound_min(0.0, hr_pot_evapotrans - hr_uptake);
 
    return  LDNDC_ERR_OK;
}

References CalcTranspirationCouvreur().

Referenced by CalcTranspiration(), and CalcTranspirationCouvreur().

Here is the call graph for this function:

Here is the caller graph for this function:

event_flood()

lerr_t ldndc::WatercycleDNDC::event_flood ( )

private

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

• surface water table
• bund height

{
    lerr_t  rc = m_eventflood.solve();
    if ( rc != LDNDC_ERR_OK)
    {
        KLOGERROR("Flooding event handler not successful!");
        return  rc;
    }
 
    /* max_percolation is gradually decreased after end of flooding event to inital value */
    double const maximum_percolation = m_eventflood.get_maximum_percolation();
    if ( cbm::is_valid( maximum_percolation))
    {
        if ( cbm::flt_greater_zero( maximum_percolation))
        {
            max_percolation = maximum_percolation * substep_hours / 24.0;
        }
    }
    else
    {
        /* time rate for gradual recovery of max_percolation (called once per time step, i.e., substeps_per_timestep substeps) */
        double const time_rate( substep_day_fraction() * double( substeps_per_timestep) * 0.1);
        max_percolation -= (max_percolation - WaterFlowSaturated( m_soillayers->soil_layer_cnt()-1, 1.0)) * time_rate;
    }
 
    return  LDNDC_ERR_OK;
}

References event_flood(), and substep_day_fraction().

Referenced by event_flood(), and solve().

Here is the call graph for this function:

Here is the caller graph for this function:

get_leaf_water()

double ldndc::WatercycleDNDC::get_leaf_water ( )

private

Collects leaf water from all canopy layers.

Parameters

[in]	None
[out]	None

Returns

Leaf water

{
    if ( m_veg->size() > 0u)
    {
        return wc_.wc_fl.sum();
    }
    return 0.0;
}

References get_leaf_water().

Referenced by balance_check(), CalcInterception(), and get_leaf_water().

Here is the call graph for this function:

Here is the caller graph for this function:

hours_per_time_step()

double ldndc::WatercycleDNDC::hours_per_time_step ( )

private

number of hours per time step

{
    return substep_hours * substeps_per_timestep;
}

References hours_per_time_step().

Referenced by CalcRaining(), CalcTranspiration(), and hours_per_time_step().

Here is the call graph for this function:

Here is the caller graph for this function:

rainfall_intensity()

double ldndc::WatercycleDNDC::rainfall_intensity ( )

private

rainfall intensity, i.e., maximum rainfall/irrigation per hour [m h-1]

{
    if ( rainfall_intensity_cache < 0.0)
    {
        if ( m_iokcomm->get_input_class< input_class_climate_t >())
        {
            rainfall_intensity_cache = m_iokcomm->get_input_class< climate::input_class_climate_t >()->station_info()->rainfall_intensity * cbm::M_IN_MM;
        }
        else
        {
            rainfall_intensity_cache = ldndc::climate::climate_info_defaults.rainfall_intensity * cbm::M_IN_MM;
        }
    }
    return rainfall_intensity_cache;
}

References rainfall_intensity(), and rainfall_intensity_cache.

Referenced by CalcIrrigation(), CalcRaining(), and rainfall_intensity().

Here is the call graph for this function:

Here is the caller graph for this function:

SoilWaterEvaporation()

void ldndc::WatercycleDNDC::SoilWaterEvaporation ( size_t _sl )

private

Parameters

[in]	_sl	Soil layer
[out]	None

Returns

None

{
    //percolation and direct evaporation from the (upper) soil
    //if there is snow cover, evaporation from soil will be set to 0, and only sublimation is allowed.
    double const wc_min_evaporation(m_sipar->WCDNDC_EVALIM_FRAC_WCMIN() * wc_.wc_wp_sl[_sl]);
    if ( cbm::flt_less( sc_.depth_sl[_sl] - sc_.h_sl[_sl], m_sipar->EVALIM()) &&
         cbm::flt_greater( wc_.wc_sl[_sl], wc_min_evaporation) &&
         cbm::flt_greater_zero( hr_pot_evapotrans))
    {
        //limiting factor for soil water infiltration
        double const f_clay( 1.0 + m_sipar->SLOPE_CLAYF() * std::min(1.0 - sc_.fcorg_sl[_sl] / cbm::CCORG, sc_.clay_sl[_sl]));
        double const f_water( cbm::bound( 0.0,
                                          pow( (wc_.wc_sl[_sl] - wc_min_evaporation) / wc_.wc_fc_sl[_sl], f_clay),
                                          1.0));
        double const h_rest( std::min( sc_.h_sl[_sl], m_sipar->EVALIM() - ( sc_.depth_sl[_sl] - sc_.h_sl[_sl])));
        double const f_depth( h_rest / m_sipar->EVALIM() * std::max(0.0,
                                       1.0 - std::min((double)m_sipar->EVALIM(),
                                                      sc_.depth_sl[_sl] - ( 0.5 * sc_.h_sl[_sl]))
                                       / m_sipar->EVALIM()));
        double const soilevaporation_max( cbm::bound_max( f_water * wc_.wc_sl[_sl] * sc_.h_sl[_sl], hr_pot_evapotrans));
 
        double const soilevaporation( cbm::bound( 0.0,
                                      hr_pot_evapotrans * f_depth * f_water,
                                      soilevaporation_max));
 
        //update of the current soil layer water content
        wc_.wc_sl[_sl] -= soilevaporation / sc_.h_sl[_sl];
        wc_.accumulated_soilevaporation += soilevaporation;
        hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - soilevaporation);
    }
}

References SoilWaterEvaporation().

Referenced by CalcSoilEvapoPercolation(), CalcTranspiration(), and SoilWaterEvaporation().

Here is the call graph for this function:

Here is the caller graph for this function:

substep_day_fraction()

double ldndc::WatercycleDNDC::substep_day_fraction ( )

private

Length of a substep as fraction of a day.

Length of a substep as fraction of a day (used to convert daily rates to rates per substep)

{
    return (substep_hours / 24.0);
}

References substep_day_fraction().

Referenced by CalcSurfaceFlux(), CalcTranspiration(), event_flood(), substep_day_fraction(), WatercycleDNDC_preferential_flow(), and WaterFlowAboveFieldCapacity().

Here is the call graph for this function:

Here is the caller graph for this function:

WatercycleDNDC_preferential_flow()

void ldndc::WatercycleDNDC::WatercycleDNDC_preferential_flow ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    //BY_PASSF is defined as the fraction of surface water that passes per hour (0.0 by default)
    if ( cbm::flt_greater_zero(m_sipar->BY_PASSF()) &&
         cbm::flt_greater_zero( wc_.surface_water))
    {
 
        double water_def_total( 0.0);
        for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
        {
            //water flux that flow through macro pores can maximally store
            water_def_total += cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]);
        }
 
        //duration of timestep (hr-1)
        double const fhr( substep_day_fraction());
 
 
        //fraction of surface water that is put into bypass flow (m timestep-1)
        double const bypass_fraction( cbm::bound_max(m_sipar->BY_PASSF() * fhr, 0.99));
        double const bypass_water( wc_.surface_water * bypass_fraction);
        wc_.surface_water -= bypass_water;
        wc_.accumulated_infiltration += bypass_water;
 
        if ( cbm::flt_greater_equal( water_def_total, bypass_water))
        {
 
            for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
            {
                double const water_def( cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
                double const water_add( water_def / water_def_total * bypass_water);
 
                wc_.wc_sl[sl] += water_add / sc_.h_sl[sl];
            }
        }
        else
        {
            for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
            {
                double const water_def( cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
                wc_.wc_sl[sl] += water_def / sc_.h_sl[sl];
            }
            /* add surplus to percolation out of last soil layer */
            wc_.accumulated_percolation += bypass_water - water_def_total;
            wc_.accumulated_waterflux_sl[m_soillayers->soil_layer_cnt()-1] += bypass_water - water_def_total;
        }
    }
}

References substep_day_fraction(), and WatercycleDNDC_preferential_flow().

Referenced by CalcTranspiration(), solve(), and WatercycleDNDC_preferential_flow().

Here is the call graph for this function:

Here is the caller graph for this function:

WaterFlowAboveFieldCapacity()

double ldndc::WatercycleDNDC::WaterFlowAboveFieldCapacity ( size_t _sl, double const & _fact_impedance )

private

factor accounting for different water travelling characteristics in litter and mineral soil

{
    double const sks( WaterFlowSaturated( _sl, _fact_impedance));
    double const fsl( get_fsl( _sl));
    double const fhr( substep_day_fraction());

    double slope(( _sl < m_soillayers->soil_layers_in_litter_cnt()) ? m_sipar->SLOPE_FF() : m_sipar->SLOPE_MS());
 
    double const travel_time( 1.0 - log10( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * fhr));
    double const travel_time_factor( cbm::flt_greater_zero( travel_time) ? (1.0 - exp(-1.0 / travel_time)) : 1.0);
    return cbm::bound_max( cbm::sqr(1.0 - (wc_.wc_fc_sl[_sl] / ((wc_.wc_sl[_sl] + wc_.ice_sl[_sl]) * slope)))
                           * wc_.wc_sl[_sl] * sc_.h_sl[_sl] * fsl * travel_time_factor * _fact_impedance,
                           sks);
}

References get_fsl(), substep_day_fraction(), and WaterFlowAboveFieldCapacity().

Referenced by CalcSoilEvapoPercolation(), CalcTranspiration(), and WaterFlowAboveFieldCapacity().

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

rainfall_intensity_cache

double ldndc::WatercycleDNDC::rainfall_intensity_cache

private

rainfall intensity (constant during a run), cached on first use (negative: not yet read) [m h-1]

Referenced by rainfall_intensity().