%load up the precip from the LGM/4XCO_2/ Jolocene and Hosing experiments -- all on the common grid runs and look at the annual mean shifts as
%a function of latitude

%use this one to visualize the optimization
clear all; close all;
IT=5; %number of interations of the optimization to perform
LW=2.1; FSz=24; FN='Arial';
 LWA=2.5; LWM = 3.5; trans=.2; LWC=3.7; LW2=2.3; LW3 =2; LWV=1.5;
gray=[0.4 0.4 0.4];
green=[0.1 0.4 0.1];
lgreen= [127 255 0]/255;
dblue=[ 0 0 .5];
 cmap=colormap('jet'); NC=length(cmap);
 
cols = cmap(roundup( NC*(1:IT)/IT),:); %discrete color for each iteration
cbar_ytic=(1:IT)/IT;
flat=-20:.01:20;
%load the data and compose a single matrix of all the experiments --
%P_perturb and P_pi for the contro run
load('data\PRECIP_LGM_common', 'lat', 'lon', 'P_PI', 'P_LGM','MODS')
[NM NT NN N]=size(P_LGM);
do_lgm=1:NM;
P_perturb=P_LGM;
P_pi=P_PI;
load('data\PRECIP_4X_common', 'lat', 'lon', 'P_PI', 'P_4X','MODS')
[NM NT NN N]=size(P_4X);
start=7+1; %start the model count with the end of the LGM
do_4x=start:(NM+start-1);
P_perturb(start:NM+start-1,1:12,1:NN,1:N)=P_4X;
P_pi(start:NM+start-1,1:12,1:NN,1:N)=P_PI;

load('data\PRECIP_HOS_common', 'lat', 'lon', 'P_PI', 'P_HOS','MODS')
[NM NT NN N]=size(P_HOS);
start=23+1; %start the model count with the end of the LGM
do_hos=start:(NM+start-1);
P_perturb(start:NM+start-1,1:12,1:NN,1:N)=P_HOS;
P_pi(start:NM+start-1,1:12,1:NN,1:N)=P_PI;


%%%%%%%%%%%%%%%%%%%%%%% and the holocene runs
load('data\P_HOL', 'P_PI', 'P_HOL', 'lat', 'lon', 'MODS')
[NM NT NN N]=size(P_HOL);
start=30+1; %start the model count with the end of the LGM
do_hol=start:(NM+start-1);
P_perturb(start:NM+start-1,1:12,1:NN,1:N)=P_HOL;
P_pi(start:NM+start-1,1:12,1:NN,1:N)=P_PI;

[NM NT NN N]=size(P_pi);
mlat=lat; %just cause the coast.mat maps have a variable called lat

%I usually just pic a single mode
for nm=31:31;%start+NM-1;
lat=mlat; %in case the coast.mat lat is loaded
    
PC=squeeze(mean(squeeze(P_pi(nm,:,:,:))));
PN=squeeze(mean(squeeze(P_perturb(nm,:,:,:))));

