function [lon] = rain_gpcc( ref, Yr, lon_w, lon_e, lat_s, lat_n, RainDatasetName, GCN_I, RainName, Rain250Name )

RainName = RainName;
Rain250Name = Rain250Name ;

L = (1:yeardays(Yr)); % array of days in the year

% Load longitude array: GPCC "lon";
lon = double(ncread(RainDatasetName,'lon')); 
% determine which indices correspond to dimension 1:
ind1 = find(lon>=lon_w & lon<=lon_e);

% Do the same for lat:
lat = double(ncread(RainDatasetName,'lat'));
ind2 = find(lat>=lat_s & lat<=lat_n);

% Clip lat and lon to their specified range:
lat = lat(ind2);
lon = lon(ind1);

%%  interpolate to new spatial resolution 
 
 %Global Curve Number (gcn) - clip lon and lat to limits; assume ARCI, ARCII, ARCIII have same extents
%load longitude array; Global Curve Number (gcn)
lon_gcn = double(ncread(GCN_I,'lon'));
% load latitude array; Global Curve Number (gcn)
lat_gcn = double(ncread(GCN_I,'lat'));
% determine which indices correspond to dimension 1: Global Curve Number (gcn)
ind1_gcn = find(lon_gcn>= lon_w & lon_gcn <= lon_e);
% determine which indices correspond to dimension 2: Global Curve Number (gcn)
ind2_gcn = find(lat_gcn>= lat_s & lat_gcn <= lat_n);
% Clip lat and lon to their specified range
lat_gcn = lat_gcn(ind2_gcn);
lon_gcn = lon_gcn(ind1_gcn);

%% create NC file - rainfal dataset resolution 

% create variable to hold rainfall at original spatial resolution
nccreate(RainName, 'Rain', 'Dimensions',{'time', (Inf), 'longitude',length(lon), 'latitude',length(lat)});
ncwriteatt(RainName, 'Rain', "units", ("mm/day"));

% create variable to hold rainfall at original spatial resolution
% 5 day rainfall summation 
nccreate(RainName, 'RainSum', 'Dimensions',{'time', (Inf), 'longitude',length(lon), 'latitude',length(lat)});
ncwriteatt(RainName, 'RainSum', "units", ("mm/5-days"));

% create variable to hold synthesised rainfall at original spatial resolution 
% annual total of synthesised rainy days 
nccreate(RainName, 'NumSynRainyDays', 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)});
ncwriteatt(RainName, 'NumSynRainyDays', "units", ("days/year"));

nccreate(RainName, 'NumUnsynRainyDays', 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)});
ncwriteatt(RainName, 'NumUnsynRainyDays', "units", ("days/year"));

nccreate(RainName, 'RainDay240', 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)});
ncwriteatt(RainName, 'RainDay240', "units", ("mm/day"));

% create variable to hold synthesised rainfall at original spatial resolution 
% annual rainfall
nccreate(RainName, 'AnnualSynRain', 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)});
ncwriteatt(RainName, 'AnnualSynRain', "units", ("mm/year"));

% create variable to hold synthesised rainfall at original spatial resolution 
% annual unsynthesised rainfall
nccreate(RainName, 'AnnualUnsynRain', 'Dimensions',{'longitude',length(lon), 'latitude',length(lat)});
ncwriteatt(RainName, 'AnnualUnsynRain', "units", ("mm/year"));

nccreate(RainName, 'longitude', 'Dimensions',{'longitude', length(lon)});
nccreate(RainName, 'latitude', 'Dimensions',{'latitude', length(lat)});
nccreate(RainName, 'time', 'Dimensions',{'time', Inf}); % create time dimension
ncwrite(RainName, 'longitude', lon); % write longitude dimension
ncwrite(RainName, 'latitude', lat); % write latitude dimension
ncwrite(RainName, 'time', int16(L)); % write time

%% NC file for high resolution synthesised data
%create variable to hold adjusted rainfall
nccreate(Rain250Name, 'Rain250', 'Dimensions',{'longitude',length(lon_gcn), 'latitude',length(lat_gcn),'time', (Inf)},'DataType', 'int16');
ncwriteatt(Rain250Name, 'Rain250', "units", ("mm/day"));

% create variable to hold synthesised rainfall at original spatial resolution
nccreate(Rain250Name, 'RainSum250', 'Dimensions',{'longitude',length(lon_gcn), 'latitude',length(lat_gcn),'time', (Inf)},'DataType', 'int16');
ncwriteatt(Rain250Name, 'RainSum250', "units", ("mm/5-day"));

% create variable to hold synthesised rainfall at original spatial resolution
nccreate(Rain250Name, 'RainAnnual250', 'Dimensions',{'longitude',length(lon_gcn), 'latitude',length(lat_gcn)});
ncwriteatt(Rain250Name, 'RainAnnual250', "units", ("mm/year"));


nccreate(Rain250Name, 'longitude', 'Dimensions',{'longitude', length(lon_gcn)});
nccreate(Rain250Name, 'latitude', 'Dimensions',{'latitude', length(lat_gcn)});
nccreate(Rain250Name, 'time', 'Dimensions',{'time', Inf}); % create time dimension

