Skip to content

fr_crop_model module

simpest.models.fr_crop_model

Daily crop growth step and disease damage mechanisms.

This module computes one day of crop growth together with the four disease damage mechanisms that couple the epidemiological state to crop performance:

  • light stealers reduce intercepted radiation,
  • RUE reducers lower radiation use efficiency,
  • assimilate sappers drain fixed carbon, and
  • senescence accelerators shorten the green-canopy duration.

Crop growth is computed through one of two branches:

  • Internal growth model. When no external crop-model series is supplied, the canopy, biomass, and yield are simulated from thermal time using logistic light-interception curves and a radiation-use-efficiency biomass model.
  • External growth model. When a daily crop-model series is supplied, the attainable light interception, biomass, and yield are taken from that series, and the damage mechanisms are applied to the daily increments to obtain the actual (disease-limited) trajectories.

In the external branch, day_after_sowing is incremented each simulated day so that the runner's maturity/safety stop and the calibration objective's "is-planted" logic operate correctly, and growing_degree_days is carried from the external series. These fields therefore reflect the simulated crop calendar rather than being left unset.

run(input_, parameters, output, output1)

Compute one daily crop growth step.

The four disease damage mechanisms are derived first from the current severity, then crop growth is advanced through either the internal growth model (logistic light interception plus a radiation-use-efficiency biomass model) or the external crop-model branch, depending on whether a daily crop-model series is attached to input_. In both branches the attainable (disease-free) and actual (disease-limited) light interception, biomass, and yield are written to output1.

Parameters:

Name Type Description Default
input_ InputsDaily

Today's daily inputs, including the optional external crop-model series.

required
parameters Parameters

Model parameters, including the crop and disease groups.

required
output Outputs

Previous day's output state, used as the starting state.

required
output1 Outputs

Today's output state, updated in place.

required

Returns:

Type Description
None

This function updates output1 in place.

