%FUN version 2.0
%This is an update to the basic FUN model code of Fisher et al. (2010) that
%predicts nitrogen (N)uptake, the carbon (C) costs of N acquisition, and
%the amount of C for growth. The code here is meant to be run offline and
%can be run at timesteps ranging from annual to daily. However, FUN can be
%coupled and intergrated into global models (e.g., JULES).
%In the model, there are four strategies to take up N: (1)passive uptake of N
%throught ET stream (no cost) (2) resorption of leaf N (3) active uptake of
%soil N (4) fixation. The model caluclates the C cost of each strategy and
%allocates C to the cheapest source of N until N demand is met or the pool
%is exhausted.
%----updated by Edward Brzostek (edbrzost@indiana.edu)---
%Citation: Fisher J.B., Sitch S., Malhi Y., Fisher R.A., Huntingford C. &
%Tan S.Y. (2010). Carbon cost of plant nitrogen acquisition: A mechanistic,
%globally applicable model of plant nitrogen uptake, retranslocation, and
%fixation. Global Biogeochemical Cycles, 24.
%The main improvements to FUN code are:
%(1) At the annual timestep, the C cost of resorption and active uptake is now iteratively
%coded to become more costly for every unit N taken from the leaf N or soil N pool,
%and (2) at all timesteps, the solubility and concentration of soil N
%in solution for passive uptake decreases as soil water depth (sd) declines.
%Previously at the annual scale, a single C cost for resoprtion or active uptake
%ed to resorption efficiencies or active uptake approaching 90-100% of available pool.
%Similarly, low soil water depth resulted in a highly concentrated pool of available N
%in soil solution that was easily exhausted by passive uptake to meet plant
%N demand.
%(1) MODEL INPUTS
%NPP = net primary production(kg C/m2/timestep)
%plantcn = C to N ratio of live standing biomass (kg C/kg N)
%ET = evapotranpiration (km/timestep)
%sd= soil water depth (km); measure of the amount of total soil volume that
%the plants can take up water from (km)
%rootbiomass=amount of C in roots (kg C/m2)
%leafN = amount of N in leaves (kg N/m2)
%soilN = amount of available soil N (kg N/m2)
%soilT = soil temperature](degrees C)
%fixer = logical yes or no if ecosystem has fixation
%(2) COST PARAMETERS
%Fixation cost parameters (cost_fixN) based off of Houlton et al. 2008
s_FIX = -5;
a_FIX = -3.62;
b_FIX = 0.27;
c_FIX = 25.15;
%Active uptake cost parameters (cost_active)
kN=1;
kC=1;
%Resorb cost parameters (cost_resorb)
kR=0.01;
%(3)MODEL CODE.
%Note: At monthly to annual timesteps, the iterative cost function of
%resoprtion and active uptake should be used. This iterative approach
%divides the leafN pool or soilN pool by the number of iterations to allow
%the cost of uptake to increase per unit N taken out of each pool.
%At daily timesteps, this is not necessary given
%the magnitude of daily N demand vs. leaf N or soil N. In addition, this iterative
%approach is computionally costly at daily timesteps. The code has three
%rounds to allow each strategy to occur if the soilN pool or leafN pool is
%exhausted.
%Below, the user can specify number of iterations to perform. Default is set
%to one for daily time steps.
iterations=1;
%preallocate leafN and soilN matrices to calculate iterative cost functions
%At daily time step this step will be skipped.
placeholder=nan(length(NPP),iterations);
leafN=horzcat(leafN,placeholder(:,2:iterations)); %leafN pool is depleted according to number of times N is taken out
soilN=horzcat(soilN,placeholder(:,2:iterations)); %soilN pool is depleted according to number of times N is taken out
soildraw=nan(length(NPP),iterations); %Number of times N is taken up from soilN pool (n-1)
soildraw(:,1)=1;
leafdraw=nan(length(NPP),iterations); %Number of times N is taken up from soilN pool (n-1)
leafdraw(:,1)=1;
%PREALLOCATION OF OUTPUT MATRICES FOR ROUND ONE
Npassive=nan(length(NPP),1); %(kg N/m2/timestep)amount of N taken up through passive uptake
soilNpassive=nan(length(NPP),1);%(kg N/m2) amount of N available for passive uptake
Ndemand=nan(length(NPP),1); %(kg N/m2/timestep)N needed to support NPP
cost_fixN=nan(length(NPP),1);%(kg C/m2/timestep) cost of N fixation
cost_active=nan(length(NPP),iterations);%(kg C/m2/timestep) cost of active N uptake
cost_resorb=nan(length(NPP),iterations);%(kg C/m2/timestep) cost of resorption of leaf N
cost_acq=nan(length(NPP),iterations);%(kg C/m2/timestep) minimun of cost (1) fixation,
%(2) active, and (3)resorb.
uptake_strategy=nan(length(NPP),iterations);%(n=1:3) indicates pathway
Cgrowth=nan(length(NPP),iterations);%(kg C/m2/timestep) C available for NPP and growth
Cgrowth(:,1)=NPP(:,1);
Cacq=nan(length(NPP),iterations);%(kg C/m2/timestep) C expended on N uptake
Cacq(:,1)=0;
Nacq=nan(length(NPP),iterations);%(kg N/m2/timestep) Total N taken up
Nacq(:,1)=0;
Cavailable=nan(length(NPP),iterations);%C available to expend to growth or N uptake
rsoilN=nan(length(NPP),1); %soilN remaining after passive uptake
Ndeficit=nan(length(NPP),iterations);%N uptake needed to meet N demand
final_leafN=nan(length(NPP),iterations); %N remaining in leaf pool
final_soilN=nan(length(NPP),iterations);%N remaining in soil pool
%MODEL CODE FOR ROUND 1
for i=1:length(NPP)
Ndemand(i,1)=NPP(i)/plantcn(i); %Calculate N demand
%Calculates Npassive with Ndemand and soilN as constraints
%the relationship between soilN and soil water depth is modeled as a first
%order rate equation with the constant (kp) equal to 0.001
kp=0.001; %constant in equation
soilNpassive(i)=soilN(i)*(1-exp(-kp*sd(i)));
Npassive(i,1)=min(soilN(i)*(ET(i)/sd(i)),Ndemand(i,1));
if Npassive(i)>soilN(i);
Npassive(i)=soilN(i);
end
rsoilN(i,1)=soilN(i)-Npassive(i);
%Calculates the cost for each uptake strategy and includes or excludes
%fixation based upon presence of N fixers.
cost_fixN(i)=s_FIX*exp(a_FIX+b_FIX*soilT(i)*(1-0.5*soilT(i)/c_FIX))-2*s_FIX;
if fixer(i)==0; % exclude fixation for non fixers
cost_fixN(i)=9999999; % random high number
end
%As you deplete leaf N or soil N pool the cost rises
for j=1:iterations
Cavailable(i,j)=NPP(i)-nansum(Cacq(i,1:j));
leafN(i,j)=((iterations-(leafdraw(i,j)-1))/iterations)*leafN(i,1);
cost_resorb(i,j)=kR/leafN(i,j);
soilN(i,j)=((iterations-(soildraw(i,j)-1))/iterations)*rsoilN(i,1);
cost_active(i,j)=(kN/soilN(i,j))*(kC/rootbiomass(i));
if Npassive(i)==Ndemand(i); % truncate uptake if Ndemand is met
cost_active(i,j)= 9999999; % random high number
cost_resorb(i,j)=9999999; % random high number
end
%Calculates which N strategy have the lowest cost
if cost_resorb(i,j)0;
Nacq(i,j)=min(Ndeficit(i,j),(Cavailable(i,j)-Npassive(i)*plantcn(i))/(cost_acq(i,j)+plantcn(i)));
Cacq(i,j)=(cost_acq(i,j)*Cavailable(i,j)-Npassive(i)*plantcn(i)*cost_acq(i,j))/(cost_acq(i,j)+plantcn(i));
Cgrowth(i,j)=(plantcn(i)*Cavailable(i,j)+Npassive(i)*plantcn(i)*cost_acq(i,j))/(cost_acq(i,j)+plantcn(i));
if uptake_strategy(i,j)==2
final_leafN(i,j)=leafN(i,j);
final_soilN(i,j)=soilN(i,j)-Nacq(i,j);
if uptake_strategy(i,j)==1
final_soilN(i,j)=soilN(i,j);
final_leafN(i,j)=leafN(i,j)-Nacq(i,j);
end
end
end
if uptake_strategy(i,j)==2 && Nacq(i,j)>soilN(i,j)-((iterations-(soildraw(i,j)))/iterations)*rsoilN(i,1) ;
Nacq(i,j)=min(Ndeficit(i,j),(soilN(i,j)-((iterations-(soildraw(i,j)))/iterations)*rsoilN(i,1)));
Cacq(i,j)=Nacq(i,j)*cost_acq(i,j);
Cgrowth(i,j)=Cgrowth(i,j)-Cacq(i,j);
final_soilN(i,j)=soilN(i,j)-Nacq(i,j);
final_leafN(i,j)=leafN(i,j);
else
if uptake_strategy(i,j)==1 && Nacq(i,j)>leafN(i,j)-((iterations-(leafdraw(i,j)))/iterations)*leafN(i,1) ;
Nacq(i,j)=min(Ndeficit(i,j),(leafN(i,j)-((iterations-(leafdraw(i,j)))/iterations)*leafN(i,1)));
Cacq(i,j)=Nacq(i,j)*cost_acq(i,j);
Cgrowth(i,j)=Cgrowth(i,j)-Cacq(i,j);
final_leafN(i,j)=leafN(i,j)-Nacq(i,j);
final_soilN(i,j)=soilN(i,j);
end
end
end
%if N demand is not met and either the soilN pool or leafN pool is exhausted switch to next cheapest N source
Ndeficit(i,j)=Ndemand(i,1)-nansum(Nacq(i,1:j));
end
end
%PREALLOCATION OF OUTPUT MATRICES FOR ROUND TWO
%If N demand is not met and either the leafN pool or soilN pool is
%exhausted, then uptake moves towards the next cheapest source until demand
%is met.
%Adjusts NPP to reflect C expended for growth in Round One.
NPP_2=nan(length(NPP),1);
for i=1:length(NPP);
NPP_2(i,1)=NPP(i,1)-nansum(Cacq(i,1:iterations));
end
leafN_2=horzcat(final_leafN(:,end),placeholder(:,2:iterations));
soilN_2=horzcat(final_soilN(:,end),placeholder(:,2:iterations));
soildraw_2=nan(length(NPP),iterations);
soildraw_2(:,1)=1;
leafdraw_2=nan(length(NPP),iterations);
leafdraw_2(:,1)=1;
cost_active_2=nan(length(NPP),iterations);%(kg C/m2/timestep) cost of active N uptake
cost_resorb_2=nan(length(NPP),iterations);%(kg C/m2/timestep) cost of resorption of leaf N
cost_acq_2=nan(length(NPP),iterations);%(kg C/m2/timestep) minimun of cost (1) fixation,
%(2) active, and (3)resorb.
uptake_strategy_2=nan(length(NPP),iterations);%(n=1:3) indicates pathway
Cgrowth_2=nan(length(NPP),iterations);%(kg C/m2/timestep) C available for NPP and growth
Cgrowth_2(:,1)=NPP_2(:,1);
Cacq_2=nan(length(NPP),iterations);%(kg C/m2/timestep) C expended on N uptake
Cacq_2(:,1)=0;
Nacq_2=nan(length(NPP),iterations);%(kg N/m2/timestep) Total N taken up
Nacq_2(:,1)=0;
Cavailable_2=nan(length(NPP),iterations);%C available to expend to growth or N uptake
Ndeficit_2=nan(length(NPP),iterations);%N uptake needed to meet N demand
final_leafN_2=nan(length(NPP),iterations); %N remaining in leaf pool
final_soilN_2=nan(length(NPP),iterations);%N remaining in soil pool
%MODEL CODE FOR ROUND 2 (Note: same as Round One except passive uptake is
%not included as it already occurred)
for i=1:length(NPP)
if Ndeficit(i,iterations)==0 || Ndeficit(i,iterations)<0 || Cgrowth(i,iterations)<0 || Cgrowth(i,iterations)==0;
continue
else
for j=1:iterations
Cavailable_2(i,j)=NPP_2(i)-nansum(Cacq_2(i,1:j));
leafN_2(i,j)=((iterations-(leafdraw_2(i,j)-1))/iterations)*final_leafN(i,end);
cost_resorb_2(i,j)=kR/leafN_2(i,j);
soilN_2(i,j)=((iterations-(soildraw_2(i,j)-1))/iterations)*final_soilN(i,end);
cost_active_2(i,j)=(kN/soilN_2(i,j))*(kC/rootbiomass(i));
%Calculates which N strategy have the lowest cost
if cost_resorb_2(i,j)0;
Nacq_2(i,j)=min(Ndeficit_2(i,j),(Cavailable_2(i,j)-Npassive(i)*plantcn(i))/(cost_acq_2(i,j)+plantcn(i)));
Cacq_2(i,j)=(cost_acq_2(i,j)*Cavailable_2(i,j)-Npassive(i)*plantcn(i)*cost_acq_2(i,j))/(cost_acq_2(i,j)+plantcn(i));
Cgrowth_2(i,j)=(plantcn(i)*Cavailable_2(i,j)+Npassive(i)*plantcn(i)*cost_acq_2(i,j))/(cost_acq_2(i,j)+plantcn(i));
if uptake_strategy_2(i,j)==2
final_leafN_2(i,j)=leafN_2(i,j);
final_soilN_2(i,j)=soilN_2(i,j)-Nacq_2(i,j);
if uptake_strategy_2(i,j)==1
final_soilN_2(i,j)=soilN_2(i,j);
final_leafN_2(i,j)=leafN_2(i,j)-Nacq_2(i,j);
end
end
end
if uptake_strategy_2(i,j)==2 && Nacq_2(i,j)>soilN_2(i,j)-((iterations-(soildraw_2(i,j)))/iterations)*soilN_2(i,1) ;
Nacq_2(i,j)=min(Ndeficit_2(i,j),(soilN_2(i,j)-((iterations-(soildraw_2(i,j)))/iterations)*soilN_2(i,1)));
Cacq_2(i,j)=Nacq_2(i,j)*cost_acq_2(i,j);
Cgrowth_2(i,j)=Cgrowth_2(i,j)-Cacq_2(i,j);
final_soilN_2(i,j)=soilN_2(i,j)-Nacq_2(i,j);
final_leafN_2(i,j)=leafN_2(i,j);
else
if uptake_strategy_2(i,j)==1 && Nacq_2(i,j)>leafN_2(i,j)-((iterations-(leafdraw_2(i,j)))/iterations)*leafN_2(i,1) ;
Nacq_2(i,j)=min(Ndeficit_2(i,j),(leafN_2(i,j)-((iterations-(leafdraw_2(i,j)))/iterations)*leafN_2(i,1)));
Cacq_2(i,j)=Nacq_2(i,j)*cost_acq_2(i,j);
Cgrowth_2(i,j)=Cgrowth_2(i,j)-Cacq_2(i,j);
final_leafN_2(i,j)=leafN_2(i,j)-Nacq_2(i,j);
final_soilN_2(i,j)=soilN_2(i,j);
end
end
end
Ndeficit_2(i,j)=Ndemand(i,1)-(nansum(Nacq(i,1:end))+nansum(Nacq_2(i,1:j)));
end
end
end
%PREALLOCATION OF OUTPUT MATRICES FOR ROUND THREE
%If N demand is not met and either the leafN pool or soilN pool is
%exhausted, then uptake moves towards the next cheapest source until demand
%is met.
%Adjusts NPP to reflect C expended for growth in Round One and Two.
NPP_3=nan(length(NPP),1);
for i=1:length(NPP);
NPP_3(i,1)=NPP_2(i,1)-nansum(Cacq_2(i,1:iterations));
end
leafN_3=horzcat(final_leafN(:,end),placeholder(:,2:iterations));
soilN_3=horzcat(final_soilN(:,end),placeholder(:,2:iterations));
soildraw_3=nan(length(NPP),iterations);
soildraw_3(:,1)=1;
leafdraw_3=nan(length(NPP),iterations);
leafdraw_3(:,1)=1;
cost_active_3=nan(length(NPP),iterations);%(kg C/m2/timestep) cost of active N uptake
cost_resorb_3=nan(length(NPP),iterations);%(kg C/m2/timestep) cost of resorption of leaf N
cost_acq_3=nan(length(NPP),iterations);%(kg C/m2/timestep) minimun of cost (1) fixation,
%(2) active, and (3)resorb.
uptake_strategy_3=nan(length(NPP),iterations);%(n=1:3) indicates pathway
Cgrowth_3=nan(length(NPP),iterations);%(kg C/m2/timestep) C available for NPP and growth
Cgrowth_3(:,1)=Cgrowth_2(:,end);
Cacq_3=nan(length(NPP),iterations);%(kg C/m2/timestep) C expended on N uptake
Cacq_3(:,1)=0;
Nacq_3=nan(length(NPP),iterations);%(kg N/m2/timestep) Total N taken up
Nacq_3(:,1)=0;
Cavailable_3=nan(length(NPP),iterations);%C available to expend to growth or N uptake
Ndeficit_3=nan(length(NPP),iterations);%N uptake needed to meet N demand
final_leafN_3=nan(length(NPP),iterations); %N remaining in leaf pool
final_soilN_3=nan(length(NPP),iterations);%N remaining in soil pool
%MODEL CODE FOR ROUND 3 (Note: same as Round One except passive uptake is
%not included as it already occurred)
for i=1:length(NPP)
if Ndeficit_2(i,iterations)==0 || Ndeficit_2(i,iterations)<0 || Cgrowth_2(i,iterations)<0 || Cgrowth_2(i,iterations)==0;
continue
else
for j=1:iterations
Cavailable_3(i,j)=NPP_3(i)-nansum(Cacq_3(i,1:j));
leafN_3(i,j)=((iterations-(leafdraw_3(i,j)-1))/iterations)*final_leafN_2(i,end);
cost_resorb_3(i,j)=kR/leafN_3(i,j);
soilN_3(i,j)=((iterations-(soildraw_3(i,j)-1))/iterations)*final_soilN_2(i,end);
cost_active_3(i,j)=(kN/soilN_3(i,j))*(kC/rootbiomass(i));
%Calculates which N strategy have the lowest cost
if cost_resorb_3(i,j)0;
Nacq_3(i,j)=min(Ndeficit_3(i,j),(Cavailable_3(i,j)-Npassive(i)*plantcn(i))/(cost_acq_3(i,j)+plantcn(i)));
Cacq_3(i,j)=(cost_acq_3(i,j)*Cavailable_3(i,j)-Npassive(i)*plantcn(i)*cost_acq_3(i,j))/(cost_acq_3(i,j)+plantcn(i));
Cgrowth_3(i,j)=(plantcn(i)*Cavailable_3(i,j)+Npassive(i)*plantcn(i)*cost_acq_3(i,j))/(cost_acq_3(i,j)+plantcn(i));
if uptake_strategy_2(i,j)==2
final_leafN_3(i,j)=leafN_3(i,j);
final_soilN_3(i,j)=soilN_3(i,j)-Nacq_3(i,j);
if uptake_strategy_3(i,j)==1
final_soilN_3(i,j)=soilN_3(i,j);
final_leafN_3(i,j)=leafN_3(i,j)-Nacq_3(i,j);
end
end
end
if uptake_strategy_3(i,j)==2 && Nacq_3(i,j)>soilN_3(i,j)-((iterations-(soildraw_3(i,j)))/iterations)*soilN_3(i,1) ;
Nacq_3(i,j)=min(Ndeficit_3(i,j),(soilN_3(i,j)-((iterations-(soildraw_3(i,j)))/iterations)*soilN_3(i,1)));
Cacq_3(i,j)=Nacq_3(i,j)*cost_acq_3(i,j);
Cgrowth_3(i,j)=Cgrowth_3(i,j)-Cacq_3(i,j);
final_soilN_3(i,j)=soilN_3(i,j)-Nacq_3(i,j);
final_leafN_3(i,j)=leafN_2(i,j);
else
if uptake_strategy_3(i,j)==1 && Nacq_3(i,j)>leafN_3(i,j)-((iterations-(leafdraw_3(i,j)))/iterations)*leafN_3(i,1) ;
Nacq_3(i,j)=min(Ndeficit_3(i,j),(leafN_3(i,j)-((iterations-(leafdraw_3(i,j)))/iterations)*leafN_3(i,1)));
Cacq_3(i,j)=Nacq_3(i,j)*cost_acq_3(i,j);
Cgrowth_2(i,j)=Cgrowth_2(i,j)-Cacq_3(i,j);
final_leafN_3(i,j)=leafN_3(i,j)-Nacq_3(i,j);
final_soilN_3(i,j)=soilN_3(i,j);
end
end
end
Ndeficit_3(i,j)=Ndemand(i,1)-(nansum(Nacq(i,1:iterations))+nansum(Nacq_2(i,1:iterations))+nansum(Nacq_3(i,1:j)));
end
end
end
%OUTPUT SUMMARY
%Preallocate output matrices
total_Nacq=nan(length(NPP),1); %Total N uptake from all strategies
total_Cacq=nan(length(NPP),1); %Total C expended for N uptake
total_Cgrowth=nan(length(NPP),1); %Total C available for growth
total_Nacq_active=nan(length(NPP),1); %Total N uptake using active strategy
total_Nacq_resorb=nan(length(NPP),1); %Total N uptake using resorb strategy
total_Nacq_fix=nan(length(NPP),1); %Total N uptake using fixation strategy
total_Cacq_active=nan(length(NPP),1); %Total C expended for active uptake
total_Cacq_resorb=nan(length(NPP),1); %Total C expended for resorption
total_Cacq_fix=nan(length(NPP),1); %Total C expended for fixation
final_N_deficit=nan(length(NPP),1); %Final N deficit
leaf_N_remaining=nan(length(NPP),1); %N remaining in leafN pool
soil_N_remaining=nan(length(NPP),1); %N remaining in soilN pool
%Loop below calculates output summary for each timestep and number of iterations
for i=1:length(NPP)
total_Nacq(i,1)=nansum(Nacq(i,:))+nansum(Nacq_2(i,:))+nansum(Nacq_3(i,:));
total_Cacq(i,1)=nansum(Cacq(i,:))+nansum(Cacq_2(i,:))+nansum(Cacq_3(i,:));
total_Cgrowth(i,1)=Cgrowth_3(i,end);
final_N_deficit(i,1)=nanmin(nanmin(Ndeficit(i),Ndeficit_2(i)),Ndeficit_3(i));
leaf_N_remaining(i,1)=leafN(i,1)-total_Nacq_resorb(i,1);
soil_N_remaining(i,1)=soilN(i,1)-total_Nacq_active(i,1);
rindx=find(uptake_strategy(i,:)==1);
rindx2=find(uptake_strategy_2(i,:)==1);
rindx3=find(uptake_strategy_3(i,:)==1);
total_Nacq_resorb(i,1)=nansum(Nacq(i,rindx))+nansum(Nacq_2(i,rindx2))+nansum(Nacq_3(i,rindx3));
total_Cacq_resorb(i,1)=nansum(Cacq(i,rindx))+nansum(Cacq_2(i,rindx2))+nansum(Cacq_3(i,rindx3));
aindx=find(uptake_strategy(i,:)==2);
aindx2=find(uptake_strategy_2(i,:)==2);
aindx3=find(uptake_strategy_3(i,:)==2);
total_Nacq_active(i,1)=nansum(Nacq(i,aindx))+nansum(Nacq_2(i,aindx2))+nansum(Nacq_3(i,aindx3));
total_Cacq_active(i,1)=nansum(Cacq(i,aindx))+nansum(Cacq_2(i,aindx2))+nansum(Cacq_3(i,aindx3));
findx=find(uptake_strategy(i,:)==3);
findx2=find(uptake_strategy_2(i,:)==3);
findx3=find(uptake_strategy_3(i,:)==3);
total_Nacq_fix(i,1)=nansum(Nacq(i,findx))+nansum(Nacq_2(i,findx2))+nansum(Nacq_3(i,findx3));
total_Cacq_fix(i,1)=nansum(Cacq(i,findx))+nansum(Cacq_2(i,findx2))+nansum(Cacq_3(i,findx3));
leaf_N_remaining(i,1)=leafN(i,1)-total_Nacq_resorb(i,1);
soil_N_remaining(i,1)=soilN(i,1)-total_Nacq_active(i,1);
end
%Clear intermediate variables
clear rindx rindx2 rindx3 aindx aindx2 aindx3 rindx3 findx findx2 findx3
clear Cacq Cacq_2 Cacq_3 Cavailable Cavailable_2 Cavailable_3 Cgrowth Cgrowth_2 Cgrowth_3
clear NPP_2 NPP_3 Ndeficit Ndeficit_2 Ndeficit_3 a_FIX b_FIX c_FIX
clear cost_acq cost_active cost_resorb cost_fixN final_leafN final_soilN
clear cost_acq_2 cost_active_2 cost_resorb_2 cost_fixN_2 final_leafN_2 final_soilN_2
clear cost_acq_3 cost_active_3 cost_resorb_3 cost_fixN_3 final_leafN_3 final_soilN_3
clear i j leafN_2 leafN_3 soilN_2 soilN_3 leafdraw leafdraw_2 leafdraw_3
clear soildraw soildraw_2 soildraw_3 placeholder s_FIX uptake_strategy uptake_strategy_2 uptake_strategy_3
clear kC kR kp kN iterations rsoilN Nacq Nacq_2 Nacq_3