%  era5u10.m  make a SAM index from 10m zonal wind, or 100m when I get it
%  sstsam5.m  Program to composite sst on SAM based on monthly data
% this one uses ERA52m Temperature
% April 2020 got latest data from psl.noaa.gov/datagriddeddata.noaa.ersse.v5.html  See references there
%  based on ssteof.m  Make eofs from SST in various basins
%  program to use sst.mnmean_ersst.nc Smith NOAA web site
%data starts in 1854  montly means 2-degree resolution 88N-88S  0-358
%  Dennis L. Hartmann  January 2016  **1944 months in todays sst.mnmean.nc V3b
%close all
clear all
%  PLAN 
%  make a data set of just North Pacific, or Atlantic data
%  a two dimensional array with basin points in one direction and time in the other
%  Compute EOFs and PCs of this
%  Regress PCs on the whole map of SST to show patterns  Would help to have land mask I think
%  Suggest loading topo map, interpolating to SST grid, and using elevation aas land-sea mask

ilag=0         %  We're going to look at the SST regression before the anomaly develops
icor= 0            % this is to make correlation maps instead of regression maps  if = 1
jfilt = 0
detrend = 0
gmremove = 0
montrend = 0
ivmx=3
fcut = 72.;
yrstart=1959;
yrbeg= 1880;
yrbeg= 1950;
yrbeg= 1854;
yrbeg= 1900
yrbeg= 1959;
yrend= 2022;
lon1=150.; lon2=290.;

   yskip= (yrbeg-yrstart)*12
it1=1+yskip
% sub select season you want

%           dmo1= 5; dm02 = 9
           dmo1= 1; dm02 = 12
%	   dmo1= 11; dm02=3
		       dmo2=dm02;
if dm02 < dmo1
	   dmo2=dm02+12;
end

%  OK we are ready to make our array, using some prescribed lat long zones and our land mask

%  below for equatorial ENSO
xlat1 = 10.;  xlat2 = -10.;
xlon1 = 120;  xlon2 = 280.;
%  below for tropical Pacific
xlat1 = 30.;  xlat2 = -30.;
xlon1 = 110;  xlon2 = 280.;
% Below for Global Pacific ST
xlat1 = 87.;  xlat2 = -87.;
xlon1 = 1;  xlon2 = 285.0;
%  below for all Tropics 
xlat1 = 20.;  xlat2 = -20.;
xlon1 = 1.;  xlon2 = 357.;
% Below for full pacific N of 65S
xlat1 = 65.;  xlat2 = -65.;
xlon1 = 120;  xlon2 = 255.;
% Below for Rachel's Atlantic
xlat1 = 65.;  xlat2 = -65.;
xlon1 = 230;  xlon2 = 320.;
% Below for Atlantic SST
xlat1 = 80.;  xlat2 = -65.;
xlon1 = 285;  xlon2 = 357.0;
expname='Atlantic';
% Below for PDO region  N of 20N
xlat1 = 65.;  xlat2 = 20.;
xlon1 = 120;  xlon2 = 260.;
expname='PDORegion';
% Below for Global SST
xlat1 = 70.;  xlat2 = -70.;
xlon1 = 0;  xlon2 = 359.75;
expname='Global';


iymx=43;
%yint=32
%delda = yint * 12
%iymx=1

coast = load( 'coast.mat');

%filenme =['/home/disk/eos10/dennis/ERAUW.d/sst.mnmean.v4.nc']  % Through December of 2017 Version 4
%filenme =['/home/disk/eos10/dennis/ERAUW.d/sst.mnmeanV5.nc']  % Through January of 2021 Version 5
%filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_2mTemp.nc']  % Through March of 2022  ERA5 monthly
%filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_SST.nc']  % Through March of 2022  ERA5 monthly
%filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_10mU.nc']  % Through March of 2022  ERA5 monthly
filenme =['/home/disk/eos10/dennis/ERAUW.d/ERA5_1959_u10.nc']  % Through March of 2022  ERA5 monthly
 %  Lats go down from North Pole at 2.0 degrees  Longs are 0-358
