% PLAN to have different moisture scalars for grass and trees

% relative drying rate (rdr) algorithm
rdr = (1 + SMparams(:,4)) ./ ((1 + (SMparams(:,4).* (last_soilm./SMparams(:,2)).^SMparams(:,5))));
rdr(ppt(:,month) > PET) = 1;                % rdr is 1 if there is more PPT than PET

current_ppt = ppt(:,month);
io_freeze = airt(:,month) < 0;              % average temperure below zero
io_melt   = airt(:,month) >= 0;             % average temperure above zero

fid = last_pack + ppt(:,month);             % add this month's ppt to last month's snowpack
pack(io_freeze) = fid(io_freeze);           % snowpack gets last month's snowpack and current ppt
current_ppt(io_freeze) = 0;                 % current ppt not available for plants (snow)

fid = current_ppt + last_pack;              % add last month's snowpack to this month's ppt
current_ppt(io_melt) = fid(io_melt);        % current ppt increases
pack(io_melt) = 0;                          % snowpack is melted

% begin estimating evapotranspiration and calculate new soil moisture

% if PET exceeds PPT, EET is limited
eeta = current_ppt + ((PET - current_ppt) .* rdr);
eetb = current_ppt + last_soilm - SMparams(:,1);

% EET becomes the smallest of EETa and EETb
io = eeta > eetb;  EET(io) = eetb(io);
io = eeta <= eetb; EET(io) = eeta(io);

% if PPT exceeds PET, EET is not limited
io = current_ppt >= PET;
fid = PET;
EET(io) = fid(io);

this_soilm = last_soilm + current_ppt - EET;

io = this_soilm > SMparams(:,2); % if estimated soil moisture exceeds field capacity, runoff will occur
fid = SMparams(:,2); this_soilm(io) = fid(io);

soilm = this_soilm;  

io = PET > 0;
warning off MATLAB:divideByZero
fid = 0.5 + (0.5.*(EET./PET));
NPPmoist_temp(io) = fid(io);

bgratio = (last_soilm - SMparams(:,1) + current_ppt) ./ PET;
warning on MATLAB:divideByZero

% different bgratio values (too little water is limiting, too much also)
% water limited:
io = (bgratio >= 0 & bgratio <1);
fid = 0.1 + (0.9.*bgratio);
bgmoist_temp(io) = fid(io);

% enough water
io = (bgratio >= 1 & bgratio <= 2);
bgmoist_temp(io) = 1;

% too much water
io = (bgratio > 2 & bgratio <= 30);
fid = (1+1/28) - ((0.5/28).*bgratio);
bgmoist_temp(io) = fid(io);

% way too much water
io = bgratio > 30;
bgmoist_temp(io) = 0.5;

% set up moisture factors for NPP calculation and BG run
% in case PET is zero
io = PET <= 0;
NPPmoist_temp(io) = NPPmoistpret(io);
bgmoist_temp(io)  = bgmoistpret(io);

NPPmoist(:,month) = NPPmoist_temp;
bgmoist(:,month)  = bgmoist_temp;

% set snow pack, soil moisture, bgmoistpret, NPPmoistpret for next step
last_pack    = pack;
last_soilm   = soilm;
bgmoistpret  = bgmoist(:,month);
NPPmoistpret = NPPmoist(:,month);   