Select Git revision
Marco Matthies authored
The code seems to work for crop plants but currently doesn't yet work for non-crop plants such as "permanent set-aside", which grow to enourmous height and do not get reset at the end of the year. I think this is due to different treatment of the growth phase transitions for these plants.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
almass.jl 33.71 KiB
### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
###
### This file is responsible for managing the crop growth modules.
###
#FIXME Looking at fields.pdf in the Persefone output, crop maturity is not recognised
# properly, so most crops (except maize) are never planted. At the same time, maize never
# seems to start growing.
#TODO write tests to compare our output to the output of the original ALMaSS algorithm
# Some functions translated from ALMaSS for use in Persefone.jl
#
# ALMaSS is licensed as follows:
#
# Copyright (c) 2011, Christopher John Topping, Aarhus University
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are permitted provided
# that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list of conditions and the
# following disclaimer.
# Redistributions in binary form must reproduce the above copyright notice, this list of conditions and
# the following disclaimer in the documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
# FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
# BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
module ALMaSS
## The two input files required by ALMaSS. Must be present in the crop directory.
const CROPFILE = "crop_data_general.csv"
const GROWTHFILE = "crop_growth_curves.csv"
using Persefone:
@rand,
@u_str,
AnnualDate,
Management,
cm,
SimulationModel,
fertiliser,
maxtemp,
mintemp,
year,
dayofyear
import Persefone:
stepagent!,
croptype,
cropname,
cropheight,
cropcover,
cropyield,
sow!,
harvest!,
isharvestable,
bounds
using Dates: Date, month, monthday
using CSV: CSV
using Unitful: unit
# We can't use Length{Float64} as that is a Union type
const Tlength = typeof(1.0cm)
"""
GrowthPhase
ALMaSS crop growth curves are split into five phases, triggered by
seasonal dates or agricultural events.
"""
@enum GrowthPhase janfirst sow marchfirst harvest1 harvest2
"""
CropCurveParams
The values in this struct define one crop growth curve.
"""
struct CropCurveParams
#TODO add Unitful
curveID::Int
highnutrients::Bool
GDD::Dict{GrowthPhase,Vector{Float64}}
LAItotal::Dict{GrowthPhase, Vector{Float64}}
LAIgreen::Dict{GrowthPhase, Vector{Float64}}
height::Dict{GrowthPhase, Vector{typeof(1.0cm)}}
end
"""
CropType
The type struct for all crops. Currently follows the crop growth model as
implemented in ALMaSS.
"""
struct CropType
#FIXME this needs thinking about. The sowing and harvest dates belong in the farm model,
# not here. Also, we need to harmonise crops across the crop growth models.
name::String
group::String
is_c4_plant::Bool # false means it is a C3 plant
minsowdate::Union{Missing,AnnualDate}
maxsowdate::Union{Missing,AnnualDate}
minharvestdate::Union{Missing,AnnualDate}
maxharvestdate::Union{Missing,AnnualDate}
mingrowthtemp::Union{Missing,Float64}
highnutrientgrowth::Union{Missing,CropCurveParams}
lownutrientgrowth::Union{Missing,CropCurveParams}
end
cropname(ct::CropType) = ct.name
"""
CropState
The state data for an ALMaSS vegetation point calculation. Usually
part of a `FarmPlot`.
"""
@kwdef mutable struct CropState
# Class in original ALMaSS code:
# `VegElement` from `Landscape/Elements.h`, line 601
croptype::CropType
events::Vector{Management} = Management[]
mature::Bool = false #TODO how do we determine this?
# biomass::Float64 #XXX I need to figure out how to calculate this
# height::Length{Float64}
# growingdegreedays::Float64
phase::GrowthPhase
vegddegs::Float64 = 0 # Vegetation growing degree days since sowing
ddegs::Float64 = 0 # Growing degree days in current phase
yddegs::Float64 = 0 # Growing degree days in current phase for previous day LAgreen::Float64 = 0
LAtotal::Float64 = 0
oldLAtotal::Float64 = 0
veg_height::Tlength = 0.0cm
newgrowth::Float64 = 0
veg_cover::Float64 = 0
veg_biomass::Float64 = 0
green_biomass::Float64 = 0
dead_biomass::Float64 = 0
growth_scaler::Float64 = 1.0 # TODO: where is this set?
forced_phase_shift::Bool = false # TODO: what does this do?
# TODO: meaning of force_* fields ?
force_growth::Bool = false
force_LAgreen::Float64 = 0
force_LAtotal::Float64 = 0
force_veg_height::Tlength = 0.0cm
end
croptype(cs::CropState) = cs.croptype
cropname(cs::CropState) = cropname(croptype(cs))
cropheight(cs::CropState) = cs.veg_height
cropcover(cs::CropState) = cs.veg_cover
cropyield(cs::CropState) = cs.veg_biomass # TODO: correct? units? dry or wet?
isharvestable(cs::CropState) = cs.mature
# Constant in original ALMaSS code:
# `EL_VEG_START_LAIT` from `Landscape/Elements.cpp`, line 238-239
const VEG_START_LAIT::Float64 = 1.08
# Constant in original ALMaSS code:
# `EL_VEG_HEIGHTSCALE` from `Landscape/Elements.cpp`, line 242-243
const VEG_HEIGHTSCALE::Tlength = 16.0cm
# Constant in original ALMaSS code:
# `cfg_beer_law_extinction_coef` from `Landscape/Elements.cpp`, line 139
# TODO: units
const BEER_LAW_EXTINCTION_COEF::Float64 = 0.6
# TODO: units
"Temperature to solar conversion factor for C3 plants."
const temperature_to_solar_conversion_c3 = [
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # -30°C to -21°C
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # -20°C to -11°C
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # -10°C to -1°C
0, 0, 0, 0, 0, 0.28, 0.56, 0.84, 1.12, 1.4, # 0°C to 9°C
1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, # 10°C to 19°C
1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.26, 1.12, 0.98, 0.84, # 20°C to 29°C
0.7, 0.56, 0.42, 0.28, 0.14, 0, 0, 0, 0, 0, # 30°C to 39°C
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # 40°C to 49°C
0 # 50°C
]
# TODO: units
"Temperature to solar conversion factor for C4 plants."
const temperature_to_solar_conversion_c4 = [
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # -30°C to -21°C
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # -20°C to -11°C
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # -10°C to -1°C
0, 0, 0, 0, 0, 0, 0, 0, 0.242857, 0.485714, # 0°C to 9°C
0.728571, 0.971429, 1.214286, 1.457143, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, # 10°C to 19°C
1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, # 20°C to 29°C
1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.53, 1.36, 1.19, 1.02, # 30°C to 39°C
0.85, 0.68, 0.51, 0.34, 0.17, 0, 0, 0, 0, 0, # 40°C to 49°C
0 # 50°C
]
"""
solar_conversion_c3(temperature)
Solar conversion factor (no units) for C3 plants.
"""
function solar_conversion_c3(temperature)
if temperature < -30 || temperature > 50
return zero(eltype(temperature_to_solar_conversion_c3))
end
idx = Int(floor(0.5 + temperature)) + 30 + 1
return temperature_to_solar_conversion_c3[idx]
end
"""
solar_conversion_c3(temperature)
Solar conversion factor (no units) for C4 plants.
"""
function solar_conversion_c4(temperature)
if temperature < -30 || temperature > 50
return zero(eltype(temperature_to_solar_conversion_c4))
end
idx = Int(floor(0.5 + temperature)) + 30 + 1
return temperature_to_solar_conversion_c4[idx]
end
# Below is the data for insolation in MJ per m2 for Ødum fitted to 1961-1990
# by Olesen (1991)
# TODO: unit is u"MJ/m^2"
const global_insolation = [
0.91, 0.921585, 0.935138, 0.950669, 0.968188,
0.987706, 1.009231, 1.032774, 1.058346, 1.085956,
1.115614, 1.14733, 1.181113, 1.216974, 1.254921,
1.294964, 1.337112, 1.381373, 1.427757, 1.476271,
1.526923, 1.57972, 1.634671, 1.691781, 1.751057,
1.812504, 1.876128, 1.941934, 2.009926, 2.080107,
2.15248, 2.227048, 2.303812, 2.382774, 2.463933,
2.547289, 2.63284, 2.720585, 2.81052, 2.902642,
2.996945, 3.093425, 3.192073, 3.292884, 3.395847,
3.500954, 3.608194, 3.717555, 3.829024, 3.942587,
4.05823, 4.175935, 4.295687, 4.417465, 4.541252,
4.667024, 4.794762, 4.924441, 5.056036, 5.189522,
5.324873, 5.462059, 5.601051, 5.741819, 5.884329,
6.02855, 6.174446, 6.321981, 6.471118, 6.621819,
6.774044, 6.927753, 7.082902, 7.239449, 7.397349,
7.556555, 7.717022, 7.8787, 8.041541, 8.205493,
8.370505, 8.536524, 8.703496, 8.871365, 9.040077,
9.209573, 9.379796, 9.550686, 9.722183, 9.894227,
10.06675, 10.2397, 10.41301, 10.58661, 10.76044,
10.93443, 11.10852, 11.28264, 11.45671, 11.63068,
11.80448, 11.97802, 12.15125, 12.3241, 12.49648,
12.66834, 12.8396, 13.01019, 13.18004, 13.34907,
13.51721, 13.6844, 13.85056, 14.01561, 14.17949,
14.34212, 14.50344, 14.66337, 14.82184, 14.97877,
15.13411, 15.28778, 15.43971, 15.58983, 15.73807,
15.88437, 16.02866, 16.17087, 16.31093, 16.44879,
16.58439, 16.71764, 16.8485, 16.97691, 17.1028,
17.22612, 17.3468, 17.46479, 17.58005, 17.6925,
17.80211, 17.90881, 18.01256, 18.11332, 18.21102,
18.30564, 18.39712, 18.48542, 18.5705, 18.65233,
18.73086, 18.80605, 18.87788, 18.94632, 19.01132,
19.07286, 19.13092, 19.18546, 19.23647, 19.28393,
19.3278, 19.36809, 19.40475, 19.4378, 19.4672,
19.49296, 19.51506, 19.53349, 19.54825, 19.55934,
19.56676, 19.5705, 19.57058, 19.56698, 19.55973,
19.54883, 19.53428, 19.51611, 19.49432, 19.46893,
19.43996, 19.40743, 19.37136, 19.33177, 19.28868,
19.24213, 19.19214, 19.13873, 19.08195, 19.02183,
18.95839, 18.89168, 18.82173, 18.74858, 18.67228,
18.59286, 18.51036, 18.42484, 18.33634, 18.2449,
18.15058, 18.05342, 17.95348, 17.8508, 17.74545,
17.63746, 17.52691, 17.41384, 17.29832, 17.18039, 17.06012, 16.93757, 16.8128, 16.68586, 16.55683,
16.42576, 16.29271, 16.15775, 16.02094, 15.88235,
15.74203, 15.60006, 15.4565, 15.3114, 15.16485,
15.0169, 14.86762, 14.71707, 14.56532, 14.41243,
14.25848, 14.10352, 13.94762, 13.79084, 13.63325,
13.47491, 13.31588, 13.15624, 12.99604, 12.83534,
12.6742, 12.5127, 12.35088, 12.18881, 12.02654,
11.86415, 11.70167, 11.53918, 11.37673, 11.21437,
11.05216, 10.89016, 10.72841, 10.56697, 10.40589,
10.24522, 10.08502, 9.925323, 9.766186, 9.607653,
9.449772, 9.292586, 9.13614, 8.980476, 8.825635,
8.67166, 8.518588, 8.366459, 8.215311, 8.06518,
7.916101, 7.768108, 7.621236, 7.475515, 7.330978,
7.187655, 7.045573, 6.904763, 6.765249, 6.62706,
6.490218, 6.35475, 6.220676, 6.088021, 5.956804,
5.827045, 5.698765, 5.571981, 5.44671, 5.32297,
5.200776, 5.080142, 4.961083, 4.843613, 4.727743,
4.613485, 4.50085, 4.38985, 4.280492, 4.172788,
4.066743, 3.962368, 3.859668, 3.758651, 3.659322,
3.561687, 3.465753, 3.371522, 3.279, 3.18819,
3.099097, 3.011723, 2.926071, 2.842145, 2.759946,
2.679477, 2.600739, 2.523735, 2.448466, 2.374933,
2.303139, 2.233084, 2.16477, 2.098198, 2.033369,
1.970285, 1.908946, 1.849355, 1.791512, 1.735419,
1.681077, 1.628489, 1.577656, 1.528581, 1.481264,
1.43571, 1.391919, 1.349896, 1.309643, 1.271164,
1.234461, 1.199539, 1.166401, 1.135053, 1.105497,
1.07774, 1.051786, 1.02764, 1.005309, 0.984798,
0.966113, 0.94926, 0.934248, 0.921081, 0.909769,
0.900317, 0.892735, 0.88703, 0.88321, 0.881284,
0.881261, 0.883149, 0.886958, 0.892696, 0.900374,
0.900374, # duplicated last datapoint in case the year has 366 days
]
"""
Base.tryparse(type, str)
Extend `tryparse` to allow parsing GrowthPhase values.
(Needed to read in the CSV parameter file.)
"""
function Base.tryparse(::Type{GrowthPhase}, str::String)
str == "janfirst" ? janfirst :
str == "sow" ? sow :
str == "marchfirst" ? marchfirst :
str == "harvest1" ? harvest1 :
str == "harvest2" ? harvest2 :
nothing
end
"""
buildgrowthcurve(data)
Convert a list of rows from the crop growth data into a CropCurveParams object.
"""
function buildgrowthcurve(data::Vector{CSV.Row})
isempty(data) && return missing
GDD = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
harvest2=>Vector{Float64}())
LAItotal = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
harvest2=>Vector{Float64}())
LAIgreen = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
harvest2=>Vector{Float64}())
height = Dict(janfirst=>Vector{Tlength}(), sow=>Vector{Tlength}(),
marchfirst=>Vector{Tlength}(), harvest1=>Vector{Tlength}(),
harvest2=>Vector{Tlength}())
for e in data
append!(GDD[e.growth_phase], e.GDD) append!(LAItotal[e.growth_phase], e.LAI_total)
append!(LAIgreen[e.growth_phase], e.LAI_green)
append!(height[e.growth_phase], e.height * cm) # assumes `height` is in units of `cm`
end
CropCurveParams(data[1].curve_id, data[1].nutrient_status=="high",
GDD, LAItotal, LAIgreen, height)
end
"""
readcropparameters(cropdirectory)
Parse a CSV file containing the required parameter values for each crop
(as produced from the original ALMaSS file by `convert_almass_data.py`).
"""
function readcropparameters(cropdirectory::String)
@debug "Reading crop parameters"
cropdata = CSV.File(joinpath(cropdirectory, CROPFILE), missingstring="NA",
types=[String,AnnualDate,AnnualDate,AnnualDate,AnnualDate,Float64,String])
growthdata = CSV.File(joinpath(cropdirectory, GROWTHFILE), missingstring="NA",
types=[Int,String,String,GrowthPhase,String,
Float64,Float64,Float64,Float64])
croptypes = Dict{String,CropType}()
for crop in cropdata
cropgrowthdata = growthdata |> filter(x -> !ismissing(x.crop_name) &&
x.crop_name == crop.name)
highnuts = buildgrowthcurve(cropgrowthdata |>
filter(x -> x.nutrient_status=="high"))
lownuts = buildgrowthcurve(cropgrowthdata |>
filter(x -> x.nutrient_status=="low"))
# TODO: set crop group temporarily until there is a column in
# the csv file
crop_group = "CROP_GROUP_NOT_SET"
# TODO: it would be better to save this in the parameter file
# (Note that this matches the current ALMaSS code though,
# which also only hardcodes maize as a C4 crop)
is_c4_plant = occursin("maize", lowercase(crop.name))
croptypes[crop.name] = CropType(crop.name, crop_group, is_c4_plant, crop.minsowdate,
crop.maxsowdate, crop.minharvestdate, crop.maxharvestdate,
crop.mingrowthtemp, highnuts, lownuts)
end
croptypes
end
"""
setphase!(cropstate, phase)
Set the growth `phase` of an ALMaSS `cropstate`.
"""
function setphase!(cs::CropState, phase::GrowthPhase, model::SimulationModel)
# Function in original ALMaSS code:
# `void VegElement::SetGrowthPhase(int a_phase)` from `Landscape/Elements.cpp`, line 2098
if phase == sow
cs.vegddegs = 0
elseif phase == harvest1
cs.vegddegs = -1
end
if phase == janfirst
cs.forced_phase_shift = false
# TODO: Original comments from ALMaSS below, but do they make sense?
#
# If it is the first growth phase of the year then we might
# cause some unnecessary hops if e.g. our biomass is 0 and we
# suddenly jump up to 20 cm. To stop this from happening we
# check here and if our settings are lower than the targets we
# do nothing.
if isvalidstart(cs, phase)
if startvalue_veg_height(cs, phase) < cs.veg_height
# we are better off with the numbers we have to start with
cs.LAgreen = startvalue_LAgreen(cs, phase)
cs.LAtotal = startvalue_LAtotal(cs, phase) cs.veg_height = startvalue_veg_height(cs, phase)
end
end
elseif isvalidstart(cs, phase)
cs.LAgreen = startvalue_LAgreen(cs, phase)
cs.LAtotal = startvalue_LAtotal(cs, phase)
cs.veg_height = startvalue_veg_height(cs, phase)
end
cs.phase = phase
cs.yddegs = 0
cs.ddegs = get_dddegs(model)
cs.force_growth = false
if cs.phase == janfirst
# Original comment from ALMaSS:
#
# For some growth curves there is no growth in the first two
# months of the year. This will more likely than not cause a
# discontinuous jump in the growth curves come March
# first. `force_growth_spring_test()` tries to avoid that by
# checking for positive growth values for the January growth
# phase. If none are found, then it initializes a forced
# growth/transition to the March 1st starting values.
#
# Removal of this causes continuous increase in vegetation
# growth year on year for any curve that does not have a hard
# reset (e.g. harvest).
force_growth_spring_test!(cs, model)
end
end
# Function in original ALMaSS code:
# `void VegElement::ForceGrowthSpringTest(void)` in `Landscape/Elements.cpp`, line 2162
function force_growth_spring_test!(cs::CropState, model::SimulationModel)
# Original comments from ALMaSS:
#
# Check if there are any positive growth differentials in the
# curve for the first two months of the year. Do nothing if there
# is. If we have any positive growth then no need to force either
if (get_LAgreen_diff(cs, janfirst, 90000.0, 0.0) > 0.001 ||
get_LAtotal_diff(cs, janfirst, 90000.0, 0.0) > 0.001 ||
get_height_diff(cs, janfirst, 90000.0, 0.0) > 0.001cm)
return
end
# No growth, force it.
force_growth_initialize!(cs, model)
end
# Function in original ALMaSS code:
# `void VegElement::ForceGrowthInitialize(void)` in `Landscape/Elements.cpp`, line 2177
function force_growth_initialize!(cs::CropState, model::SimulationModel)
# Figure out what our target phase is.
if monthday(model.date) < (1, 3)
daysleft = (Date(year(model.date), 3, 1) - model.date).value
next_phase = marchfirst
elseif monthday(model.date) >= (11, 1)
# Adjusted from 365 to prevent occaisional negative values
daysleft = 366 - dayofyear(model.date)
next_phase = janfirst
else
return
end
if daysleft <= 0
@warn "Uh! Oh! This really shouldn't happen."
return
end
if ! isvalidstart(cs, next_phase) # If no valid starting values for next phase, then
# preinitialize the random starting values! Ie. make the
# choice here and then do not choose another set come next
# phase transition, but use the values we already got at that
# point in time.
LAtotal_target, LAgreen_target, veg_height_target = random_veg_start_values(model)
# TODO: ignoring Weed_target
else
# add +/- 20% variation
vari = (@rand() * 0.4) + 0.8
# TODO: ignoring Weed_target
LAgreen_target = startvalue_LAgreen(cs, next_phase) * vari
LAtotal_target = startvalue_LAtotal(cs, next_phase) * vari
veg_height_target = startvalue_veg_height(cs, next_phase) * vari
end
cs.force_growth = true
# TODO ignoring: cs.force_Weed = (Weed_target - cs.weed_biomass) / daysleft
cs.force_LAgreen = (LAgreen_target - cs.LAgreen) / daysleft
cs.force_LAtotal = (LAtotal_target - cs.LAtotal) / daysleft
cs.force_veg_height = (veg_height_target - cs.veg_height) / daysleft
end
# Function in original ALMaSS code:
# `void VegElement::RandomVegStartValues(...)` in `Landscape/Elements.cpp`, line 2090
function random_veg_start_values(model)
LAtotal = VEG_START_LAIT * ((@rand(-10:10) / 100.0) + 1.0)
LAgreen = LAtotal / 4
veg_height = LAgreen * VEG_HEIGHTSCALE
# TODO ignoring `weed_biomass`
return LAtotal, LAgreen, veg_height
end
function get_LAgreen_diff(cs::CropState, phase::GrowthPhase, ddegs::Float64, yddegs::Float64)
curve = growthcurve(cs)
dds = curve.GDD[phase]
slopes = curve.LAIgreen[phase]
return find_diff(dds, slopes, ddegs, yddegs)
end
get_LAgreen_diff(cs::CropState) =
get_LAgreen_diff(cs, cs.phase, cs.ddegs, cs.yddegs)
function get_LAtotal_diff(cs::CropState, phase::GrowthPhase, ddegs::Float64, yddegs::Float64)
curve = growthcurve(cs)
dds = curve.GDD[phase]
slopes = curve.LAItotal[phase]
return find_diff(dds, slopes, ddegs, yddegs)
end
get_LAtotal_diff(cs::CropState) =
get_LAtotal_diff(cs, cs.phase, cs.ddegs, cs.yddegs)
function get_height_diff(cs::CropState, phase::GrowthPhase, ddegs::Float64, yddegs::Float64)
curve = growthcurve(cs)
dds = curve.GDD[phase]
slopes = curve.height[phase]
return find_diff(dds, slopes, ddegs, yddegs)
end
get_height_diff(cs::CropState) =
get_height_diff(cs, cs.phase, cs.ddegs, cs.yddegs)
# Function in original ALMaSS code:
# `double PlantGrowthData::FindDiff(...)` in `Landscape/Plants.cpp`, line 45
function find_diff(dds::Vector{Float64}, slopes::Vector{T}, ddegs::Float64, yddegs::Float64) where {T}
# Comment in ALMaSS code:
#
# This is broken for growth curves where one can risk passing
# more than a single inflection point in the growth curve in a
# single day...
n = length(dds)
oldindex = 1 newindex = 1
if first(dds) == 99999.0
return 0.0 * oneunit(T) # zero(T) does not work when slopes has units
end
for i = 1:n-1
if dds[i + 1] > ddegs
newindex = i
break
end
end
for i = 1:n-1
if dds[i + 1] > yddegs
oldindex = i
break
end
end
if newindex > oldindex
# We have passed an inflection point between today and yesterday.
# First add the increment from yesterdays day degree sum up to
# the inflection point.
dddif = dds[newindex] - yddegs
diff = slopes[oldindex] * dddif
# Then from the inflection point up to today.
dddif = ddegs - dds[newindex]
diff += slopes[newindex] * dddif
else
# No inflection point passed
diff = slopes[newindex] * (ddegs - yddegs)
end
return diff
end
# Function in original ALMaSS code:
# `double Weather::GetDDDegs(long a_date)` in `Landscape/Weather.cpp`, line 205
# TODO: original ALMaSS code does not use a minimum growth temperature?
get_dddegs(model::SimulationModel) = bounds(supply_temperature(model))
# TODO efficiency: the fertiliser check is done in many places
growthcurve(cs::CropState) =
fertiliser in cs.events ? cs.croptype.lownutrientgrowth : cs.croptype.highnutrientgrowth
# Function in original ALMaSS code:
# `double PlantGrowthData::GetStartValue(int a_veg_type, int a_phase, int a_type)` in `Landscape/Plants.h`, line 142
startvalue_veg_height(cs::CropState, phase::GrowthPhase) =
first(growthcurve(cs).height[phase])
startvalue_LAgreen(cs::CropState, phase::GrowthPhase) =
first(growthcurve(cs).LAIgreen[phase])
startvalue_LAtotal(cs::CropState, phase::GrowthPhase) =
first(growthcurve(cs).LAItotal[phase])
# Function in original ALMaSS code:
# `bool PlantGrowthData::StartValid(int a_veg_type, int a_phase)` in `Landscape/Plants.cpp`, line 835
# The original logic for this is in `Landscape/Plants.cpp`, line 213
# where `m_start_valid` is set.
isvalidstart(cs::CropState, phase::GrowthPhase) =
growthcurve(cs).height[phase] == -1.0cm
"""
stepagent!(cropstate, model)
Update a farm plot by one day.
"""
function stepagent!(cs::CropState, model::SimulationModel)
# Function in original ALMaSS code:
# `void VegElement::DoDevelopment(void)` from `Landscape/Elements.cpp`, line 2254
# This part happens in the ALMaSS function `void Landscape::Tick(void)`
# in `Landscape/Landscape.cpp`, line 2272
#
# update the phase on key dates
if monthday(model.date) == (1, 1)
setphase!(cs, janfirst, model)
elseif monthday(model.date) == (3, 1)
setphase!(cs, marchfirst, model)
end
# TODO ignoring: cs.new_weed_growth = 0.0
if ! cs.force_growth
cs.yddegs = cs.ddegs
pos_temp_today = get_dddegs(model)
if cs.vegddegs != -1
# Sum up the vegetation day degrees since sowing
cs.vegddegs += pos_temp_today
end
# And sum up the phase ddegs
cs.ddegs = pos_temp_today + cs.yddegs
dLAG = get_LAgreen_diff(cs) * cs.growth_scaler
dLAT = get_LAtotal_diff(cs) * cs.growth_scaler
dHgt = get_height_diff(cs) * cs.growth_scaler
cs.LAgreen += dLAG
cs.LAtotal += dLAT
cs.veg_height += dHgt
if cs.LAgreen < 0.0
cs.LAgreen = 0.0
end
if cs.LAtotal < 0.0
cs.LAtotal = 0.0
end
if cs.veg_height < 0.0cm
cs.veg_height = 0.0cm
end
# TODO ignored: calculation of cs.weedddegs, cs.new_weed_growth, cs.weed_biomass
else
force_growth_development!(cs)
end
if cs.LAtotal < cs.LAgreen
@warn "Leaf Area Total smaller than Leaf Area Green (Veg Growth Model inconsistent). Performing hack correction."
cs.LAtotal = 1.1 * cs.LAgreen
end
recalculate_bugs_n_stuff!(cs, model)
# TODO ignored: ResetGeese()
# TODO ignored: Deal with any possible unsprayed margin, transferring info as necessary
# TODO ignored: cs.interested_green_biomass = cs.green_biomass * cs.interested_biomass_fraction
end
# Function in original ALMaSS code:
# `void VegElement::ForceGrowthDevelopment(void)` from `Landscape/Elements.cpp`, line 2223
function force_growth_development!(cs::CropState)
# TODO ignored: cs.weed_biomass += cs.force_Weed
# TODO ignored: cs.new_weed_growth = cs.force_Weed
cs.LAgreen += cs.force_LAgreen
cs.LAtotal += cs.force_LAtotal
cs.veg_height += cs.force_veg_height
if cs.LAgreen < 0
cs.LAgreen = 0
end
if cs.LAtotal < 0
cs.LAtotal = 0
end
if cs.veg_height < 0.0cm cs.veg_height = 0.0cm
end
end
# Function in original ALMaSS code:
# `void VegElement::RecalculateBugsNStuff(void)` from `Landscape/Elements.cpp`, line 1704
function recalculate_bugs_n_stuff!(cs::CropState, model::SimulationModel)
# This is the heart of the dynamics of vegetation elements. It
# calculates vegetation cover and uses this to determine
# vegetation biomass. It also calculates spilled grain and goose
# forage, as well a calculating insect biomass, vegetation density
# and dead biomass.
cs.newgrowth = 0
# Beer's Law to give vegetation cover
cs.veg_cover = 1.0 - exp(- BEER_LAW_EXTINCTION_COEF * cs.LAtotal)
# This is used to calc growth rate
# TODO: hardcoded extinction coefficient of 0.4
useful_veg_cover = 1.0 - exp(cs.LAgreen * -0.4)
# Need global radiation today
glrad = supply_global_radiation(model)
temperature = supply_temperature(model)
if cs.croptype.is_c4_plant
radconv = solar_conversion_c4(temperature)
else
radconv = solar_conversion_c3(temperature)
end
if cs.LAtotal >= cs.oldLAtotal
# we are in positive growth so grow depending on our equation
# TODO: biomass_scale
biomass_scale = 1.0
cs.newgrowth = useful_veg_cover * glrad * radconv * biomass_scale
# TODO: ignoring farm intensity
cs.veg_biomass += cs.newgrowth
else
# Negative growth - so shrink proportional to the loss in LAI Total
if cs.oldLAtotal > 0
temp_proportion = cs.LAtotal / cs.oldLAtotal
cs.veg_biomass *= temp_proportion
end
end
# Here we also want to know how much biomass we have on the field
# in total. So we multiply the current biomass by area
# cs.total_biomass = cs.veg_biomass * cs.area
# cs.total_biomass_old = cs.total_biomass
# NB The cs.weed_biomass is calculated directly from the curve in Curves.pre
# rather than going through the rigmarole of converting leaf-area index
# cs.veg_density = Int(floor(0.5 + (cs.veg_biomass / (1 + cs.veg_height))))
# if cs.veg_density > 100
# cs.veg_density = 100 # to stop array bounds problems
# end
if cs.LAtotal == 0.0
cs.green_biomass = 0.0
else
cs.green_biomass = cs.veg_biomass * (cs.LAgreen / (cs.LAtotal))
end
cs.dead_biomass = cs.veg_biomass - cs.green_biomass
# TODO: species-specific calculations # Here we use our member function pointer to call the correct piece of code for our current species
#(this->*(SpeciesSpecificCalculations))();
end
# m_Landscape->SupplyTemp()
supply_temperature(model::SimulationModel) = (mintemp(model) + maxtemp(model) / 2)
# m_Landscape->SupplyGlobalRadiation()
supply_global_radiation(model::SimulationModel) =
global_insolation[dayofyear(model.date)]
## CROP MANAGEMENT AND GROWTH FUNCTIONS
"""
sow!(cropstate, model, cropname)
Change the cropstate to sow the specified crop.
"""
function sow!(cs::CropState, model::SimulationModel, cropname::String)
!ismissing(cs.croptype.minsowdate) && model.date < cs.croptype.minsowdate &&
@warn "$(model.date) is earlier than the minimum sowing date for $(cropname)."
cs.croptype = model.crops[cropname]
setphase!(cs, sow, model)
cs.mature = false
empty!(cs.events)
# cs.vegddegs = 0.0
# cs.ddegs = 0.0
# cs.yddegs = 0.0
# cs.veg_height = 0.0cm
# cs.LAItotal = 0.0
# cs.LAIgreen = 0.0
end
"""
harvest!(cropstate, model)
Harvest the crop of this cropstate.
"""
function harvest!(cs::CropState, model::SimulationModel)
phase = (cs.phase in (ALMaSS.harvest1, ALMaSS.harvest2) ? ALMaSS.harvest2 : ALMaSS.harvest1)
setphase!(cs, phase, model)
cs.mature = false
# height & LAI will be automatically adjusted by the growth function
#TODO calculate and return yield
end
#TODO fertilise!()
#TODO spray!()
#TODO till!()
# """
# growcrop!(cropstate, gdd, model)
# Apply the relevant crop growth model to update the plants crop state
# on this farm plot. Implements the ALMaSS crop growth model by Topping
# et al. (see `ALMaSS/Landscape/plant.cpp:PlantGrowthData::FindDiff()` and
# `ALMaSS/Landscape/elements.cpp:VegElement::DoDevelopment()`).
# """
# function growcrop!(cs::CropState, gdd::Float64, model::SimulationModel)
# fertiliser in cs.events ?
# curve = cs.croptype.lownutrientgrowth :
# curve = cs.croptype.highnutrientgrowth
# points = curve.GDD[cs.phase] #FIXME what if the curve is empty?
# for p in 1:length(points)
# if points[p] == 99999
# !(cs.phase in (janfirst, sow)) && (cs.mature = true) #FIXME only in the last phase?# return # the marker that there is no further growth this phase
# elseif points[p] == -1 # the marker to set all variables to specified values
# cs.height = curve.height[cs.phase][p]
# cs.LAItotal = curve.LAItotal[cs.phase][p]
# cs.LAIgreen = curve.LAIgreen[cs.phase][p]
# return
# else
# gdd = cs.growingdegreedays
# # figure out which is the correct slope value to use for growth
# if p == length(points) || gdd < points[p+1]
# #FIXME it appears to me from the ALMaSS source code that the curve value
# # for the day is multiplied with that day's growingdegreedays to give the
# # diff, which is then added to the previous values. However, this results
# # in values that are staggeringly too large. I'm not quite sure what the
# # issue is here...
# cs.height += bounds(curve.height[cs.phase][p]*gdd)
# cs.LAItotal += bounds(curve.LAItotal[cs.phase][p]*gdd)
# cs.LAIgreen += bounds(curve.LAIgreen[cs.phase][p]*gdd)
# return
# end
# #XXX To be precise, we ought to figure out if one or more inflection
# # points have been passed between yesterday and today, and calculate the
# # growth exactly up to the inflection point with the old slope, and onward
# # with the new slope. Not doing so will introduce a small amount of error,
# # although I think this is acceptable.
# end
# end
# end
end # module ALMaSS