diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5b7ce0062ef1f0004ad86815c06c5b1fd5dd198d..039086baf9ea6f9fc08215e6d8f7478aad02c6a9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,10 +10,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
*Aim: 3 species, 2 crop growth models, farm model, GAEC scenarios, experimental analysis*
-## [0.6.0] - unreleased
-
*Plan: decouple CairoMakie (#81), fix & test ALMaSS, set up first experiments*
+## [0.6.0]
+
+**This point release re-implements the ALMaSS crop model**
+
+
### Added
- `crop.cropdirectory` parameter specifies folder in which all crop data files for the
@@ -28,6 +31,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- `simulate()` and `initialise()` now take a `params` keyword argument that can be used to
override parameters from other input sources
+- The ALMaSS crop data config file
+ `data/crops/almass/crop_data_general.csv` now has extra columns for
+ `is_c4_plant`, `sowingdensity`, and `phase_before_harvest`
+
### Deprecated
### Removed
@@ -37,6 +44,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
### Fixed
+- The implementation of the ALMaSS vegetation model in Persefone has
+ been completely rewritten, hopefully more faithfully reproducing the
+ logic in ALMaSS. The resulting plant heights are now more realistic
+ and do not produce the extreme plant heights seen previously (which
+ was due to an erroneous interpretation of the ALMaSS growth curves).
+
+
---
## [0.5.5] - 09-08-2024
diff --git a/Project.toml b/Project.toml
index b34a8e3bb74ffe2a863e273a054352ef359724ac..da7164d8dd4153c9eb27a504f91d966a84799468 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
name = "Persefone"
uuid = "039acd1d-2a07-4b33-b082-83a1ff0fd136"
authors = ["Daniel Vedder <daniel.vedder@idiv.de>"]
-version = "0.5.5"
+version = "0.6.0"
[deps]
ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
diff --git a/data/crops/almass/crop_data_general.csv b/data/crops/almass/crop_data_general.csv
index c779dec5fa35c5777bb5363b88a81dbd8556bddc..1b34e04465790bf212b5622cdcc6d8ff1a651f0a 100644
--- a/data/crops/almass/crop_data_general.csv
+++ b/data/crops/almass/crop_data_general.csv
@@ -1,28 +1,28 @@
-name,minsowdate,maxsowdate,minharvestdate,maxharvestdate,mingrowthtemp,group
-"spring rape",NA,NA,NA,NA,NA,"grain"
-"winter rape","20 August","25 August",NA,NA,NA,"grain"
-"winter wheat","15 October","31 October",NA,NA,0,"grain"
-"spring wheat",NA,NA,NA,NA,NA,"grain"
-"winter barley","15 September","30 September",NA,NA,0,"grain"
-"spring barley","1 March","10 April",NA,NA,0,"grain"
-"undersown spring barley",NA,NA,NA,NA,0,"grain"
-"winter rye","23 September","15 October",NA,NA,NA,"grain"
-"triticale","25 September","10 October",NA,NA,NA,"grain"
-"oats",NA,NA,NA,NA,NA,"grain"
-"maize","15 April","30 April",NA,NA,8,"grain"
-"potatoes",NA,NA,NA,NA,4,"root"
-"carrots",NA,NA,NA,NA,NA,"root"
-"beet","15 March","10 May",NA,NA,NA,"root"
-"sunflower","25 March","15 April",NA,NA,NA,"other"
-"lucerne",NA,NA,NA,NA,NA,"legumes"
-"peas/beans","15 February","15 March",NA,NA,5,"legumes"
-"silage clover/grass",NA,NA,NA,NA,NA,"legumes"
-"fodder/clover",NA,NA,NA,NA,NA,"legumes"
-"lawn",NA,NA,NA,NA,NA,"grass"
-"permanent grassland (grazed)",NA,NA,NA,NA,NA,"grass"
-"permanent grassland (seeded)",NA,NA,NA,NA,NA,"grass"
-"permanent grassland (low yield)",NA,NA,NA,NA,NA,"grass"
-"permanent set-aside",NA,NA,NA,NA,NA,"semi-natural"
-"natural grass",NA,NA,NA,NA,NA,"semi-natural"
-"no growth",NA,NA,NA,NA,NA,"semi-natural"
-"heath",NA,NA,NA,NA,NA,"semi-natural"
+name,minsowdate,maxsowdate,minharvestdate,maxharvestdate,mingrowthtemp,group,biomass_scale,is_c4_plant,sowingdensity,phase_before_harvest
+spring rape,NA,NA,NA,NA,NA,grain,1.071,false,75,marchfirst
+winter rape,20 August,25 August,NA,NA,NA,grain,1.071,false,75,marchfirst
+winter wheat,15 October,31 October,NA,NA,0,grain,1.0,false,400,marchfirst
+spring wheat,NA,NA,NA,NA,NA,grain,1.0,false,475,marchfirst
+winter barley,15 September,30 September,NA,NA,0,grain,0.857,false,340,marchfirst
+spring barley,1 March,10 April,NA,NA,0,grain,0.857,false,360,marchfirst
+undersown spring barley,NA,NA,NA,NA,0,grain,0.857,false,NA,marchfirst
+winter rye,23 September,15 October,NA,NA,NA,grain,0.857,false,290,marchfirst
+triticale,25 September,10 October,NA,NA,NA,grain,1.0,false,320,marchfirst
+oats,NA,NA,NA,NA,NA,grain,0.857,false,385,marchfirst
+maize,15 April,30 April,NA,NA,8,grain,1.0,true,13,sow
+potatoes,NA,NA,NA,NA,4,root,0.857,false,5,marchfirst
+carrots,NA,NA,NA,NA,NA,root,0.7857,false,160,marchfirst
+beet,15 March,10 May,NA,NA,NA,root,0.857,false,16,marchfirst
+sunflower,25 March,15 April,NA,NA,NA,other,1.0,false,8,marchfirst
+lucerne,NA,NA,NA,NA,NA,legumes,1.2,false,NA,marchfirst
+peas/beans,15 February,15 March,NA,NA,5,legumes,0.857,false,35/80,marchfirst
+silage clover/grass,NA,NA,NA,NA,NA,legumes,1.1,false,NA,marchfirst
+fodder/clover,NA,NA,NA,NA,NA,legumes,1.2,false,NA,marchfirst
+lawn,NA,NA,NA,NA,NA,grass,0.5,false,NA,marchfirst
+permanent grassland (grazed),NA,NA,NA,NA,NA,grass,1.1,false,NA,NA
+permanent grassland (seeded),NA,NA,NA,NA,NA,grass,1.1,false,NA,NA
+permanent grassland (low yield),NA,NA,NA,NA,NA,grass,1.0,false,NA,NA
+permanent set-aside,NA,NA,NA,NA,NA,semi-natural,0.7857,false,NA,NA
+natural grass,NA,NA,NA,NA,NA,semi-natural,0.567,false,NA,NA
+no growth,NA,NA,NA,NA,NA,semi-natural,0.0,false,NA,NA
+heath,NA,NA,NA,NA,NA,semi-natural,0.567,false,NA,NA
diff --git a/src/crop/almass.jl b/src/crop/almass.jl
index a3d47c238c953b703ad3dccd72b6eac9d0efec7d..59b94b9f0cf6ace3c5ebd6869dc10aa56e6b3542 100644
--- a/src/crop/almass.jl
+++ b/src/crop/almass.jl
@@ -3,6 +3,31 @@
### 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.
@@ -16,14 +41,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 +67,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 +90,24 @@ 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
+
+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
"""
@@ -74,6 +121,8 @@ struct CropType
# 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}
@@ -81,6 +130,21 @@ struct CropType
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
@@ -91,25 +155,197 @@ 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[]
+
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?
+ 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?
-isharvestable(cs::CropState) = cs.mature
+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
+
+# 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)
@@ -142,9 +378,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)
@@ -163,8 +399,11 @@ Parse a CSV file containing the required parameter values for each crop
"""
function readcropparameters(cropdirectory::String)
@debug "Reading crop parameters"
+ # 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])
+ 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])
@@ -176,41 +415,428 @@ function readcropparameters(cropdirectory::String)
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"
- croptypes[crop.name] = CropType(crop.name, crop_group, crop.minsowdate, crop.maxsowdate,
- crop.minharvestdate, crop.maxharvestdate,
- crop.mingrowthtemp, highnuts, lownuts)
+ 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)
- 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
+
+ # 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
- # 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)
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
@@ -223,13 +849,17 @@ 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
- cs.mature = false
- cs.events = Vector{Management}()
+ 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
"""
@@ -238,62 +868,14 @@ 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
- cs.mature = false
+ 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
- #TODO calculate and return yield
+ return cs.veg_biomass # TODO: units ok of cs.veg_biomass ?
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
diff --git a/src/crop/cropmodels.jl b/src/crop/cropmodels.jl
index 3d84e98f7d6e87967d76163cae07d3028549dfc5..92722a11be435ccfe5d75a1900331bf45a6142cd 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}()
+ croptype = model.crops["natural grass"],
+ phase = phase,
)
elseif @param(crop.cropmodel) == "simple"
cs = SimpleCrop.CropState(
diff --git a/src/crop/simplecrop.jl b/src/crop/simplecrop.jl
index 348d9448561f693d1e1f951d33b9b7781ea2187a..bc4f11dc1c9cd4b52d18274ff626af7e73cddabd 100644
--- a/src/crop/simplecrop.jl
+++ b/src/crop/simplecrop.jl
@@ -43,7 +43,7 @@ cropyield(cs::CropState) = 0.0 # TODO: units?
isharvestable(cs::CropState) = true
"""
- stepagent!(farmplot, model)
+ stepagent!(cropstate, model)
Update a crop state by one day.
"""
diff --git a/src/farm/farm.jl b/src/farm/farm.jl
index 3bc3109e4efc086672a26681f147c3deb78065ce..c60cb93db7b39d423f388131bdd9a3f5ececa103 100644
--- a/src/farm/farm.jl
+++ b/src/farm/farm.jl
@@ -70,7 +70,7 @@ function stepagent!(farmer::BasicFarmer, model::SimulationModel)
@warn "minsowdate and/or maxsowdate is `missing` for crop \"$(ctype.name)\", not assigning sowdates"
end
end
- elseif cropname(field) == "no growth" && model.date == farmer.sowdates[f]
+ elseif cropname(field) == "no growth" && haskey(farmer.sowdates, f) && model.date == farmer.sowdates[f]
# if a field has been harvested, check if the next crop can be sown
@sow(farmer.croprotations[f][1])
end
diff --git a/test/crop_tests.jl b/test/crop_tests.jl
index 106f3854133a6d55e805551e822f1d4b092d2bfc..b412030f5b822b968c616f76524c96174d54d6d7 100644
--- a/test/crop_tests.jl
+++ b/test/crop_tests.jl
@@ -19,15 +19,15 @@ const TESTPARAM_SIMPLECROP = joinpath(pkgdir(Persefone), "test", "test_parameter
end
@testset "Submodule ALMaSS" begin
- ct = Ps.ALMaSS.CropType("olive tree", "no_crop_group",
- missing, missing, missing, missing, missing, missing, missing)
+ ct = Ps.ALMaSS.CropType("olive tree", "no_crop_group", false, missing,
+ missing, missing, missing, missing, missing, missing, missing, 1.0)
id = 0
pixels = [(0, 0)]
farmer = 0
mature = false
events = Ps.Management[]
- fp = FarmPlot(id, pixels, farmer,
- Ps.ALMaSS.CropState(ct, Ps.ALMaSS.janfirst, 0.0, 0.0m, 0.0, 0.0, mature, events))
+ force_growth = false
+ 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/cropparams-almass/crop_data_general.csv b/test/cropparams-almass/crop_data_general.csv
deleted file mode 100644
index c779dec5fa35c5777bb5363b88a81dbd8556bddc..0000000000000000000000000000000000000000
--- a/test/cropparams-almass/crop_data_general.csv
+++ /dev/null
@@ -1,28 +0,0 @@
-name,minsowdate,maxsowdate,minharvestdate,maxharvestdate,mingrowthtemp,group
-"spring rape",NA,NA,NA,NA,NA,"grain"
-"winter rape","20 August","25 August",NA,NA,NA,"grain"
-"winter wheat","15 October","31 October",NA,NA,0,"grain"
-"spring wheat",NA,NA,NA,NA,NA,"grain"
-"winter barley","15 September","30 September",NA,NA,0,"grain"
-"spring barley","1 March","10 April",NA,NA,0,"grain"
-"undersown spring barley",NA,NA,NA,NA,0,"grain"
-"winter rye","23 September","15 October",NA,NA,NA,"grain"
-"triticale","25 September","10 October",NA,NA,NA,"grain"
-"oats",NA,NA,NA,NA,NA,"grain"
-"maize","15 April","30 April",NA,NA,8,"grain"
-"potatoes",NA,NA,NA,NA,4,"root"
-"carrots",NA,NA,NA,NA,NA,"root"
-"beet","15 March","10 May",NA,NA,NA,"root"
-"sunflower","25 March","15 April",NA,NA,NA,"other"
-"lucerne",NA,NA,NA,NA,NA,"legumes"
-"peas/beans","15 February","15 March",NA,NA,5,"legumes"
-"silage clover/grass",NA,NA,NA,NA,NA,"legumes"
-"fodder/clover",NA,NA,NA,NA,NA,"legumes"
-"lawn",NA,NA,NA,NA,NA,"grass"
-"permanent grassland (grazed)",NA,NA,NA,NA,NA,"grass"
-"permanent grassland (seeded)",NA,NA,NA,NA,NA,"grass"
-"permanent grassland (low yield)",NA,NA,NA,NA,NA,"grass"
-"permanent set-aside",NA,NA,NA,NA,NA,"semi-natural"
-"natural grass",NA,NA,NA,NA,NA,"semi-natural"
-"no growth",NA,NA,NA,NA,NA,"semi-natural"
-"heath",NA,NA,NA,NA,NA,"semi-natural"
diff --git a/test/nature_tests.jl b/test/nature_tests.jl
index d4b723e5716c7128786213340a6e20c84978ff23..0a786cabd2cd81eacd45fd12d9930ee561d57566 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}()
+ croptype = model.crops["winter wheat"],
+ phase = Ps.ALMaSS.janfirst,
)
)
push!(model.farmplots, fp)
diff --git a/test/test_parameters.toml b/test/test_parameters.toml
index 47c5a08e094537da28b0d519f4776e1228f8c04b..a8a614e780e7ee1ea9368de9b0db168358a88207 100644
--- a/test/test_parameters.toml
+++ b/test/test_parameters.toml
@@ -41,4 +41,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
[crop]
cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "test/cropparams-almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
diff --git a/test/test_parameters_almass.toml b/test/test_parameters_almass.toml
index fa520c207a6980016d2595ed5d16c78f76137af4..4eaae831355219d74fe90b3fe8ebc1e217c8d900 100644
--- a/test/test_parameters_almass.toml
+++ b/test/test_parameters_almass.toml
@@ -41,4 +41,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
[crop]
cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "test/cropparams-almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
diff --git a/test/test_parameters_simplecrop.toml b/test/test_parameters_simplecrop.toml
index 43bd55087732774a792cb39d098b7e3abcfce43c..b5cd61e6f69f60f1268db06588e9d86547aec52a 100644
--- a/test/test_parameters_simplecrop.toml
+++ b/test/test_parameters_simplecrop.toml
@@ -41,4 +41,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
[crop]
cropmodel = "simple" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "test/cropparams-almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
diff --git a/test/zzz-cropparams-almass/crop_data_general.csv b/test/zzz-cropparams-almass/crop_data_general.csv
new file mode 100644
index 0000000000000000000000000000000000000000..9af47222ca7e54057816d197135aa966266a2b1f
--- /dev/null
+++ b/test/zzz-cropparams-almass/crop_data_general.csv
@@ -0,0 +1,28 @@
+name,minsowdate,maxsowdate,minharvestdate,maxharvestdate,mingrowthtemp,group,biomass_scale,is_c4_plant,sowingdensity
+spring rape,NA,NA,NA,NA,NA,grain,1.071,false,75
+winter rape,20 August,25 August,NA,NA,NA,grain,1.071,false,75
+winter wheat,15 October,31 October,NA,NA,0,grain,1.0,false,400
+spring wheat,NA,NA,NA,NA,NA,grain,1.0,false,475
+winter barley,15 September,30 September,NA,NA,0,grain,0.857,false,340
+spring barley,1 March,10 April,NA,NA,0,grain,0.857,false,360
+undersown spring barley,NA,NA,NA,NA,0,grain,0.857,false,NA
+winter rye,23 September,15 October,NA,NA,NA,grain,0.857,false,290
+triticale,25 September,10 October,NA,NA,NA,grain,1.0,false,320
+oats,NA,NA,NA,NA,NA,grain,0.857,false,385
+maize,15 April,30 April,NA,NA,8,grain,1.0,true,13
+potatoes,NA,NA,NA,NA,4,root,0.857,false,5
+carrots,NA,NA,NA,NA,NA,root,0.7857,false,160
+beet,15 March,10 May,NA,NA,NA,root,0.857,false,16
+sunflower,25 March,15 April,NA,NA,NA,other,1.0,false,8
+lucerne,NA,NA,NA,NA,NA,legumes,1.2,false,NA
+peas/beans,15 February,15 March,NA,NA,5,legumes,0.857,false,35/80
+silage clover/grass,NA,NA,NA,NA,NA,legumes,1.1,false,NA
+fodder/clover,NA,NA,NA,NA,NA,legumes,1.2,false,NA
+lawn,NA,NA,NA,NA,NA,grass,0.5,false,NA
+permanent grassland (grazed),NA,NA,NA,NA,NA,grass,1.1,false,NA
+permanent grassland (seeded),NA,NA,NA,NA,NA,grass,1.1,false,NA
+permanent grassland (low yield),NA,NA,NA,NA,NA,grass,1.0,false,NA
+permanent set-aside,NA,NA,NA,NA,NA,semi-natural,0.7857,false,NA
+natural grass,NA,NA,NA,NA,NA,semi-natural,0.567,false,NA
+no growth,NA,NA,NA,NA,NA,semi-natural,0.0,false,NA
+heath,NA,NA,NA,NA,NA,semi-natural,0.567,false,NA
diff --git a/test/cropparams-almass/crop_growth_curves.csv b/test/zzz-cropparams-almass/crop_growth_curves.csv
similarity index 100%
rename from test/cropparams-almass/crop_growth_curves.csv
rename to test/zzz-cropparams-almass/crop_growth_curves.csv