% Script to build a process model for Pacific Climate

% dCHL/dt = a*AL + b*NPO + c*EPW + d*CPW -g*CHL

% or equivalent to the discretize form:
% CHL(t+1)= a*AL(t) + b*NPO(t) + c*EPW(t) + d*CPW(t) +e*CHL(t)
% e=(1-g)

% STEP 1: Generate the forcing functions for the process model
% EPW/CPW as the EOFs 1 and 2 of SSTa [15S-15N]   (epw=eastern pacific
% warming)
pac.xlim=[-266  -64];
pac.ylim=[-15 15];          %tropics
trange=[datenum(1950, 1,15) datenum(2008, 12, 15)];
tmp=rnc_Extract_NOAA_SST(pac.xlim,pac.ylim,trange,0);  %tmp=temporary
tmp = rnc_forcdStats(tmp,'SST',  [1950 2008]);
%[ano, years]=rnc_yearavg(tmp.datenum, tmp.ano);
tmp.ano=rnt_filter3d(tmp.ano,tmp.mask, 3, 'blackman'); tmp.ano(isnan(tmp.ano))=0;
in=find(tmp.month==2);
[eoft,pct,vexpt]=rnt_doEof(tmp.ano(:,:,in), tmp.mask);
% Assign the forcing functions
f.year=tmp.year(in);
f.epw=nn(pct(:,1));     %epw- 1st eof
f.cpw=nn(pct(:,2));     %cpw- 2nd eof

% AL/NPO     (Aleutian Low/N Pacific Oscilation)
pac.xlim=[-245  -110];
pac.ylim=[20 70]; %Change for latitudes
trange=[datenum(1950, 1,15) datenum(2008, 12, 15)];
tmp=rnc_Extract_NCEP_SLP(pac.xlim,pac.ylim,trange,0);
tmp = rnc_forcdStats(tmp,'SLP',  [1950 2008]);
%[ano, years]=rnc_yearavg(tmp.datenum, tmp.ano);  
%using a 3 month running avg, blackman style-weihted toward the
%point(triangle)
tmp.ano=rnt_filter3d(tmp.ano,tmp.mask, 3, 'blackman'); tmp.ano(isnan(tmp.ano))=0;
in=find(tmp.month==2);
[eofa,pca,vexpa]=rnt_doEof(tmp.ano(:,:,in), tmp.mask);
clf;rnc_map(eofa(:,:,2), tmp);
% Assign the forcing functions
f.al=nn(pca(:,1));      %al- 1st eof
f.npo=nn(pca(:,2));     %npo- 2nd eof

% STEP 2: Use model to hindcast NP(north pacific) SSTa
pac.xlim=[-245  -110];
pac.ylim=[20 70]; %Change for latitudes
trange=[datenum(1950, 1,15) datenum(2008, 12, 15)];
tmp=rnc_Extract_NOAA_SST(pac.xlim,pac.ylim,trange,0);
sstnp = rnc_forcdStats(tmp,'SST',  [1950 2008]);
%[ano, years]=rnc_yearavg(sstnp.datenum, sstnp.ano);  
sstnp.ano=rnt_filter3d(sstnp.ano,sstnp.mask, 3, 'blackman'); sstnp.ano(isnan(sstnp.ano))=0;
in=find(sstnp.month==2);

sstnp.an.time=sstnp.year(in);
sstnp.an.ano=sstnp.ano(:,:,in);

%choosing one point to compare model w/ noaa data of ssta
i=56; j=16;  i=40; j=13 
i=60; j=13;  i=60; j=25
clf; rnc_map(sstnp.mean, sstnp);
plot(sstnp.lon(i,j), sstnp.lat(i,j), '*b');

data=nn(sq(sstnp.an.ano(i,j,:)));
% y=F*beta
% F'*y=F'*F*beta; CME= F'*F;
% F'*y=CME*beta; beta=inv(CME)*F'*y;

y=data(2:end,1);
F(:,1)=f.al(1:end-1);         %forcings
F(:,2)=f.npo(1:end-1);
F(:,3)=f.epw(1:end-1);
F(:,4)=f.cpw(1:end-1);
F(:,5)=data(1:end-1);    %memory of previous timestep(like 2 years ago)

CME= F'*F;
beta=inv(CME)*F'*y;

a=beta(1); b=beta(2); c=beta(3); d=beta(4);
e=beta(5);
g=1-e; tao=1/g;

yhat=data(1);
for t=1:length(data)-1;
 yhat(t+1)= a*f.al(t) + b*f.npo(t) + c*f.epw(t) + d*f.cpw(t) +e*yhat(t);
end
yhat=nn(yhat);
red_signi(f.year, data, f.year, yhat, 1, 1500);


%Mapping the Correlations
i=40; j=13   i=40; j=25
i=60; j=13;  i=60; j=25

[I,J,K]=size(sstnp.ano);
map=zeros(I,J)*nan; CORR=map; CORR_LP=map;
A=map; B=map; C=map; D=map; E=map;
for i=1:I;
    for j=1:J
        data=nn(sq(sstnp.an.ano(i,j,:)));
        if ~isnan(data(1))
            y=data(2:end,1);
            F(:,1)=f.al(1:end-1);
            F(:,2)=f.npo(1:end-1);
            F(:,3)=f.epw(1:end-1);
            F(:,4)=f.cpw(1:end-1);
            F(:,5)=data(1:end-1);
            CME= F'*F;
            beta=inv(CME)*F'*y;

            a=beta(1); b=beta(2); c=beta(3); d=beta(4);
            e=beta(5);
            g=1-e; tao=1/g;
            A(i,j)=a; B(i,j)=b; C(i,j)=c; D(i,j)=d; E(i,j)=e;
            yhat=data(1);
            for t=1:length(data)-1;
                %yhat(t+1)= a*f.al(t) + b*f.npo(t) + c*f.epw(t) + d*f.cpw(t) +e*yhat(t);
                yhat(t+1)=  c*f.epw(t) + d*f.cpw(t) +e*yhat(t);
            end
            yhat=nn(yhat);
            corr=correlation(f.year, data, f.year, yhat, 1);
            CORR(i,j)=corr;
            CORR_LP(i,j)=correlation(f.year, lowpass(data,1/7), f.year, lowpass(yhat',1/7), 1);
        else
            map(i,j)=nan;
        end
    end
end
clf;rnc_map(map, sstnp)