% Process binary data for each dataset 

% DATASETS THAT WERE REGRIDDED IN GRADS TO FIT 2.5 x 2.5 domain
% erslp.2.5_2.5.gs (more coarse box averaging)
% 20thC.2.5_2.5.gs  (more coarse box averaging)
% kaplanslp.2.5_2.5.gs (more fine bi-linear)
% hadslp.2.5_2.5.gs (more fine bi-linear)

% (1) Reduce to 10N-10S to save space (delete data immediately)
% (2) Remove climo (based on full record) to form anomalies but keeps raw data too

% NINE DATASETS (not including ICOADS)
%CDAS: Jan 1950 Mar 2012  (A)
%ERA40: Jan 1958 Aug 2002 (E)
%20thC: Jan 1900 Dec 2010 (F)
%CFSR: Jan 1979 Mar 2012  (G)
%MERRA: Jan 1979 Feb 2012 (H)
%ERA-INT: Jan 1979 Jan 2012 (I)
%ERSLP     (B)
%HADSLP    (C)
%KaplanSLP (D)

%%%%%%%%% FIXED
% Reduce to 110E to (280) 80W (Indo - Pacific basin) 
%lat/lon index 
lat = 73;
lon = 144;
latsi = [-90:2.5:90]; %73
lonsi = [0:2.5:357.5]; %144
nw_latsi = latsi(33:41); % NEW RANGE OF LATITUDES (10N - 10S)
latp1 = 33;
latp2 = 41;
nw_lonsi = lonsi(45:113);  % NEW RANGE OF LONGITUDES (110E - 80W)
lonp1 = 45;
lonp2 = 113;

%%%%%%%%%%%%

