function SST_prediction(startmonth)

%Various info 
OBS_ENS=10;
dlon=360;
dlat=180;
nx=320;
ny=384;
year_start=1985
year_end=2009
nb_yr=year_end-year_start + 1;
O_PATH='/work/shared/nn9039k/NorCPM/Obs/SST/HADISST2/';
FREE_PATH='/work/ywang/NorESM-ME/Processed/';
ASSIM_PATH = '/work/ywang/yiguo_hindcasts/Processed/';
ASSIM_PATH0 = '/work/ywang/ana_19800115_me_hindcasts/Processed/';
Input_PATH = '/work/shared/nn9039k/NorCPM/Input/NorESM/f19_g16/HADISST2/';
grid_PATH='/work/shared/noresm/inputdata/ocn/micom/gx1v6/20101119/grid.nc';

%reading pivot
ipiv = ncgetvar([Input_PATH 'pivots_SST.nc'], 'ipiv');
jpiv = ncgetvar([Input_PATH 'pivots_SST.nc'], 'jpiv');
for j = 1:dlat
      lat(j) = 90 - j;
end
for i = 1:dlon
   lon(i) = -180 + (i - 1);
end

%spatial distance
for j = 2:dlat-1
   dx = Haversin_dist(lon(1), lat(j), lon(3), lat(j)) / 2;
   dy = Haversin_dist(lon(1), lat(j-1), lon(1), lat(j+1)) / 2;
   parea(1:dlon, j) = dx * dy;
