%compare the seasonal heating from surface fluxes and direct solar
%absorption
clear all; close all;
%load the climatologiclay data
dr =    {'C:/research/aqua_gcm/slab_2.4m/';  %directory listing dfor each run
         'C:/research/aqua_gcm/slab_6m/';
         'C:/research/aqua_gcm/slab_12m/';
         'C:/research/aqua_gcm/slab_24m/';
         'C:/research/aqua_gcm/slab_50m/'};
     
NR=5;
colormap('jet');cmap=colormap; %save the RGB colormap vector
CN=length(cmap);
figure(1); clf;
depth={'2.4';'6';'12';'24';'50'}
FSz=22; LW=2; FN='Arial'; %figure plottinf options
LW2=1.5; %for the phasing info
LW=2; LWA=2.7;
g=9.81;%gravity
cp=1007; %specific heat at constant pressurre (J/(kg *K(\)
depths=[2.4 6 12 24 50];
col=[         0         0    0.5625
         0    0.5000    1.0000
    0.4375    1.0000    0.5625
    1.0000    0.6250         0
    0.5         0         0]; % discrete colormap (stolen from jet)

for nr=1:5;  %loop over run (depths)
     
     cn=round(CN*nr/NR); %index of the given run in the colomap matrix (could make it scale as depth if you want
     load([char(dr(nr)) 'energy_climatology.mat']);
     load([char(dr(nr)) 'MSE_TEN.mat']);
     load([char(dr(nr)) 'T.mat']);
     
     %zonal average the temperature
     T=squeeze(mean(shiftdim(T,3)));
     NL=length(Pfull); NN=length(lat);
     %then take the fundamental amplitude and phase
     cross=40; %crossover latitude poleward of which to average over
     doz=find(Pfull>00); %troposphere
     %doz2=find(Pfull > 10 & Pfull <300); %upper troposphere
     doz2=find(Pfull > 800 & Pfull <1000); %upper troposphere
     do=find(lat>cross); %domain of latitudes
     [coef, phase] = fund2_3d(T(:,:,do));
     PHASE=((2*pi-phase)*365/(2*pi)+15)-175; %day of maximum relative to solstice
     
     Phalf=mean([Pfull(1:end-1) Pfull(2:end)]')';
     Pdiff=[Phalf;1000]-[0;Phalf]; %pressure differenc by layer
     W=zeros(NL,1);W(doz)=Pdiff(doz); W=W/sum(W); %normalized weighting vector for mass integrals over the troposphere
     W2=zeros(NL,1);W2(doz2)=Pdiff(doz2); W2=W2/sum(W2); %normalized weighting vector for mass integrals over the troposphere
     
     for mn=1:12;
         field=squeeze(T(mn,:,:));
         TV(mn,1:NN)= W'*field; %mass weighted vertical average of T for each month and latitude 
         TV2(mn,1:NN)= W2'*field; %mass weighted vertical average of T for each month and latitude 
     end
     
     color_trans;
     
     %make the vertical structure of the temperature
     avg=ones(length(doz),1)*mean(coef(doz,:)); %column average seasonal cycle 
     
     
%      MAX=1.3*max(max(coef(doz,:)));
    [LAT,LEV]=meshgrid(lat(do),Pfull(doz)); %mesh for vertical cross sections
%      figure(10+nr)%;subplot(1,2,1);
%      pcolor(LAT, -LEV, coef(doz,:)./avg); colormap(cmap);
%      caxis([.8 1.2]);
%      colorbar;
     
          %color_dhoe3_fine;
%      figure(20+nr)
%      colormap('jet');
%      pcolor(LAT, -LEV, coef(doz,:)); colorbar 
%      ytic=-1000:100:-300;
%      yticl=-ytic;
%      xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
%      ylabel('Pressure level (hPa)', 'FontSize', FSz, 'FontName', FN)
%      shading interp
%      set(gca,'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTick', ytic, 'YTickLabel', yticl); 
%      title([char(depth(nr,:)) ' Meters'], 'FontSize', FSz,'FontName', FN) 
     
     
%      figure(30+nr)
%      %ytic=-1000:100:-300;
%      %yticl=-ytic;
%      colormap('jet');
%      pcolor(LAT, -LEV, PHASE(doz,:)); caxis([45 150]); colorbar;%caxis([45 115]);colorbar 
%      xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
%      ylabel('Pressure level (hPa)', 'FontSize', FSz, 'FontName', FN)
%      shading interp
%      set(gca,'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTick', ytic, 'YTickLabel', yticl); 
%      
     %load([char(dr(n)) 'T_clim.mat'], 'TS','T','lev');
     %nl = find2(lev, 700); %index for the vertical level you want
     
     SHF = squeeze(mean(shiftdim(FS -US + DS,2))); %the surface heat flux upward from the ocean (not counting the solar that passes through the atmosphere

     NN=length(lat); N=length(lon); %dimensionality
     
     %define the absorbed radiation
     SWABS = squeeze(mean(shiftdim(DT - UT -DS + US,2))); %SW radiation absorbed within the atmosphere
     ZDT=squeeze(mean(shiftdim(DT,2)));
     ZOLR=squeeze(mean(shiftdim(OLR,2)));
     ZTETEN3=squeeze(mean(shiftdim(TTEN,2)));
     ZTETEN=MSE_TEN/2.03; %the 2 month difference was divided by thirty days when it sohould have been divided by 365/6 (this corrects that)
     ZDSETEN=DSE_TEN/2.03;
     ZTETEN2=(sum(Pdiff(doz))*100*cp/g) * ( TV([2:12 1],:) - TV([12 1:11],:)) / (24*3600*365/6);
     %ZTEDIV=squeeze(mean(shiftdim(TEDIV,2)));
     ZTEDIV=ZTETEN+ZOLR-SWABS-SHF; %as a residual
     

     %now define the seasonal amplitude
     SF=2^.5;
     for nn=1:NN;
        ind =(sqrt(12/11))*squeeze(ZDT(:,nn)-squeeze(mean(ZDT(:,nn))))/std(squeeze(ZDT(:,nn))); %normalized index to give regression values (accounts for small sample size of std
         S_SWABS(nn) = SF*mean(ind.*squeeze(SWABS(:,nn)));
         S_SHF(nn) = SF*mean(ind.*squeeze(SHF(:,nn)));
         S_TEDIV(nn) = SF*mean(ind.*squeeze(ZTEDIV(:,nn)));
        S_OLR(nn) = -SF*mean(ind.*squeeze(ZOLR(:,nn)));
        S_TETEN(nn) = -SF*mean(ind.*squeeze(ZTETEN(:,nn)));
        %S_TETEN2(nn) = -SF*mean(ind.*squeeze(ZTETEN2(:,nn)));
        %S_TETEN(nn) = -SF*mean(ind.*squeeze(MSE_TEN(:,nn)));
     end
     

     
        
        %average poleward of a certaon latitude
        x=sind(lat); xdo=x(do); %sine of latitude (for area weighting over the polar cap
        norm=trapz(xdo,ones(size(xdo))); %normaliztion for averages over this region
        for mn=1:12;
        RM=1; %remove annual mean? 1 for yes and 0 for no    
        R_SWABS(mn)=trapz(xdo, SWABS(mn,do)-RM*mean(SWABS(:,do)))/norm; %monthly average over the polar cap
        R_SHF(mn)=trapz(xdo, SHF(mn,do)-RM*mean(SHF(:,do)))/norm; %monthly average over the polar cap
        R_TEDIV(mn)=trapz(xdo, ZTEDIV(mn,do)-RM*mean(ZTEDIV(:,do)))/norm; %monthly average over the polar cap
        R_OLR(mn)=trapz(xdo, ZOLR(mn,do)-RM*mean(ZOLR(:,do)))/norm; %monthly average over the polar cap
        R_DT(mn)=trapz(xdo, ZDT(mn,do)-RM*mean(ZDT(:,do)))/norm; %monthly average over the polar cap
        R_TETEN(mn)=trapz(xdo, ZTETEN(mn,do)-RM*mean(ZTETEN(:,do)))/norm; %monthly average over the polar cap
        R_TV(mn)=trapz(xdo, TV(mn,do)-RM*mean(TV(:,do)))/norm;
        R_TV2(mn)=trapz(xdo, TV2(mn,do)-RM*mean(TV2(:,do)))/norm;
        end
        
        %now crank out the phases of the terms
        SWABS_ph(nr)=fund_phase(R_SWABS); %day of the peak
        DT_ph(nr)=fund_phase(R_DT); %day of the peak
        SHF_ph(nr)=fund_phase(R_SHF); %day of the peak
        TETEN_ph(nr)=fund_phase(R_TETEN); %day of the peak
        TV_ph(nr)=fund_phase(R_TV); %day of the peak
        TV2_ph(nr)=fund_phase(R_TV2); %day of the peak
        OLR_ph(nr)=fund_phase(R_OLR); %day of the peak
        TEDIV_ph(nr)=fund_phase(-R_TEDIV); %day of the peak
        HEAT_ph(nr)=fund_phase(R_SHF+R_SWABS);
        

        SWABS_amp(nr)=fund(R_SWABS); %day of the peak
        DT_amp(nr)=fund(R_DT); %day of the peak
        SHF_amp(nr)=fund(R_SHF); %day of the peak
        TETEN_amp(nr)=fund(R_TETEN); %day of the peak
        TV_amp(nr)=fund(R_TV); %day of the peak
        TV2_amp(nr)=fund(R_TV2); %day of the peak
        OLR_amp(nr)=fund(R_OLR); %day of the peak
        TEDIV_amp(nr)=fund(-R_TEDIV); %day of the peak
        HEAT_amp(nr)=fund(R_SHF+R_SWABS);
        
        TS=R_SHF+R_SWABS;
        ind=(2^.5*(12/11)^.5)*(TS-mean(TS))./std(TS); %dot and average to get the amplitude in phase with the insolation
        SWABS_AMP(nr)=mean(ind.*R_SWABS);
        SHF_AMP(nr)=mean(ind.*R_SHF);
        TEDIV_AMP(nr)=mean(ind.*-R_TEDIV);
        OLR_AMP(nr)=mean(ind.*R_OLR);
        TETEN_AMP(nr)= mean(ind.*R_TETEN);

end

%the phasing of the terms in each run
figure(7);clf;
loc=1:5; %location for each run on the phase plot
REF=monthday(6)+21; %reference date
MS=18;
  col=[0 .8 1];
for nr=1:5;
    hold on;plot(HEAT_ph(nr)-REF, loc(nr), 'bd', 'MarkerSize', MS,'LineWidth', 3, 'MarkerFaceColor', 'r');
    hold on;plot(TV_ph(nr)-REF, loc(nr), 'kd', 'MarkerSize', MS, 'MarkerFaceColor', 'k');
    hold on;plot(TETEN_ph(nr)-REF, loc(nr), 'cd', 'MarkerSize', MS, 'MarkerFaceColor', col, 'MarkerEdgeColor', col);
    hold on;plot(OLR_ph(nr)-REF, loc(nr), 'gd', 'MarkerSize', MS, 'MarkerFaceColor', 'g');
    hold on;plot(TEDIV_ph(nr)-REF, loc(nr), 'md', 'MarkerSize', MS, 'MarkerFaceColor', 'm');
    %hold on;plot(TV2_ph(nr)-REF, loc(nr), 'kd', 'MarkerSize', MS, 'MarkerFaceColor', 'w');
end






figure(7);clf;
loc=1:5; %location for each run on the phase plot
REF=monthday(6)+21; %reference date
MS=18;
for nr=1:5;
    hold on;plot(HEAT_ph(nr)-REF, loc(nr), 'bd', 'MarkerSize',MS,'LineWidth', 3, 'MarkerFaceColor', 'r');
    hold on;plot(TV_ph(nr)-REF, loc(nr), 'kd', 'MarkerSize', MS, 'MarkerFaceColor', 'k');
    hold on;plot(TETEN_ph(nr)-REF, loc(nr), 'cd', 'MarkerSize', MS, 'MarkerFaceColor', 'c');
    hold on;plot(OLR_ph(nr)-REF, loc(nr), 'gd', 'MarkerSize', MS, 'MarkerFaceColor', 'g');
    hold on;plot(TEDIV_ph(nr)-REF, loc(nr), 'md', 'MarkerSize', MS, 'MarkerFaceColor', 'm');
    %hold on;plot(TV2_ph(nr)-REF, loc(nr), 'kd', 'MarkerSize', MS, 'MarkerFaceColor', 'w');
end

figure(7); hold on;
legend('Atmospheric Heating', 'Tropospheric Temperature', 'Column Tendency', 'OLR', 'Heat Transport div.');%,'Stratospheric temp.')
ymn=0.5; ymx=7.5;
plot([0 0], [ymn ymx], 'r--', 'LineWidth', 1.5); hold on;
plot(91.2*[1 1], [ymn ymx], 'k--', 'LineWidth', 1.5); hold on;
xlabel('Days past the summer solstice', 'FontSize', FSz, 'FontName', FN);
ylabel('Mixed Layer Depth (meters)', 'FontSize', FSz, 'FontName', FN);
title({'Seasonal Phasing of energy fluxes and'; 'temperature as a function of mixed layer depth'}, 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTick', 1:7,'YTickLabel', {[num2str(2.4) '  ']; [num2str(6) '  ']; [num2str(12) '  ']; [num2str(24) '  ']; [num2str(50) '  ' ]; 'NH Obs.  '; 'SH Obs.  '})
axis([-60 100 ymn ymx])
box on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% AMPLITUDES
figure(8); clf;
for nr=1:5;
    ref=SWABS_AMP(nr)+SHF_AMP(nr); %reference everything to the total amplitude of heatinf
    hold on;plot(SWABS_AMP(nr)/ref, loc(nr), 'rd', 'MarkerSize', MS, 'MarkerFaceColor', 'r');
     hold on;plot(SHF_AMP(nr)/ref, loc(nr), 'bd', 'MarkerSize', MS, 'MarkerFaceColor', 'b');
   
     if nr == 4;
    hold on;plot((-TETEN_AMP(nr)/ref)+.01, loc(nr), 'cd', 'MarkerSize', MS, 'MarkerFaceColor', col, 'MarkerEdgeColor', col);
     else
            hold on;plot((-TETEN_AMP(nr)/ref), loc(nr), 'cd', 'MarkerSize', MS, 'MarkerFaceColor', col, 'MarkerEdgeColor', col);
     end
    hold on;plot(-OLR_AMP(nr)/ref, loc(nr), 'gd', 'MarkerSize', MS, 'MarkerFaceColor', 'g');
    hold on;plot(-TEDIV_AMP(nr)/ref, loc(nr), 'md', 'MarkerSize', MS, 'MarkerFaceColor', 'm');
    %hold on;plot(TV2_ph(nr)-ref, loc(nr), 'kd', 'MarkerSize', MS, 'MarkerFaceColor', 'w');
end

%plot the relative amplitudes of the terms (relative to the total heatsing
figure(8); hold on;
legend('SWABS', 'Surface Heating', 'Column Tendency', 'OLR', 'Heat Transport div.');%,'Stratospheric temp.')
xmn=-1; xmx=1.65;
plot([0 0], [ymn ymx], 'r--', 'LineWidth', 1.5); hold on;
%plot(91.2*[1 1], [ymn ymx], 'k--', 'LineWidth', 1.5); hold on;
xlabel('Seasonal amplitude of enrgy flux in phase with total heating', 'FontSize', FSz, 'FontName', FN);
ylabel('Mixed Layer Depth (meters)', 'FontSize', FSz, 'FontName', FN);
title({'Seasonal Amplitude of energy fluxes and'; 'normalized by magnitude of total heating'}, 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTick', 1:7,'YTickLabel', {[num2str(2.4) '  ']; [num2str(6) '  ']; [num2str(12) '  ']; [num2str(24) '  ']; [num2str(50) '  ' ]; 'NH Obs.  '; 'SH Obs.  '})
axis([xmn xmx ymn ymx]); box on;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% observations
 load C:\Users\thedhoe\Desktop\CMIP_seasonal_lap\obs\data\seasonal_heat2.mat

load C:\Users\thedhoe\Desktop\CMIP_seasonal_lap\obs\data\MASK.mat
load('C:\research\aqua_depth\figures\energy_flux_phases\T_trop_obs.mat','readme','tlat','tlon','plev', 'VT'); %the vertically averaged temperature

%interpolate VT to the grids everything else is on
[TLON,TLAT]=meshgrid(tlon, tlat);
[LON, LAT]=meshgrid(lon,lat);[NN, N]=size(LON);

for mn=1:12;
    vt(mn,1:NN,1:N)=interp2(TLON, TLAT, squeeze(VT(mn,:,:)), LON, LAT);
end

NN=length(lat);
ABS=DT-UT-DS+US;
%compute the annual means first
MAX=250; MIN=0;
FSz=25; FN='Arial';
%sphcolor2jet(mean(FS(:,:,:)),lat,lon,1); caxis([MIN MAX]); title('Annual mean surface heating (W m^-^2)', 'FontSize', FSz, 'FontName', FN); set(gca, 'FontSize', FSz, 'FontName', FN)
%sphcolor2jet(mean(ABS(:,:,:)),lat,lon,2); caxis([MIN MAX]); title('Annual mean atmospheric solar absorption (W m^-^2)', 'FontSize', FSz, 'FontName', FN); set(gca, 'FontSize', FSz, 'FontName', FN)

%zonal means of ABS and FS
LW=2; ymx=1.1;
x=sind(lat);


%%%%%%%%%%%%%%%            now the seasonal cycle   %%%%%%%%%%%%%%%%%%
%define zonal mean anomalies from the annual mean
ZABS=squeeze(mean(shiftdim(ABS,2))) - ones(12,1)*mean(squeeze(mean(shiftdim(ABS,2))));
ZFS=squeeze(mean(shiftdim(FS,2))) - ones(12,1)*mean(squeeze(mean(shiftdim(FS,2))));
ZTETEN=squeeze(mean(shiftdim(TETEN,2))) - ones(12,1)*mean(squeeze(mean(shiftdim(TETEN,2))));
ZTEDIV=squeeze(mean(shiftdim(TEDIV,2))) - ones(12,1)*mean(squeeze(mean(shiftdim(TEDIV,2))));
ZVT=squeeze(mean(shiftdim(vt,2))) - ones(12,1)*mean(squeeze(mean(shiftdim(vt,2))));
ZOLR=squeeze(mean(shiftdim(OLR,2))) - ones(12,1)*mean(squeeze(mean(shiftdim(OLR,2))));
ZDT=squeeze(mean(shiftdim(DT,2))) - ones(12,1)*mean(squeeze(mean(shiftdim(DT,2))));
min_lat=cross; 
max_lat=90;

%average over this domain - first creat monthly time series
donn=find(lat>min_lat & lat < max_lat); %latitude indexes within the domain
doss=find(lat>-max_lat & lat<= -min_lat); %southern hemisphere
Rx=sind(lat(donn)); NORM=trapz(Rx,ones(size(Rx))); %the sin of latitude vector over the region and any area weighting normalization
Rs=sind(lat(doss));
SNORM =trapz(Rs,ones(size(Rs))); 


%monthly averages over the domain
NR_SHF(1:12) = trapz(Rx, squeeze(ZFS(:,donn))')/NORM; % average surface heat flux anomaly (from annual mean) over the region (W m^-2)
NR_SWABS(1:12) = trapz(Rx, squeeze(ZABS(:,donn))')/NORM; % average direct SW absorption anomaly (from annual mean) over the region (W m^-2)
NR_TEDIV(1:12) = trapz(Rx, squeeze(ZTEDIV(:,donn))')/NORM; % average direct SW absorption anomaly (from annual mean) over the region (W m^-2)
NR_TETEN(1:12) = trapz(Rx, squeeze(ZTETEN(:,donn))')/NORM; % average direct SW absorption anomaly (from annual mean) over the region (W m^-2)
NR_OLR(1:12) = trapz(Rx, squeeze(ZOLR(:,donn))')/NORM; % average direct SW absorption anomaly (from annual mean) over the region (W m^-2)
NR_VT(1:12) = trapz(Rx, squeeze(ZVT(:,donn))')/NORM ; %vertically averaged temperature

SR_SHF(1:12) = trapz(Rs, squeeze(ZFS(:,doss))')/NORM; % average surface heat flux anomaly (from annual mean) over the region (W m^-2)
SR_SWABS(1:12) = trapz(Rs, squeeze(ZABS(:,doss))')/NORM; % average direct SW absorption anomaly (from annual mean) over the region (W m^-2)
SR_TEDIV(1:12) = trapz(Rs, squeeze(ZTEDIV(:,doss))')/NORM; % average direct SW absorption anomaly (from annual mean) over the region (W m^-2)
SR_TETEN(1:12) = trapz(Rs, squeeze(ZTETEN(:,doss))')/NORM; % average direct SW absorption anomaly (from annual mean) over the region (W m^-2)
SR_OLR(1:12) = trapz(Rs, squeeze(ZOLR(:,doss))')/NORM; % average direct SW absorption anomaly (from annual mean) over the region (W m^-2)
SR_VT(1:12) = trapz(Rs, squeeze(ZVT(:,doss))')/NORM;

        
        TS=NR_SHF+NR_SWABS;
        Nind=(2^.5*(12/11)^.5)*(TS-mean(TS))./std(TS); %dot and average to get the amplitude in phase with the insolation
        
        TS=SR_SHF+SR_SWABS;
        Sind=(2^.5*(12/11)^.5)*(TS-mean(TS))./std(TS); %dot and average to get the amplitude in phase with the insolation
        
        
        SWABS_AMP_O=[ mean(Nind.*NR_SWABS) mean(Sind.*SR_SWABS)];
        SHF_AMP_O=[mean(Nind.*NR_SHF) mean(Sind.*SR_SHF) ];
        TEDIV_AMP_O=[mean(Nind.*NR_TEDIV) mean(Sind.*SR_TEDIV)];
        OLR_AMP_O=[mean(Nind.*NR_OLR) mean(Sind.*SR_OLR)] ;
        TETEN_AMP_O= [mean(Nind.*NR_TETEN) mean(Sind.*SR_TETEN)];
        VT_AMP_O= [mean(Nind.*NR_VT) mean(Sind.*SR_VT)];
        
        SWABS_PH_O=[ fund_phase(NR_SWABS) fund_phase(SR_SWABS)];
        SHF_PH_O=[fund_phase(NR_SHF) fund_phase(SR_SHF) ];
        TEDIV_PH_O=[fund_phase(NR_TEDIV) fund_phase(SR_TEDIV)];
        OLR_PH_O=[fund_phase(NR_OLR) fund_phase(SR_OLR)] ;
        TETEN_PH_O= [fund_phase(NR_TETEN) fund_phase(SR_TETEN)];
        VT_PH_O= [fund_phase(NR_VT) fund_phase(SR_VT)];
        HEAT_PH_O=[fund_phase(NR_SHF+NR_SWABS) fund_phase(SR_SHF+SR_SWABS)] ;
        
        LOC=6:7; %yaxis location of the observations on Figure 8
        figure(7); hold on;
        OREF=[monthday(6)+21 monthday(12)+21]; %soltice dates in each hemisphere
        
        DATA=[HEAT_PH_O; VT_PH_O; TETEN_PH_O; OLR_PH_O; TEDIV_PH_O]- (OREF'*[1  1   1  1 1])';
        DATA(find(DATA<-180)) = DATA(find(DATA<-180))+365; %add a year to the ones that go through day 0
        
        for nr=1:2;
    hold on;plot(DATA(1,nr), LOC(nr), 'bd', 'MarkerSize', MS,'LineWidth',3, 'MarkerFaceColor', 'r');
    hold on;plot(DATA(2,nr), LOC(nr), 'kd', 'MarkerSize', MS, 'MarkerFaceColor', 'k');
    hold on;plot(DATA(3,nr), LOC(nr), 'cd', 'MarkerSize', MS, 'MarkerFaceColor', 'c');
    hold on;plot(DATA(4,nr), LOC(nr), 'gd', 'MarkerSize', MS, 'MarkerFaceColor', 'g');
    hold on;plot(DATA(5,nr), LOC(nr), 'md', 'MarkerSize', MS, 'MarkerFaceColor', 'm');
    %hold on;plot(TV2_ph(nr)-REF, loc(nr), 'kd', 'MarkerSize', MS, 'MarkerFaceColor', 'w');
        end

        figure(8);hold on;

        for nr=1:2;
    ref=SWABS_AMP_O(nr)+SHF_AMP_O(nr); %reference everything to the total amplitude of heatinf
    hold on;plot(SWABS_AMP_O(nr)/ref, LOC(nr), 'rd', 'MarkerSize', MS, 'MarkerFaceColor', 'r');
    hold on;plot(SHF_AMP_O(nr)/ref, LOC(nr), 'bd', 'MarkerSize', MS, 'MarkerFaceColor', 'b');
    hold on;plot((-TETEN_AMP_O(nr)/ref), LOC(nr), 'cd', 'MarkerSize', MS, 'MarkerFaceColor', col, 'MarkerEdgeColor', col);
   hold on; plot(-OLR_AMP_O(nr)/ref, LOC(nr), 'gd', 'MarkerSize', MS, 'MarkerFaceColor', 'g');
    hold on;plot(-TEDIV_AMP_O(nr)/ref, LOC(nr), 'md', 'MarkerSize', MS, 'MarkerFaceColor', 'm');
        end