%work out the change in p_cent and AHT_EQ in the LGM hosing experiments
clear all ; close all;
dir='C:\Users\thedhoe\Desktop\hosing\cc\'
fil='abruptc_2member_19yr.nc'; %the control climnatology
FIL=[dir fil];
ncdisp(FIL);
lat=ncread(FIL, 'lat', [1], [Inf]); NN=length(lat);
lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
lev=ncread(FIL, 'lev', [1], [Inf]); NL=length(lev);
time=ncread(FIL, 'time', [1], [Inf]); NT=length(time);

daysperyear=365;
year=rounddown(time/daysperyear);
month=1+round((time-year*daysperyear)/(365/12)); %seesm to start in December



fields={'PRECC'; %1 convenctive precipitation in m/s
       'PRECL'; %2 large scale precipitation
       'LHFLX'; %3 latent heat flux
       'SHFLX'; %4 sensible heat flux
       'FLNT'; %5 OLR
       'FSNT'; %6 ASR-- net SW at TOA
       'FSNS'; %7 net shortwave at surface
       'FLNS'; %8 shortwave at surface
       'TS'; %9 surface temperaryre
       'LANDFRAC';}; %10 Landfraction
   for n=1:length(fields);
       for mn=1:12;
       do=find(month == mn); %indexes for each month
       twod(n,mn,1:NN,1:N) = double(mean(squeeze(permute( ncread(FIL, char(fields(n,:)), [1 1 do(1)], [Inf Inf Inf], [1 1 12] ), [3 2 1])))); %monthly mean of the given two d variable  
       end
   end
   
   P=squeeze(twod(1,:,:,:) + twod(2,:,:,:));
   LH = squeeze(twod(3,:,:,:));
   SH = squeeze(twod(4,:,:,:));
   OLR = squeeze(twod(5,:,:,:));
   ASR = squeeze(twod(6,:,:,:));
   NS = squeeze(twod(7,:,:,:));
   LNS = squeeze(twod(8,:,:,:));
   TS = squeeze(twod(9,:,:,:));
   LC_p=squeeze(mean(P));
   Lthresh=0.3; %threshold below which to call something ocean domain
   SHF= LH + SH + LNS - NS; %surface heat flux to the atmosphere 
   %check out the global mean of quantities that should be zero
   ZSHF=squeeze(mean(shiftdim(SHF,2))); ASHF=squeeze(mean(ZSHF));
   ZASR=squeeze(mean(shiftdim(ASR,2))); AASR=squeeze(mean(ZASR));
   ZOLR=squeeze(mean(shiftdim(OLR,2))); AOLR=squeeze(mean(ZOLR));
   ZP=squeeze(mean(shiftdim(P,2))); AP=squeeze(mean(ZP)); LC_P=AP;
   ZTS=squeeze(mean(shiftdim(TS,2))); LC_TS=squeeze(mean(ZTS));
   LMASK=squeeze(mean(squeeze(twod(10,:,:,:))));
   LMASK(find(LMASK >=Lthresh))=NaN;
   LMASK(find(isnan(LMASK) == 0)) =1; %NaN and ones ocean mask
   SST=squeeze(mean(TS)).*LMASK;
   LC_SST = squeeze(nanmean(SST'));
   
   x=sind(lat);
   norm=trapz(x, ones(size(x))); %normalization for global integrals, should be 2
   
   GSHF=trapz(x,ASHF)/norm
   GNET=trapz(x,AASR-AOLR)/norm
   
   %now precipitation centroid
   flat=-20:0.05:20;
   for mn=1:12;
       fl=interp1(lat, squeeze(ZP(mn,:)), flat);
       cdf=cumsum(fl)/sum(fl);
       pcent(mn)=flat(find2(cdf,0.5));
   end
   
    fl=interp1(lat, AP, flat);
       cdf=cumsum(fl)/sum(fl);
       PCENT=flat(find2(cdf,0.5));
   
       %work on the heat transport
       mlat=[-90; (lat(1:end-1)+lat(2:end))/2;90];
       dlat=diff(mlat); %differential latitude
       ae=6.36*10^6;
       dA=(ae*dlat*pi/180).*(cosd(lat)*2*pi*ae); %differewntial area
       
       %calculate the heat transport at mlat implied by the TOA and SHF
       %imbalances
       for mn=1:12;
           mht(mn,1)=0;
           oht(mn,1)=0;
           for nn=2:NN+1;
               gnet=trapz(x, squeeze(ZASR(mn,:) - ZOLR(mn,:)))/norm;
               gshf=trapz(x, squeeze(ZSHF(mn,:)))/norm;
               mht(mn,nn)=mht(mn,nn-1) +( ZASR(mn,nn-1) - ZOLR(mn,nn-1) - gnet ).*dA(nn-1);
               oht(mn,nn)=oht(mn,nn-1) -( ZSHF(mn,nn-1) - gshf ).*dA(nn-1);
           end
       end
               
       aht=mht-oht;
       MHT=squeeze(mean(mht)); LC_MHT=MHT;
       OHT=squeeze(mean(oht)); LC_OHT=OHT;
       AHT=squeeze(mean(aht)); LC_AHT=AHT;
       
       
       LW=2; PW=10^-15;
       FSz=26; FN='Arial';
       xmn=-90; xmx=90;
       ymn=-7; ymx=-ymn;
       figure(11); clf;
       hold on; plot(mlat, PW*MHT, 'k-', 'LineWidth', LW);
       hold on; plot(mlat, PW*AHT, 'r-', 'LineWidth', LW);
       hold on; plot(mlat, PW*OHT, 'b-', 'LineWidth', LW);
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Meridional Energy Transport (PW)', 'FontSize', FSz, 'FontName', FN);
       axis([xmn xmx ymn ymx]);
       hold on; plot([xmn xmx], 0*[ymn ymx],'k--', 'LineWidth', LW);
       hold on; plot(0*[xmn xmx], [ymn ymx],'k--', 'LineWidth', LW);
       box on;
       
    ieq=find(mlat== 0);
    AHT_EQ=PW*AHT(ieq)
    PCENT; LC_PCENT=PCENT;
    
    %store what you want for comparison to the hosing results
    C_PCENT=PCENT; CAHT=AHT; COHT=OHT; CMHT=MHT; CP=AP; C_AHT_EQ=AHT_EQ;
    
    LC_P=AP;
    LC_AHT_EQ=AHT_EQ;
       
    
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   dir='C:\Users\thedhoe\Desktop\hosing\cc\'
fil='abrupt_2member_19yr.nc'; %the control climnatology
FIL=[dir fil];
ncdisp(FIL);
lat=ncread(FIL, 'lat', [1], [Inf]); NN=length(lat);
lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
lev=ncread(FIL, 'lev', [1], [Inf]); NL=length(lev);
time=ncread(FIL, 'time', [1], [Inf]); NT=length(time);

daysperyear=365;
year=rounddown(time/daysperyear);
month=1+round((time-year*daysperyear)/(365/12)); %seesm to start in December



fields={'PRECC'; %1 convenctive precipitation in m/s
       'PRECL'; %2 large scale precipitation
       'LHFLX'; %3 latent heat flux
       'SHFLX'; %4 sensible heat flux
       'FLNT'; %5 OLR
       'FSNT'; %6 ASR-- net SW at TOA
       'FSNS'; %7 net shortwave at surface
       'FLNS'; %8 shortwave at surface
       'TS'; %9 surface temperaryre
       'LANDFRAC';} %10 Landfraction
   for n=1:length(fields);
       for mn=1:12;
       do=find(month == mn); %indexes for each month
       twod(n,mn,1:NN,1:N) = double(mean(squeeze(permute( ncread(FIL, char(fields(n,:)), [1 1 do(1)], [Inf Inf Inf], [1 1 12] ), [3 2 1])))); %monthly mean of the given two d variable  
       end
   end
   
   P=squeeze(twod(1,:,:,:) + twod(2,:,:,:));
   LH = squeeze(twod(3,:,:,:));
   SH = squeeze(twod(4,:,:,:));
   OLR = squeeze(twod(5,:,:,:));
   ASR = squeeze(twod(6,:,:,:));
   NS = squeeze(twod(7,:,:,:));
   LNS = squeeze(twod(8,:,:,:));
   TS = squeeze(twod(9,:,:,:));
   LH_p=squeeze(mean(P));
   SHF= LH + SH + LNS - NS; %surface heat flux to the atmosphere 
   %check out the global mean of quantities that should be zero
   ZSHF=squeeze(mean(shiftdim(SHF,2))); ASHF=squeeze(mean(ZSHF));
   ZASR=squeeze(mean(shiftdim(ASR,2))); AASR=squeeze(mean(ZASR));
   ZOLR=squeeze(mean(shiftdim(OLR,2))); AOLR=squeeze(mean(ZOLR));
   ZP=squeeze(mean(shiftdim(P,2))); AP=squeeze(mean(ZP)); LH_P=AP;
   ZTS=squeeze(mean(shiftdim(TS,2))); LH_TS=squeeze(mean(ZTS));
   LMASK=squeeze(mean(squeeze(twod(10,:,:,:))));
   LMASK(find(LMASK >=Lthresh))=NaN;
   LMASK(find(isnan(LMASK) == 0)) =1; %NaN and ones ocean mask
   SST=squeeze(mean(TS)).*LMASK;
   LH_SST = squeeze(nanmean(SST'));
   x=sind(lat);
   norm=trapz(x, ones(size(x))); %normalization for global integrals, should be 2
   
   GSHF=trapz(x,ASHF)/norm
   GNET=trapz(x,AASR-AOLR)/norm
   
   %now precipitation centroid
   flat=-20:0.05:20;
   for mn=1:12;
       fl=interp1(lat, squeeze(ZP(mn,:)), flat);
       cdf=cumsum(fl)/sum(fl);
       pcent(mn)=flat(find2(cdf,0.5));
   end
   
    fl=interp1(lat, AP, flat);
       cdf=cumsum(fl)/sum(fl);
       PCENT=flat(find2(cdf,0.5));
   
       %work on the heat transport
       mlat=[-90; (lat(1:end-1)+lat(2:end))/2;90];
       dlat=diff(mlat); %differential latitude
       ae=6.36*10^6;
       dA=(ae*dlat*pi/180).*(cosd(lat)*2*pi*ae); %differewntial area
       
       %calculate the heat transport at mlat implied by the TOA and SHF
       %imbalances
       for mn=1:12;
           mht(mn,1)=0;
           oht(mn,1)=0;
           for nn=2:NN+1;
               gnet=trapz(x, squeeze(ZASR(mn,:) - ZOLR(mn,:)))/norm;
               gshf=trapz(x, squeeze(ZSHF(mn,:)))/norm;
               mht(mn,nn)=mht(mn,nn-1) +( ZASR(mn,nn-1) - ZOLR(mn,nn-1) - gnet ).*dA(nn-1);
               oht(mn,nn)=oht(mn,nn-1) -( ZSHF(mn,nn-1) - gshf ).*dA(nn-1);
           end
       end
               
       aht=mht-oht;
       MHT=squeeze(mean(mht)); LH_MHT=MHT;
       OHT=squeeze(mean(oht)); LH_OHT=OHT;
       AHT=squeeze(mean(aht)); LH_AHT=AHT;
       
       
       LW=2; PW=10^-15;
       FSz=26; FN='Arial';
       xmn=-90; xmx=90;
       ymn=-7; ymx=-ymn;
       figure(11); hold on;
       hold on; plot(mlat, PW*MHT, 'k--', 'LineWidth', LW);
       hold on; plot(mlat, PW*AHT, 'r--', 'LineWidth', LW);
       hold on; plot(mlat, PW*OHT, 'b--', 'LineWidth', LW);
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Meridional Energy Transport (PW)', 'FontSize', FSz, 'FontName', FN);
       axis([xmn xmx ymn ymx]);
       hold on; plot([xmn xmx], 0*[ymn ymx],'k--', 'LineWidth', LW);
       hold on; plot(0*[xmn xmx], [ymn ymx],'k--', 'LineWidth', LW);
       box on;
       
    ieq=find(mlat== 0);
    AHT_EQ=PW*AHT(ieq)
    PCENT;
    LH_P=AP;
    LH_PCENT=PCENT;
    LH_AHT_EQ=AHT_EQ;   
    break
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure(12); clf;
    LWP=1.5;
    SF=1000*3600*24; %from m/s to mm/day
    hold on; plot(lat, SF*CP, 'k-', 'LineWidth', LW);
    hold on; plot(lat, SF*AP, 'b-', 'LineWidth', LW);
  
    ymx=7; ymn=0;
    xmn=-30; xmx=-xmn;
    axis([xmn xmx ymn ymx]); box on;
    hold on; plot(C_PCENT*[1 1], [ymn ymx], 'k--', 'LineWidth',LWP)
    hold on; plot(PCENT*[1 1], [ymn ymx], 'b--', 'LineWidth',LWP)
    hold on; plot(0*[1 1], [ymn ymx], 'm--', 'LineWidth',LWP)
    legend('Modern Control', 'Modern HOSED', 'Control P_C_E_N_T', 'HOSED P_C_E_N_T', 'EQUATOR');
    
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Precipitation (mm/day)', 'FontSize', FSz, 'FontName', FN);
    
    
   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODERN
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DATY
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MEAN
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STATE
    %work out the change in p_cent and AHT_EQ in the LGM hosing experiments

dir='C:\Users\thedhoe\Desktop\hosing\cc\'
fil='fwc_3n6member_19yr.nc'; %the MODERN control climnatology
FIL=[dir fil];
ncdisp(FIL);
lat=ncread(FIL, 'lat', [1], [Inf]); NN=length(lat);
lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
lev=ncread(FIL, 'lev', [1], [Inf]); NL=length(lev);
time=ncread(FIL, 'time', [1], [Inf]); NT=length(time);

daysperyear=365;
year=rounddown(time/daysperyear);
month=1+round((time-year*daysperyear)/(365/12)); %seesm to start in December



fields={'PRECC'; %1 convenctive precipitation in m/s
       'PRECL'; %2 large scale precipitation
       'LHFLX'; %3 latent heat flux
       'SHFLX'; %4 sensible heat flux
       'FLNT'; %5 OLR
       'FSNT'; %6 ASR-- net SW at TOA
       'FSNS'; %7 net shortwave at surface
       'FLNS'; %8 shortwave at surface
       'TS'; %9 surface temperaryre
       'LANDFRAC';}; %10 Landfraction
   
   for n=1:length(fields);
       for mn=1:12;
       do=find(month == mn); %indexes for each month
       twod(n,mn,1:NN,1:N) = double(mean(squeeze(permute( ncread(FIL, char(fields(n,:)), [1 1 do(1)], [Inf Inf Inf], [1 1 12] ), [3 2 1])))); %monthly mean of the given two d variable  
       end
   end
   
   P=squeeze(twod(1,:,:,:) + twod(2,:,:,:));
   LH = squeeze(twod(3,:,:,:));
   SH = squeeze(twod(4,:,:,:));
   OLR = squeeze(twod(5,:,:,:));
   ASR = squeeze(twod(6,:,:,:));
   NS = squeeze(twod(7,:,:,:));
   LNS = squeeze(twod(8,:,:,:));
   TS = squeeze(twod(9,:,:,:));
   
   SHF= LH + SH + LNS - NS; %surface heat flux to the atmosphere 
   %check out the global mean of quantities that should be zero
   ZSHF=squeeze(mean(shiftdim(SHF,2))); ASHF=squeeze(mean(ZSHF));
   ZASR=squeeze(mean(shiftdim(ASR,2))); AASR=squeeze(mean(ZASR));
   ZOLR=squeeze(mean(shiftdim(OLR,2))); AOLR=squeeze(mean(ZOLR));
   ZP=squeeze(mean(shiftdim(P,2))); AP=squeeze(mean(ZP)); MC_P=AP;
   ZTS=squeeze(mean(shiftdim(TS,2))); MC_TS=squeeze(mean(ZTS));
   LMASK=squeeze(mean(squeeze(twod(10,:,:,:))));
   LMASK(find(LMASK >=Lthresh))=NaN;
   LMASK(find(isnan(LMASK) == 0)) =1; %NaN and ones ocean mask
   SST=squeeze(mean(TS)).*LMASK;
   MC_SST = squeeze(nanmean(SST'));
   x=sind(lat);
   norm=trapz(x, ones(size(x))); %normalization for global integrals, should be 2
   
   GSHF=trapz(x,ASHF)/norm
   GNET=trapz(x,AASR-AOLR)/norm
   
   %now precipitation centroid
   flat=-20:0.05:20;
   for mn=1:12;
       fl=interp1(lat, squeeze(ZP(mn,:)), flat);
       cdf=cumsum(fl)/sum(fl);
       pcent(mn)=flat(find2(cdf,0.5));
   end
   
    fl=interp1(lat, AP, flat);
       cdf=cumsum(fl)/sum(fl);
       PCENT=flat(find2(cdf,0.5));
   
       %work on the heat transport
       mlat=[-90; (lat(1:end-1)+lat(2:end))/2;90];
       dlat=diff(mlat); %differential latitude
       ae=6.36*10^6;
       dA=(ae*dlat*pi/180).*(cosd(lat)*2*pi*ae); %differewntial area
       
       %calculate the heat transport at mlat implied by the TOA and SHF
       %imbalances
       for mn=1:12;
           mht(mn,1)=0;
           oht(mn,1)=0;
           for nn=2:NN+1;
               gnet=trapz(x, squeeze(ZASR(mn,:) - ZOLR(mn,:)))/norm;
               gshf=trapz(x, squeeze(ZSHF(mn,:)))/norm;
               mht(mn,nn)=mht(mn,nn-1) +( ZASR(mn,nn-1) - ZOLR(mn,nn-1) - gnet ).*dA(nn-1);
               oht(mn,nn)=oht(mn,nn-1) -( ZSHF(mn,nn-1) - gshf ).*dA(nn-1);
           end
       end
               
       aht=mht-oht;
       MHT=squeeze(mean(mht)); MC_MHT=MHT;
       OHT=squeeze(mean(oht)); MC_OHT=OHT;
       AHT=squeeze(mean(aht)); MC_AHT=AHT;
       
       
       LW=2; PW=10^-15;
       FSz=26; FN='Arial';
       xmn=-90; xmx=90;
       ymn=-7; ymx=-ymn;
       figure(1); clf;
       hold on; plot(mlat, PW*MHT, 'k-', 'LineWidth', LW);
       hold on; plot(mlat, PW*AHT, 'r-', 'LineWidth', LW);
       hold on; plot(mlat, PW*OHT, 'b-', 'LineWidth', LW);
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Meridional Energy Transport (PW)', 'FontSize', FSz, 'FontName', FN);
       axis([xmn xmx ymn ymx]);
       hold on; plot([xmn xmx], 0*[ymn ymx],'k--', 'LineWidth', LW);
       hold on; plot(0*[xmn xmx], [ymn ymx],'k--', 'LineWidth', LW);
       box on;
       
    ieq=find(mlat== 0);
    AHT_EQ=PW*AHT(ieq)
    C_OHT_EQ=PW*OHT(ieq);
    C_MHT_EQ=PW*MHT(ieq);
    
    MC_PCENT=PCENT;
    MC_AHT_EQ=AHT_EQ;
    %store what you want for comparison to the hosing results
    C_PCENT=PCENT; CAHT=AHT; COHT=OHT; CMHT=MHT; CP=AP; C_AHT_EQ=AHT_EQ;
    
    MC_P=AP;
       
    
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   dir='C:\Users\thedhoe\Desktop\hosing\cc\'
fil='fw_3n6member_19yr.nc'; %the control climnatology
FIL=[dir fil];
ncdisp(FIL);
lat=ncread(FIL, 'lat', [1], [Inf]); NN=length(lat);
lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
lev=ncread(FIL, 'lev', [1], [Inf]); NL=length(lev);
time=ncread(FIL, 'time', [1], [Inf]); NT=length(time);

daysperyear=365;
year=rounddown(time/daysperyear);
month=1+round((time-year*daysperyear)/(365/12)); %seesm to start in December



fields={'PRECC'; %1 convenctive precipitation in m/s
       'PRECL'; %2 large scale precipitation
       'LHFLX'; %3 latent heat flux
       'SHFLX'; %4 sensible heat flux
       'FLNT'; %5 OLR
       'FSNT'; %6 ASR-- net SW at TOA
       'FSNS'; %7 net shortwave at surface
       'FLNS'; %8 shortwave at surface
       'TS'; %9 surface temperaryre
       'LANDFRAC';}; %10 Landfraction
   
   for n=1:length(fields);
       for mn=1:12;
       do=find(month == mn); %indexes for each month
       twod(n,mn,1:NN,1:N) = double(mean(squeeze(permute( ncread(FIL, char(fields(n,:)), [1 1 do(1)], [Inf Inf Inf], [1 1 12] ), [3 2 1])))); %monthly mean of the given two d variable  
       end
   end
   
   P=squeeze(twod(1,:,:,:) + twod(2,:,:,:));
   LH = squeeze(twod(3,:,:,:));
   SH = squeeze(twod(4,:,:,:));
   OLR = squeeze(twod(5,:,:,:));
   ASR = squeeze(twod(6,:,:,:));
   NS = squeeze(twod(7,:,:,:));
   LNS = squeeze(twod(8,:,:,:));
   TS = squeeze(twod(9,:,:,:));
   
   SHF= LH + SH + LNS - NS; %surface heat flux to the atmosphere 
   %check out the global mean of quantities that should be zero
   ZSHF=squeeze(mean(shiftdim(SHF,2))); ASHF=squeeze(mean(ZSHF));
   ZASR=squeeze(mean(shiftdim(ASR,2))); AASR=squeeze(mean(ZASR));
   ZOLR=squeeze(mean(shiftdim(OLR,2))); AOLR=squeeze(mean(ZOLR));
   ZP=squeeze(mean(shiftdim(P,2))); AP=squeeze(mean(ZP)); MH_P=AP;
   ZTS=squeeze(mean(shiftdim(TS,2))); MH_TS=squeeze(mean(ZTS));
      LMASK=squeeze(mean(squeeze(twod(10,:,:,:))));
   LMASK(find(LMASK >=Lthresh))=NaN;
   LMASK(find(isnan(LMASK) == 0)) =1; %NaN and ones ocean mask
   SST=squeeze(mean(TS)).*LMASK;
   MH_SST = squeeze(nanmean(SST'));
   
   x=sind(lat);
   norm=trapz(x, ones(size(x))); %normalization for global integrals, should be 2
   
   GSHF=trapz(x,ASHF)/norm
   GNET=trapz(x,AASR-AOLR)/norm
   
   %now precipitation centroid
   flat=-20:0.05:20;
   for mn=1:12;
       fl=interp1(lat, squeeze(ZP(mn,:)), flat);
       cdf=cumsum(fl)/sum(fl);
       pcent(mn)=flat(find2(cdf,0.5));
   end
   
    fl=interp1(lat, AP, flat);
       cdf=cumsum(fl)/sum(fl);
       PCENT=flat(find2(cdf,0.5));
   
       %work on the heat transport
       mlat=[-90; (lat(1:end-1)+lat(2:end))/2;90];
       dlat=diff(mlat); %differential latitude
       ae=6.36*10^6;
       dA=(ae*dlat*pi/180).*(cosd(lat)*2*pi*ae); %differewntial area
       
       %calculate the heat transport at mlat implied by the TOA and SHF
       %imbalances
       for mn=1:12;
           mht(mn,1)=0;
           oht(mn,1)=0;
           for nn=2:NN+1;
               gnet=trapz(x, squeeze(ZASR(mn,:) - ZOLR(mn,:)))/norm;
               gshf=trapz(x, squeeze(ZSHF(mn,:)))/norm;
               mht(mn,nn)=mht(mn,nn-1) +( ZASR(mn,nn-1) - ZOLR(mn,nn-1) - gnet ).*dA(nn-1);
               oht(mn,nn)=oht(mn,nn-1) -( ZSHF(mn,nn-1) - gshf ).*dA(nn-1);
           end
       end
               
       aht=mht-oht;
       MHT=squeeze(mean(mht)); MH_MHT=MHT;
       OHT=squeeze(mean(oht)); MH_OHT=OHT;
       AHT=squeeze(mean(aht)); MH_AHT=AHT;
       
       
       LW=2; PW=10^-15;
       FSz=26; FN='Arial';
       xmn=-90; xmx=90;
       ymn=-7; ymx=-ymn;
       figure(1); hold on;
       hold on; plot(mlat, PW*MHT, 'k--', 'LineWidth', LW);
       hold on; plot(mlat, PW*AHT, 'r--', 'LineWidth', LW);
       hold on; plot(mlat, PW*OHT, 'b--', 'LineWidth', LW);
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Meridional Energy Transport (PW)', 'FontSize', FSz, 'FontName', FN);
       axis([xmn xmx ymn ymx]);
       hold on; plot([xmn xmx], 0*[ymn ymx],'k--', 'LineWidth', LW);
       hold on; plot(0*[xmn xmx], [ymn ymx],'k--', 'LineWidth', LW);
       box on;
       
    ieq=find(mlat== 0);
    AHT_EQ=PW*AHT(ieq)
    MH_PCENT=PCENT;
    MH_AHT_EQ=AHT_EQ;
    OHT_EQ=PW*OHT(ieq);
    MHT_EQ=PW*MHT(ieq);
    MH_P=AP;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure(2); clf;
    LWP=1.5;
    SF=1000*3600*24; %from m/s to mm/day
    hold on; plot(lat, SF*CP, 'k-', 'LineWidth', LW);
    hold on; plot(lat, SF*AP, 'b-', 'LineWidth', LW);
  
    ymx=7; ymn=0;
    xmn=-30; xmx=-xmn;
    axis([xmn xmx ymn ymx]); box on;
    hold on; plot(C_PCENT*[1 1], [ymn ymx], 'k--', 'LineWidth',LWP)
    hold on; plot(PCENT*[1 1], [ymn ymx], 'b--', 'LineWidth',LWP)
    hold on; plot(0*[1 1], [ymn ymx], 'm--', 'LineWidth',LWP)
    legend('Modern Control', 'Modern HOSED', 'Control P_C_E_N_T', 'HOSED P_C_E_N_T', 'EQUATOR');
    
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Precipitation (mm/day)', 'FontSize', FSz, 'FontName', FN);
       
       %%% plot up the surface temperatures -- make sure these are the runs
       %%% you think they are
       figure(5); clf;
       subplot(1,2,1);
       hold on; plot(lat, LC_TS, 'b-', 'LineWidth', LW);
       hold on; plot(lat, LH_TS, 'b--', 'LineWidth', LW);
       hold on; plot(lat, MC_TS, 'r-', 'LineWidth', LW);
       hold on; plot(lat, MH_TS, 'r--', 'LineWidth', LW);
       legend('LGM control', 'LGM hosed', 'Modern control', 'Modern Hosed');
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2); box on;
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Surface Temperature (K)', 'FontSize', FSz, 'FontName', FN);
       
              subplot(1,2,2);
       hold on; plot(lat, LC_SST, 'b-', 'LineWidth', LW);
       hold on; plot(lat, LH_SST, 'b--', 'LineWidth', LW);
       hold on; plot(lat, MC_SST, 'r-', 'LineWidth', LW);
       hold on; plot(lat, MH_SST, 'r--', 'LineWidth', LW);
       legend('LGM control', 'LGM hosed', 'Modern control', 'Modern Hosed');
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2); box on;
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Surface Temperature (K)', 'FontSize', FSz, 'FontName', FN);
       
       LW=2.5;
       LWP=1.2;
       ymx=7; ymn=0;
       figure(6); clf;
       hold on; plot(lat, SF*LC_P, 'b-', 'LineWidth', LW);
       hold on; plot(lat, SF*LH_P, 'b--', 'LineWidth', LW);
       hold on; plot(lat, SF*MC_P, 'r-', 'LineWidth', LW);
       hold on; plot(lat, SF*MH_P, 'r--', 'LineWidth', LW);
       legend('LGM control', 'LGM hosed', 'Modern control', 'Modern Hosed');
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2); box on;
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Precipitation (mm/day)', 'FontSize', FSz, 'FontName', FN);
       hold on; plot(LC_PCENT*[1 1], [ymn ymx], 'b-', 'LineWidth', LWP);
       hold on; plot(LH_PCENT*[1 1], [ymn ymx], 'b--', 'LineWidth', LWP);
       hold on; plot(MC_PCENT*[1 1], [ymn ymx], 'r-', 'LineWidth', LWP);
       hold on; plot(MH_PCENT*[1 1], [ymn ymx], 'r--', 'LineWidth', LWP);
       hold on; plot(0*[1 1], [ymn ymx], 'm-', 'LineWidth', 2);
       axis([-25 25 ymn ymx])
       
       figure(7); clf;
       hold on; plot(mlat, PW*(LH_MHT-LC_MHT), 'k-', 'LineWidth', LW);
       hold on; plot(mlat, PW*(MH_MHT-MC_MHT), 'k--', 'LineWidth', LW);
       
       hold on; plot(mlat, PW*(LH_OHT-LC_OHT), 'b-', 'LineWidth', LW);
       hold on; plot(mlat, PW*(LH_AHT-LC_AHT), 'r-', 'LineWidth', LW);
       
       hold on; plot(mlat, PW*(MH_OHT-MC_OHT), 'b--', 'LineWidth', LW);
       hold on; plot(mlat, PW*(MH_AHT-MC_AHT), 'r--', 'LineWidth', LW);
       legend('LGM MHT', 'Modern MHT', 'LGM OHT', 'LGM AHT')
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTickLabel', ['']); box on;
       xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
       ylabel('Change in meridional heat transport (Hosed-control -- in PW )', 'FontSize', FSz, 'FontName', FN);
       ymn=-0.4; ymx=-ymn;
       xmn=-90; xmx=90;
       axis([xmn xmx ymn ymx]);
       hold on; plot([xmn xmx], 0*[ymn ymx], 'm-', 'LineWidth', 2);
       hold on; plot(0*[xmn xmx], [ymn ymx], 'm-', 'LineWidth', 2);
       
       %% and a scatter plot of AHT_EQ and P_CENT
       MS=16; MLW=3;
       figure(8);clf;
       hold on; plot(LC_AHT_EQ, LC_PCENT, 'bd', 'MarkerSize', MS, 'LineWidth', MLW, 'MarkerFaceColor', 'b');
       hold on; plot(LH_AHT_EQ, LH_PCENT, 'bd', 'MarkerSize', MS, 'LineWidth', MLW, 'MarkerFaceColor', 'w');
       
       hold on; plot(MC_AHT_EQ, MC_PCENT, 'rd', 'MarkerSize', MS, 'LineWidth', MLW, 'MarkerFaceColor', 'r');
       hold on; plot(MH_AHT_EQ, MH_PCENT, 'rd', 'MarkerSize', MS, 'LineWidth', MLW, 'MarkerFaceColor', 'w');
       xmn=-0.5; xmx=-xmn;
       ymn=-2; ymx=-ymn;
       
       axis([xmn xmx ymn ymx]); box on;
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTickLabel', [''])
       hold on; plot([xmn xmx], 0*[ymn ymx], 'k-', 'LineWidth', 2)
       hold on; plot(0*[xmn xmx], [ymn ymx], 'k-', 'LineWidth', 2)
       xlabel('AHT_E_Q (PW)', 'FontSize', FSz, 'FontName', FN);
       ylabel('P_C_E_N_T (^o)', 'FontSize', FSz, 'FontName', FN);
       legend('LGM control', 'LGM hosed', 'Modern control', 'Modern Hosed')
       
       X = [LC_AHT_EQ LH_AHT_EQ MC_AHT_EQ MH_AHT_EQ];
       Y = [LC_PCENT LH_PCENT MC_PCENT MH_PCENT];
       [a b ci R]=lreg3(X,Y)
       
       %compose the SST gradient across the equator
       cross=30;
       don = find(lat>0 & lat<=cross);
       dos = find(lat<0 & lat >=-cross);
       LC_GSST= mean(LC_SST(don)) - mean(LC_SST(dos));
       LH_GSST= mean(LH_SST(don)) - mean(LH_SST(dos));
       MC_GSST= mean(MC_SST(don)) - mean(MC_SST(dos));
       MH_GSST= mean(MH_SST(don)) - mean(MH_SST(dos));
       
       figure(9);clf;
       hold on; plot(LC_GSST, LC_PCENT, 'bd', 'MarkerSize', MS, 'LineWidth', MLW, 'MarkerFaceColor', 'b');
       hold on; plot(LH_GSST, LH_PCENT, 'bd', 'MarkerSize', MS, 'LineWidth', MLW, 'MarkerFaceColor', 'w');
       
       hold on; plot(MC_GSST, MC_PCENT, 'rd', 'MarkerSize', MS, 'LineWidth', MLW, 'MarkerFaceColor', 'r');
       hold on; plot(MH_GSST, MH_PCENT, 'rd', 'MarkerSize', MS, 'LineWidth', MLW, 'MarkerFaceColor', 'w');
       xmn=-0.5; xmx=-xmn;
       ymn=-2; ymx=-ymn;
       
       axis([xmn xmx ymn ymx]); box on;
       set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);%, 'YTickLabel', [''])
       hold on; plot([xmn xmx], 0*[ymn ymx], 'k-', 'LineWidth', 2)
       hold on; plot(0*[xmn xmx], [ymn ymx], 'k-', 'LineWidth', 2)
       xlabel('Hemispheric difference in Tropical SST  (NH-SH in K)', 'FontSize', FSz, 'FontName', FN);
       ylabel('P_C_E_N_T (^o)', 'FontSize', FSz, 'FontName', FN);
       legend('LGM control', 'LGM hosed', 'Modern control', 'Modern Hosed')
       
       X = [LC_GSST LH_GSST MC_GSST MH_GSST];
       Y = [LC_PCENT LH_PCENT MC_PCENT MH_PCENT];
       [a b ci R]=lreg3(X,Y)
       
        
       