end
parea(1:dlon, 1) = parea(1:dlon, 2);
parea(1:dlon,180) = parea(1:dlon,179);
plon=repmat(lon', 1, dlat);
plat=repmat(lat, dlon, 1);

% Nino 3.4 [5S-5N, 170W-120W]
latmin=86
latmax=96
lonmin=11
lonmax=61
paream=nansum(nansum(parea(lonmin:lonmax,latmin:latmax)));

if exist(['SST_prediction_' num2str(startmonth,'%2.2d') '.mat'])
  load(['SST_prediction_' num2str(startmonth,'%2.2d') '.mat']);
else
  cnt = 1;

  fsst = zeros(dlon,dlat,nb_yr*4);
  asst = zeros(dlon,dlat,nb_yr*4);
  asst0 = zeros(dlon,dlat,nb_yr*4);
  for yr = year_start : year_end
   yr
      for mm = 1:3:12
         yr1=yr + floor((mm+startmonth-1)./12);
    	 mm1=mod((mm+startmonth-1),12)+1;

         %Data
         d_sst = squeeze(nanmean(ncgetvar([O_PATH num2str(yr1) '_' num2str(mm1,'%2.2d') '.nc'],'sst'),4)) - 273.15;
         o_avg = squeeze(nanmean(ncgetvar([O_PATH 'Anomaly/SST_avg_' num2str(mm1,'%2.2d') '-1980-2010.nc'],'sst'),4)) - 273.15;
         d_sst = d_sst - o_avg;

         filename=['ocn' num2str(yr) num2str(mm,'%2.2d') '_' num2str(yr1) num2str(mm1,'%2.2d') '.nc'];

         %reading free
         f_sst = ncgetvar([FREE_PATH 'free-average' num2str(yr1) '-' num2str(mm1,'%2.2d') '.nc'], 'sst');
         avg_sst_assim = ncgetvar([Input_PATH '../Free-average' num2str(mm1,'%2.2d') '-1980-2010.nc'],'sst');
         f_sst = f_sst - avg_sst_assim;

         %reading Assim
         a_sst = ncgetvar([ASSIM_PATH  filename], 'sst');
         a_sst = a_sst - avg_sst_assim;
%         a_sic = ncgetvar([ASSIM_PATH filename], 'fice');
%         mask = find(a_sic ~= 0);
%         a_sst(mask) = nan;

         %reading Assim
         a_sst0 = ncgetvar([ASSIM_PATH0 filename], 'sst');
         a_sst0 = a_sst0 - avg_sst_assim;
%         a_sic = ncgetvar([ASSIM_PATH0 filename], 'fice');
%         mask = find(a_sic ~= 0);
%         a_sst0(mask) = nan;

      %interp_model to obs
      f_osst = zeros(dlon,dlat);
      a_osst = zeros(dlon,dlat);
      a_osst0 = zeros(dlon,dlat);
      nb_aosst = zeros(dlon,dlat);

      for i = 1:dlon
         for j = 1:dlat
            if ipiv(i,j) ~=0 & jpiv(i,j) ~=0
               f_osst(i,j) = f_osst(i,j) + f_sst(ipiv(i,j),jpiv(i,j));
               a_osst(i,j) = a_osst(i,j) + a_sst(ipiv(i,j),jpiv(i,j));
               a_osst0(i,j) = a_osst0(i,j) + a_sst0(ipiv(i,j),jpiv(i,j));
               nb_aosst(i,j) = nb_aosst(i,j) + 1;
            end
         end
      end
     fsst(:,:,cnt) = fsst(:,:,cnt) + f_osst./nb_aosst;
     asst(:,:,cnt) = asst(:,:,cnt) + a_osst./nb_aosst;
     asst0(:,:,cnt) = asst0(:,:,cnt) + a_osst0./nb_aosst;
    osst(:,:,cnt) = d_sst;
    cnt = cnt + 1;
  end
end
cnt = cnt - 1;
disp(cnt)

for i = 1:dlon
      for j = 1:dlat
         ssta.R_fd(i,j)=corr2(squeeze(detrend(fsst(i,j,:))),squeeze(detrend(osst(i,j,:))));
         ssta.R_f(i,j)=corr2(squeeze(fsst(i,j,:)),squeeze(osst(i,j,:)));
         ssta.R_ad(i,j)=corr2(squeeze(detrend(asst(i,j,:))),squeeze(detrend(osst(i,j,:))));
         ssta.R_a(i,j)=corr2(squeeze(asst(i,j,:)),squeeze(osst(i,j,:)));
         ssta.R_ad0(i,j)=corr2(squeeze(detrend(asst0(i,j,:))),squeeze(detrend(osst(i,j,:))));
         ssta.R_a0(i,j)=corr2(squeeze(asst0(i,j,:)),squeeze(osst(i,j,:)));
      end
end

for i = 1:cnt
   anino(i) = nansum(nansum(asst(lonmin:lonmax,latmin:latmax,i).*parea(lonmin:lonmax,latmin:latmax)))/paream;   
   anino0(i) = nansum(nansum(asst0(lonmin:lonmax,latmin:latmax,i).*parea(lonmin:lonmax,latmin:latmax)))/paream;
   onino(i) = nansum(nansum(osst(lonmin:lonmax,latmin:latmax,i).*parea(lonmin:lonmax,latmin:latmax)))/paream;
end
r=corr2(onino,anino)
r0=corr2(onino,anino0)

% save data
save(['SST_prediction_' num2str(startmonth,'%2.2d') '.mat'],'ssta');
end

%%%%%%%%%%%%
%%%%Plot%%%%
%%%%%%%%%%%%
close all
figure(1)
set(gcf, 'Renderer', 'opengl');
set(gcf, 'InvertHardCopy', 'off');
set(gca,'XAxisLocation','top','fontsize',16);
whitebg('w');
%m_proj('equidistant','lon',[-70 30],'lat',[36 82])
m_proj('stereographic','lon',-47,'lat',80.0,'rad',42,'rec','on')
P=m_pcolor(plon, plat, ssta.R_a)
set(P,'LineStyle','none');
hold on
title(['SSTA corr @' num2str(startmonth,'%2.2d') ' months lead time'],'FontSize', 22);
%ylabel(' ','FontSize', 28);
caxis([-1 1]);
colormap(anomisuperwide(20));
m_grid;
c=colorbar('h','FontSize', 16);
m_coast('patch',[0.5 0.5 0.5]);
%xlabel(c,'%','FontSize', 22);
print('-depsc', ['SST_prediction_' num2str(startmonth,'%2.2d_SSTA_TSA') '.eps'])

figure(2)
set(gcf, 'Renderer', 'opengl');
set(gcf, 'InvertHardCopy', 'off');
set(gca,'XAxisLocation','top','fontsize',16);
whitebg('w');
%m_proj('equidistant','lon',[-70 30],'lat',[36 82])
%m_proj('stereographic','lon',-10,'lat',65.0,'rad',20,'rec','on')
m_proj('stereographic','lon',-47,'lat',80.0,'rad',42,'rec','on')
P=m_pcolor(plon, plat, ssta.R_a0)
set(P,'LineStyle','none');
hold on
P=m_pcolor(plon-360,plat, ssta.R_a0);
set(P,'LineStyle','none');
title(['SSTA corr @' num2str(startmonth,'%2.2d') ' months lead time'],'FontSize', 22);
%ylabel(' ','FontSize', 28);
caxis([-1 1]);
colormap(anomisuperwide(20));
m_grid;
m_coast('patch',[0.5 0.5 0.5]);
c=colorbar('h','FontSize', 16);
%xlabel(c,'%','FontSize', 22);
print('-depsc', ['SST_prediction_' num2str(startmonth,'%2.2d_SSTA') '.eps'])
end
