function [RainfallAnnualExtentTotal, RunoffAnnualExtentTotal, PQconversionAnnualTotalExtent] = Q_compute( ref, Yr, lon_w, lon_e, lat_s, lat_n, F, m, lambda, RainDatasetName, Rain250Name, CN250_assigned, season_code, path_netCDF, n)

% OUTPUT data
% define location & name of NC file to hold annual data
annual = append(string(path_netCDF), string(Yr),'_annual_Part_B_',string(ref),'.nc');


% define location & name of NC file to hold daily data
daily  = append(string(path_netCDF), string(Yr),'_daily_',string(ref),'.nc');

D =int16(datenum(['1-jan-',num2str(Yr)]):datenum(['31-dec-',num2str(Yr)])); % create vector of data range

L = (1:yeardays(Yr));

% data range, years
%range = (Yr_s:Yr_e)';
%index_yr = int8(find(range == Yr));


%% clip lon and lat of rain data
% load longitude array; Global Curve Number (gcn)
lon = (ncread(Rain250Name,'longitude')); 
% load latitude array; Global Curve Number (gcn)
lat = (ncread(Rain250Name,'latitude'));

%% Output NC file, create variables and write attributes
% create variable to hold annual sum of daily runoff (mm)
nccreate(annual, ['Runoff_Annual_',num2str(Yr)], 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)}, 'DataType', 'int16');
ncwriteatt(annual, ['Runoff_Annual_',num2str(Yr)], "units", ("mm/year"));

% create variable to hold annual sum of daily rainfall (mm)
nccreate(annual, ['Rainfall_Annual_',num2str(Yr)], 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)}, 'DataType', 'int16');
ncwriteatt(annual, ['Rainfall_Annual_',num2str(Yr)], "units", ("mm/year"));

% create variable to hold annual number of runoff days (days)
nccreate(annual, ['Q_days_',num2str(Yr)], 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)}, 'DataType', 'int16');
ncwriteatt(annual, ['Q_days_',num2str(Yr)], "units", ("days/year"));

% create variable to hold mean of CN values (-) on days runoff occurs
nccreate(annual, ['CN_mean_annual_Q_days_',num2str(Yr)], 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)},'DataType', 'int16');
ncwriteatt(annual, ['CN_mean_annual_Q_days_',num2str(Yr)], "units", ("-")); 

% create variable to hold mean of CN values (-) on days runoff occurs
nccreate(annual, ['CN_mean_annual_',num2str(Yr)], 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)},'DataType', 'int16');
ncwriteatt(annual, ['CN_mean_annual_',num2str(Yr)], "units", ("-")); 

% create variable to hold mean of CN values (-) on days runoff occurs 
% (void filled on no Q days with mean annual CN values)
nccreate(annual, ['mean_annual_CN_value_Q_days_void_filled_',num2str(Yr)], 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)},'DataType', 'int16');
ncwriteatt(annual, ['mean_annual_CN_value_Q_days_void_filled_',num2str(Yr)], "units", ("-")); 


% create variable to hold conversion of rainfal to runoff at cell level
nccreate(annual, ['P_Q_conversion_',num2str(Yr)], 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)},'DataType', 'int16');
ncwriteatt(annual, ['P_Q_conversion_',num2str(Yr)], "units", ("%")); 


% create dimensions
nccreate(annual, 'longitude', 'Dimensions',{'longitude', length(lon)}); % longitude, needs DOUBLE
nccreate(annual, 'latitude', 'Dimensions',{'latitude', length(lat)}); % latitude

% Output NC file, write dimensions
ncwrite(annual, 'longitude', lon); % write longitude dimension
ncwrite(annual, 'latitude', lat); % write latitude dimension


%% clip lon and lat of rain data
%load longitude array; Global Curve Number (gcn)
lon_r = (ncread(Rain250Name,'longitude')); % needs to be DOUBLE
% load latitude array; Global Curve Number (gcn)
lat_r = (ncread(Rain250Name,'latitude'));
% determine which indices correspond to dimension 1: Global Curve Number (gcn)
ind1_r = int16(find(lon_r>= lon_w & lon_r <= lon_e)); % integer, max value 20,000
% determine which indices correspond to dimension 2: Global Curve Number (gcn)
ind2_r = int16(find(lat_r>= lat_s & lat_r <= lat_n));
% Clip lat and lon to their specified range
lat_r = lat_r(ind2_r); % needs to be DOUBLE
lon_r = lon_r(ind1_r); % needs to be DOUBLE


