% CLIMATOLOGY OF PRECIPITATION IN SOUTH AMERICA


% % CESM
% mdl_prior='cesm_lme010';
% [gcm_syr,gcm_eyr]=modelyrs(mdl_prior);
% mon_avg_o=1;mon_avg_f=12;
% p_yrs=[851:1848]-gcm_syr+1;
% state_tp = 'p';
% [X,xlat,xlon] = load_Smon(mdl_prior,state_tp,mon_avg_o,mon_avg_f,p_yrs);
% clat=xlat;clon=xlon;
% 
% % Calculate the land mask
% %varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
% %sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
% 
% %---------------------
% % Compute climatology
% %---------------------
% 
% % Just take last half of data
% 
% X2=zeros(size(X,1),12,size(X,2)/12);
% j=1;k=12;
% for i=1:size(X,2)/12
%    X2(:,:,i)=X(:,j:k);
%    j=j+12;k=k+12;
% end
% 
% % Climo
% Xclim=mean(X2,3);
% %Xann(Xann==0)=NaN; % add back fill nans
% Xclim=reshape(Xclim,length(xlat),length(xlon),12);
% 
% 
% % percent of precip that falls in a given season
% pann=sum(Xclim,3);
% pdjf=sum(Xclim(:,:,[12,1,2]),3);
% pmam=sum(Xclim(:,:,[3,4,5]),3);
% pjja=sum(Xclim(:,:,[6,7,8]),3);
% pson=sum(Xclim(:,:,[9,10,11]),3);
% 
% 
% s=load('coast');mlat=s.lat;mlon=s.long;
% cmx=0.75;
% 
% figure
% h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
% 'grid','on','MLineLocation',[xlon(1:72) (xlon(73:end)-360)],'PLineLocation',xlat,'GLineStyle','-');
% rto=pdjf./pann;
% rto(rto<0.33)=NaN;
% pcolorm(xlat,xlon,rto)
% caxis([0 cmx])
% hold on
% plotm(mlat,mlon,'color',[0 0 0])
% title(['DJF'])
% set(gca,'FontSize',14)
% colorbar
% tightmap
% %print(['/home/nsteiger/Dropbox/sa-phyda/sa_precip_climo_cesm_djf'],'-dpng','-r300')
% 
% % figure
% % h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
% % 'grid','on','MLineLocation',[xlon(1:72) (xlon(73:end)-360)],'PLineLocation',xlat,'GLineStyle','-');
% % pcolorm(xlat,xlon,pmam./pann)
% % caxis([0 cmx])
% % hold on
% % plotm(mlat,mlon,'color',[0 0 0])
% % title(['MAM'])
% % set(gca,'FontSize',14)
% % colorbar
% % tightmap
% % %print(['/home/nsteiger/Dropbox/sa-phyda/sa_precip_climo_cesm_mam'],'-dpng','-r300')
% 
% 
% figure
% h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
% 'grid','on','MLineLocation',[xlon(1:72) (xlon(73:end)-360)],'PLineLocation',xlat,'GLineStyle','-');
% rto=pjja./pann;
% rto(rto<0.33)=NaN;
% pcolorm(xlat,xlon,rto)
% caxis([0 cmx])
% hold on
% plotm(mlat,mlon,'color',[0 0 0])
% title(['JJA'])
% set(gca,'FontSize',14)
% colorbar
% tightmap
% %print(['/home/nsteiger/Dropbox/sa-phyda/sa_precip_climo_cesm_jja'],'-dpng','-r300')
% 
% 
% % figure
% % h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
% % 'grid','on','MLineLocation',[xlon(1:72) (xlon(73:end)-360)],'PLineLocation',xlat,'GLineStyle','-');
% % pcolorm(xlat,xlon,pson./pann)
% % caxis([0 cmx])
% % hold on
% % plotm(mlat,mlon,'color',[0 0 0])
% % title(['SON'])
% % set(gca,'FontSize',14)
% % colorbar
% % tightmap
% % %print(['/home/nsteiger/Dropbox/sa-phyda/sa_precip_climo_cesm_son'],'-dpng','-r300')





% LOAD OBSERVATIONAL DATA
%mon_avg_o=12;mon_avg_f=2;
%obs_dtp1='vanderS_pdsi';
%mon_avg_o=4;mon_avg_f=3;
%obs_dtp1='bearth'; % particular data type specified by input
%[O_t,xlat,xlon,x_yrs] = load_obs_S(obs_dtp1,mon_avg_o,mon_avg_f);


% GPCP
%X=ncread('/d1/nsteiger/climate-data/gpcp/gpcp_v2.3_climatology.nc','precip');
%X=permute(X,[2 1 3]);
%xlat=double(ncread('/d1/nsteiger/climate-data/gpcp/gpcp_v2.3_climatology.nc','lat'));
%xlon=double(ncread('/d1/nsteiger/climate-data/gpcp/gpcp_v2.3_climatology.nc','lon'));

% CRU
X=ncread('/d1/nsteiger/climate-data/cru_ts3/cru_ts3.23.1901.2014.pre.dat.nc','pre');
X=permute(X,[2 1 3]);
xlon = double(ncread('/d1/nsteiger/climate-data/cru_ts3/cru_ts3.23.1901.2014.pre.dat.nc','lon'));
xlat = double(ncread('/d1/nsteiger/climate-data/cru_ts3/cru_ts3.23.1901.2014.pre.dat.nc','lat'));
[X,xlon]=regrid2gcm(X,xlon);
xlat=xlat(:);xlon=xlon(:);