%xlat1=5.0;  xlat2=-5.0;
%ila1=int32(1+(88-xlat1)/2.); ila2=int32(1+(88-xlat2)/2.);     %  Specifies indices


  level='10m';
var='u';

time = ncread(filenme,'time');
timxo=length(time)
timx=length(time)-(yskip);
%year = time/365.242199 + 1800;
%time=(time-time(1));   %now in months
				lat = ncread(filenme,'latitude');
				lon = ncread(filenme,'longitude');
%				ssto = ncread(filenme,'t2m');
%				ssto = ncread(filenme,'sst');
				ssto = ncread(filenme,'u10');
expver=1;



latx = int32(length(lat));
  lonx = int32(length(lon));
% subsample quarter-degree to one degree
ids=4;  % pick every 4th data point

%timx=size(ssto,3)
%timx=timx-3  ;  % data ends with march so subtract 3 to get on integer years
timx=timx  ;  % 1959 data was selected with integer years until 2021
%sst=squeeze(ssto(1:ids:lonx,1:ids:latx,expver,1:timx));  % reduced volume of data
sst=squeeze(ssto(1:ids:lonx,1:ids:latx,1:timx));  % reduced volume of data no experver in 1959 data


lat=lat(1:ids:latx);
lon=lon(1:ids:lonx);
lonx=length(lon); latx=length(lat);
  lat = double(lat);
lonxp = lonx+1
  dlon=lon(2)-lon(1);
  lonp=cat(1,lon,lon(lonx)+dlon);
lonp= double(lonp);
%ssto2=ssto(:,:,it1:timxo);
% make even number of years by adding missing values to end
%ssta=NaN(lonx,latx,timx);    %The plus zero here has to be changed depending on what month the series ends on
%ssta(:,:,1:timx)=ssto2;sst=ssta;
%sst(:,:,1:timx)=ssto(:,:,1:timx);
%timx=timx+9    % modified to include December 2017  year = [1:timx]/12.+yrbeg;
  month = mod(1:timx,12);
year = [1:timx]/12.+yrbeg;
month(month<1)=12;
clear ssto:
clear ssto2:
clear ssta;
maxsst=max(max(max(sst)))

timx=int32(size(sst,3))

% replace missing values with NaN
sst(sst< -100)= NaN;
sst=double(sst);
maxsst=max(max(max(sst)))

% here we develop a mask for topography
load('topo.mat','topo')
  topo(topo<0.)=0.;
lont=[0:359];
  latt=[-89:90];
lontx=length(lont);
  lattx=length(latt);
%  need to make two-D arrays of lats and longs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[lontb, lattb] = ndgrid(lont, latt);
[lonb, latb] = ndgrid(lon, lat);
lontb=double(lontb);
lonb=double(lonb);
lattb=double(lattb);
latb=double(latb);
topoz=interp2(lontb',lattb',topo,lonb,latb);

