%% Computation of runoff using SCS-CN method
% MATLAB version: R2020a

%% Inputs
% global curve number dataset netCDF files: - "dry", "normal", "wet"
% global precipitation dataset annual daily series

%% Outputs
% direct runoff using SCS-CN method 
% mean CN during flood conditions

%% functions
% CN_unassigned.m
% rain_gpcc.m
% rain_sm2rain (optional)
% rain_sm2rain_syn.m (optional)
% CN_assign.m
% Q_compute.m 


%% clear and close commands
clc,clear,close all 

%% user defined variables - setup stage
%write run reference 
ref = "extent_run1_"; 

% start and end year of precipitation dataset
year_start = double(2020);
year_end = double(2020);

%% Extent section - extents limits, units = degrees
% Extent "ext_AN_1" americas (north) 
%lat_s = 43; lat_n = 53;
%lon_w = -122; lon_e = -109;

% Extent "ext_AN_2" americas (north)
%lat_s = 34; lat_n = 37;
%lon_w = -107; lon_e = -103;

% Extent "ext_AS_1" americas (south)
%lat_s = -9; lat_n = -3;
%lon_w = -41; lon_e = -37;

% Extent ref. "ext_ME_1" middle east 
%lat_s = 32; lat_n = 35;
%lon_w = 34; lon_e = 37;

% Extent "ext_SA_1" southern africa
%lat_s = -35; lat_n = -31;
%lon_w = 18; lon_e = 22;

% Extent "ext_SA_2" southern africa
%lat_s = -29; lat_n = -23;
%lon_w = 24; lon_e = 32;

% Extent "ext_SA_3" southern africa
%lat_s = -34; lat_n = -28;
%lon_w = 24; lon_e = 32;

% Extent "ext_SA_4" southern africa 
%lat_s = -33; lat_n = -28;
%lon_w = 22; lon_e = 27;

% Extent "ext_AUS_1" australasia
% lat_s = -36; lat_n = -32;
% lon_w = 116; lon_e = 119;

% Extent "ext_AUS_2" australasia 
% lat_s = -24; lat_n = -20;
% lon_w = 117; lon_e = 121;

% Extent "ext_AUS_3" australasia
% lat_s = -21; lat_n = -17;
% lon_w = 137; lon_e = 146;

% Extent "ext_test" australasia
lat_s = -19; lat_n = -17;
lon_w = 137; lon_e = 139;

%% location of input datasets                                           

