Improve ALMaSS vegetation model: harvest!(), growcrop!()

Redo ALMaSS calculation to start treating -1 and 99999.0 entries
explicitly.

src/crop/almass.jl

@@ -245,6 +245,8 @@ function sow!(cs::CropState, model::SimulationModel, cropname::String)
     cs.croptype = model.crops[cropname]
     cs.phase = ALMaSS.sow
     cs.growingdegreedays = 0.0
+    # TODO: set from crop curve params
+    #       cs.height = cs.croptype.lownutrientgrowth.height[sow]  # or highnutrientgrowth
     cs.height = 0.0cm
     cs.LAItotal = 0.0
     cs.LAIgreen = 0.0
@@ -263,7 +265,27 @@ function harvest!(cs::CropState, model::SimulationModel)
         cs.phase = ALMaSS.harvest1
     cs.growingdegreedays = 0.0
     cs.mature = false
-    # height & LAI will be automatically adjusted by the growth function
+
+    # adjust height, LAItotal, LAIgreen
+    fertiliser in cs.events ?
+        curve = cs.croptype.lownutrientgrowth :
+        curve = cs.croptype.highnutrientgrowth
+    if ismissing(curve)
+        error("Missing curve = $curve")
+    end
+    # TODO: cropname == "no growth" has curve.GDD[harvest1] == 99999.0
+    if length(curve.GDD[cs.phase]) == 0 || !(curve.GDD[cs.phase][1] == -1.0 || curve.GDD[cs.phase][1] == 99999.0)
+        error("Bad crop curve GDD\n  for phase == $(cs.phase), croptype = $(cs.croptype)):\n  $(curve.GDD)")
+    end
+    if (length(curve.height[cs.phase]) == 0
+        || length(curve.LAItotal[cs.phase]) == 0
+        || length(curve.LAIgreen[cs.phase]) == 0)
+        error("One of these has length zero:\n  $curve.height\n  $curve.LAItotal\n  $curve.LAIgreen")
+    end
+    cs.height = curve.height[cs.phase][1]
+    cs.LAItotal = curve.LAItotal[cs.phase][1]
+    cs.LAIgreen = curve.LAIgreen[cs.phase][1]
+
     #TODO calculate and return yield
 end
 
@@ -287,31 +309,30 @@ function growcrop!(cs::CropState, delta_gdd::Real, 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
-            # figure out which is the correct slope value to use for growth
-            if p == length(points) || total_gdd < points[p+1]
-                cs.height += curve.height[cs.phase][p] * delta_gdd
-                cs.LAItotal += curve.LAItotal[cs.phase][p] * delta_gdd
-                cs.LAIgreen += curve.LAIgreen[cs.phase][p] * delta_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
+    if ismissing(curve)
+         #FIXME what to do if the curve is empty?
+        error("Missing curve")
+    end
+    points = curve.GDD[cs.phase]
+    if length(points) == 1 && points[1] == 99999.0
+        !(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
     end
+
+    # TODO: calc oldidx, catch cases where the delta_gdd crosses more
+    #       than one inflection points
+    # TODO: check logic for idx
+    idx = findfirst(x -> x > total_gdd, points)
+    if isnothing(idx)
+        idx = lastindex(points)
+    elseif idx != firstindex(points)
+        idx -= 1  # we need the last index that is not > total_gdd
+    end
+    cs.height += curve.height[cs.phase][idx] * delta_gdd
+    cs.LAItotal += curve.LAItotal[cs.phase][idx] * delta_gdd
+    cs.LAIgreen += curve.LAIgreen[cs.phase][idx] * delta_gdd
+
+    if cs.phase == janfirst && length(points
 end
 
 end # module ALMaSS