% Jim Carton 7-7-2020
% This program combines several monthly river discharge data sets.  We begin with our primary data set: the Dai-Trenberth (2002) set of
%		925 rivers.  
%		1) First, we fill in missing data with climatology to produce a continuous record 1/1980-12/2019 (480 months).  
%		2) Next we replace seven Arctic rivers whose time series are maintained at Greatarcticrivers.org.  This data may have come, 
%		originally, from UNH/GRDC: ESTIMATED RIVER RUNOFF, and ultimately from a Russian source.  You will see a bit of tinkering 
%		with Pechora since the gauge closest to the delta does not span our full period of interest 1980-2019.
%		3) We then inflate a larger set of 12 Arctic rivers by 30% to bring the basin-mean discharge up to the published estimate of 0.1Sv. 
%		4) Finally we  
%
%
%
% Get the Dai-Trenberth river discharge data set for 1/1980-12/2019 (480 months)
%
[fa_full,fa_clim,fa_lon,fa_lat] = get_rivers
%
% Get the GreatarcticRivers discharge data set for seven Arctic rivers for 1/1980-12/2019 (480 months)
%
[Aa_full] = get_Arctic_rivers; % loads Aa_full with Arctic river discharge time series [480x7]
%
fa_full(:,7) = Aa_full(:,1)*1.3;% replace yenise
fa_full(:,9) = Aa_full(:,2)*1.3;% replace lena
fa_full(:,13)= Aa_full(:,3)*1.3;% replace ob
fa_full(:,19)= Aa_full(:,4)*1.3;% replace mackenzie
fa_full(:,24)= Aa_full(:,5)*1.3;% replace yukon
fa_full(1:205,31)=fa_full(1:205,31)*1.3;
fa_full(206:480,31)= Aa_full(206:480,6)*1.3*4.90/3.48; %replace part of Pechora (see below)
fa_full(:,35)= Aa_full(:,7)*1.3;% replace kolyma
fa_full(:,33)=fa_full(:,33)*1.3;% inflate Khatanga
fa_full(:,37)=fa_full(:,37)*1.3;% inflate Sev. Dvina
fa_full(:,64)=fa_full(:,64)*1.3;% inflate Indigirka
fa_full(:,97)=fa_full(:,97)*1.3;% inflate Olenek
fa_full(:,104)=fa_full(:,104)*1.3;% inflate Yana
%
% Note: D-T bases their analysis on the gauge at Oksino 141 km from delta.  I extend the time series using the record from Ust'-Tsil'ma (300+ km from delta)
%      since the discharge is lower at Ust'-Tsil'ma I multiply the U-T-based extension by the ratio of the mean transports at the two gauges, 4.90/3.48=1.41. 
% Note: We then follow the guidence of AOMIP (https://www.whoi.edu/page.do?pid=30587) to inflate the discharge of the 12 gauged rivers by 30% 
%      to account for the missing ungauged discharge.  The total should be ~0.10Sv according to Aagard and Carmack (1989).
%
%  Map river discharge estimates onto a 1x1 grid
%
gridded_discharg = zeros(480,360,720); %initialize a 1/2x1/2-deg discharge field spanning [-90-90N, 0E-360E] and [1/1980 to 12/2019]. 

for irev = 1:size(fa_full,2) % outer loop over the 925 river discharge time series
	for itime = 1:480 % itime=1 corresponds to 1/1980, itime=480 corresponds to 12/2019
% convert discharge locations [lat,lon] to gridpoint space [j,i] by rounding and adding 1 (to guarranty i,j >= 1)
	i = (round(fa_lon(irev,1)*2)+360)+1; % longitudes are given in the range [-180 to 180E] while we shift to a grid where i=1 corresponds to 0E
	j = (round(fa_lat(irev,1)*2)+180)+1;
	gridded_discharg(itime,j,i) = fa_full(itime,irev);
end;end;
%
% Load the Mankoff et al Greenland discharge data set for 267 discharge locations for the period 1/1987-12/2019 (396 months) 
%  (currently I do not include greenland discharge for 1/1980-12/1986)
%
[Ga_full,Ga_lon,Ga_lat] = get_Greenland_rivers
%
%
for irev = 1:size(Ga_full,2) % outer loop over the 267 glacial discharge time series
	for itime = 85:480 % itime=85 corresponds to 1/1987, itime=480 corresponds to 12/2019
% convert discharge locations [lat,lon] to [j,i]
	i = (round(Ga_lon(irev,1)*2)+360)+1; % longitudes are given in the range -180 to 180E while we shift to a grid where i=1 corresponds to 0E
	j = (round(Ga_lat(irev,1)*2)+180)+1;
	gridded_discharg(itime,j,i) = Ga_full(itime-84,irev); % itime=85 corresponds to 1/1987, itime=480 corresponds to 12/2019
end;end;
gridded_discharg(isnan(gridded_discharg))=0; %replace any NaNs that sneek in  with 0 -- I didn't check to see if there are any

for i = 1:720;
	for j=1:360;
		x(i)=(i-361)/2.; %this inverts line #60
		y(j)=(j-181)/2;%this inverts line #60
%		x(i)=i;
%		y(j)=j;
end;end;
for it = 1:480
	timenum(it)=it;
end;
nccreate('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','time','dimensions',{'time',480},'Format','netcdf4');
ncwrite('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','time',timenum);
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','time','standard_name','time');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','time','long_name','time');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','time','units','months since 1979-12-15 00:00:00');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','time','calendar','standard');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','time','axis','T');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','time','cell_methods','time: mean');