%figure
%contour(lon,lat,topoz',[0 0])
%  title('Interpolated Topography')


%  OK we are ready to make our array, using some prescribed lat long zones and our land mask
%xlat1 = 60.;  xlat2 = 20.;
%xlon1 = 130;  xlon2 = 255.;
dlon=lon(2)-lon(1);
dlat=lat(2)-lat(1);
ilat1= 1 + (xlat1-88.0)/dlat;  ilat2= 1 + (xlat2-88.0)/dlat;
ilon1= 1 + xlon1/dlon;   ilon2= 1 + xlon2/dlon;
ilat1=int32(ilat1);    ilat2=int32(ilat2);
ilon1=int32(ilon1);    ilon2=int32(ilon2);

%  remove the global mean at each time step
lat=double(lat);
cs=cos(lat*pi/180.);
ls=ones(lonx,1,'double');
[CS,LS]=meshgrid(cs,ls);
CSr=reshape(CS,lonx*latx,1);
sstr=reshape(sst,lonx*latx,timx);
for im=1:timx
	 weit(:,im)=isfinite(sstr(:,im));
	 mnwt(im)=squeeze(nansum(weit(:,im).*CSr,1));
	 mnsst(im)=squeeze(nansum(sstr(:,im).*CSr,1));
mnsst(im)=mnsst(im)/mnwt(im);
ssta(:,im)=sstr(:,im)-mnsst(im);
end
sstdif=sstr-ssta;
maxdif=squeeze(nanmax(nanmax(nanmax((sstdif)))));
maxdif
for im=1:timx
	 mnsst2(im)=squeeze(nansum(ssta(:,im).*CSr,1));
mnsst2(im)=mnsst2(im)/mnwt(im);

%sstb(:,im)=ssta(:,im)-mnsst2(im);
end
if gmremove == 1
  sst=reshape(ssta,lonx,latx,timx);
end
figure
plot(year,mnsst)
  title ('Global mean at each month')
  maxx=max(mnsst)*1.001;
minn=min(mnsst)*0.999;
  axis([yrbeg 2020 minn maxx]);
figure
plot(year,mnsst2)
  title ('Global mean at each month second time')
  maxx=max(mnsst2)*1.001;
minn=min(mnsst2)*1.001;
  axis([yrbeg 2020 minn maxx]);

%  Before we do anything else we want to remove the climatology, errors and weight anomalies by sqrt(cos(lat)
%  compute and remove climatology
whos sst
yrmx=int32(timx/12)
sstr = reshape(sst,lonx*latx,12,yrmx);
sstc = nanmean(sstr,3);
for j=1:12
for k=1:yrmx
sstr(:,j,k)=sstr(:,j,k)-sstc(:,j);   % removes climatology
end
end
maxsstr=max(max(max(sstr)))
%  Also should remove trend
%  OK, now climatology has been removed,  Next remove trend

if montrend ==1
%%  OK I am going to remove the trend from each month separately
t=double([1:yrmx]);
t=(t-mean(t))/std(t);
   ssts=reshape(sstr, lonx*latx*12,yrmx);
sstm=nanmean(ssts,2);
for it = 1:yrmx
ssts(:,it)=ssts(:,it)-sstm;
end
ssts(isnan(ssts))=0.0;     %  need to set anomalies to zero instead of NaN before computing trend, idiot
maxssts=max(max(max(ssts)))
corr = ssts*t'/double(yrmx);
sstrend=corr*t;
if detrend == 1
ssts=ssts-sstrend;          %  HERE is where TREND is removed
end
maxssts=max(max(max(ssts)))
sst=reshape(ssts,lonx*latx,12*yrmx);
ssts=sst;
else
t=double([1:timx]);
t=(t-mean(t))/std(t);
ssts=reshape(sstr,lonx*latx,timx);
sstm=nanmean(ssts,2);
for it = 1:timx
ssts(:,it)=ssts(:,it)-sstm;
end
ssts(isnan(ssts))=0.0;     %  need to set anomalies to zero before computing trend, idiot
maxssts=max(max(max(ssts)))
corr = ssts*t'/double(timx);
%  plot trend
figure
corll = reshape(corr,lonx,latx);
contourf(lon,lat,corll')
				title(['map of trend  ',num2str(yrbeg)])
colorbar
drawnow;
	maxcorr=max(max(abs(corr)))
sstrend=corr*t;
if detrend == 1
ssts=ssts-sstrend;          %  HERE is where TREND is removed
end
maxssts=max(max(max(ssts)))
  sst=reshape(ssts,lonx,latx,timx);
end  % end of montrend if

%  OK if this has worked we have sst(lonx*latx,timx) with seasons and trend removed and ready for something
%  Next thing would be filterin, then seasonal stratification
%  ******************************************************************************************************************
  % OK this is where we are going to start with our SAM compositing  We will want to cycle over all months starting in 1979, 
% ending in 2022  43 years
% read in SAM index 
%make the index
  % j=(lat-lat(1))/dlat + 1  dlat=lat(2)-lat(1)
  lat1=-35.; lat2=-45.; lat3=-55.; lat4=-65.;
dlon=lon(2)-lon(1);
lon1=150.; lon2=290.;
lon1=0.; lon2=359.;
lon1=lon(1); lon2=lon(length(lon));
	     k1=int32(lon1/dlon+1);
k2=int32(lon2/dlon+1);
dlat=lat(2)-lat(1);
j1=(lat1-lat(1))/dlat+1;
j2=(lat2-lat(1))/dlat+1;
j3=(lat3-lat(1))/dlat+1;
j4=(lat4-lat(1))/dlat+1;
sstz=squeeze(nanmean(sst(k1:k2,:,:),1));
x1=squeeze(nanmean(sstz(j3:j4,:),1));
x2=squeeze(nanmean(sstz(j1:j2,:),1));
  xc=x1-x2;
%  xc=x1;
max(xc)
%xc=xc-nanmean(xc);
%xc=xc/std(xc);

fileout=['SAM5-59',var,level,num2str(lon1),'-',num2str(lon2),'.mat']
  save (fileout,'lat1','lat2','lat3','lat4','xc','level','lon1','lon2')
% Need to normalized xc
    xc=xc-mean(xc);
max(xc)
% need to remove seasonal cycle from index
for i=1:12
	ann(i)=nanmean(xc(i:12:12*yrmx));
	end
	xc=reshape(xc,12,yrmx);
for i=1:yrmx
	xc(:,i)=xc(:,i)-ann';
end
xc=reshape(xc,12*yrmx,1);
max(xc)

xc=xc/std(xc,0);
%xc=detrend(squeeze(xc));
imx=yrmx*12;
sdcrit= 1.0;
%make storage arrays
n=zeros(2,1);
xm=zeros(lonx*latx,2);
yrmx=40
  im1=1; im2=12;
%  im1=5; im2=11;
if im2>12
yrmx=yrmx-1;
end
for iy=1:yrmx
  for imo=im1:im2  %  Loop on 480 months from 1979-2018
	    im=(iy-1)*12+imo;
if xc(im)>sdcrit
xm(:,1)=xm(:,1)+ssts(:,im);
n(1)=n(1)+1;
elseif xc(im)<-sdcrit
xm(:,2)=xm(:,2)+ssts(:,im);
n(2)=n(2)+1;
end % if on sdcrit
  end  %imx
end %yrmx  ******

% OK now take difference and plot
xmd=xm(:,1)/n(1)-xm(:,2)/n(2);
xmp=reshape(xmd,lonx,latx);


%  Let's plot the variance of SST, Just for fun
    sstv=nanstd(ssts,0,2);
    maxx=squeeze(max(max(abs(sstv))))
%  maxx=2.0
    sstd=reshape(sstv, lonx,latx);
	    sstdp=cat(1,sstd,sstd(1,:));
cint=maxx/10;
     conts=[cint:cint:maxx];
     maxe=maxx*0.8;
pslat1=-90.; pslat2=90.0;
FH = figure;     %  Try doing a polar stereographic projection

ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]);

%plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4);

%contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50);

[cmap]=buildcmap('wr');
colormap(cmap);
[C, H] =   contourfm(lat,lonp,sstdp',conts,'linewidth',0.5);
caxis([0 maxe]);
%,'LineWidth',0.5,'Fill','on','LineColor','black');
%h5 = clabelm(C,H); %,'Rotation',0);  %,'LabelSpacing',10);
hold on
%plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3);
geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black')
    set(gca,'FontSize',12);
colorbar('Location','SouthOutside');

ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]);

title(['STD SST ',num2str(yrbeg),'_, Season ',num2str(dmo1),'-',num2str(dm02),', ',num2str(yrbeg,4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)]);
geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black')

set(gca,'Color','none')
hold off
saveas(gca,['SSTSIG_G',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc']);
export_fig ['SSTSIG_G',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc'] -transparent

% now plot the SAM Composite difference
    maxx=squeeze(max(max(abs(xmp))))
  cint=maxx/5;
%  maxx=2.0
	    sstdp=cat(1,xmp,xmp(1,:));
lon=double(lon);
dlon=lon(2)-lon(1);
lonp=cat(1,lon,lon(lonx)+dlon);
lat=double(lat);
cint=maxx/5.;
     conts=[-maxx:cint:maxx];
maxe=maxx;

pslat1=-90.; pslat2=90.0;
FH = figure;     %  Try doing a polar stereographic projection

ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]);


[cmap]=buildcmap('bwr');
colormap(cmap);
[C, H] =   contourfm(lat,lonp,sstdp',conts,'linewidth',0.5);
caxis([-maxe maxe]);

hold on

geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black')
    set(gca,'FontSize',12);
colorbar('Location','SouthOutside');

ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]);

