/
coffee.jl
309 lines (252 loc) · 17 KB
/
coffee.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
"""
Coffee crop model
Computes the coffee crop growth and yield. This function is called from [`dynacof`](@ref) and should not be called
by the user.
# Return
Nothing, modify the DataFrame of simulation `Sim` in place. See [`dynacof`](@ref) for more details.
"""
function coffee_model!(Sim, Parameters, Met_c, i, n_i=size(Sim, 1))
# Light interception ------------------------------------------------------
# Metamodel LUE coffee:
Sim.lue[i] = Base.invokelatest(Parameters.lue, Sim, Met_c, i)
#GPP Coffee
Sim.GPP[i] = Sim.lue[i] * Sim.APAR[i]
# Maintenance respiration -------------------------------------------------
# Rm is computed at the beginning of the day on the drymass of the previous day.
# This is considered as the highest priority for the plant (to maintain its dry mass)
after_2 = i <= 2 ? 0 : 1 # used to start respiration after two days so there are some dry mass.
# Resprout (branches) wood:
Sim.Rm_Shoot[i] = after_2 * (Parameters.pa_Shoot * Sim.DM_Shoot[previous_i(i)] * Parameters.NC_Shoot *
Parameters.MRN * Parameters.Q10_Shoot^((Sim.TairCanopy[i] - Parameters.TMR) / 10.0))
# Stump and Coarse roots (perennial wood):
Sim.Rm_SCR[i] = after_2 * (Parameters.pa_SCR * Sim.DM_SCR[previous_i(i)] * Parameters.NC_SCR * Parameters.MRN *
Parameters.Q10_SCR^((Sim.TairCanopy[i] - Parameters.TMR) / 10.0))
# Fruits:
Sim.Rm_Fruit[i] = after_2 * (Parameters.pa_Fruit * Sim.DM_Fruit[previous_i(i)] * Parameters.NC_Fruit * Parameters.MRN *
Parameters.Q10_Fruit^((Sim.TairCanopy[i] - Parameters.TMR) / 10.0))
# Leaves:
Sim.Rm_Leaf[i] = after_2 * (Parameters.pa_Leaf * Sim.DM_Leaf[previous_i(i)] * Parameters.NC_Leaf * Parameters.MRN *
Parameters.Q10_Leaf^((Sim.TairCanopy[i] - Parameters.TMR) / 10.0))
# Fine roots:
Sim.Rm_FRoot[i] = after_2 * (Parameters.pa_FRoot * Sim.DM_FRoot[previous_i(i)] * Parameters.NC_FRoot * Parameters.MRN *
Parameters.Q10_FRoot^((Sim.TSoil[i] - Parameters.TMR) / 10.0))
# Total plant maintenance respiration
Sim.Rm[i] = Sim.Rm_Fruit[i] + Sim.Rm_Leaf[i] + Sim.Rm_Shoot[i] + Sim.Rm_SCR[i] + Sim.Rm_FRoot[i]
# Coffee Allocation -------------------------------------------------------
# Potential use of reserves:
Sim.Consumption_RE[i] = Parameters.kres * Sim.CM_RE[previous_i(i)]
# Supply function
Sim.Supply[i] = max(Sim.GPP[i] - Sim.Rm[i] + Sim.Consumption_RE[i], 0.0)
# If the respiration is greater than the GPP + reserves use, then take this carbon
# from mortality of each compartments' biomass equally (not for fruits or reserves):
Sim.Carbon_Lack_Mortality[i] = -min(0.0, Sim.GPP[i] - Sim.Rm[i] + Sim.Consumption_RE[i])
# 1-Resprout wood ---------------------------------------------------------
# Allocation priority 1, see Charbonnier 2012.
Sim.Alloc_Shoot[i] = Parameters.lambda_Shoot * Sim.Supply[i]
Sim.NPP_Shoot[i] = Sim.Alloc_Shoot[i] / Parameters.epsilon_Shoot
Sim.Rg_Shoot[i] = Sim.Alloc_Shoot[i] - Sim.NPP_Shoot[i]
Sim.Mnat_Shoot[i] = Sim.CM_Shoot[previous_i(i)] / Parameters.lifespan_Shoot
# Pruning
if (Sim.Plot_Age[i] >= Parameters.MeanAgePruning) & (Met_c.DOY[i] == Parameters.D_pruning)
Sim.Mprun_Shoot[i] = Sim.CM_Shoot[previous_i(i)] * Parameters.WoodPruningRate
end
Sim.Mortality_Shoot[i] = min((Sim.Mnat_Shoot[i] + Sim.Mprun_Shoot[i]), Sim.CM_Shoot[previous_i(i)])
# 2-Stump and coarse roots (perennial wood) ------------------------------
Sim.Alloc_SCR[i] = Parameters.lambda_SCR * Sim.Supply[i]
Sim.NPP_SCR[i] = Sim.Alloc_SCR[i] / Parameters.epsilon_SCR
Sim.Rg_SCR[i] = Sim.Alloc_SCR[i] - Sim.NPP_SCR[i]
Sim.Mnat_SCR[i] = Sim.CM_SCR[previous_i(i)] / Parameters.lifespan_SCR
Sim.Mortality_SCR[i] = Sim.Mnat_SCR[i]
# Ratio of number of new nodes per LAI unit as affected by canopy air temperature
# according to Drinnan & Menzel, 1995
# Source "0 Effect T on yield and vegetative growth.xlsx", sheet
# "Std20dComposWinterNodeperBr"
# NB: computed at the end of the vegetatitve growth only to have Tcan of the
# whole period already computed
# NB2 : This is the total number of productive nodes on the coffee plant, i.e. the
# number of green wood nodes that potentially carry flower buds. Green wood mass (and
# so number of nodes) are related to leaf area (new leaves appear on nodes) :
# GUTIERREZ et al. (1998)
if Met_c.DOY[i] == Parameters.DVG2
T_VG = Sim.TairCanopy[(Met_c.year.==Met_c.year[i]).&(Met_c.DOY.>=Parameters.DVG1).&(Met_c.DOY.<=Parameters.DVG2)]
T_VG = sum(T_VG) / length(T_VG)
Sim.ratioNodestoLAI[Met_c.year.>=Met_c.year[i]] .= Parameters.RNL_base * CN(T_VG)
end
# Flower Buds + Flower + Fruits -------------------------------------------
# (1) Buds induction
# Buds start appearing for the very first time from 5500 dd. After that,
# they appear every "Parameters.F_Tffb" degree days until flowering starts
if Sim.BudInitPeriod[i]
Sim.Budinit[i] = (Parameters.a_bud + Parameters.b_bud * (Sim.PAR_Trans_Tree[i] / Parameters.FPAR)) *
Sim.LAI[i-1] * Sim.ratioNodestoLAI[i-1] * Sim.DegreeDays_Tcan[i]
# NB: Number of nodes= Sim.LAI[i-1] * Sim.ratioNodestoLAI[i-1]
Sim.Bud_available[i] = Sim.Budinit[i]
end
# NB: number of fruits ~1200 / year / coffee tree, source : Castro-Tanzi et al. (2014)
# Sim%>%group_by(Plot_Age)%>%summarise(N_Flowers= sum(BudBreak))
# (2) Cumulative degree days experienced by each bud cohort :
dd_i = cumsum(Sim.DegreeDays_Tcan[i:-1:previous_i(i, 1000)])
# (3) Find the window where buds are under dormancy (find the dormant cohorts)
# Bud develops during F_buds1 (840) degree days after initiation, so they cannot
# be dormant less than F_buds1 before i. But they can stay under dormancy until
# F_buds2 dd maximum, so they cannot be older than F_buds2 dd before i.
OldestDormancy = i - (maximum(findall(dd_i .< Parameters.F_buds2)) - 1)
YoungestDormancy = i - (maximum(findall(dd_i .< Parameters.F_buds1)) - 1)
# Idem above (reduce the days computed, F_buds2 is ~300 days and F_buds1 ~80-100 days)
# (4) Test if the condition of minimum required rain for budbreak is met, and if
# not, which cohort first met the condition (starting from younger to older cohorts):
CumRain = cumsum(Met_c.Rain[YoungestDormancy:-1:OldestDormancy])
# (5) Compute the period were all cohorts have encountered all conditions to break
# dormancy :
DormancyBreakPeriod = OldestDormancy:(YoungestDormancy-sum(CumRain .< Parameters.F_rain))
# (6) Temperature effect on bud phenology
Sim.Temp_cor_Bud[i] = Base.invokelatest(Parameters.Bud_T_correction)(Sim.Tleaf_Coffee[i])
# (7) Bud dormancy break, Source, Drinnan 1992 and Rodriguez et al., 2011 eq. 13
Sim.pbreak[i] = 1.0 / (1.0 + exp(Parameters.a_p + Parameters.b_p * Sim.PSIL[i]))
# (8) Compute the number of buds that effectively break dormancy in each cohort:
Sim.BudBreak_cohort[DormancyBreakPeriod] .=
map(min, Sim.Bud_available[DormancyBreakPeriod],
Sim.Budinit[DormancyBreakPeriod] .* Sim.pbreak[i] .* Sim.Temp_cor_Bud[DormancyBreakPeriod])
# NB 1: cannot exceed the number of buds of each cohort
# NB 2: using Budinit and not Bud_available because pbreak is fitted on total bud cohort
# (9) Remove buds that did break dormancy from the pool of dormant buds
Sim.Bud_available[DormancyBreakPeriod] = Sim.Bud_available[DormancyBreakPeriod] .- Sim.BudBreak_cohort[DormancyBreakPeriod]
# (10) Sum the buds that break dormancy from each cohort to compute the total number
# of buds that break dormancy on day i :
Sim.BudBreak[i] = min(sum(Sim.BudBreak_cohort[DormancyBreakPeriod]), Parameters.Max_Bud_Break)
# Rodriguez et al. state that the maximum number of buds that may break dormancy
# during each dormancy-terminating episode was set to 12 (see Table 1).
# Fruits :
FruitingPeriod = i .- findall(dd_i .< Parameters.F_over) .+ 1
# NB : Fruits that are older than the FruitingPeriod are overripped
# Demand from each fruits cohort present on the coffee tree (not overriped),
# same as Demand_Fruit but keeping each value :
demand_distribution = logistic_deriv.(dd_i[1:length(FruitingPeriod)], Parameters.u_log, Parameters.s_log) .*
[0.0; Base.diff(dd_i[1:length(FruitingPeriod)])]
# NB: we use diff because the values are not evenly distributed (it is not grided, e.g. not 1 by 1 increment)
demand_distribution[demand_distribution.==Inf] .= 0.0
Demand_Fruit_Cohort_Period = Sim.BudBreak[FruitingPeriod] .* Parameters.DE_opt .* demand_distribution
# Total C demand of the fruits :
Sim.Demand_Fruit[i] = sum(Demand_Fruit_Cohort_Period)
# C supply to Fruits (i.e. what is left from Supply after removing the consumption
# by previous compartments and Rm):
Sim.Supply_Fruit[i] = Sim.Supply[i] - Sim.Alloc_Shoot[i] - Sim.Alloc_SCR[i]
# Total C allocation to all fruits on day i :
Sim.Alloc_Fruit[i] = min(Sim.Demand_Fruit[i], Sim.Supply_Fruit[i])
# Allocation to each cohort, relative to each cohort demand :
if Sim.Demand_Fruit[i] > 0.0
Rel_DE = Demand_Fruit_Cohort_Period ./ Sim.Demand_Fruit[i]
else
Rel_DE = 0.0
end
Sim.Alloc_Fruit_Cohort[FruitingPeriod] .= Sim.Alloc_Fruit[i] .* Rel_DE
Sim.NPP_Fruit_Cohort[FruitingPeriod] .= Sim.Alloc_Fruit_Cohort[FruitingPeriod] ./ Parameters.epsilon_Fruit
Sim.CM_Fruit_Cohort[FruitingPeriod] .= Sim.CM_Fruit_Cohort[FruitingPeriod] .+ Sim.NPP_Fruit_Cohort[FruitingPeriod]
Sim.DM_Fruit_Cohort[FruitingPeriod] .= Sim.CM_Fruit_Cohort[FruitingPeriod] ./ Parameters.CC_Fruit
# Using CM_Fruit_Cohort_remain to keep track of the fruit mass that is created, but it is updated by removing the overriped
# fruits then, so the overriped fruits can only be removed once.
Sim.CM_Fruit_Cohort_remain[FruitingPeriod] .= Sim.CM_Fruit_Cohort_remain[FruitingPeriod] .+ Sim.NPP_Fruit_Cohort[FruitingPeriod]
# Overriped fruits that fall onto the ground (= to mass of the cohort that overripe) :
overriped_day = max(minimum(FruitingPeriod) - 1, 1)
Sim.Overriped_Fruit[i] = sum(Sim.CM_Fruit_Cohort_remain[previous_i.(overriped_day, 0:10)])
Sim.CM_Fruit_Cohort_remain[previous_i.(overriped_day, 0:10)] .= 0.0
# Sim.Overriped_Fruit[i]= Sim.CM_Fruit_Cohort[minimum(FruitingPeriod)-1.0] * Parameters.epsilon_Fruit
# Duration of the maturation of each cohort born in the ith day (in days):
Sim.Maturation_duration[FruitingPeriod] .= 1:length(FruitingPeriod)
# Sucrose content of each cohort:
Sim.SC[FruitingPeriod] .= Sucrose_cont_perc.(Sim.Maturation_duration[FruitingPeriod], Parameters.S_a, Parameters.S_b,
Parameters.S_x0, Parameters.S_y0)
# Sucrose mass of each cohort
Sim.SM[FruitingPeriod] .= Sim.DM_Fruit_Cohort[FruitingPeriod] .* Sim.SC[FruitingPeriod]
# Harvest maturity:
Sim.Harvest_Maturity_Pot[i] = sum(Sim.SM[FruitingPeriod]) / sum(Sim.DM_Fruit_Cohort[FruitingPeriod] .* ((Parameters.S_y0 .+ Parameters.S_a) ./ 100.0))
# NB : here harvest maturity is computed as the average maturity of the cohorts, because
# all cohorts present in the Coffea are within the 'FruitingPeriod' window.
# It could be computed as the percentage of cohorts that are fully mature (Pezzopane
# et al. 2012 say at 221 days after flowering)
# Optimal sucrose concentration around 8.8% of the dry mass
Sim.NPP_Fruit[i] = Sim.Alloc_Fruit[i] / Parameters.epsilon_Fruit
Sim.Rg_Fruit[i] = Sim.Alloc_Fruit[i] - Sim.NPP_Fruit[i]
# Harvest. Made one day only for now (TODO: make it a period of harvest)
if Parameters.harvest == "quantity"
is_harvest = (Sim.Plot_Age[i] >= Parameters.ageMaturity) &
all(Sim.NPP_Fruit[previous_i.(i, 0:10)] .< Sim.Overriped_Fruit[previous_i.(i, 0:10)]) &
(Sim.CM_Fruit[previous_i(i)] > Parameters.Min_Fruit_CM)
# Made as soon as the fruit dry mass is decreasing for 10 consecutive days.
# This condition is met when fruit overriping is more important than fruit NPP
# for 10 days.
# This option is the best one when fruit maturation is not well known or when the
# harvest is made throughout several days or weeks with the assumption that fruits
# are harvested when mature.
elseif Parameters.harvest == "quality"
is_harvest = (Sim.Plot_Age[i] >= Parameters.ageMaturity) &
(mean(Sim.Harvest_Maturity_Pot[previous_i.(i, 0:9)]) < mean(Sim.Harvest_Maturity_Pot[previous_i.(i, 10:19)]))
# Made as soon as the overall fruit maturation is optimal (all fruits are mature)
else
is_harvest = false
end
if is_harvest
# Save the date of harvest:
Sim.Date_harvest[i] = Met_c.DOY[i]
Sim.Harvest_Fruit[i] = Sim.CM_Fruit[i-1] + Sim.NPP_Fruit[i] - Sim.Overriped_Fruit[i]
Sim.Harvest_Maturity[i] = Sim.Harvest_Maturity_Pot[i]
Sim.CM_Fruit[i-1] = 0.0
Sim.NPP_Fruit[i] = 0.0
Sim.Overriped_Fruit[i] = 0.0
Sim.CM_Fruit_Cohort .= zeros(n_i)
Sim.CM_Fruit_Cohort_remain .= zeros(n_i)
# RV: could harvest mature fruits only (To do).
else
Sim.Harvest_Fruit[i] = 0.0
end
# Leaves ------------------------------------------------------------------
Sim.Supply_Leaf[i] = Parameters.lambda_Leaf_remain * (Sim.Supply[i] - Sim.Alloc_Fruit[i] - Sim.Alloc_Shoot[i] - Sim.Alloc_SCR[i])
Sim.Alloc_Leaf[i] = min(Parameters.DELM * (Parameters.Stocking_Coffee / 10000.0),
Sim.Supply_Leaf[i])
Sim.NPP_Leaf[i] = Sim.Alloc_Leaf[i] / Parameters.epsilon_Leaf
Sim.Rg_Leaf[i] = Sim.Alloc_Leaf[i] - Sim.NPP_Leaf[i]
Sim.Mnat_Leaf[i] = Sim.CM_Leaf[previous_i(i)] / Parameters.lifespan_Leaf
Sim.NPP_RE[i] = Sim.NPP_RE[i] + (Sim.Supply_Leaf[i] - Sim.Alloc_Leaf[i])
Sim.M_ALS[i] = after_2 * max(0.0, Sim.CM_Leaf[previous_i(i)] * Sim.ALS[i])
if (Sim.Plot_Age[i] >= Parameters.MeanAgePruning) & (Met_c.DOY[i] == Parameters.D_pruning)
Sim.Mprun_Leaf[i] = Sim.CM_Leaf[previous_i(i)] * Parameters.LeafPruningRate
else
Sim.Mprun_Leaf[i] = 0.0
end
Sim.Mortality_Leaf[i] = Sim.Mnat_Leaf[i] + Sim.Mprun_Leaf[i] + Sim.M_ALS[i]
# Fine Roots --------------------------------------------------------------
Sim.Supply_FRoot[i] = Parameters.lambda_FRoot_remain * (Sim.Supply[i] - Sim.Alloc_Fruit[i] - Sim.Alloc_Shoot[i] - Sim.Alloc_SCR[i])
Sim.Alloc_FRoot[i] = max(0.0, min(Sim.Alloc_Leaf[i], Sim.Supply_FRoot[i]))
Sim.NPP_FRoot[i] = Sim.Alloc_FRoot[i] / Parameters.epsilon_FRoot
Sim.Rg_FRoot[i] = Sim.Alloc_FRoot[i] - Sim.NPP_FRoot[i]
Sim.NPP_RE[i] = Sim.NPP_RE[i] + (Sim.Supply_FRoot[i] - Sim.Alloc_FRoot[i])
Sim.Mnat_FRoot[i] = Sim.CM_FRoot[previous_i(i)] / Parameters.lifespan_FRoot
Sim.Mprun_FRoot[i] = Parameters.m_FRoot * Sim.Mprun_Leaf[i]
Sim.Mortality_FRoot[i] = Sim.Mnat_FRoot[i] + Sim.Mprun_FRoot[i]
# Biomass -----------------------------------------------------------------
CM_tot = Sim.CM_Leaf[previous_i(i)] + Sim.CM_Shoot[previous_i(i)] + Sim.CM_SCR[previous_i(i)] + Sim.CM_FRoot[previous_i(i)]
Sim.CM_Leaf[i] = Sim.CM_Leaf[previous_i(i)] + Sim.NPP_Leaf[i] - Sim.Mortality_Leaf[i] -
Sim.Carbon_Lack_Mortality[i] * Sim.CM_Leaf[previous_i(i)] / CM_tot
Sim.CM_Shoot[i] = Sim.CM_Shoot[previous_i(i)] + Sim.NPP_Shoot[i] - Sim.Mortality_Shoot[i] -
Sim.Carbon_Lack_Mortality[i] * Sim.CM_Shoot[previous_i(i)] / CM_tot
Sim.CM_Fruit[i] = Sim.CM_Fruit[previous_i(i)] + Sim.NPP_Fruit[i] - Sim.Overriped_Fruit[i]
Sim.CM_SCR[i] = Sim.CM_SCR[previous_i(i)] + Sim.NPP_SCR[i] - Sim.Mortality_SCR[i] -
Sim.Carbon_Lack_Mortality[i] * Sim.CM_SCR[previous_i(i)] / CM_tot
Sim.CM_FRoot[i] = Sim.CM_FRoot[previous_i(i)] + Sim.NPP_FRoot[i] - Sim.Mortality_FRoot[i] -
Sim.Carbon_Lack_Mortality[i] * Sim.CM_FRoot[previous_i(i)] / CM_tot
Sim.CM_RE[i] = Sim.CM_RE[previous_i(i)] + Sim.NPP_RE[i] - Sim.Consumption_RE[i]
Sim.DM_Leaf[i] = Sim.CM_Leaf[i] / Parameters.CC_Leaf
Sim.DM_Shoot[i] = Sim.CM_Shoot[i] / Parameters.CC_Shoot
Sim.DM_Fruit[i] = Sim.CM_Fruit[i] / Parameters.CC_Fruit
Sim.DM_SCR[i] = Sim.CM_SCR[i] / Parameters.CC_SCR
Sim.DM_FRoot[i] = Sim.CM_FRoot[i] / Parameters.CC_FRoots
Sim.DM_RE[i] = Sim.CM_RE[i] / Parameters.CC_SCR
# Total Respiration and NPP -----------------------------------------------
Sim.Rg[i] = Sim.Rg_Fruit[i] + Sim.Rg_Leaf[i] + Sim.Rg_Shoot[i] + Sim.Rg_SCR[i] + Sim.Rg_FRoot[i]
Sim.Ra[i] = Sim.Rm[i] + Sim.Rg[i]
Sim.NPP[i] = Sim.NPP_Shoot[i] + Sim.NPP_SCR[i] + Sim.NPP_Fruit[i] + Sim.NPP_Leaf[i] + Sim.NPP_FRoot[i]
# Compute LAI for the day after based on the CM of this day ---------------
# CM is in gC m-2soil, so use C content to transform in dry mass
Sim.LAI[min(i + 1, n_i)] = Sim.CM_Leaf[i] * Parameters.SLA / 1000.0 / Parameters.CC_Leaf
Sim.LAIplot[min(i + 1, n_i)] = Sim.LAIplot[min(i + 1, n_i)] + Sim.LAI[min(i + 1, n_i)]
end