nccreate('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','x','dimensions',{'x',720});
ncwrite('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','x',x);
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','x','standard_name','longitude');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','x','long_name','longitude');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','x','units','degrees_E');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','x','axis','X');

nccreate('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','y','dimensions',{'y',360});
ncwrite('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','y',y);
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','y','standard_name','latitude');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','y','long_name','latitude');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','y','units','degrees_N');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','y','axis','Y');

%nccreate('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','lat','dimensions',{'x','y'});
%ncwrite('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','lat',lat);
%nccreate('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','lon','dimensions',{'x','y'});
%ncwrite('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','lon',lon);

nccreate('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','discharge','Dimensions',{'time','y','x'},'ChunkSize',[480 360 1]);
ncwrite('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','discharge',gridded_discharg);
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','discharge','standard_name','discharge');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','discharge','long_name','river_discharge');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','discharge','units','m^3/s');
%ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','discharge','_FillValue','-1.e+20f');
%ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','discharge','missing_value','-1.e+20f');
ncwriteatt('/aosc/horse/carton/rivers/gridded_discharge_mo_360x720x480.nc','discharge','cell_methods','time: mean');


function [Ga_full,Ga_lon,Ga_lat] = get_Greenland_rivers
%
% download Greenland discharge data.  To check file content: ncdisp('/aosc/horse/carton/rivers/greenland/gate.nc');
%
Ga_lat = ncread('/aosc/horse/carton/rivers/greenland/gate.nc','mean_lat'); % mean lat of discharge gates [267x1]
Ga_lon = ncread('/aosc/horse/carton/rivers/greenland/gate.nc','mean_lon'); % mean lon of discharge gates [267x1]
jday = ncread('/aosc/horse/carton/rivers/greenland/gate.nc','time')+juliandate(1986,04,15); % days since 1986-04-15 converted to JDAY [445x1]
discharge = ncread('/aosc/horse/carton/rivers/greenland/gate.nc','discharge')*31.7; % land_ice_volume_flow_rate_due_to_calving_and_ice_front_melting [converted to m3/s] [445x267]
%
%  Greenland discharge is NOT sampled uniformly in time.  Here I do a spline interpolation onto a  monthly grid for the period  1/1987 through 12/2019
%
ind = 0
for iyear = 1987:2019
	for imonth = 1:12
		ind = ind+1;
		timevec(ind) = juliandate(iyear,imonth,15); % timevec is a monthly time count beginning 1/15/1987 for 33 years [1x396]
end;end;
Ga_full = interp1(double(jday),discharge,timevec,'spline'); % Greenland discharge beginning 1/1987 [396x267]
% Helheim glacier: 172; jakobshavn glacier: 133; For jakobshavn wikipedia reports an average discharge rate of ~1100 m3/s which is roughly consistent

