From c8538b0aef1c0946c07b6e09828d44ccefaf2bcd Mon Sep 17 00:00:00 2001
From: Marco Matthies <71844+marcom@users.noreply.github.com>
Date: Tue, 15 Oct 2024 18:07:46 +0200
Subject: [PATCH] Add some extra scaffolding for the ALMaSS crop model

---
 src/crop/almass.jl     | 62 ++++++++++++++++++++++++++++++++++++++++--
 src/crop/cropmodels.jl |  2 +-
 test/crop_tests.jl     |  5 ++--
 test/nature_tests.jl   |  2 +-
 4 files changed, 65 insertions(+), 6 deletions(-)

diff --git a/src/crop/almass.jl b/src/crop/almass.jl
index a3d47c2..c400682 100644
--- a/src/crop/almass.jl
+++ b/src/crop/almass.jl
@@ -74,6 +74,7 @@ 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
     minsowdate::Union{Missing,AnnualDate}
     maxsowdate::Union{Missing,AnnualDate}
     minharvestdate::Union{Missing,AnnualDate}
@@ -102,6 +103,7 @@ mutable struct CropState
     #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
 end
 
 croptype(cs::CropState) = cs.croptype
@@ -111,6 +113,58 @@ cropcover(cs::CropState) = 0.0  # TODO: related to LAItotal, LAIgreen?
 cropyield(cs::CropState) = 0.0  # TODO: units? needs biomass?
 isharvestable(cs::CropState) = cs.mature
 
+"Temperature to solar conversion factor for C3 plants."
+const temperature_to_solar_conversion_c3::Vector{Float64} = [
+    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
+]
+
+"Temperature to solar conversion factor for C4 plants."
+const temperature_to_solar_conversion_c4::Vector{Float64} = [
+    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
+]
+
+"""
+    solar_conversion_c3(temperature)
+
+Solar conversion factor (no units) for C3 plants.
+"""
+function solar_conversion_c3(temperature)
+    if temperature < -30 || temperature > 50
+        return 0.0
+    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 0.0
+    end
+    idx = Int(floor(0.5 + temperature)) + 30 + 1
+    return temperature_to_solar_conversion_c4[idx]
+end
+
 """
     Base.tryparse(type, str)
 
@@ -179,8 +233,12 @@ function readcropparameters(cropdirectory::String)
         # 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,
+        # TODO: it would be better to save this in the parameter file
+        # (Note that this matches the current ALMaSS code though,
+        # which also only hardcodes maize as a C4 crop)
+        is_c4_plant = occursin("maize", lowercase(crop.name))
+        croptypes[crop.name] = CropType(crop.name, crop_group, is_c4_plant, crop.minsowdate,
+                                        crop.maxsowdate, crop.minharvestdate, crop.maxharvestdate,
                                         crop.mingrowthtemp, highnuts, lownuts)
     end
     croptypes
diff --git a/src/crop/cropmodels.jl b/src/crop/cropmodels.jl
index 3d84e98..3480da2 100644
--- a/src/crop/cropmodels.jl
+++ b/src/crop/cropmodels.jl
@@ -67,7 +67,7 @@ function makecropstate(model::SimulationModel)
         cs = ALMaSS.CropState(
                 model.crops["natural grass"],
                 phase,
-                0.0, 0.0cm, 0.0, 0.0, false, Vector{Management}()
+                0.0, 0.0cm, 0.0, 0.0, false, Vector{Management}(), false
         )
     elseif @param(crop.cropmodel) == "simple"
         cs = SimpleCrop.CropState(
diff --git a/test/crop_tests.jl b/test/crop_tests.jl
index 106f385..7182ce0 100644
--- a/test/crop_tests.jl
+++ b/test/crop_tests.jl
@@ -19,15 +19,16 @@ const TESTPARAM_SIMPLECROP = joinpath(pkgdir(Persefone), "test", "test_parameter
 end
 
 @testset "Submodule ALMaSS" begin
-    ct = Ps.ALMaSS.CropType("olive tree", "no_crop_group",
+    ct = Ps.ALMaSS.CropType("olive tree", "no_crop_group", false,
                             missing, missing, missing, missing, missing, missing, missing)
     id = 0
     pixels = [(0, 0)]
     farmer = 0
     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))
+                  Ps.ALMaSS.CropState(ct, Ps.ALMaSS.janfirst, 0.0, 0.0m, 0.0, 0.0, mature, events, force_growth))
     @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 d4b723e..a392287 100644
--- a/test/nature_tests.jl
+++ b/test/nature_tests.jl
@@ -52,7 +52,7 @@ end) # end eval
         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}()
+            0.0, 0.0m, 0.0, 0.0, false, Vector{Ps.Management}(), false
         )
     )
     push!(model.farmplots, fp)
-- 
GitLab