%% clip lon and lat of Curve Number data
%load longitude array; Global Curve Number (gcn)
lon_cn = ncread(CN250_assigned,'longitude'); % needs to be DOUBLE
% load latitude array; Global Curve Number (gcn)
lat_cn = ncread(CN250_assigned,'latitude'); % needs to be DOUBLE
% determine which indices correspond to dimension 1: Global Curve Number (gcn)
ind1_cn = int16(find(lon_cn>= lon_w & lon_cn <= lon_e));
% determine which indices correspond to dimension 2: Global Curve Number (gcn)
ind2_cn = int16(find(lat_cn>= lat_s & lat_cn <= lat_n));
% Clip lat and lon to their specified range
lat_cn = lat_cn(ind2_cn);
lon_cn = lon_cn(ind1_cn);

%% Output NC file, create variables and write attributes

% create variable to hold P daily values (mm/day)
nccreate(daily, 'P', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time', (Inf)},'DataType', 'int16');
ncwriteatt(daily, 'P', "units", ("mm/day")); 

% create variable to hold Q daily values (mm/day)
% holds total Q without cut-off applied
nccreate(daily, 'Q_all', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time', (Inf)},'DataType', 'int16');
ncwriteatt(daily, 'Q_all', "units", ("mm/day")); 

% create variable to hold runoff, daily values, mm/day
nccreate(daily, 'Q', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time', (Inf)},'DataType', 'int16');
ncwriteatt(daily, 'Q', "units", ("mm/day")); 

% create variable to hold daily 'S' values (mm)
nccreate(daily, 'S', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time', (Inf)},'DataType', 'int16');
ncwriteatt(daily, 'S', "units", ("mm")); 

% create variable to hold daily CN values
nccreate(daily, 'CN', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time', (Inf)},'DataType', 'int16');
ncwriteatt(daily, 'CN', "units", ("-")); 

% create variable to hold daily lambda values
nccreate(daily, 'lambda', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time', (Inf)},'DataType', 'int16');
ncwriteatt(daily, 'lambda', "units", ("-")); 

% create variable to hold daily Ia (mm)values
nccreate(daily, 'Ia', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time', (Inf)},'DataType', 'int16');
ncwriteatt(daily, 'Ia', "units", ("mm")); 

% create variable to hold days when P > Ia (or P > lambda x S)
% grt = "greater than"
nccreate(daily, 'P_grt_Ia', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time', (Inf)},'DataType', 'int16');
ncwriteatt(daily, 'P_grt_Ia', "units", ("mm")); 

% create variable to hold CN values for runoff producing days(dimensionless)
nccreate(daily, 'CN_Q_days', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r),'time'},'DataType', 'int16');
ncwriteatt(daily, 'CN_Q_days', "units", ("dimensionless"));

% create variable to hold annual sum of daily runoff (mm)
nccreate(daily, 'AnnualRunoff', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r)},'DataType', 'int16');
ncwriteatt(daily, 'AnnualRunoff', "units", ("mm/year")); 

% create variable to hold annual sum of daily rainfall (mm)
nccreate(daily, 'AnnualRainfall', 'Dimensions',{'longitude',length(lon_r), 'latitude',length(lat_r)},'DataType', 'int16');
ncwriteatt(daily, 'AnnualRainfall', "units", ("mm/year")); 

% create dimensions
nccreate(daily, 'longitude', 'Dimensions',{'longitude', length(lon_r)}); % longitude, needs DOUBLE
nccreate(daily, 'latitude', 'Dimensions',{'latitude', length(lat_r)}); % latitude
nccreate(daily, 'time', 'Dimensions',{'time', Inf}); % create time dimension

% Output NC file, write dimensions
ncwrite(daily, 'longitude', lon_r); % write longitude dimension
ncwrite(daily, 'latitude', lat_r); % write latitude dimension
ncwrite(daily, 'time', int16(L)); % write time

%% 
% lambda array
arraylambda = lambda*ones(length(ind1_cn), length(ind2_cn)); % lambda 

array25400 = int16(25400*ones(length(ind1_cn), length(ind2_cn))); % array size of CN array with value of 25400
array254 = int16(254*ones(length(ind1_cn), length(ind2_cn))); % array size of CN array with value of 254
array1 = int16(ones(length(ind1_cn), length(ind2_cn))); % array size of CN array with value of 1

array_neg999 = int16(ones(length(lon_cn), length(lat_cn))*-999); % 

CN_no_data = int16(ones(length(lon_cn), length(lat_cn))*-127); % 
CN_Q_days = int16(ones(length(lon_cn), length(lat_cn))*-999); % to hold CN values for runoff days only
CN_Q_days_annual = int16(NaN(length(lon_cn), length(lat_cn))); % to hold CN values for runoff days only annual
CN_annual_total = int16(zeros(length(lon_cn), length(lat_cn))); % to hold annual total of CN values
MeanAnnualCNvalueQdaysVoidFilled = int16(zeros(length(lon_cn), length(lat_cn))); % to hold the annual mean CN values on Q days ......
% ...with void cells (no runoff in the year) with the annual mean CN value 
CN_annual_mean_voids = int16(zeros(length(lon_cn), length(lat_cn))); % void values (cells with zero annual runoff)

Q = int16(zeros(length(lon_r), length(lat_r))); % same DATATYPE as create in NC file
Qannual = int16(zeros(length(lon_r), length(lat_r))); % to hold annual sum of daily runoff
Q_zeros = int16(zeros(length(lon_r), length(lat_r))); %

Q_days_count = int16(zeros(length(lon_cn), length(lat_cn))); % to hold number of runoff days
Q_days_total = int16(zeros(length(lon_cn), length(lat_cn))); % to hold number of runoff days in year
Q_days_1day = int16(ones(length(lon_cn), length(lat_cn))); % array to count 1 day of runoff
Qsum_daily = int16(NaN(365,1));
P = int16(zeros(length(lon_r), length(lat_r)));
Pannual = int16(zeros(length(lon_r), length(lat_r))); % same DATATYPE as create in NC file

% i = 1:yeardays(Yr)
for i=1:yeardays(Yr)
    
    
    
    %% Precipitation
    
    P = int16(ncread(Rain250Name, 'Rain250', [ind1_r(1), ind2_r(1), i], [length(ind1_r), length(ind2_r), 1]));
    
   %ncwrite(daily, 'P', P, [1 1 i]); % write daily P to NC file 
    
    Pannual = Pannual + P;
    
    %% Curve Number
    
    CN = int16(ncread(CN250_assigned, 'CurveNumber', [ind1_cn(1), ind2_cn(1), i], [length(ind1_cn), length(ind2_cn), 1])); % read
    CN(CN==127) = -999; % set no data
    
    %ncwrite(daily, 'CN', CN, [1 1 i]); % write daily CN NC file    
    
    CN_annual_total = CN_annual_total + CN; % tally CN annual total
     
    %% Mask
    
    % create mask based on Curve Number no data cells
    
    mask = int32(find(CN == -999));
    
    
    %% S
    
    S_i = array25400./CN; % first stage of computing S
    S_i(S_i<0) = -999;
    S = S_i - array254; % compute S, in mm, ref Equation 4b (Soulis and Valiantzas, 2012)  % assume lowest CN of 15 then Smax is 1,439
    S(S<0) = -999;
    %ncwrite(daily, 'S', S, [1 1 i]); % write daily values of 'S' to NC file
    
       
    %% Initial abstraction 
    % Ia (mm) = lambda x S
        
    S = double(S);                          % change S from integer to double
        
    lambdaS = int16(arraylambda.*S);        % calculate Initial Abstraction
    lambdaS(lambdaS < 0) = -999;            % set no data 
    
    %ncwrite(daily, 'Ia', lambdaS, [1 1 i]); % write daily Ia to NC file
    
    S = int16(S);                           % change S back to integer
     
    %% P > Ia
    
    P_grt_Ia = Q_zeros;
    P_grt_Ia(P > lambdaS) = 1; % lambdaS = Ia
    P_grt_Ia(mask) = -999;
    
    %ncwrite(daily, 'P_grt_Ia', P_grt_Ia, [1 1 i]); % write daily (P > Ia) to NC file
    
    %% Q 'all' runoff 
    
    Q_all_num = (P - lambdaS).^2;   % compute numerator of runoff equation Q
    
    Q_all_den = (P - lambdaS + S);  % compute denominator of runoff equation Q
    
    % below is not all actual runoff as Q only occurs if P > Ia
    Q_all = Q_all_num./Q_all_den;   % compute runoff equation
    
    Q_all(mask) = -999; % set mask cells to no data
    
    %ncwrite(daily, 'Q_all', Q_all, [1 1 i]); % write daily total Q values, no cut-off
    
    %% Q runoff
    
    % find cells precipitation > than initial abstraction (lambda x S)
    Q_yes = int32(find(P > lambdaS)); 
    
    Q = Q_zeros;
    
    Q(Q_yes) = Q_all(Q_yes); % Q
    
    Q(mask) = -999; 
    
    % use line below is daily record of variable is required
    %ncwrite(daily, 'Q', Q, [1 1 i]); % write daily runoff Q to NC file
    
    % annual runoff 
    Qannual = Qannual + Q; % record update annual sum of runoff
    Qannual(mask) = -999;    % mask no data cells based on CN data
       
    %% CN for flow days/cells
    
    CN_no_Q = int32(find(Q <= 0));                     % find cells no Q and nodata cells
    
    CN_Q_days = CN;
    
    CN_Q_days(CN_no_Q) = NaN;                          % populate CN values for cells with runoff only
    
    CN_Q_days_annual = CN_Q_days_annual + CN_Q_days;   % add CN on Q days to tally
     
    %ncwrite(daily, 'CN_Q_days', CN_Q_days, [1 1 i]);  % write CN on runoff days to NC file
    
    
    %% count Q days 
        
    Q_count = Q_zeros;
    
    Q_count(Q>0) = 1;
    
    Q_days_total = Q_days_total + Q_count;
       
    %% count "rainy days" 
    
    %rainy_days_count = Q_zeros;
    %rainy_days_count(P>=1) = 1;
    
    
end

%% write annual runoff to NC file
Qannual(Qannual<0) = nan;


ncwrite(annual, ['Runoff_Annual_',num2str(Yr)], Qannual, [1 1]); % write annual runoff to NC file

RunoffAnnualExtentTotal = sum(Qannual, 'all', 'omitnan'); % single value for entire extent


%% write annual rainfall to NC file

Pannual(Pannual<0) = nan;

ncwrite(annual, ['Rainfall_Annual_',num2str(Yr)], Pannual, [1 1]); % write annual runoff to NC file

RainfallAnnualExtentTotal = sum(Pannual, 'all', 'omitnan'); % single value for entire extent

%% find cell rainfall-runoff conversion, and write to NC file

PQconversion = int16((double(Qannual) ./ double(Pannual))*100);
ncwrite(annual, ['P_Q_conversion_',num2str(Yr)], PQconversion); 

%% Total Extent Conversion of Rainfall to Runoff (%)
PQconversionAnnualTotalExtent = (RunoffAnnualExtentTotal / RainfallAnnualExtentTotal)*100; 

%% write total annual number of runoff days

Q_days_total(mask) = -999; % mask water bodies

ncwrite(annual, ['Q_days_',num2str(Yr)], Q_days_total, [1 1]); % write to NC file


%% calculate annual mean CN (all days, runoff and no runoff days)

CN_annual_total(mask) = -999; % mask area of water bodies

CN_mean_annual = CN_annual_total ./ yeardays(Yr); % find mean

CN_mean_annual(mask)= -999; % apply mask water bodies 

% write mean annual runoff CN values (all days)
ncwrite(annual, ['CN_mean_annual_',num2str(Yr)], CN_mean_annual, [1 1]);

%% calculate annual mean CN only on runoff days

CN_Q_days_mean_annual = CN_Q_days_annual ./ Q_days_total;

CN_Q_days_mean_annual(mask) = -999; % mask water bodies 

% write mean annual runoff CN values for runoff days
% this array contains voids (cells where CN is zero as no runff. Q, occurs
% throughout the entire year)
ncwrite(annual, ['CN_mean_annual_Q_days_',num2str(Yr)], CN_Q_days_mean_annual, [1 1]); 

%% create CN annual mean Q days (void filled)
% for cells with no CN value (no Q in year) fill with mean CN annual value
% find cells with no CN value
mask_zero_CN = int32(find(CN_Q_days_mean_annual == 0));

CN_annual_mean_voids(mask_zero_CN) = CN_mean_annual(mask_zero_CN);

MeanAnnualCNvalueQdaysVoidFilled  = CN_Q_days_mean_annual + CN_annual_mean_voids;

ncwrite(annual, ['mean_annual_CN_value_Q_days_void_filled_',num2str(Yr)], MeanAnnualCNvalueQdaysVoidFilled, [1 1]); 


%% reference for Q calculation
% Cronshey, R.G., Roberts, R.T., Miller, N., 1985. Urban Hydrology for Small Watersheds (Tr-55 Rev.). US Department of Agriculture, Soil Conservation Service, Engineering Division.Cronshey, R.G., Roberts, R.T., Miller, N., 1985. Urban Hydrology for Small Watersheds (Tr-55 Rev.). US Department of Agriculture, Soil Conservation Service, Engineering Division.
% Soulis, K.X., Valiantzas, J.D., 2012. SCS-CN parameter determination using rainfall-runoff data in heterogeneous watersheds - the two-CN system approach. Hydrol. earth Syst. Sci. https://doi.org/10.5194/hess-16-1001-2012
% Silveira, L., Charbonnier, F., & Genta, J. L. (2000). L’humidité antérieure des sols dans la méthode “curve number.” 
% Hydrological Sciences Journal. https://doi.org/10.1080/02626660009492302



end

% created by robert g. delaney (2023)
% This work is licensed under the Creative Commons Attribution-
% NonCommercial-ShareAlike 4.0 International License. 
% To view a copy of this license, visit 
% http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a 
% letter to Creative Commons, PO Box 1866, 
% Mountain View, CA 94042, USA.