Something went wrong on our end
-
Marco Matthies authoredMarco Matthies authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
almass.jl 34.32 KiB
### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
###
### This file is responsible for managing the crop growth modules.
###
# 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.
#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
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,
AbstractCropState,
AbstractCropType,
AnnualDate,
Management,
cm,
SimulationModel,
SoilType,
fertiliser,
maxtemp,
mintemp,
year,
dayofyear
import Persefone:
stepagent!,
croptype,
cropname,
cropheight,
cropcover,
cropyield,
makecropstate,
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
function Base.show(io::IO, ::MIME"text/plain", curve::CropCurveParams)
println(io, "CropCurveParams")
println(io, " curveID = ", curve.curveID)
println(io, " highnutrients = ", curve.highnutrients)
println(io, "\nGDD:")
show(io, MIME"text/plain"(), curve.GDD)
println(io, "\n\nLAItotal:")
show(io, MIME"text/plain"(), curve.LAItotal)
println(io, "\n\nLAIgreen:")
show(io, MIME"text/plain"(), curve.LAIgreen)
println(io, "\n\nheight:")
show(io, MIME"text/plain"(), curve.height)
return nothing
end
"""
CropType
The type struct for all crops. Currently follows the crop growth model as
implemented in ALMaSS.
"""
struct CropType <: AbstractCropType
#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
phase_before_harvest::Union{Missing,GrowthPhase} # the phase where the crop can be harvested
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}
biomass_scale::Float64
end
function Base.show(io::IO, ::MIME"text/plain", ct::CropType) println(io, "CropType")
println(io, "=========")
for f in filter(x -> ! (x in (:highnutrientgrowth, :lownutrientgrowth)), fieldnames(typeof(ct)))
println(io, " $f = ", string(getfield(ct, f)))
end
println(io, "\n\nlownutrientgrowth:")
println(io, "--------------------")
show(io, MIME"text/plain"(), ct.lownutrientgrowth)
println(io, "\n\nhighnutrientgrowth:")
println(io, "--------------------")
show(io, MIME"text/plain"(), ct.highnutrientgrowth)
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 <: AbstractCropState
# Class in original ALMaSS code:
# `VegElement` from `Landscape/Elements.h`, line 601
croptype::CropType
events::Vector{Management} = Management[]
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?
function isharvestable(cs::CropState)
phase_before_harvest = croptype(cs).phase_before_harvest
if ismissing(phase_before_harvest) || cs.phase != phase_before_harvest
return false
end
curve = growthcurve(cs)
gdd = curve.GDD[phase_before_harvest]
return cs.ddegs >= last(gdd)
end
function makecropstate(croptype::CropType, model::SimulationModel, soiltype::SoilType)
phase = month(model.date) < 3 ? janfirst : marchfirst
return CropState(; croptype, phase)
end
# 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 for ALMaSS crop model"
# TODO: the second-last column (sowingdensity) uses `String` as
# type. This is because the entry for "peas/beans" has a value of
# "35/80" for the sowingdensity.
cropdata = CSV.File(joinpath(cropdirectory, CROPFILE), missingstring="NA",
types=[String,AnnualDate,AnnualDate,AnnualDate,AnnualDate,Float64,String,Float64,Bool,String,GrowthPhase])
growthdata = CSV.File(joinpath(cropdirectory, GROWTHFILE), missingstring="NA",
types=[Int,String,String,GrowthPhase,String,
Float64,Float64,Float64,Float64])
croptypes = Dict{String,AbstractCropType}() 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"))
croptypes[crop.name] = CropType(crop.name, crop.group, crop.is_c4_plant, crop.phase_before_harvest,
crop.minsowdate, crop.maxsowdate, crop.minharvestdate, crop.maxharvestdate,
crop.mingrowthtemp, highnuts, lownuts, crop.biomass_scale)
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...
# Note: the code is slightly different to the original ALMaSS
# code, which I think is caused by ALMaSS not storing the -1 entry
# in the dds array.
if first(dds) == 99999.0
return 0.0 * oneunit(T) # zero(T) does not work when slopes has units
end
function get_index(array::Tarr, threshold::eltype(Tarr)) where {Tarr <: DenseArray}
idx = findfirst(x -> x > threshold, array)
if isnothing(idx)
idx = lastindex(array) + 1
end
idx -= 1
if idx < firstindex(array)
error("Index too small, idx = $idx")
end
return idx
end
newindex = get_index(dds, ddegs)
oldindex = get_index(dds, yddegs)
dds_newindex = dds[newindex] == -1.0 ? zero(eltype(dds)) : dds[newindex]
slopes_oldindex = dds[oldindex] == -1.0 ? zero(eltype(slopes)) : slopes[oldindex]
slopes_newindex = dds[newindex] == -1.0 ? zero(eltype(slopes)) : slopes[newindex]
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
function growthcurve(cs::CropState)
if ismissing(cs.croptype.lownutrientgrowth) && ismissing(cs.croptype.highnutrientgrowth)
error("No growth curves available for cropstate:\n $cs")
end
if fertiliser in cs.events
curve = cs.croptype.highnutrientgrowth
if ismissing(curve)
@warn "fertiliser used, but highnutrientgrowth curve is missing. Using lownutrientgrowth curve."
curve = cs.croptype.lownutrientgrowth
end
else
curve = cs.croptype.lownutrientgrowth
if ismissing(curve)
# Note: This warning is commented out, as this is quite
# common for e.g. "permanent grassland"
# @warn "fertiliser not used, but lownutrientgrowth curve is missing. Using highnutrientgrowth curve."
curve = cs.croptype.highnutrientgrowth
end
end
return curve
end
# 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
# Note: calls `VegElement::RecalculateBugsNStuff()` in the original ALMaSS code
recalculate_vegetation!(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_vegetation!(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
biomass_scale = cs.croptype.biomass_scale
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)
# TODO: is this correct to remove all events? this would
# e.g. remove fertilising event?
#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)
# height & LAI will be automatically adjusted by the growth function
return cs.veg_biomass # TODO: units ok of cs.veg_biomass ?
end
#TODO fertilise!()
#TODO spray!()
#TODO till!()
end # module ALMaSS