function nGR = aGR3(T,P,Du, parms)
% aGR3 is a model for max algae growth rate, depending on environmental input
%   (these can be vectors or matrices, but must be the same size).  The
%   thermal response in growth is a Briere curve.  A type II nutrient
%   uptake is assumed, and a death rate which grows with temperature as
%   well as offset from an optimal pH
%
%   This is a cleaned-up and commented version of the original code,
%   written by James Powell, USU, in May 2023
%
%   Input variables:
%       T    = temperature, degrees C
%       P    = pH
%       Du   = nutrient concentration (dust)
%
%   Parameters controlling available nutrient response to pH and T
%   these were fit using log normal variance and a logistic response (see
%   the script Nut_model.m)
%
%        KP = additive correction for low dust (1.405)
%        a = max nutrient content for dust (multiplier) (.9425)
%        mu, sig2 = mean and variance for log-normal nutrient model
%        beta(1) = half sat constant for logisttic (-9.7098)
%        beta(2) = multiplier for pH in exponent (+1.5483)
%        beta(3) = Arrhenius constant 1/(T+273) (54.122)
%
% Parameters for the rate curves are in parms:
%
%   parms(1:4) are b_j controlling thermal rate curve (Briere)
%       parms(1) = b1 = base rate
%       parms(2) = b2 = estimated upper threshold in T
%       parms(3) = b3 = adjustment of upper threshold with available R
%   parms(4:6) are c_j controlling type II nutrient uptake 
%       parms(4) = K = base half-sat value for type II uptake (Thomas)
%       parms(5) = c2 = exp change of half-sat with pH (relative to 7.3)
%       parms(6) = c3 = exp change of half-sat with T (relative to 23)
%   parms(7:10) are d_j controlling death rate from Thomas + pH response
%       parms(7) = d0 = base death rate
%       parms(8) = d1 = proportionality of exp. death rate increase 
%       parms(9) = d2 = exp change with T (relative to 23)
%       parms(10) = d3 = exponent of death response to pH (relative to 7.3)
%
%   NOTE: the roles of c2 and c3, as well as d2 and d3, are switched from
%   what appears in the final published paper!
%
a=.9425; KP=1.405; 
mu=-0.092667; sig2=0.11497; % parameters for log normal correction
beta = [-9.7098, 1.5483, 54.122];
b1=parms(1); b2=parms(2); b3=parms(3); 
K=abs(parms(4)); 
c2=parms(5); c3=parms(6); 
d0=parms(7); d1=parms(8); d2=parms(9); d3=parms(10); 

Ts=T-23;  % for arguments with shifted temperatures
% Model for nutrient assimilation -> growth
    % available nutrient depends on pH:
    R = exp(mu+.5*sig2)*(  KP+(a*Du)./( 1+exp(beta(1)+beta(2)*P-beta(3)./(T+273)) )  ); % includes log normal correction
    % growth rate depends on Type II uptake and temp/pH:
    gR=abs(b1)*max(0,T-5).*T.*sqrt(max(0,1-T./(abs(b2)+b3*R))).^1.*R./(K*exp(c2*(P-7.3)+c3*Ts)+R); 

% Model for death rate
    % death rate with temperature from Thomas + pH
    dRT=abs(d0)+abs(d1)*exp(abs(d2)*Ts-d3*(P-7.3).^2);

nGR=gR-dRT;  % net growth rate = increase - death
    
end

