%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)<cost_active(i,j) && cost_resorb(i,j)< cost_fixN(i);
                 cost_acq(i,j)=cost_resorb(i,j);
                 uptake_strategy(i,j)=1;
                 else  
                 if cost_active(i,j)<cost_resorb(i,j) && cost_active(i,j)< cost_fixN(i);
                 cost_acq(i,j)=cost_active(i,j);
                 uptake_strategy(i,j)=2;
                 else
                 if cost_fixN(i)<cost_resorb(i,j) && cost_fixN(i)< cost_active(i,j);
                 cost_acq(i,j)=cost_fixN(i);
                 uptake_strategy(i,j)=3;
                 end
                 end
                 end
                 if uptake_strategy(i,j)==3;
                    soildraw(i,j+1)=soildraw(i,j);
                    leafdraw(i,j+1)=leafdraw(i,j);
                 else
                 if uptake_strategy(i,j)==1;
                     soildraw(i,j+1)=soildraw(i,j);
                     leafdraw(i,j+1)=leafdraw(i,j)+1;
                 else   
                 if uptake_strategy(i,j)==2;
                     leafdraw(i,j+1)=leafdraw(i,j);
                     soildraw(i,j+1)=soildraw(i,j)+1;
                 else    
                 end
                 end
                 end
        %This code below simultaneously solves eq 6a-6d from the Fisher et al. 
        %(2010) to optimize C allocated to growth and N uptake. N uptake is
        %truncated if Ndemand is met or leafN pool or soilN pool is
        %exhausted.  When using the iterative cost functions N uptake is
        %truncated when the portion of the leafN or soilN pool in that
        %iteration is exhausted.
            Ndeficit(i,j)=Ndemand(i,1)-nansum(Nacq(i,1:j));
            if Ndeficit(i,j)<0 || Ndeficit(i,j)==0;       
            Nacq(i,j)=0;
            Cacq(i,j)=Nacq(i,j)*cost_acq(i,j);
            Cgrowth(i,j)=Cgrowth(i,j)-Cacq(i,j);
            else  
            if Ndeficit(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)<cost_active_2(i,j) && cost_resorb_2(i,j)< cost_fixN(i);
                 cost_acq_2(i,j)=cost_resorb(i,j);
                 uptake_strategy_2(i,j)=1;
                 else  
                 if cost_active_2(i,j)<cost_resorb_2(i,j) && cost_active_2(i,j)< cost_fixN(i);
                 cost_acq_2(i,j)=cost_active_2(i,j);
                 uptake_strategy_2(i,j)=2;
                 else
                 if cost_fixN(i)<cost_resorb_2(i,j) && cost_fixN(i)< cost_active_2(i,j);
                 cost_acq_2(i,j)=cost_fixN(i);
                 uptake_strategy_2(i,j)=3;
                 end
                 end
                 end
                 if uptake_strategy_2(i,j)==3;
                    soildraw_2(i,j+1)=soildraw_2(i,j);
                    leafdraw_2(i,j+1)=leafdraw_2(i,j);
                 else
                 if uptake_strategy_2(i,j)==1;
                     soildraw_2(i,j+1)=soildraw_2(i,j);
                     leafdraw_2(i,j+1)=leafdraw_2(i,j)+1;
                 else   
                 if uptake_strategy_2(i,j)==2;
                     leafdraw_2(i,j+1)=leafdraw_2(i,j);
                     soildraw_2(i,j+1)=soildraw_2(i,j)+1;
                 else    
                 end
                 end
                 end
        %This code below simultaneously solves eq 6a-6d from the Fisher et al. 
        %(2010) to optimize C allocated to growth and N uptake. 
            Ndeficit_2(i,j)=Ndemand(i,1)-(nansum(Nacq(i,1:iterations))+nansum(Nacq_2(i,1:j)));
            if Ndeficit_2(i,j)<0 || Ndeficit_2(i,j)==0;       
            Nacq_2(i,j)=0;
            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);
            else  
            if Ndeficit_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)<cost_active_3(i,j) && cost_resorb_3(i,j)< cost_fixN(i);
                 cost_acq_3(i,j)=cost_resorb(i,j);
                 uptake_strategy_3(i,j)=1;
                 else  
                 if cost_active_3(i,j)<cost_resorb_3(i,j) && cost_active_3(i,j)< cost_fixN(i);
                 cost_acq_3(i,j)=cost_active_3(i,j);
                 uptake_strategy_3(i,j)=2;
                 else
                 if cost_fixN(i)<cost_resorb_3(i,j) && cost_fixN(i)< cost_active_3(i,j);
                 cost_acq_3(i,j)=cost_fixN(i);
                 uptake_strategy_3(i,j)=3;
                 end
                 end
                 end
                 if uptake_strategy_3(i,j)==3;
                    soildraw_3(i,j+1)=soildraw_3(i,j);
                    leafdraw_3(i,j+1)=leafdraw_3(i,j);
                 else
                 if uptake_strategy_3(i,j)==1;
                     soildraw_3(i,j+1)=soildraw_3(i,j);
                     leafdraw_3(i,j+1)=leafdraw_3(i,j)+1;
                 else   
                 if uptake_strategy_3(i,j)==2;
                     leafdraw_3(i,j+1)=leafdraw_3(i,j);
                     soildraw_3(i,j+1)=soildraw_3(i,j)+1;
                 else    
                 end
                 end
                 end
        %This code below simultaneously solves eq 6a-6d from the Fisher et al. 
        %(2010) to optimize C allocated to growth and N uptake. 
            Ndeficit_3(i,j)=Ndemand(i,1)-(nansum(Nacq(i,1:iterations))+nansum(Nacq_2(i,1:iterations))+nansum(Nacq_3(i,1:j)));
            if Ndeficit_3(i,j)<0 || Ndeficit_3(i,j)==0;       
            Nacq_3(i,j)=0;
            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);
            else  
            if Ndeficit_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
       
        