% folder path of global curve number dataset
path_GCN = ('D:\documents\PhD\OneDrive - Lancaster University\GCN250\');
% location of netCDF files holding global curve numbers
GCN_I   = ([path_GCN,'GCN250_I.nc']); % dry condition
GCN_II  = ([path_GCN,'GCN250_II.nc']); % normal condition
GCN_III = ([path_GCN,'GCN250_III.nc']); % wet condition 

% folder path of precipitation dataset 
path_rain = ('D:\documents\PhD\OneDrive - Lancaster University\GPCC_full_data_daily_product_2022\');
% precipitation netCDF file minus year and extension
rain_name = ('full_data_daily_v2022_10_');

% folder path of output files 
path_netCDF = ('D:\documents\PhD\OneDrive - Lancaster University\netCDF\test\'); 

% end of user set-up 

%% create counter - used to place results in output file
n = 1; 

%% create table to hold results 
fid = append(ref,'outputs.txt');
fopen(fid, 'a+');
T = ["year", "synthesis F", "synthesis m", "season code", "lambda", "RainfallAnnualExtentTotal", "RunoffAnnualExtentTotal", "PQconversionAnnualTotalExtent"];
writematrix(T, fid, 'Delimiter', ',' );


%% loop to process annual precipitation series 
% start and end year of precipitation dataset
for Yr = year_start:year_end

    Yr % output to update user on progress
% path to precipitation dataset - file name minus year
RainDatasetName = [path_rain, rain_name, num2str(Yr),'.nc'];

%% location of outputs 

CN250_unassigned = ([path_netCDF, 'CN250_unassigned_', num2str(Yr),'.nc']);
CN250_assigned = ([path_netCDF, 'CN250_assigned_', num2str(Yr),'.nc']);

% define location & name of NC file to holding NC file of unadjusted (input) rain data
RainName = ([path_netCDF, 'Rain_',num2str(Yr),'.nc']);

% define location & name of NC file to holding NC file of adjusted (output) rain data
Rain250Name = ([path_netCDF, 'Rain250_',num2str(Yr),'.nc']);


% rainfall synthesis parameters
for F_count = 2:2 % target rainy days total factor
  F = F_count*5; % increase interval step by multiplier
    
for m_count = 3:3 % top tier target rainy days factor
    m = m_count * 20; %increase interval step by multipler
    

% SCS_CN runoff parameters
for lambda_count = 4:4 % lambda, the ratio of initial abstraction (mm) to S (mm)
    lambda = lambda_count / 20; 


%% paths inputs

% SCS-CN runoff equation season code 
% 1 for dormant season; 2 for growing season
for season_code = 1:1

% increase counter by one
n = n +1; 

%% Unassigned Curve Number extent - call function
Yr % output to update user on progress
" processing of 'CN_unassigned.m' function" % output to update user on progress
CN_unassigned(ref, Yr, lon_w, lon_e, lat_s, lat_n, RainDatasetName, path_GCN, GCN_I, GCN_II, GCN_III, path_netCDF, CN250_unassigned, CN250_assigned, RainName, Rain250Name, F, m, lambda, season_code);

%% select rainfall dataset with or without synthesise

%function for synthesising SM2RAIN-ASCAT rainfall dataset
%rain_sm2rain_syn( ref, Yr, lon_w, lon_e, lat_s, lat_n, F, m, RainDatasetName, GCN_I, RainName, Rain250Name );
% function for non-synthesised SM2RAIN-ASCAT rainfall dataset
%rain_sm2rain(ref, Yr, lon_w, lon_e, lat_s, lat_n, RainDatasetName, GCN_I, RainName, Rain250Name);
% function for GPCC rainfall dataset
Yr % output to update user on progress
" processing of 'rain_gpcc.m' function" % output to update user on progress
rain_gpcc(ref, Yr, lon_w, lon_e, lat_s, lat_n, RainDatasetName, GCN_I, RainName, Rain250Name);

%% assign Curve Number - call function
Yr % output to update user on progress
" processing of 'CN_assign.m' function"  % output to update user on progress
CN_assign( ref, Yr, lon_w, lon_e, lat_s, lat_n, Rain250Name, CN250_unassigned, CN250_assigned, season_code);


%% compute Runoff Q - call function
Yr % output to update user on progress
" processing of 'Q_compute.m' function"  % output to update user on progress
[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);



%% add results to a table
T(n,:) = [Yr, F, m, season_code, lambda, RainfallAnnualExtentTotal, RunoffAnnualExtentTotal, PQconversionAnnualTotalExtent];
writematrix(T, fid, 'Delimiter', ',' );

end

end

end

end

end

% write text file with run details 
T  = table({'reference'; 'longitude west'; 'longitude east'; 'latitude south'; 'latitude north'; ...
   'rainfall dataset'; 'gcn1'; 'gcn2'; 'gcn3';...
    'RainName'; 'Rain250Name'; 'CN250 unassigned'; 'CN250 assigned'; 'path netCDF'},...
    [string(ref); string(lon_w); string(lon_e); string(lat_s); string(lat_n);...
    string(RainDatasetName); string(GCN_I); string(GCN_II); string(GCN_III);...
    string(RainName); string(Rain250Name); string(CN250_unassigned); string(CN250_assigned); string(path_netCDF)]);

% write "setup"  details to text file 
writetable(T, append(ref,'setup.txt'), 'Delimiter', ',' ); 

% close files
fclose('all'); 

%% dataset references
% curve number dataset
% Jaafar, H. H., Ahmad, F. A., and El Beyrouthy, N.: GCN250, new global gridded curve numbers for hydrologic modeling and design, Sci. Data, 6, 145, https://doi.org/10.1038/s41597-019-0155-x, 2019.
% precipitation dataset
% Ziese, Markus; Rauthe-Schöch, Armin; Becker, Andreas; Finger, Peter; Rustemeier, Elke; Hänsel, Stephanie; Schneider, U.: GPCC Full Data Daily Version 2022 at 1.0°, https://doi.org/10.5676/DWD_GPCC/FD_D_V2022_100, 2022.


% created 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.