%work out the change in p_cent and AHT_EQ in the LGM hosing experiments
clear all ; close all;
LW=2; FSz=26; FN='Arial'; PW=10^(-15);

dir='C:\Users\thedhoe\Desktop\hosing\cc\'
fil='fwc_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}
   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));
   ZTS=squeeze(mean(shiftdim(TS,2))); ATS=squeeze(mean(ZTS));
   LMASK=squeeze(mean(squeeze(twod(10,:,:,:))));
   LMASK(find(LMASK >=0.3))=NaN;
   LMASK(find(isnan(LMASK) == 0)) =1; %NaN and ones ocean mask
   SST=squeeze(mean(TS)).*LMASK;
   MC_SST = squeeze(mean(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
   
   %and hemispheric averages
   donh=find(lat>=0);
   dosh=find(lat<=0);
   HASR = [trapz(x(donh),AASR(donh))/trapz(x(donh), ones(size(donh)))      trapz(x(dosh),AASR(dosh))/trapz(x(dosh), ones(size(dosh)))]
   HOLR = [trapz(x(donh),AOLR(donh))/trapz(x(donh), ones(size(donh)))      trapz(x(dosh),AOLR(dosh))/trapz(x(dosh), ones(size(dosh)))]
   %convert to the associated heat transport with the hemispheric asymetry
   %of OLR and ASR
   ae=6.36*10^6;
   ASR_brack=PW*(2*pi*ae^2)*(HASR(2)-HASR(1))/2
   OLR_brack=PW*(2*pi*ae^2)*(HOLR(1)-HOLR(2))/2
   
   
   figure(19);clf;
   subplot(1,2,1);
   hold on; plot(lat(donh), ATS(donh), 'r', 'LineWidth', LW);
   hold on; plot(abs(lat(dosh)), ATS(dosh), 'b', 'LineWidth', LW)
   legend('NH', 'SH'); 
   set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);
   xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
   ylabel('Surface Temperature (K)', 'FontSize', FSz, 'FontName', FN);
   title('Surface Temperature', 'FontSize', FSz, 'FontName', FN)
   
   subplot(1,2,2);
   hold on; plot(lat(donh), AOLR(donh), 'r', 'LineWidth', LW);
   hold on; plot(abs(lat(dosh)), AOLR(dosh), 'b', 'LineWidth', LW)
   legend('NH', 'SH'); 
   set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2);
   xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
   ylabel('Outgoinhg Longwave Radiation (W m^-^2)', 'FontSize', FSz, 'FontName', FN);
   title('OLR', 'FontSize', FSz, 'FontName', FN)
   
   
   
   %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));
       OHT=squeeze(mean(oht));
       AHT=squeeze(mean(aht));
       
       
    
       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);
    
    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;
    
       
    
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   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
       
   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));
   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));
       OHT=squeeze(mean(oht));
       AHT=squeeze(mean(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)
    PCENT
    OHT_EQ=PW*OHT(ieq);
    MHT_EQ=PW*MHT(ieq);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    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);
       
   
   