Source code in simpest/models/fr_crop_model.py
def run(input_: InputsDaily, parameters: Parameters,
        output: Outputs, output1: Outputs) -> None:
    """Compute one daily crop growth step.

    The four disease damage mechanisms are derived first from the current
    severity, then crop growth is advanced through either the internal growth
    model (logistic light interception plus a radiation-use-efficiency biomass
    model) or the external crop-model branch, depending on whether a daily
    crop-model series is attached to ``input_``. In both branches the attainable
    (disease-free) and actual (disease-limited) light interception, biomass, and
    yield are written to ``output1``.

    Args:
        input_ (InputsDaily): Today's daily inputs, including the optional
            external crop-model series.
        parameters (Parameters): Model parameters, including the crop and disease
            groups.
        output (Outputs): Previous day's output state, used as the starting state.
        output1 (Outputs): Today's output state, updated in place.

    Returns:
        None: This function updates ``output1`` in place.
    """
    pc = parameters.par_crop
    pd = parameters.par_disease

    # -----------------------------------------------------------------------
    # Damage mechanisms – always computed before selecting model branch
    # -----------------------------------------------------------------------
    severity = output.disease.disease_severity

    output1.disease.damage_mechanisms.light_stealers = (
        1.0 - _light_stealers_fn(severity, pd.light_stealer_damage)
    )
    output1.disease.damage_mechanisms.rue_reducers = (
        1.0 - _rue_reduction_fn(severity, pd.rue_reducer_damage)
    )
    output1.disease.damage_mechanisms.assimilate_sappers = (
        _assimilate_sappers_fn(severity, pd.assimilate_sappers_damage)
    )
    output1.disease.damage_mechanisms.senescence_accelerators = (
        _senescence_accelerator_fn(severity, pd.senescence_accelerator_damage)
    )

    dm = output1.disease.damage_mechanisms
    cmd = input_.crop_model_data

    # -----------------------------------------------------------------------
    # Branch A: internal crop growth model (no external crop-model series)
    # -----------------------------------------------------------------------
    if cmd is None or len(cmd.f_int) == 0:
        t_ave = (input_.tmin + input_.tmax) / 2.0
        t_func = t_response(t_ave, pc.tbase_crop, pc.topt_crop, pc.tmax_crop)

        # Accumulate GDD
        output1.crop.growing_degree_days = (
            output.crop.growing_degree_days
            + t_func * (pc.topt_crop - pc.tbase_crop)
        )

        if output.crop.cycle_completion_percentage <= 100.0:
            gdd = output1.crop.growing_degree_days

            # Phenological code (1 = vegetative, 2 = reproductive)
            output1.crop.pheno_code = _pheno_code(
                gdd, pc.flowering_start / 100.0 * pc.cycle_length
            )

            # Attainable light interception
            f_int_att, senescence_started = _f_int_compute(
                pc.cycle_length, pc.slope_growth, pc.half_int_growth,
                pc.slope_senescence, pc.half_int_senescence, gdd
            )
            output1.crop.light_interception_attainable = f_int_att
            output1.crop.senescence_started = senescence_started

            # Save peak fInt at the onset of senescence
            if senescence_started and output1.crop.f_int_peak == 0.0:
                output1.crop.f_int_peak = output.crop.light_interception_attainable

            # Shift senescence half-point due to senescence accelerators
            half_int_sen_shifted = (
                pc.half_int_senescence - dm.senescence_accelerators * 100.0
            )

            # Actual light interception (with disease pressure)
            f_int_act, _ = _f_int_compute(
                pc.cycle_length, pc.slope_growth, pc.half_int_growth,
                pc.slope_senescence, half_int_sen_shifted, gdd
            )
            output1.crop.light_interception_actual = f_int_act

            # Apply light stealers
            output1.crop.light_interception_actual -= (
                output1.crop.light_interception_actual * dm.light_stealers
            )
            if output1.crop.light_interception_actual < 0.0:
                output1.crop.light_interception_actual = 0.0

            # Potential biomass accumulation
            carbon_rate_pot = _carbon_rate(
                pc.radiation_use_efficiency, input_.rad, t_func,
                output1.crop.light_interception_attainable
            ) * 10.0  # g m⁻² → kg ha⁻¹
            output1.crop.agb_attainable = output.crop.agb_attainable + carbon_rate_pot

            # Actual biomass accumulation
            carbon_rate_act = _carbon_rate(
                pc.radiation_use_efficiency
                - pc.radiation_use_efficiency * dm.rue_reducers,
                input_.rad, t_func,
                output1.crop.light_interception_actual
            ) * 10.0
            carbon_rate_act -= dm.assimilate_sappers
            output1.crop.agb_actual = output.crop.agb_actual + carbon_rate_act
            if output1.crop.agb_actual < 0.0:
                output1.crop.agb_actual = 0.0

            # Potential yield
            output1.crop.yield_attainable = output.crop.yield_attainable + _yield_rate(
                output1.crop.pheno_code, carbon_rate_pot, pc.partitioning_maximum
            )

            # Actual yield
            output1.crop.yield_actual = output.crop.yield_actual + _yield_rate(
                output1.crop.pheno_code, carbon_rate_act, pc.partitioning_maximum
            )
            if output1.crop.yield_actual < 0.0:
                output1.crop.yield_actual = 0.0

            # Days after sowing
            output1.crop.day_after_sowing = output.crop.day_after_sowing + 1

            # Cycle completion (clamped at 100 %)
            output1.crop.cycle_completion_percentage = min(
                gdd / pc.cycle_length * 100.0, 100.0
            )
        else:
            # Crop cycle complete → reset
            output1.crop = CropOutputs()
            output.crop = CropOutputs()

    # -----------------------------------------------------------------------
    # Branch B: External crop model supplies daily f_int, AGB, yield, GDD
    # -----------------------------------------------------------------------
    else:
        # Normalise date key (InputsDaily.date is datetime; dict keys are date)
        today: date = (
            input_.date.date() if isinstance(input_.date, datetime) else input_.date
        )

        if today in cmd.f_int:
            # --- Attainable values from external model ---
            output1.crop.light_interception_attainable = cmd.f_int[today]
            output1.crop.agb_attainable = cmd.agb[today]
            output1.crop.yield_attainable = cmd.yield_.get(today, 0.0)

            # --- Actual light interception ---
            f_int_att = output1.crop.light_interception_attainable
            f_int_act = (f_int_att - f_int_att * dm.light_stealers - dm.senescence_accelerators)
            if f_int_act < 0.0:
                f_int_act = 0.0
            output1.crop.light_interception_actual = f_int_act

            # Senescence flag: senescence starts when yield appears
            if cmd.yield_.get(today, 0.0) > 0.0:
                output1.crop.senescence_started = True

            # Peak fInt across entire season
            output1.crop.f_int_peak = max(cmd.f_int.values())

            # --- Daily rates (attainable) ---
            prev_day = today - timedelta(days=1)
            if prev_day in cmd.agb:
                pot_agb_rate = cmd.agb[today] - cmd.agb[prev_day]
                pot_yield_rate = (cmd.yield_.get(today, 0.0) - cmd.yield_.get(prev_day, 0.0))
            else:  # first day of season
                pot_agb_rate = cmd.agb[today]
                pot_yield_rate = cmd.yield_.get(today, 0.0)

            # --- Actual rates ---
            if f_int_att > 0.0:
                damage_ratio = (f_int_att - f_int_act) / f_int_att
                act_agb_rate = pot_agb_rate - pot_agb_rate * damage_ratio
                act_yield_rate = pot_yield_rate - pot_yield_rate * damage_ratio
            else:
                # First day, no light interception yet
                act_agb_rate = output.crop.agb_actual
                act_yield_rate = output.crop.yield_actual

            # Clamp negative rates before applying further reductions
            if act_agb_rate < 0.0:
                act_agb_rate = 0.0
            if act_yield_rate < 0.0:
                act_yield_rate = 0.0

            # Apply RUE reducers and assimilate sappers
            act_agb_rate = (
                act_agb_rate
                - act_agb_rate * dm.rue_reducers
                - dm.assimilate_sappers
            )
            act_yield_rate = (
                act_yield_rate
                - act_yield_rate * dm.rue_reducers
                - dm.assimilate_sappers
            )

            # Update state variables
            output1.crop.agb_actual = output.crop.agb_actual + act_agb_rate
            output1.crop.yield_actual = output.crop.yield_actual + act_yield_rate
            if output1.crop.agb_actual < 0.0:
                output1.crop.agb_actual = 0.0
            if output1.crop.yield_actual < 0.0:
                output1.crop.yield_actual = 0.0

            # Cycle completion from external model
            output1.crop.cycle_completion_percentage = cmd.cycle_percentage.get(today, 0.0)

            # Track the crop calendar: increment days-after-sowing so the
            # maturity/safety stop and the calibration objective's is-planted
            # logic operate correctly (see the module docstring).
            output1.crop.day_after_sowing = output.crop.day_after_sowing + 1

            # Carry thermal time from the external series. This is reporting-only
            # when use_gdd=False and drives cycle completion when use_gdd=True.
            output1.crop.growing_degree_days = cmd.gdd.get(today, 0.0)

        else:
            # Date not covered by external model → crop harvested, reset
            output1.crop = CropOutputs()
            output.crop = CropOutputs()