%calculate the precipitation centroid
fP=interp1(lat, nanmean(PC'),flat);
cdf=cumsum(fP)/sum(fP);
PCENT(nm)=flat(find2(cdf,0.5));

fP=interp1(lat, nanmean(PN'),flat);
cdf=cumsum(fP)/sum(fP);
LPCENT(nm)=flat(find2(cdf,0.5));

sp=0.5;
ilat=-20:sp:20; iNN=length(ilat);
ilat_ext=-30:sp:30; iNN2=length(ilat_ext);
ilat_ext2=-45:sp:45; iNN3=length(ilat_ext);
[LON LAT]=meshgrid(lon,lat);
[ILON ILAT]=meshgrid(lon, ilat);
[ILON_ext ILAT_ext]=meshgrid(lon, ilat_ext);
[ILON_ext2 ILAT_ext2]=meshgrid(lon, ilat_ext2);

iPN=interp2(LON, LAT, PN,ILON,ILAT);
iPN_ext=interp2(LON, LAT, PN,ILON_ext,ILAT_ext);
iPN_ext2=interp2(LON, LAT, PN,ILON_ext2,ILAT_ext2);
iZPN=nanmean(iPN'); % the perturbed precip on the interpolation grid
iZPN_ext=nanmean(iPN_ext'); % the perturbed precip on the interpolation grid
iZPN_ext2=nanmean(iPN_ext2'); % the perturbed precip on the interpolation grid

iPC=interp2(LON, LAT, PC,ILON,ILAT);
iPC_ext=interp2(LON, LAT, PC,ILON_ext,ILAT_ext);
iPC_ext2=interp2(LON, LAT, PC,ILON_ext2,ILAT_ext2);
iZPC=nanmean(iPC'); % the perturbed precip on the interpolation grid
iZPC_ext=nanmean(iPC_ext'); % the perturbed precip on the interpolation grid
iZPC_ext2=nanmean(iPC_ext2'); % the perturbed precip on the interpolation grid

 DELP(nm,1:iNN2, 1:N) = iPN_ext-iPC_ext; 
 shift=-2:0.02:2; NS=length(shift);
 int=.80:.008:1.15; NI=length(int);
 con=.85:.002:1.1; NC=length(con);
 %start with zonal mean shifts, expansion and intesnification considered
 %collectively
 
 %start with initial guesses from performing each shift independently
 figure(30); clf; plot(ilat, iZPN, 'k-', 'LineWidth', LWC);
  figure(31); clf; plot(ilat, iZPN, 'k-', 'LineWidth', LWC);
   figure(32); clf; plot(ilat, iZPN, 'k-', 'LineWidth', LWC);
 figure(1); clf;
 
 %%%%%%%%%%%%%%%%%%%%%%% SHIFTING MODE %%%%%%%%%%%%%%%%%%%%%%%
err_ref = (mean((iZPN-iZPC).^2))^.5; %the error in the unchanged precipitation
 sh_old=0;
 con_old=1;
 int_old=1;

 tic
 for i=1:IT; %iterations of the error minimization
     
     
 
    for ns=1:NS;
         ZP_CI = int_old*interp1(con_old*(ilat_ext2), iZPC_ext2, ilat_ext); %the contracted and intensified precipitation from the previous iteration
         ZP_S = interp1((ilat_ext+shift(ns)), ZP_CI, ilat);
        err_s(ns) = mean(squeeze(mean(((iZPN-ZP_S).^2))))^.5;
    end
 
    sh_new=shift(find(err_s == min(err_s)));
    SH_I(nm,i)=sh_new; %save for posteriority
    figure(1); hold on;
    subplot(1,3,1); plot(shift, err_s/err_ref, '-', 'Color', cols(i,1:3),'LineWidth', LW);
    hold on; plot(sh_new*[1 1], [min(err_s) max(err_s)]/err_ref, '--', 'Color', cols(i,1:3),'LineWidth', LWV);
    
    %visualize the initial and optimaly shifted precip fo each iteration
    ZP_S_new=interp1((ilat_ext+sh_new), ZP_CI, ilat); % the optimally shifted precip for this iteration
    figure(30); 
    hold on; plot(ilat_ext, ZP_CI, '-', 'Color', cols(i,1:3),'LineWidth', LW);
    hold on; plot(ilat, ZP_S_new, '--', 'Color', cols(i,1:3),'LineWidth', LW);
    
  %%%%%%%%%%%%%%%%%%%%%%% CONTRACTION MODE %%%%%%%%%%%%%%%%%%%%%%%
    for nc=1:NC;
         ZP_SI = int_old*interp1((ilat_ext2+sh_old), iZPC_ext2, ilat_ext); %the shifted and intensified precipitation from the previous iteration
         ZP_CON = interp1((ilat_ext*con(nc)), ZP_SI, ilat);
         err_c(nc) = mean(squeeze(mean(((iZPN-ZP_CON).^2))))^.5;
    end
 
    con_new=con(find(err_c == min(err_c)));
    CON_I(nm,i)=con_new; %save for posteriority
    figure(1); hold on;
    subplot(1,3,2); plot(con, err_c/err_ref, '-','Color', cols(i,1:3), 'LineWidth', LW);
    hold on; plot(con_new*[1 1], [min(err_c) max(err_c)]/err_ref, '--', 'Color', cols(i,1:3),'LineWidth', LWV)
 
    %visualize the initial and optimaly shifted precip fo each iteration
    ZP_C_new=interp1((con_new*ilat_ext), ZP_SI, ilat); % the optimally shifted precip for this iteration
    figure(31); 
    hold on; plot(ilat_ext, ZP_SI, '-', 'Color', cols(i,1:3),'LineWidth', LW);
    hold on; plot(ilat, ZP_C_new, '--', 'Color', cols(i,1:3),'LineWidth', LW);
    
 %%%%%%%%%%%%%%%%%%%%%%% INTENSIFICATION MODE %%%%%%%%%%%%%%%%%%%%%%%
    for ni=1:NI;
         ZP_SC = interp1(con_old*(ilat_ext2+sh_old), iZPC_ext2, ilat_ext);  %the shifted and contracted precipitation from the previous iteration 
         ZP_INT = int(ni)*interp1((ilat_ext), ZP_SC, ilat);
         err_i(ni) = mean(squeeze(mean(((iZPN-ZP_INT).^2))))^.5;
    end
 
    int_new=int(find(err_i== min(err_i)));
    INT_I(nm,i)=int_new; %save for posteriority
    figure(1); hold on;
    subplot(1,3,3); plot(int, err_i/err_ref, 'Color', cols(i,1:3), 'LineWidth', LW);
    hold on; plot(int_new*[1 1], [min(err_i) max(err_i)]/err_ref, '--','Color', cols(i,1:3), 'LineWidth', LWV)
 

    %visualize the initial and optimaly shifted precip fo each iteration
    ZP_I_new=int_new*interp1((ilat_ext), ZP_SC, ilat); % the optimally intensified precip for this iteration
    figure(32); 
    hold on; plot(ilat_ext, ZP_SC, '-', 'Color', cols(i,1:3),'LineWidth', LW);
    hold on; plot(ilat, ZP_I_new, '--', 'Color', cols(i,1:3),'LineWidth', LW);
    
 %%%%%%%%%%%%%%%%% prep for next interation
 inc=0.5; %the increment of the optimal change to step forward by
 sh_old=sh_old+inc*(sh_new-sh_old);
 con_old=con_old+inc*(con_new-con_old);
 int_old=int_old+inc*(int_new-int_old);    
        
end
             

%graphical stuff
YMX=0.9*max([err_c err_s err_i]/err_ref);
YMN=0.9*min([err_c err_s err_i]/err_ref);
figure(1); hold on;
subplot(1,3,1); hold on; plot([min(shift) max(shift)], [1 1], 'k--', 'LineWidth', LW);
xlabel('shift (^o latitude)', 'FontSize', FSz, 'FontName', FN);
ylabel('Root mean squared precipitation error over the tropics (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('Shift', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
axis([min(shift) max(shift) YMN YMX]);

subplot(1,3,2); hold on; plot([min(con) max(con)], [1 1], 'k--', 'LineWidth', LW);
xlabel('Contraction (scalar multiplier)', 'FontSize', FSz, 'FontName', FN);
ylabel('Root mean squared precipitation error over the tropics (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('Contraction', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
axis([min(con) max(con) YMN YMX]);

subplot(1,3,3); hold on; plot([min(int) max(int)], [1 1], 'k--', 'LineWidth', LW);
xlabel('Intensification (scalar multiplier)', 'FontSize', FSz, 'FontName', FN);
ylabel('Root mean squared precipitation error over the tropics (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('Intesnification', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
axis([min(int) max(int) YMN YMX]);
H=colorbar('peer', gca);
set(get(H,'YLabel'),'string', 'Iteration Number', 'FontSize', FSz, 'FontName', FN);
colorbar('YTick', cbar_ytic, 'YTickLabel', 1:IT, 'LineWidth', LWA);

figure(32); hold on; 
colormap('jet')
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
ylabel('Zonal mean precipitation (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('Intesnification', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
H=colorbar('peer', gca);
set(get(H,'YLabel'),'string', 'Iteration Number', 'FontSize', FSz, 'FontName', FN);
colorbar('YTick', cbar_ytic, 'YTickLabel', 1:IT, 'LineWidth', LWA);
legend('TARGET', 'Initial guess', 'Optimized Intensification')

figure(31); hold on; 
colormap('jet')
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
ylabel('Zonal mean precipitation (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('Contraction', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
H=colorbar('peer', gca);
set(get(H,'YLabel'),'string', 'Iteration Number', 'FontSize', FSz, 'FontName', FN);
colorbar('YTick', cbar_ytic, 'YTickLabel', 1:IT, 'LineWidth', LWA);
legend('TARGET', 'Initial guess', 'Optimized Contraction')

figure(30); hold on; 
colormap('jet')
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
ylabel('Zonal mean precipitation (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('Shift', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
H=colorbar('peer', gca);
set(get(H,'YLabel'),'string', 'Iteration Number', 'FontSize', FSz, 'FontName', FN);
colorbar('YTick', cbar_ytic, 'YTickLabel', 1:IT, 'LineWidth', LWA);
legend('TARGET', 'Initial guess', 'Optimized Shift')

%percent of change explained by all three, and excluding on mechanism
ZDP=iZPN-iZPC;
ZALL_err = int_old*interp1(con_old*(ilat_ext+sh_old), iZPC_ext , ilat) - iZPN; 
ZSHIFT_err = int_old*interp1(con_old*(ilat_ext), iZPC_ext , ilat) - iZPN;
ZCON_err = int_old*interp1((ilat_ext+sh_old), iZPC_ext , ilat) - iZPN;
ZINT_err = interp1(con_old*(ilat_ext+sh_old), iZPC_ext , ilat) - iZPN;

ZPER(nm,1:4) = [1- ( (mean((ZALL_err).^2))^.5) / ((mean((ZDP).^2))^.5);
               1- ( (mean((ZSHIFT_err).^2))^.5) / ((mean((ZDP).^2))^.5);
               1- ( (mean((ZCON_err).^2))^.5) / ((mean((ZDP).^2))^.5);
               1- ( (mean((ZINT_err).^2))^.5) / ((mean((ZDP).^2))^.5)];
           
           
%ALTERNATIVELY  -- error reduction resulting from a single operation from
%the means state
ZSHIFT_err2 = interp1((ilat_ext+sh_old), iZPC_ext , ilat) - iZPN;
ZCON_err2 = interp1(con_old*(ilat_ext), iZPC_ext , ilat) - iZPN;
ZINT_err2 = int_old*interp1((ilat_ext), iZPC_ext , ilat) - iZPN;

ZPER_ASCI2(nm,1:4) = [ 1-rms(ZALL_err)/ rms(ZDP);
                       1-rms(ZSHIFT_err2)/ rms(ZDP);
                       1-rms(ZCON_err2)/ rms(ZDP);
                       1-rms(ZINT_err2)/ rms(ZDP)];
            
 %lastly, the best answer from the first iteration
ZSHIFT_err3 = interp1((ilat_ext+SH_I(nm,1)), iZPC_ext , ilat) - iZPN;
ZCON_err3 = interp1(CON_I(nm,1)*(ilat_ext), iZPC_ext , ilat) - iZPN;
ZINT_err3 = INT_I(nm,1)*interp1((ilat_ext), iZPC_ext , ilat) - iZPN;

ZPER_ASCI3(nm,1:4) = [ 1-rms(ZALL_err)/ rms(ZDP);
                       1-rms(ZSHIFT_err3)/ rms(ZDP);
                       1-rms(ZCON_err3)/ rms(ZDP);
                       1-rms(ZINT_err3)/ rms(ZDP)];
toc

ZPER_ASCI(nm, 1:4) = [ZPER(nm,1);
                      ZPER(nm,1)-ZPER(nm,2);
                      ZPER(nm,1)-ZPER(nm,3);
                      ZPER(nm,1)-ZPER(nm,4)];

zonal_percentages = [nm NM round(100*ZPER_ASCI(nm,1:4));
                     nm NM round(100*ZPER_ASCI2(nm,1:4));
                     nm NM round(100*ZPER_ASCI3(nm,1:4));]           
           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2D
           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           %%%%%%%%%%%%%%%%%%% start from initial guesses above
           
%start with initial guesses from performing each shift independently

 
 %%%%%%%%%%%%%%%%%%%%%%% SHIFTING MODE %%%%%%%%%%%%%%%%%%%%%%%
 err_ref = mean(squeeze(mean((iPN-iPC).^2)))^.5; %the error in the unchanged precipitation
tic 
figure(71); clf;

 for i=1:IT; %iterations of the error minimization

     
 
    for ns=1:NS;
         P_CI = int_old*interp2(ILON_ext2, con_old*(ILAT_ext2), iPC_ext2, ILON_ext, ILAT_ext); %the contracted and intensified precipitation from the previous iteration
         P_S = interp2(ILON_ext, (ILAT_ext+shift(ns)), P_CI, ILON, ILAT);
        err_s(ns) = mean(squeeze(mean(((iPN-P_S).^2))))^.5;
    end
 
    sh_new=shift(find(err_s == min(err_s)));
    SH_I_T(nm,i)=sh_new; %save for posteriority
    figure(71); hold on;
    subplot(1,3,1); plot(shift, err_s/err_ref, '-', 'Color', cols(i,1:3),'LineWidth', LW);
    hold on; plot(sh_new*[1 1], [min(err_s) max(err_s)]/err_ref, '--', 'Color', cols(i,1:3),'LineWidth', LWV);
    
    
  %%%%%%%%%%%%%%%%%%%%%%% CONTRACTION MODE %%%%%%%%%%%%%%%%%%%%%%%
    for nc=1:NC;
         P_SI = int_old*interp2(ILON_ext2,(ILAT_ext2+sh_old), iPC_ext2, ILON_ext, ILAT_ext); %the shifted and intensified precipitation from the previous iteration
         P_CON = interp2(ILON_ext, (ILAT_ext*con(nc)), P_SI, ILON, ILAT);
         err_c(nc) = mean(squeeze(mean(((iPN-P_CON).^2))))^.5;
    end
 
    con_new=con(find(err_c == min(err_c)));
    CON_I_T(nm,i)=con_new; %save for posteriority
    figure(71); hold on;
    subplot(1,3,2); plot(con, err_c/err_ref, '-','Color', cols(i,1:3), 'LineWidth', LW);
    hold on; plot(con_new*[1 1], [min(err_c) max(err_c)]/err_ref, '--', 'Color', cols(i,1:3),'LineWidth', LWV)

 %%%%%%%%%%%%%%%%%%%%%%% INTENSIFICATION MODE %%%%%%%%%%%%%%%%%%%%%%%
    for ni=1:NI;
         P_SC = interp2(ILON_ext2, con_old*(ILAT_ext2+sh_old), iPC_ext2, ILON_ext, ILAT_ext);  %the shifted and contracted precipitation from the previous iteration 
         P_INT = int(ni)*interp2(ILON_ext,(ILAT_ext), P_SC, ILON, ILAT);
         err_i(ni) = mean(squeeze(mean(((iPN-P_INT).^2))))^.5;
    end
 
    int_new=int(find(err_i== min(err_i)));
    INT_I_T(nm,i)=int_new;
    figure(71); hold on;
    subplot(1,3,3); plot(int, err_i/err_ref, 'Color', cols(i,1:3), 'LineWidth', LW);
    hold on; plot(int_new*[1 1], [min(err_i) max(err_i)]/err_ref, '--','Color', cols(i,1:3), 'LineWidth', LWV)
 

 %%%%%%%%%%%%%%%%% prep for next interation
 %the increment of the optimal change to step forward by
 sh_old=sh_old+inc*(sh_new-sh_old);
 con_old=con_old+inc*(con_new-con_old);
 int_old=int_old+inc*(int_new-int_old);   
        
       
        
end
             

%graphical stuff
YMX=0.9*max([err_c err_s err_i]/err_ref);
YMN=0.9*min([err_c err_s err_i]/err_ref);
figure(71); hold on;
colormap('jet');
subplot(1,3,1); hold on; plot([min(shift) max(shift)], [1 1], 'k--', 'LineWidth', LW);
hold on; plot(SH_I(nm,IT)*[1 1], [YMN YMX], 'k--', 'LineWidth', LW);
xlabel('shift (^o latitude)', 'FontSize', FSz, 'FontName', FN);
ylabel('Root mean squared precipitation error over the tropics (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('2D Shift', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
axis([min(shift) max(shift) YMN YMX]);

subplot(1,3,2); hold on; plot([min(con) max(con)], [1 1], 'k--', 'LineWidth', LW);
hold on; plot(CON_I(nm,IT)*[1 1], [YMN YMX], 'k--', 'LineWidth', LW);
xlabel('Contraction (scalar multiplier)', 'FontSize', FSz, 'FontName', FN);
ylabel('Root mean squared precipitation error over the tropics (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('2D Contraction', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
axis([min(con) max(con) YMN YMX]);

subplot(1,3,3); hold on; plot([min(int) max(int)], [1 1], 'k--', 'LineWidth', LW);
hold on; plot(INT_I(nm,IT)*[1 1], [YMN YMX], 'k--', 'LineWidth', LW);
xlabel('Intensification (scalar multiplier)', 'FontSize', FSz, 'FontName', FN);
ylabel('Root mean squared precipitation error over the tropics (mm/day)', 'FontSize', FSz, 'FontName', FN);
title('2D Intensification', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN,'LineWidth', LWA); box on ; grid on; 
axis([min(int) max(int) YMN YMX]);
H=colorbar('peer', gca);
set(get(H,'YLabel'),'string', 'Iteration Number', 'FontSize', FSz, 'FontName', FN);
colorbar('YTick', cbar_ytic, 'YTickLabel', 1:IT, 'LineWidth', LWA);



%percent of change explained by all three, and excluding on mechanism
DP=iPN-iPC;
ALL_err = int_old*interp2(ILON_ext, con_old*(ILAT_ext+sh_old), iPC_ext , ILON, ILAT) - iPN; 
SHIFT_err = int_old*interp2(ILON_ext, con_old*(ILAT_ext), iPC_ext , ILON, ILAT) - iPN;
CON_err = int_old*interp2(ILON_ext, (ILAT_ext+sh_old), iPC_ext , ILON, ILAT) - iPN;
INT_err = interp2(ILON_ext, con_old*(ILAT_ext+sh_old), iPC_ext , ILON, ILAT) - iPN;

TPER(nm,1:4) = [1- ( (mean(squeeze(mean((ALL_err).^2))))^.5) /err_ref ;
               1- ( mean(squeeze(mean((SHIFT_err).^2) ))^.5) / err_ref;
               1- ( (mean(squeeze(mean((CON_err).^2))))^.5) / err_ref;
               1- ( mean(squeeze(mean((INT_err).^2)))^.5) / err_ref];
           
TPER_ASCI(nm, 1:4) = [TPER(nm,1);
                      TPER(nm,1)-TPER(nm,2);
                      TPER(nm,1)-TPER(nm,3);
                      TPER(nm,1)-TPER(nm,4)];
                  

                  
 %ALTERNATIVELY -- try just perturb the single mode off the mean state
 SHIFT_err2 = interp2(ILON_ext, (ILAT_ext+sh_old), iPC_ext , ILON, ILAT) - iPN;
CON_err2 = interp2(ILON_ext, con_old*(ILAT_ext), iPC_ext , ILON, ILAT) - iPN;
INT_err2 = int_old*interp2(ILON_ext, (ILAT_ext), iPC_ext , ILON, ILAT) - iPN;

TPER_ASCI2(nm,1:4) = [1- ( (mean(squeeze(mean((ALL_err).^2))))^.5) /err_ref ;
               1- ( mean(squeeze(mean((SHIFT_err2).^2) ))^.5) / err_ref;
               1- ( (mean(squeeze(mean((CON_err2).^2))))^.5) / err_ref;
               1- ( mean(squeeze(mean((INT_err2).^2)))^.5) / err_ref];
           
 %LASTLY -- the answer from the first iteration
 SHIFT_err3 = interp2(ILON_ext, (ILAT_ext+SH_I_T(nm,1)), iPC_ext , ILON, ILAT) - iPN;
CON_err3 = interp2(ILON_ext, CON_I_T(nm,1)*(ILAT_ext), iPC_ext , ILON, ILAT) - iPN;
INT_err3 = INT_I_T(nm,1)*interp2(ILON_ext, (ILAT_ext), iPC_ext , ILON, ILAT) - iPN;

TPER_ASCI3(nm,1:4) = [1- ( (mean(squeeze(mean((ALL_err).^2))))^.5) /err_ref ;
               1- ( mean(squeeze(mean((SHIFT_err3).^2) ))^.5) / err_ref;
               1- ( (mean(squeeze(mean((CON_err3).^2))))^.5) / err_ref;
               1- ( mean(squeeze(mean((INT_err3).^2)))^.5) / err_ref];
                  
DP = iPN - iPC;
DP_ALL = int_old*interp2(ILON_ext, con_old*(ILAT_ext+sh_old), iPC_ext , ILON, ILAT) - iPC;
DP_SHIFT = interp2(ILON_ext, (ILAT_ext+sh_old), iPC_ext , ILON, ILAT) - iPC;
DP_CON = interp2(ILON_ext, con_old*(ILAT_ext), iPC_ext , ILON, ILAT) - iPC;
DP_INT = int_old*interp2(ILON_ext, (ILAT_ext), iPC_ext , ILON, ILAT) - iPC;

DEL_P_EXT(nm, 1:length(ilat_ext), 1:N) = iPN_ext-iPC_ext;

DELP_P_ASCI(nm, 1:length(ilat_ext), 1:N) = int_old*interp2(ILON_ext2, con_old*(ILAT_ext2+sh_old), iPC_ext2 , ILON_ext, ILAT_ext) - iPC_ext;

%and for the individual modess
DEL_P_S(nm, 1:length(ilat_ext), 1:N) = interp2(ILON_ext2, (ILAT_ext2+sh_old), iPC_ext2 , ILON_ext, ILAT_ext) - iPC_ext;
DEL_P_C(nm, 1:length(ilat_ext), 1:N) = interp2(ILON_ext2, con_old*(ILAT_ext2), iPC_ext2 , ILON_ext, ILAT_ext) - iPC_ext;
DEL_P_I(nm, 1:length(ilat_ext), 1:N) = int_old*interp2(ILON_ext2, (ILAT_ext2), iPC_ext2 , ILON_ext, ILAT_ext) - iPC_ext;

MXLAT=35;
coast=load('coast');
load coast;
colortrans;
long(find(long<=0)) = long(find(long<=0))+360;
long(find(long <0.5 | long>359)) = NaN;
lat(find(isnan(long) == 1))  = NaN;
figure(61); clf; pcolor(ILON, ILAT, iPN-iPC); colormap(cmap);colorbar; shading interp;
hold on; plot(long,lat,'Color', gray, 'LineWidth', 2);
cfrac=0.8;
CMAX=max(squeeze(max(abs(iPN-iPC))));
caxis(cfrac*[-CMAX CMAX]);

axis([0 360 -MXLAT MXLAT]); box on; 
conts=[4:4:30];
AMP=6; dnum=4; %fraction for the optimal change plot and number of contours to display within that fraction
dcont=[-CMAX:CMAX/dnum:-CMAX/dnum   CMAX/dnum:CMAX/dnum:CMAX]; %will be scaled up by 1/dfrac to put it on the same colorbar
figure(61); hold on; contour(ILON, ILAT,iPC, conts,'k-', 'LineWidth', LW);
figure(61); hold on; contour(ILON, ILAT,AMP*DP_ALL, dcont,'-', 'LineWidth', LW);colorbar;
title(['Total % = ' num2str(round(100*TPER_ASCI(nm,1)))], 'FontSize', FSz, 'FontName', FN);

figure(62); clf; pcolor(ILON, ILAT, iPN-iPC); colormap(cmap);colorbar; shading interp;
caxis(cfrac*[-CMAX CMAX]);
hold on; plot(long,lat,'Color', gray, 'LineWidth', 2);
figure(62); hold on; contour(ILON, ILAT,AMP*DP_SHIFT, dcont,'-', 'LineWidth', LW);colorbar;
figure(62); hold on; contour(ILON, ILAT,iPC, conts,'k-', 'LineWidth', LW);
title(['SHIFT % = ' num2str(round(100*TPER_ASCI(nm,2)))], 'FontSize', FSz, 'FontName', FN);

figure(63); clf; pcolor(ILON, ILAT, iPN-iPC); colormap(cmap);colorbar; shading interp;
caxis(cfrac*[-CMAX CMAX]);
hold on; plot(long,lat,'Color', gray, 'LineWidth', 2);
figure(63); hold on; contour(ILON, ILAT,AMP*DP_CON, dcont,'-', 'LineWidth', LW);colorbar;
figure(63); hold on; contour(ILON, ILAT,iPC, conts,'k-', 'LineWidth', LW);
title(['CONTRACTION % = ' num2str(round(100*TPER_ASCI(nm,3)))], 'FontSize', FSz, 'FontName', FN);

figure(64); clf; pcolor(ILON, ILAT, iPN-iPC); colormap(cmap);colorbar; shading interp;
caxis(cfrac*[-CMAX CMAX]);
hold on; plot(long,lat,'Color', gray, 'LineWidth', 2);
figure(64); hold on; contour(ILON, ILAT,AMP*DP_INT, dcont,'-', 'LineWidth', LW);colorbar;
figure(64); hold on; contour(ILON, ILAT,iPC, conts,'k-', 'LineWidth', LW);
title(['INTENSIFICATION % = ' num2str(round(100*TPER_ASCI(nm,4)))], 'FontSize', FSz, 'FontName', FN);

twod_percentages = [nm NM round(100*TPER_ASCI(nm,1:4));
                    nm NM round(100*TPER_ASCI2(nm,1:4));
                    nm NM round(100*TPER_ASCI3(nm,1:4))] 
toc
end