ncwrite(Rain250Name, 'longitude', lon_gcn); % write longitude dimension
ncwrite(Rain250Name, 'latitude', lat_gcn); % write latitude dimension
ncwrite(Rain250Name, 'time', int16(L)); % write time



%% 
%N = maxNumCompThreads;
%LASTN = maxNumCompThreads(2);
% N = maxNumCompThreads;

% read netCDF file holding rainfall data
% GCPP "precip", read "lon", "lat", "time"
start = [ind1(1) ind2(1) 1]; % GPCC 
stride = [length(ind1) length(ind2) yeardays(Yr) ]; %GPCC
    
% GPCC "precip"; read [lon lat time]
P = (ncread(RainDatasetName, 'precip', start, stride)); % lon - lat - time

P(P<0) = nan;

%P = permute(P, [2 3 1]);
    
[Pm,Pn,Pr] = size(P); % obtain size of rain array 

%% analysis of original rainfall dataset, prior to synthesis

NumberUnsynRainyDays3D = zeros(size(P));
NumberUnsynRainyDays3D(P>0) = 1;
NumberUnsynRainyDays2D = sum(NumberUnsynRainyDays3D, 3, 'omitnan');

%% sort rainfall 

Zeros3D = zeros(size(P), 'double');
Ones2D =  ones([Pm Pn], 'int16');

P_OriginalSum2D = sum(P,3, 'native', 'omitnan'); % find array of annual rainfall, before re-allocation

% write unsynthesised rainfall in original spatial resolution to nc file
ncwrite(RainName, 'Rain', permute(P,[ 3 1 2]) ); 

% write synthesised annual rainfall in original spatial resolution to nc file
ncwrite(RainName, 'AnnualUnsynRain', P_OriginalSum2D);

% write total number of synthesised rainy days in the year to NC file
ncwrite(RainName, 'NumUnsynRainyDays', NumberUnsynRainyDays2D);

%% 5 day summation of unsynthesised rainfall 
    
 %Rain = ncread(RainDataName,'Rainfall',start,stride); 
RainUnsynSum = P; % find rainfall summation
Rain6 = movsum(RainUnsynSum,[5 0], 'omitnan'); % computes sum 5 trailing values AND value of cell
RainUnsynSum = int16((Rain6) - double(P)); % subtract rain value for value in cell to leave the sum of rain for 5 days preceding

% write synthesised rainfall in original spatial resolution to nc file
ncwrite(RainName, 'RainSum', permute(RainUnsynSum, [3 1 2]));

% syn rain for day 240 only
%RainSynDay240 = zeros(size(Ones2D));
RainSynDay240 = RainUnsynSum(:,:,240);
ncwrite(RainName, 'RainDay240', RainSynDay240);
 
 %% RainSum (rs) - clip lon and lat to limits
% Load longitude array 
lon_rs = double(ncread(RainName,'longitude'));
% Load latitude array 
lat_rs = double(ncread(RainName,'latitude'));
% determine which indices correspond to dimension 1:
ind1_rs = find(lon_rs >= lon_w & lon_rs <= lon_e);
% determine which indices correspond to dimension 2:
ind2_rs = find(lat_rs >= lat_s & lat_rs <= lat_n);
% Clip lat and lon to their specified range: 
lat_rs = lat_rs(ind2_rs);
lon_rs = lon_rs(ind1_rs);

% create meshgrid
%[Lat_rs, Lon_rs] = meshgrid(lat_rs, lon_rs);
%Lat_rs = fliplr(Lat_rs); % arrange array to that of Global Curve Number

%% create template arrays

P250Annual = int16(zeros(length(lon_gcn), length(lat_gcn)));

 
 %% write 5 day rain summation in NC file at resolution of Global Curve Number product
[XX, YY] = meshgrid(lat_rs,lon_rs); % meshgrid of array low resolution
[XXq, YYq] = meshgrid(lat_gcn, lon_gcn); % meshgrid of array of high resolution


for i=1:yeardays(Yr)
    
   
    XY = ncread(RainName, 'Rain', [i, ind1_rs(1), ind2_rs(1)], [1, length(ind1_rs), length(ind2_rs)]);
    P250 = int16(interp2(XX,YY,permute(XY,[2 3 1]),XXq,YYq,'nearest'));
    ncwrite(Rain250Name, 'Rain250', P250, [1 1 i] );
    P250Annual = P250Annual + P250;
    XY = ncread(RainName, 'RainSum', [i, ind1_rs(1), ind2_rs(1)], [1, length(ind1_rs), length(ind2_rs)]);
    SUM250 = int16(interp2(XX,YY,permute(XY,[2 3 1]),XXq,YYq,'nearest'));
    ncwrite(Rain250Name, 'RainSum250', SUM250, [1 1 i] );
end

ncwrite(Rain250Name, 'RainAnnual250', P250Annual, [1 1] );

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.