% LOAD CDAS
tim = 747
inityr = 1950;
finyr = 2012;
clim1 = 1950;
clim2 = 2012;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.CDAS1.jan1950-mar2012.dat');
B = fread(fid,'real*4');
fclose(fid);
%B(B>999.)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoA = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoA;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(eqslpoA);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoA);
latd = dom(1)*dom(2);
slp = reshape(eqslpoA,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpA = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpA = reshape(anomslpA,[dom(1) dom(2) tim]);

clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA anomslpA

% LOAD ERSLP
tim = 1176
inityr = 1900;
finyr = 1997;
clim1 = 1900;
clim2 = 1997;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.erslp.2.5x2.5.1900-1997.dat');
B = fread(fid,'real*4');
fclose(fid);
B(B==32767)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoB = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoB;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(eqslpoB);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoB);
latd = dom(1)*dom(2);
slp = reshape(eqslpoB,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpB = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpB = reshape(anomslpB,[dom(1) dom(2) tim]);



clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA anomslpA eqslpoB anomslpB


% LOAD HADSLP
tim = 1260
inityr = 1900;
finyr = 2004;
clim1 = 1900;
clim2 = 2004;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.hadslp2.2.5x2.5.1900-2004.dat');
B = fread(fid,'real*4');
fclose(fid);
B(B==32767)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoC = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoC;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(eqslpoC);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoC);
latd = dom(1)*dom(2);
slp = reshape(eqslpoC,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpC = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpC = reshape(anomslpC,[dom(1) dom(2) tim]);


clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA anomslpA eqslpoB anomslpB eqslpoC anomslpC


% LOAD KaplanSLP 
tim = 1219
inityr = 1900;
finyr = 2001;
clim1 = 1900;
clim2 = 2001;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.kaplan.2.5x2.5.1900-jul2001.dat');
B = fread(fid,'real*4');
fclose(fid);
B(B==32767)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoD = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoD;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(eqslpoD);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoD);
latd = dom(1)*dom(2);
slp = reshape(eqslpoD,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpD = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpD = reshape(anomslpD,[dom(1) dom(2) tim]);





clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA anomslpA eqslpoB anomslpB eqslpoC anomslpC ...
    eqslpoD anomslpD

% LOAD ERA40 
tim = 536
inityr = 1958;
finyr = 2002;
clim1 = 1958;
clim2 = 2002;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.ERA40.1958-aug2002.dat');
B = fread(fid,'real*4');
fclose(fid);
%B(B==32767)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoE = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoE;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(eqslpoE);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoE);
latd = dom(1)*dom(2);
slp = reshape(eqslpoE,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpE = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpE = reshape(anomslpE,[dom(1) dom(2) tim]);




clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA eqslpoB eqslpoC eqslpoD eqslpoE ...
    anomslpA anomslpB anomslpC anomslpD anomslpE
    
% LOAD 20thC
tim = 1332
inityr = 1900;
finyr = 2010;
clim1 = 1900;
clim2 = 2010;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.20thC.2.5x2.5.1900-2010.dat');
B = fread(fid,'real*4');
fclose(fid);
%B(B==32767)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoF = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoF;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(datfi);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoF);
latd = dom(1)*dom(2);
slp = reshape(eqslpoF,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpF = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpF = reshape(anomslpF,[dom(1) dom(2) tim]);





clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA eqslpoB eqslpoC eqslpoD eqslpoE eqslpoF ...
    anomslpA anomslpB anomslpC anomslpD anomslpE anomslpF

% LOAD CFSR
tim = 399
inityr = 1979;
finyr = 2012;
clim1 = 1979;
clim2 = 2012;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.CFSR.jan1979-mar2012.dat');
B = fread(fid,'real*4');
fclose(fid);
%B(B==32767)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoG = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoG;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(datfi);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoG);
latd = dom(1)*dom(2);
slp = reshape(eqslpoG,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpG = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpG = reshape(anomslpG,[dom(1) dom(2) tim]);



clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA eqslpoB eqslpoC eqslpoD eqslpoE eqslpoF eqslpoG ...
    anomslpA anomslpB anomslpC anomslpD anomslpE anomslpF anomslpG

% LOAD MERRA
tim = 398
inityr = 1979;
finyr = 2012;
clim1 = 1979;
clim2 = 2012;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.merra.jan1979-feb2012.dat');
B = fread(fid,'real*4');
fclose(fid);
%B(B==32767)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoH = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoH;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(datfi);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoH);
latd = dom(1)*dom(2);
slp = reshape(eqslpoH,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpH = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpH = reshape(anomslpH,[dom(1) dom(2) tim]);




clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA eqslpoB eqslpoC eqslpoD eqslpoE eqslpoF eqslpoG ...
    eqslpoH ...
    anomslpA anomslpB anomslpC anomslpD anomslpE anomslpF anomslpG anomslpH

% LOAD ERA-INT
tim = 397
inityr = 1979;
finyr = 2012;
clim1 = 1979;
clim2 = 2012;
noyrz = (finyr - inityr)+1;
fid = fopen('MSLP.eraint.jan1979-jan2012.dat');
B = fread(fid,'real*4');
fclose(fid);
%B(B==32767)=NaN;
dataslp = reshape(B,lat*lon,tim); 
clear B
% Reduce to 10N-10S
redslp = reshape(dataslp,[144 73 tim]);  
eqslpoI = redslp(lonp1:lonp2,latp1:latp2,:);

% Calculate climo on full record 
datfi = eqslpoI;
mthidx = [1:12];
mx = repmat(mthidx,1,noyrz);
[x a j] = size(datfi);
%j = size(datfi);
mx1 = mx(1:j);
yrs = [inityr:inityr+noyrz-1];
b = repmat(yrs,12,1);
yx= b(:);
yx1 = yx(1:j)';

% Reshape for analysis
dom = size(eqslpoI);
latd = dom(1)*dom(2);
slp = reshape(eqslpoI,latd,tim);

% Remove seasonal cycle based on reduced lat/lon
tonslp = [ ];
for q = 1:12
mon = find((mx1==q)&(yx1>=clim1)&(yx1<=clim2));
onslp = mean2(slp(:,mon)');
tonslp = [tonslp onslp'];
end
%plot tonslp (to look at climo)
clim = repmat(tonslp,1,noyrz);   
clim = clim(:,1:length(mx1));
% subtract from raw
anomslpI = slp - clim;
clear clim 
% Reshape again to isolate Lat/Lon for domain
anomslpI = reshape(anomslpI,[dom(1) dom(2) tim]);


clearvars -except lat lon latsi lonsi nw_latsi latp1 latp2 nw_lonsi ...
    lonp1 lonp2 eqslpoA eqslpoB eqslpoC eqslpoD eqslpoE eqslpoF eqslpoG ...
    eqslpoH eqslpoI anomslpI ...
    anomslpA anomslpB anomslpC anomslpD anomslpE anomslpF anomslpG anomslpH



