diff --git a/src/crop/almass.jl b/src/crop/almass.jl index 45de95a127d0f23d11379fcbd2a67756eb6fd0d6..bc560c1f87a9f353494a994e8fd1fe9f2595148c 100644 --- a/src/crop/almass.jl +++ b/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