%---------------------
% Compute climatology
%---------------------

% Just take last half of data
% monthly indices
r_o=1901;r_f=2014;
ovrlprng=rng_intr([r_o,r_f],[1950,2014]);
mon_r=r_o:(1/12):(r_f+11/12);
irm1_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm1_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

X1=X(:,:,irm1_a:irm1_b);
X2=zeros(size(X,1),size(X,2),12,length(irm1_a:irm1_b)/12);
j=1;k=12;
for i=1:length(irm1_a:irm1_b)/12
   X2(:,:,:,i)=X1(:,:,j:k);
   j=j+12;k=k+12;
end

% Climo
Xclim=mean(X2,4);
%Xann(Xann==0)=NaN; % add back fill nans


% percent of precip that falls in a given season
pann=sum(Xclim,3);
pdjf=sum(Xclim(:,:,[12,1,2]),3);
pmam=sum(Xclim(:,:,[3,4,5]),3);
pjja=sum(Xclim(:,:,[6,7,8]),3);
pson=sum(Xclim(:,:,[9,10,11]),3);


% s=load('coast');mlat=s.lat;mlon=s.long;
% cmx=0.75;
% 
% figure
% %h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
% h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
% 'grid','on','MLineLocation',[clon(1:72) (clon(73:end)-360)],'PLineLocation',clat,'GLineStyle','-');
% pcolorm(xlat,xlon,pdjf./pann)
% caxis([0 cmx])
% hold on
% plotm(mlat,mlon,'color',[0 0 0])
% title(['DJF'])
% set(gca,'FontSize',14)
% colorbar
% tightmap
% %print(['/home/nsteiger/Dropbox/sa-phyda/sa_precip_climo_djf'],'-dpng','-r300')
% 
% % figure
% % h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
% % pcolorm(xlat,xlon,pmam./pann)
% % caxis([0 cmx])
% % hold on
% % plotm(mlat,mlon,'color',[0 0 0])
% % title(['MAM'])
% % set(gca,'FontSize',14)
% % colorbar
% % tightmap
% % %print(['/home/nsteiger/Dropbox/sa-phyda/sa_precip_climo_mam'],'-dpng','-r300')
% 
% 
% figure
% %h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
% h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
% 'grid','on','MLineLocation',[clon(1:72) (clon(73:end)-360)],'PLineLocation',clat,'GLineStyle','-');
% pcolorm(xlat,xlon,pjja./pann)
% caxis([0 cmx])
% hold on
% plotm(mlat,mlon,'color',[0 0 0])
% title(['JJA'])
% set(gca,'FontSize',14)
% colorbar
% tightmap
% %print(['/home/nsteiger/Dropbox/sa-phyda/sa_precip_climo_jja'],'-dpng','-r300')
% 
% 
% % figure
% % h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
% % pcolorm(xlat,xlon,pson./pann)
% % caxis([0 cmx])
% % hold on
% % plotm(mlat,mlon,'color',[0 0 0])
% % title(['SON'])
% % set(gca,'FontSize',14)
% % colorbar
% % tightmap
% % %print(['/home/nsteiger/Dropbox/sa-phyda/sa_precip_climo_son'],'-dpng','-r300')



%--------------------------------
% PLOT THE WHOLE GLOBE

s=load('coast');mlat=s.lat;mlon=s.long;
cmx=0.75;

figure
h=axesm('MapProjection','eqdcylin');%,'MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
%h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
%'grid','on','MLineLocation',[clon(1:72) (clon(73:end)-360)],'PLineLocation',clat,'GLineStyle','-');
pcolorm(xlat,xlon,pdjf./pann)
caxis([0 cmx])
hold on
plotm(mlat,mlon,'color',[0 0 0])
title(['DJF'])
set(gca,'FontSize',14)
colorbar
tightmap
print(['/home/nsteiger/Dropbox/sa-phyda/glbl_precip_climo_djf'],'-dpng','-r300')

figure
h=axesm('MapProjection','eqdcylin');
%h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
pcolorm(xlat,xlon,pmam./pann)
caxis([0 cmx])
hold on
plotm(mlat,mlon,'color',[0 0 0])
title(['MAM'])
set(gca,'FontSize',14)
colorbar
tightmap
print(['/home/nsteiger/Dropbox/sa-phyda/glbl_precip_climo_mam'],'-dpng','-r300')


figure
h=axesm('MapProjection','eqdcylin');
%h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
%h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
%'grid','on','MLineLocation',[clon(1:72) (clon(73:end)-360)],'PLineLocation',clat,'GLineStyle','-');
pcolorm(xlat,xlon,pjja./pann)
caxis([0 cmx])
hold on
plotm(mlat,mlon,'color',[0 0 0])
title(['JJA'])
set(gca,'FontSize',14)
colorbar
tightmap
print(['/home/nsteiger/Dropbox/sa-phyda/glbl_precip_climo_jja'],'-dpng','-r300')


figure
h=axesm('MapProjection','eqdcylin');
%h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
pcolorm(xlat,xlon,pson./pann)
caxis([0 cmx])
hold on
plotm(mlat,mlon,'color',[0 0 0])
title(['SON'])
set(gca,'FontSize',14)
colorbar
tightmap
print(['/home/nsteiger/Dropbox/sa-phyda/glbl_precip_climo_son'],'-dpng','-r300')
