title(['SST SAM',num2str(yrbeg),'_, Season ',num2str(im1),'-',num2str(im2),', ',num2str(yrbeg,4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)]);
geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black')

set(gca,'Color','none')
hold off
saveas(gca,['SSTSAM_G',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc']);
export_fig ['SSTSAM_G',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc'] -transparent;



pause


pslat1=-90.; pslat2=90.0;
FH = figure;     %  Try doing a polar stereographic projection

ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]);

%plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4);

%contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50);

[cmap]=buildcmap('wr');
colormap(cmap);
[C, H] =   contourfm(lat,lonp,sstdp',conts,'linewidth',0.5);
caxis([0 maxe]);
%,'LineWidth',0.5,'Fill','on','LineColor','black');
%h5 = clabelm(C,H); %,'Rotation',0);  %,'LabelSpacing',10);
hold on
%plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3);
geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black')
    set(gca,'FontSize',12);
colorbar('Location','SouthOutside');

ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]);

title(['STD SST ',num2str(yrbeg),'_, Season ',num2str(dmo1),'-',num2str(dm02),', ',num2str(yrbeg,4),'Filt=',num2str(jfilt),num2str(icor),' Lag=',num2str(ilag)]);
geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black')

set(gca,'Color','none')
hold off
saveas(gca,['SSTSIG_G',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc']);
export_fig ['SSTSIG_G',num2str(yrbeg),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc'] -transparent



%  *******   ********    *************   *********  *********
if jfilt == 1
% Next let's low-pass filter the data
N=8
cut=2.0/fcut;    % make cut at 2/Tscale , so that cut is at Tscale  30 day cut here
[b,a]=butter(N,cut);
[h,w]=freqz(b,a,128);
r=abs(h);
rlow=r.*r;
f=w/pi/2.;
%figure
%plot(f,r,':r',f,rlow)
%title(['Butterworth Response: N= ',num2str(N),' cut= ',num2str(cut)])

%xlabel('Cycles per time step')
%ylabel('Response Function')
%drawnow

%  OK, now we will try to filter the data sst

for i = 1:lonx*latx
sst90(i,:)=filtfilt(b,a,ssts(i,:));
end

sst=sst90;
  clear sst90;
	clear ssts;
  else     %  this else is applied for jfilt not= 1 i.e.0
    sst=ssts;
	clear ssts;
end  %  end of filter if loop *******************************

%  NOW seasonal subsetting


		       delmo=dmo2-dmo1;
if dmo2 > 12
	   yrmx=yrmx-1;
end
	   im1=1; im2=dmo2-dmo1+1;
imo1=dmo1; imo2=dmo2;

	     spcdim=int32(lonx*latx);
timdim=int32(yrmx*(delmo+1))
  z=zeros(spcdim,timdim,'double');
  years=zeros(timdim,1,'double');
	   for iy = 1:yrmx
  z(:,im1:im2)=sst(:,imo1:imo2);
  years(im1:im2)=year(imo1:imo2);
im1=im2+1;  im2=im1+delmo;
	      imo1=imo1+12;  imo2=imo2+12;
end

	%  OK seasonal subset is in z(lonx*latx,??)
	% Next do sqrt(cos) weighting of anomalies
		    sqc=sqrt(cos(lat*(pi/180.)));
			     timxz=size(z,2)
			     sst=reshape(z,lonx,latx,timxz);
			     for il = 1:latx
			     sst(:,il,:)=sst(:,il,:)*sqc(il);
end

%  OK NOW I think we are ready to select regional subsets for EOF Analysis
mask = zeros(lonx,latx);
ind=1
			     for i=ilon1:ilon2
			     for j=ilat1:ilat2
			     if topoz(i,j) < 1.0e-5
			     ss(ind,:)=sst(i,j,:);
mask(i,j)=1.0;
ind=ind+1;
end
end
end
ind
if ilag == 0
figure
contour(lon,lat,mask')
hold on
contour(lon,lat,topoz',[0 0],'Color','b')
title('selection mask')
maxss=max(max(max(ss)))
ss(isnan(ss))=0.0;
maxsss=max(max(max(ss)))
end
%  now ss includes the region I want
%  Hey, before we do SVD, let's estimate the autocorrelation of the data
	NT = size(ss,2)
	aa = diag(ss*ss')/NT;   % this is variance at each point
        ab = diag(ss(:,1:NT-1)*ss(:,2:NT)')/(NT-1);   % one lag covariance
        ac=ab./aa;
        autt = nanmean(ac.*aa)/nanmean(aa);
	dof=NT*(1-autt*autt)/(1+autt*autt)
			     [U S V] = svd(ss);
			     sd=diag(S);
%  Get the sign right for 1900 to be consisten as El Nino
if yrbeg == 1900
V(:,1)=-V(:,1);
end

sd=sd.*sd/double(ind);
sd=sd/squeeze(sum(sd))*100;
	factor=sqrt(2./dof);
	sde=sd*factor;
	neofs = size(ss,1);
	indx=[1:neofs];
figure
	kk=12;
	errorbar(indx(1:kk),sd(1:kk),sde(1:kk),'LineWidth',1.5)
%	     plot(sd(1:10))
title('Eigenvalue Spectrum','FontSize',16)
ylabel('Percent of Variance','FontSize',16)
xlabel('Mode Number','FontSize',16)
set(gca,'FontSize',16)
set(gca,'Color','none')
saveas(gca,['EigenSpec_',expname,num2str(xlat1),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc']);
export_fig ['EigenSpec_',expname,num2str(xlat1),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc'] -transparent



%  OK, time for regression
% loop on eigenvalues
%year=([1:timxz])/12. + 1900.;

for iv = 1:ivmx      %  loop on eigenmodes  *** *****   ****   ***   ***   ***

pc=V(:,iv);
pc=(pc-mean(pc))/std(pc);
if yrbeg == 1900
if xlat2 == -70.

	pc=-pc;

end
end
if yrbeg == 1900
if xlat2 == -30.
if iv < 3
	pc=-pc;
end
end
end
if yrbeg == 1900
if xlat2 == -30.
if iv == 3
	pc=-pc;
end
end
end
if yrbeg == 1900
if xlat2 == -30
if iv ==3
	pc=-pc;
end
end
end
if yrbeg == 1880
if xlat2 == -20
if iv == 1
	pc=-pc;
end
end
end


if yrbeg == 1900
if xlat2 == -20
if iv < 5
	pc=-pc;
end
end
end
if yrbeg == 1880
if xlat2 == 20
if iv == 2
	pc=-pc;
end
end
end
if yrbeg == 1880
if detrend == 0
if iv < 3
	pc=-pc;
end
end
end

if xlat2 == -30
if yrbeg == 1900
if iv <3
pc=-pc;
end
 else
if iv == 1
	pc=-pc;
end
end
end
if xlat2 == -70.
	if iv < 1
	pc=-pc;
end
if iv > 2
	pc=-pc;
end
end

if iv == 1
if yrbeg == 1900
pc=pc;
end
if yrbeg == 1900
pc=pc;
end
end
if iv == 2
if yrbeg == 1900
pc=-pc;
end
end
if yrbeg == 1900
pc=-pc;
end
if xlat1 == 30.
if iv < 3
	pc=-pc;
end
end
if xlat2 == -70.
if iv == 2
	pc=-pc;
end
end

if ilag == 0
% Let's add some code to filter the pc time series
Wn=1.0/36.;
nf = 7;
[b,a] = butter(nf,Wn);
[H,W] = freqz(b,a,100);
%figure
%plot(W,H);
pcf = filtfilt(b,a,pc);
figure
ykmx=length(years);
ykmx5=ykmx-11;     % how do I knnow how many to cut off??  I am assuming data through Jan. 2015
	plot(years(1:ykmx5),pc(1:ykmx5),'b',years(1:ykmx5),pcf(1:ykmx5),'r');
hold on
	plot(years(1:ykmx5),pcf(1:ykmx5),'r','LineWidth',1.5);

title(['G_1900_V(:,1)',' PC#',num2str(iv)])
	set(gca,'FontSize',14);
	ylabel('PC Amplitude');
	xlabel('Year')
	axis([yrbeg 2020 -3 4]);

pbaspect([2 1 1]);

set(gca,'Color','none')
		  saveas(gca,['PC_G_',num2str(yrbeg),'_TimeSeries_',num2str(iv),'_',num2str(yrbeg),'.epsc']);
export_fig ['PC_G_',num2str(yrbeg),'_TimeSeries_',num2str(iv),'_',num2str(yrbeg),'.epsc'] -transparent
figure
ykmxm=ykmx-72;
	plot(years(ykmxm:ykmx5),pc(ykmxm:ykmx5));
		  title(['G_',num2str(yrbeg),'_V(:,1)',' PC#',num2str(iv)])
	set(gca,'FontSize',20);
	ylabel('PC Amplitude');
	xlabel('Year')
	axis([2010 2016 -3 4]);
pbaspect([2 1 1]);

end  %  end of plot if ilag

if ilag == 0
	pc=pc';
	fileout=['PDO_',expname,'_',num2str(yrbeg),'_NoTrend',num2str(iv),'.mat']
	save (fileout,'pc','year');


%  OK we can now compute the power spectra of the principle components, if we want
  m=256;  n=m/2;  m2m=int32(n/3); m2=m/2;
  h = spectrum.welch('hamming',m,50);
P1 = pwelch(pc,m,n,m);
f=(0:m2)/m;
f=f*12;
figure
plot(f(1:m2m),P1(1:m2m),'Color','r','LineWidth',1.5)
legend('Power Spectrum')
	title(['Power Spectrum','-',num2str(m),'-',num2str(yrbeg,4),' ',num2str(iv),])
xlabel('Frequency in Cycles per Year','FontSize',14)
ylabel('Power','FontSize',14)
set(gca,'FontSize',14)
%xtick=[0:0.2:2];
%ytick=[0:1:3];

%set(gca,'Xtick',xtick);
%set(gca,'Ytick',ytick);

set(gca,'Color','none')
saveas(gca,['PDO_',expname,num2str(yrbeg,4),'_Power_',num2str(iv),'_',num2str(m),'_',num2str(yrbeg),'.epsc'])
export_fig ['PDO_',expname,num2str(yrbeg,4),'_Power_',num2str(iv),'_',num2str(m),'_',num2str(yrbeg),'.epsc'] -transparent
end  %  of plot ilag if

%  regression
sstr = reshape(sst,lonx*latx,timxz);
if icor == 1    % make regression maps instead of regression maps
% remove standard deviation from sstr
sstrm=nanmean(sstr,2);
sstrd=nanstd(sstr,1,2);
for i = 1:timxz
sstra(:,i)=(sstr(:,i)-sstrm)./sstrd;
end
sstr=sstra;
end   %  end of correlation loop
		  %  OK, here we have an option to do lagged regression of SST on iTself
							NT=double(timxz);
if ilag == 0
regr=sstr*pc'/NT;
elseif ilag > 0
clear sstrs;
clear pcs;
NTM = NT -ilag;
sstrs(:,1:NTM)=sstr(:,ilag+1:NT);
pcs(1:NTM)=pc(1:NTM);
regr=sstrs*pcs'/NTM;
elseif ilag < 0
clear sstrs;
clear pcs;
jlag = abs(ilag)
NTM = NT -jlag;
sstrs(:,1:NTM)=sstr(:,1:NTM);
pcs(1:NTM)=pc(jlag+1:NT);
regr=sstrs*pcs'/NTM;
end  %  end of LAG decision loop

reg = reshape(regr,lonx,latx);

if iv == 1
  maxx=max(max(abs(reg)))*1.05
cint=maxx/10.
maxx=1.3
cint=0.1
  contm=[-maxx:cint:-cint,cint:cint:maxx];
%  contm=[-maxx:cint:maxx];
maxe=maxx*0.7;
conts=contm;
if icor == 1
	conts = [-1:0.1:1.0];
		  maxe=0.7;
end
end   % end of icor==1 loop
%  maxx=max(max(abs(reg)))
[cmap]=buildcmap('bwr');
colormap(cmap)
%need to expand plotted arrays
regp=cat(1,reg,reg(1,:));

%  Marc,
%  The code that makes the figures is below.  I think you can save either the data
%  or the located data on the plot.  Maybe just saving the data is easier.
%  I guess you just need the lat, long and value
%  Thanks for taking a look at this.

pslat1=-90.; pslat2=90.0;
FH =figure;     %  Try doing a polar stereographic projection

ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]);

%plotm(coast.lat,coast.long,'Color','k','LineWidth',0.4);

%contourm(lat,lonp,regp',contm,'linecolor','magenta','LineWidth',0.50);

[cmap]=buildcmap('bwr');
colormap(cmap);

[C, H] =   contourfm(lat,lonp,regp',conts,'linewidth',0.5);
[C, H] =   contourm(lat,lonp,regp',[0. 0.],'linewidth',2.0);
caxis([-maxe*1.04 maxe]);
%,'LineWidth',0.5,'Fill','on','LineColor','black');
%h5 = clabelm(C,H); %,'Rotation',0);  %,'LabelSpacing',10);
hold on
%plotm(coast.lat,coast.long,'Color',[0.5 0.5 0.5],'LineWidth',0.3);
geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black')

    set(gca,'FontSize',12);

 
ax = axesm('MapProjection','robinson','Grid','on','Frame','on','FLineWidth',1.0,'MapLatLimit',[pslat1 pslat2],'MapLonLimit',[30 390]);

title(['PDO Regression',expname,', Season ',num2str(dmo1),'-',num2str(dm02),' mode=',num2str(iv),', ',num2str(yrbeg,4),'Filt=',num2str(jfilt),num2str(icor),num2str(gmremove),' Lag=',num2str(ilag)]);
geoshow(coast.lat,coast.long,'DisplayType','polygon','FaceColor','black')

set(gca,'Color','none')
hold off
saveas(gca,['PDO_',expname,'_',num2str(iv),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc']);
export_fig ['PDO_',expname,'_',num2str(iv),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.epsc'] -transparent

R = georasterref('RasterSize',size(regp'),'Latitudelimits',[-90 90],'LongitudeLimits',[0 360]);
geotiffwrite(['PDO_Atlantic_65N-65S_1900_',num2str(iv),'_',num2str(jfilt),'_',num2str(yrbeg,4),'_',num2str(dmo1),'-',num2str(dm02),'-',num2str(icor),num2str(gmremove),num2str(ilag),'.tif'],regp',R);




end   %  end of eigenvalue loop