From 262d27bd9c5a57838573d9e70d8e13bf5f88ec55 Mon Sep 17 00:00:00 2001
From: Marco Matthies <71844+marcom@users.noreply.github.com>
Date: Thu, 17 Oct 2024 01:07:53 +0200
Subject: [PATCH] Merge port of ALMaSS vegetation code
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.
---
src/crop/almass.jl | 709 +++++++++++++++++++++++++++++++++++------
src/crop/cropmodels.jl | 5 +-
test/crop_tests.jl | 3 +-
test/nature_tests.jl | 4 +-
4 files changed, 616 insertions(+), 105 deletions(-)
diff --git a/src/crop/almass.jl b/src/crop/almass.jl
index c400682..1903cbb 100644
--- a/src/crop/almass.jl
+++ b/src/crop/almass.jl
@@ -9,6 +9,30 @@
#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.
@@ -16,14 +40,17 @@ const CROPFILE = "crop_data_general.csv"
const GROWTHFILE = "crop_growth_curves.csv"
using Persefone:
+ @rand,
+ @u_str,
AnnualDate,
Management,
- Length,
cm,
SimulationModel,
fertiliser,
maxtemp,
- mintemp
+ mintemp,
+ year,
+ dayofyear
import Persefone:
stepagent!,
@@ -39,6 +66,10 @@ import Persefone:
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
@@ -58,9 +89,9 @@ struct CropCurveParams
curveID::Int
highnutrients::Bool
GDD::Dict{GrowthPhase,Vector{Float64}}
- LAItotal::Dict{GrowthPhase,Vector{Float64}}
- LAIgreen::Dict{GrowthPhase,Vector{Float64}}
- height::Dict{GrowthPhase,Vector{Length{Float64}}}
+ LAItotal::Dict{GrowthPhase, Vector{Float64}}
+ LAIgreen::Dict{GrowthPhase, Vector{Float64}}
+ height::Dict{GrowthPhase, Vector{typeof(1.0cm)}}
end
"""
@@ -92,29 +123,64 @@ cropname(ct::CropType) = ct.name
The state data for an ALMaSS vegetation point calculation. Usually
part of a `FarmPlot`.
"""
-mutable struct CropState
- #TODO add Unitful
+@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
- growingdegreedays::Float64
- height::Length{Float64}
- LAItotal::Float64
- LAIgreen::Float64
- #biomass::Float64 #XXX I need to figure out how to calculate this
- mature::Bool #TODO how do we determine this?
- events::Vector{Management} #FIXME does this do anything?
- force_growth::Bool # true means the crop is not sowed but grows by itself, e.g. grass
+ 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.height
-cropcover(cs::CropState) = 0.0 # TODO: related to LAItotal, LAIgreen?
-cropyield(cs::CropState) = 0.0 # TODO: units? needs biomass?
+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::Vector{Float64} = [
+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
@@ -126,17 +192,18 @@ const temperature_to_solar_conversion_c3::Vector{Float64} = [
0 # 50°C
]
+# TODO: units
"Temperature to solar conversion factor for C4 plants."
-const temperature_to_solar_conversion_c4::Vector{Float64} = [
+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, # 0°C to 9°C
- 0, 0, 0, 0, 0, 0, 0, 0, 0.242857, 0.485714, # 10°C to 19°C
- 0.728571, 0.971429, 1.214286, 1.457143, 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.7, 1.7, 1.7, 1.7, # 30°C to 39°C
- 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.53, 1.36, 1.19, 1.02, # 40°C to 49°C
- 0.85, 0.68, 0.51, 0.34, 0.17, 0, 0, 0, 0, 0, # 50°C
- 0
+ 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
]
"""
@@ -146,7 +213,7 @@ Solar conversion factor (no units) for C3 plants.
"""
function solar_conversion_c3(temperature)
if temperature < -30 || temperature > 50
- return 0.0
+ return zero(eltype(temperature_to_solar_conversion_c3))
end
idx = Int(floor(0.5 + temperature)) + 30 + 1
return temperature_to_solar_conversion_c3[idx]
@@ -159,12 +226,92 @@ Solar conversion factor (no units) for C4 plants.
"""
function solar_conversion_c4(temperature)
if temperature < -30 || temperature > 50
- return 0.0
+ 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)
@@ -196,9 +343,9 @@ function buildgrowthcurve(data::Vector{CSV.Row})
LAIgreen = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
harvest2=>Vector{Float64}())
- height = Dict(janfirst=>Vector{Length{Float64}}(), sow=>Vector{Length{Float64}}(),
- marchfirst=>Vector{Length{Float64}}(), harvest1=>Vector{Length{Float64}}(),
- harvest2=>Vector{Length{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)
@@ -244,31 +391,394 @@ function readcropparameters(cropdirectory::String)
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)
- cs.phase = ALMaSS.janfirst
- cs.growingdegreedays = 0.0
+ setphase!(cs, janfirst, model)
elseif monthday(model.date) == (3, 1)
- cs.phase = ALMaSS.marchfirst
- cs.growingdegreedays = 0.0
+ 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
- # update growing degree days
- # if no crop-specific base temperature is given, default to 5°C
- # (https://www.eea.europa.eu/publications/europes-changing-climate-hazards-1/heat-and-cold/heat-and-cold-2014-mean)
- basetemp = cs.croptype.mingrowthtemp
- ismissing(basetemp) && (basetemp = 5.0)
- gdd = (maxtemp(model)+mintemp(model))/2 - basetemp
- gdd > 0 && (cs.growingdegreedays += gdd)
- # update crop growth
- growcrop!(cs, gdd, model)
+ 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
@@ -281,13 +791,16 @@ 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]
- cs.phase = ALMaSS.sow
- cs.growingdegreedays = 0.0
- cs.height = 0.0cm
- cs.LAItotal = 0.0
- cs.LAIgreen = 0.0
+ setphase!(cs, sow, model)
cs.mature = false
- cs.events = Vector{Management}()
+ 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
"""
@@ -296,10 +809,8 @@ end
Harvest the crop of this cropstate.
"""
function harvest!(cs::CropState, model::SimulationModel)
- cs.phase in [ALMaSS.harvest1, ALMaSS.harvest2] ?
- cs.phase = ALMaSS.harvest2 :
- cs.phase = ALMaSS.harvest1
- cs.growingdegreedays = 0.0
+ 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
@@ -309,49 +820,51 @@ end
#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
+
+# """
+# 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
diff --git a/src/crop/cropmodels.jl b/src/crop/cropmodels.jl
index 3480da2..92722a1 100644
--- a/src/crop/cropmodels.jl
+++ b/src/crop/cropmodels.jl
@@ -65,9 +65,8 @@ function makecropstate(model::SimulationModel)
if @param(crop.cropmodel) == "almass"
phase = (month(model.date) < 3 ? ALMaSS.janfirst : ALMaSS.marchfirst)
cs = ALMaSS.CropState(
- model.crops["natural grass"],
- phase,
- 0.0, 0.0cm, 0.0, 0.0, false, Vector{Management}(), false
+ croptype = model.crops["natural grass"],
+ phase = phase,
)
elseif @param(crop.cropmodel) == "simple"
cs = SimpleCrop.CropState(
diff --git a/test/crop_tests.jl b/test/crop_tests.jl
index 7182ce0..38ae55a 100644
--- a/test/crop_tests.jl
+++ b/test/crop_tests.jl
@@ -27,8 +27,7 @@ end
mature = false
events = Ps.Management[]
force_growth = false
- fp = FarmPlot(id, pixels, farmer,
- Ps.ALMaSS.CropState(ct, Ps.ALMaSS.janfirst, 0.0, 0.0m, 0.0, 0.0, mature, events, force_growth))
+ fp = FarmPlot(id, pixels, farmer, Ps.ALMaSS.CropState(croptype=ct, phase=Ps.ALMaSS.janfirst))
@test fp isa FarmPlot
@test fp isa FarmPlot{Ps.ALMaSS.CropState}
@test croptype(fp) isa Ps.ALMaSS.CropType
diff --git a/test/nature_tests.jl b/test/nature_tests.jl
index a392287..0a786ca 100644
--- a/test/nature_tests.jl
+++ b/test/nature_tests.jl
@@ -51,8 +51,8 @@ end) # end eval
fp = Ps.FarmPlot(
1, [(6,6)], 1,
Ps.ALMaSS.CropState(
- model.crops["winter wheat"], Ps.ALMaSS.janfirst,
- 0.0, 0.0m, 0.0, 0.0, false, Vector{Ps.Management}(), false
+ croptype = model.crops["winter wheat"],
+ phase = Ps.ALMaSS.janfirst,
)
)
push!(model.farmplots, fp)
--
GitLab