end; 
% - - - - - 
%
function [fa_full,fa_clim,fa_lon,fa_lat] = get_rivers
%
% Load Dai-Trenberth river discharge data.  To check file content: ncdisp('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc');
%
fa_lon = ncread('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc','lon_mou'); % lon of river mouth [925x1]
fa_lat = ncread('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc','lat_mou'); % lat of river mouth [925x1]
flow = ncread('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc','FLOW'); % monthly mean transport for each river gauge. Begins 1/1/1900 [925x1428] 
ratio_m2s = ncread('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc','ratio_m2s'); % ratio of transport at river mouth and at gauge station [925x1]
riv_name = ncread('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc','riv_name'); % river name [30x925]
ocn_name = ncread('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc','ocn_name'); % name of ocean into which river discharges
time = ncread('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc','time'); % 'time as YYYYMM' [1428x1] [961 is 1/1/1980]
%
% Compute the monthly climatology for each river.  We'll use this to fill in missing values
%
fa_trans = flow(:,961:1428); %load the 40-yr monthly discharge 1/1980-12/2018
fa1 = fa_trans';%put flow into matrix fa with time in rows, rivers in columns [468,925]
for irev = 1:size(fa1,2); % correct for the fact that the measurements are made at gauge, not river mouth
	fa(:,irev)=fa1(:,irev)*ratio_m2s(irev);
end

%%% check: S = nansum(fa,2); plot(S);  hold on; plot(movmean(S,12)) %The total river discharge as reported by Dai and Trenberth (2002) is 1.18Sv
clear fa_trans fa1
%
% 
%
fa_clim = zeros(12,size(fa,2)); %initialize the climatology
for irev = 1:size(fa,2); % loop over 925 rivers and reshape each so that there are 12 columns for the 12  months
	far = reshape(fa(:,irev),12,[]);
	faa = nanmean(far,2); %compute the climotological monthly average ignoring NANs
	fa_clim(:,irev) = faa; %load it into fa_clim
end
clear far faa
%
% build the 40-yr long river discharge array fa_full by beginning with monthly climatology and inserting observed values
%
fa_full = repmat(fa_clim,40,1); % initialize fa_full with the monthly climatology 1/1980-12/2019 [480,925]
for irev = 1:size(fa,2)
	for itime = 1:468
		if ~isnan(fa(itime,irev)) 
			fa_full(itime,irev)=fa(itime,irev); % insert observed values
end;end;end;


%%% check: S = nansum(fa_full,2); hold on; plot(S); hold on; plot(movmean(S,12));  % The Global sum river discharge as reported by Dai and Trenberth (2002) is 1.18Sv and this gives about 0.8Sv
%%% check: the Ob: plot(fa_full(:,13)), lena: plot(fa_full(:,9))

end; % end function get_rivers
% - - - - - 
%
function [Aa_full] = get_Arctic_rivers
%
% Load seven monthly Arctic river discharge time series from greatarcticrivers.org.  
% I converted the daily data to monthly and filled gaps in excel.  So this is from my excel files.
% To check file content: ncdisp('/aosc/horse/carton/rivers/dai-tren/coastal-stns-Vol-monthly.updated-May2019.nc');
%
[yenise] = textread('/aosc/horse/carton/rivers/arctic_rivers/yenise-1980-2019.txt','%f'); %#7
[lena] = textread('/aosc/horse/carton/rivers/arctic_rivers/lena-1980-2019.txt','%f');% #9
[ob] = load('/aosc/horse/carton/rivers/arctic_rivers/ob-1980-2019.txt','%f');% #13
[mackenzie] = load('/aosc/horse/carton/rivers/arctic_rivers/mackenzie-1980-2019.txt','%f'); % #19
[yukon] = load('/aosc/horse/carton/rivers/arctic_rivers/yukon-1980-2019.txt','%f'); % #24
[pechora] = load('/aosc/horse/carton/rivers/arctic_rivers/pechora-1980-2019.txt','%f'); % #31
[kolyma] = load('/aosc/horse/carton/rivers/arctic_rivers/kolyma-1980-2019.txt','%f'); % #35

Aa_full = zeros(480,7);
Aa_full(:,1) = yenise(:);
Aa_full(:,2) = lena(:);
Aa_full(:,3) = ob(:);
Aa_full(:,4) = mackenzie(:);
Aa_full(:,5) = yukon(:);
Aa_full(:,6) = pechora(:);
Aa_full(:,7) = kolyma(:);

end; % end get_Arctic_rivers




