%%%%%%%%%%%%%%%%%%%%%%
%% Compile Trizaxus Data %%
%%%%%%%%%%%%%%%%%%%%%%
% Compile Profile data from Triaxus
% each deployment as seprate legs which have a variety of profiles

% MLD computing algorithm from Talle & holte 2008 
%     http://journals.ametsoc.org/doi/abs/10.1175/2009JTECHO543.1

%clear all

% Triaxus data
 TRIpath = '/Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus/Triaxus/processing/in2020_v08/cap/cappro/avg';

% Add utility paths% add orther relevant paths
addpath(genpath('/Users/roh029/Desktop/Matlab/SOLACE'))
addpath(genpath('/Users/roh029/Desktop/Matlab/Utilities'))

% load
 cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
 load TRI_MLD
 
 
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
%load correctons

cd /Users/roh029/Desktop/Data.nosync/SOLACE/CHL
load('ChlFlo_Ctd_Corrections.mat')

%% Profiles on CTD casts
% get discrete and sensored
% find the resolution of sensors, smooth and interpolate to say every
% meter. Add discrete samples. 

cd(TRIpath)
File = dir('*g.mat'); file={File(:).name};

% pick min dpeth
mind=25;

% 
% % % Pick Best MLD by Eye  --- good up to 30
%   eyetest(i) = 4;
%   
%   does the mld vary over the leg?
%   unimld(i) =  true;
%   
%   Assign MLd and range
%   tri_mld(i) = tri_mld_all(i,eyetest(i));
%   tri_min(i) = min(tri_mld_all(i,:)); 
%   tri_max(i) = max(tri_mld_all(i,:)); 

%
 
% flag:17, 16, 15, 14,13, 12, 11,9,4
for i = 20:30 %1:length(file) %[1:3,5:8,10,18:length(file)] % 4:19 for SOTS. % 20:30 for SOuth
  cd(TRIpath)
  fprintf(['\n Building Deployment: ' num2str(i)])
  load(file{i});
  
  cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
if i < 20
    load('poc_corrections_sots.mat')
else
    load('poc_corrections_south.mat')
end

% Build profiles for full leg
    
    % Pressure (dbar)
    Pres = avgStruct.sensor.pressure.data; 
    
    % Wetlabs Chl conc (mg/m3)
    Chl = avgStruct.sensor.chlorophyll.data;   
    Chl = chl_ctd_cf_south(1) + chl_ctd_cf_south(2).*Chl;  % Correct based on ctd - shipboard correction

    % PAR (uE/m^2/sec)
    Par = avgStruct.sensor.par.data;   
    
    % Transmisometer (%)
    Trans = avgStruct.sensor.transmissometer.data; 
    Trans = Trans - trans_offset;
    c = -(1/.25) .* log(Trans ./100);    % Calibrated for SAZ at SOTS using same Wetlabs CStar beam
                                       % transmissometer (25 cm path
                                       % length). Units of m^-1
    cp  = c - c_sw_ctd;                % Contirbution from particulate equals total attenutation minus csw;  Units of m^-1
                                       % Using schallenbergMultiYearObservationsFluorescence2019
                                       % relationship, calibrated to SOTS,
                                       % but with 2 transfomations
                                       
    POC_sol = cf_cp_poc_0.*cp;         % POC using SOLACe measurements

    % Chl:Car
    ChlPOC = Chl./(POC_sol.*12);
     
    % Temperature (C) - Averaged across both sensors
    Temp = nanmean([avgStruct.sensor.temperature.data; ...
                                             avgStruct.sensor.temperature_2.data]);
                                         
    % Salinity (PSU) - Averaged across both sensors
    Sal = nanmean([avgStruct.sensor.salinity.data; ...
                                            avgStruct.sensor.salinity_2.data]);
                                        
    % oxygen (uM or mmol/m3) - Averaged across both sensors
    O2 =  nanmean([avgStruct.sensor.oxygen.data; ...
                                           avgStruct.sensor.oxygen_2.data]);
                                         
   
    
    
    %% sepreate into discrete profiles
    diffp = diff(Pres);
    diffp(diffp>.1) =1; diffp(diffp<-.1) = -1;
    
    des = find(diffp==-1);
    asc = find(diffp==1);
    
    di = find(diff(des)~=1); 
    DES = des(di);
    ai = find(diff(asc)~=1);
    ASC = asc(ai);
    
    prof_ind =sort([ASC,DES,des(end),asc(end)]);
   
    % number of profiles
      profnums(i) = length(prof_ind); % ups and downs i belive

%
% %     %% plot
%             yyaxis left
%             plot(Pres)
%             for i= 1:length(prof_ind); xline(prof_ind(i),'r','linewidth',2); end
%             set(gca,'ydir','rev')
%             yyaxis right
%             plot(diffp)
%             
    
    
    %% Average each profile
    % reset full leg averages
      tri_mld_p = [];  tri_time_utc_p = [];
      
      % Make dummies for each profile leg
      PRES   = nan(length(prof_ind)-2,200);
      TEMP   = nan(length(prof_ind)-2,200);
      SAL    = nan(length(prof_ind)-2,200);
     
      % Loop through and average all profies in leg
      clearvars tri_time_loc_p tri_time_utc_p tri_lat_p tri_lon_p tri_mld_p data_dis_ts
      for p = 1:length(prof_ind)-2  
        pind = (prof_ind(p)+1:prof_ind(p+1)-1); 
        %Only compute if section is 100m deep
        if ~isempty(pind)
        if pind(end) - pind(1) > 100
        pres = Pres(pind); [pres, dirind] = sort(pres); 
        temp = Temp(pind);      temp = temp(dirind);        temp=movmean(temp,5);
        sal  = Sal(pind);       sal = sal(dirind);          sal=movmean(sal,5);
        chl  = Chl(pind);       chl = chl(dirind);          chl=movmean(chl,10,'omitnan');
        poc  = POC_sol(pind);   poc = poc(dirind);          poc=movmean(poc,10,'omitnan');
        chlc = ChlPOC(pind);    chlc = chlc(dirind);        chlc=movmean(chlc,10,'omitnan');
        
        % remove repeat pressure and remove first in pair
        dp_ind = find(diff(pres)==0);
        pres(dp_ind) = [];  pres_dcm =pres;
        temp(dp_ind) = [];
        sal(dp_ind) = [];
        chl(dp_ind) = [];
        poc(dp_ind) = [];
        chlc(dp_ind) = [];
     
    % Find MLd of each profile
        % MLD 
            % Remove min dpeth
            pres(1:mind-1) = [];    temp(1:mind-1) = [];    sal(1:mind-1) = [];

            % remove nans
            nanind = any([isnan(pres);isnan(temp);isnan(sal)],1);
            pres(nanind) = [];      temp(nanind) = [];      sal(nanind) = [];

            % Find MLD
            findmld
            % using all methods
            tri_mld_p(p) = mean([gtmldp,mixedtp, mixedsp, mixeddp, mldepthptmpp, mldepthdensp, gdmldp]);
            % using just density thresh and grad
            tri_mld_p(p) = mean([gdmldp,mldepthdensp]);
            data_dis_ts(13,2,p) = tri_mld_p(p);
            
       % Meta Data
          tri_time_utc_p(p) =  nanmean(avgStruct.time(pind));
          tri_time_loc_p(p) =  tri_time_utc_p(p) +11/24;
          tri_lat_p(p)      =  nanmean(avgStruct.latitude(pind));
          tri_lon_p(p)      =  nanmean(avgStruct.longitude(pind));    
       
    % FIND DCM 
        prof = chl;   
        [max_m, max_i] = max(prof);
        % Assign metrics if theres a DCM
        if pres_dcm(squeeze(max_i))>tri_mld_p(p) && max_m> 1.1.*nanmean(prof(pres_dcm<tri_mld_p(p))) %  % IS DCM BELOW MLD
        dcm_thr =  nanmean(prof(pres_dcm<tri_mld_p(p))) + .5.* abs(max_m-nanmean(prof(pres_dcm<tri_mld_p(p))));  % 50% of difference between MLD and max value 
        DCM = find(prof>dcm_thr) ; 
        DCM( pres_dcm(DCM) < (tri_mld_p(p)-30)) = [] ;  % Dont condsider values 15 m above MLD
        data_dis_ts(4,1,p) = squeeze(max_m); % Size of DCM
        data_dis_ts(4,2,p) = pres_dcm(squeeze(max_i)); % Depth of DCM
        data_dis_ts(4,3,p) = pres_dcm(DCM(1));         % Top of DCM
        data_dis_ts(4,4,p) = pres_dcm(DCM(end));       % Base of DCM
        data_dis_ts(4,5,p) = pres_dcm(DCM(end)) - pres_dcm(DCM(1));    % Width of DCM (m)
        data_dis_ts(4,6,p) = nanmean(prof(DCM)).*(pres_dcm(DCM(end)) - pres_dcm(DCM(1))); % integrate of DCM (mg/m2)
        data_dis_ts(4,7,p) = 100.*(max_m-nanmean(prof(pres_dcm<tri_mld_p(p))))./nanmean(prof(pres_dcm<tri_mld_p(p))); 
        data_dis_ts(4,8,p) = 100.*(nanmean(prof(DCM))-nanmean(prof(pres_dcm<tri_mld_p(p))))./nanmean(prof(pres_dcm<tri_mld_p(p))); 
        else
        data_dis_ts(4,1:8,p) = nan;
        end
       
       
    
     % Find DBM 
        prof = poc;        
        [max_m, max_i] = max(prof);
        % Assign metrics if theres a DCM
        if pres_dcm(squeeze(max_i))>tri_mld_p(p)  && max_m> 1.1.*nanmean(prof(pres_dcm<tri_mld_p(p)))%  % IS DCM BELOW MLD
        dcm_thr =  nanmean(prof(pres_dcm<tri_mld_p(p))) + .5.* abs(nanmean(prof(pres_dcm<tri_mld_p(p)))- max_m);  % 50% of difference between MLD and max value 
        DBM = find(prof>dcm_thr) ; 
        DBM( pres_dcm(DBM) < (tri_mld_p(p)-30)) = [] ;  % Dont condsider values 15 m above MLD
        data_dis_ts(11,3,p) = pres_dcm(DBM(1));         % Top of DCM
        data_dis_ts(11,4,p) = pres_dcm(DBM(end));       % Base of DCM
        data_dis_ts(11,5,p) = pres_dcm(DBM(end)) - pres_dcm(DBM(1));    % Width of DCM (m)
        data_dis_ts(11,6,p) = nanmean(prof(DBM)).*(pres_dcm(DBM(end)) - pres_dcm(DBM(1))); % integrate of DCM (mg/m2)
        data_dis_ts(11,7,p) = 100.*(max_m-nanmean(prof(pres_dcm<tri_mld_p(p))))./nanmean(prof(pres_dcm<tri_mld_p(p))); 
        else
        data_dis_ts(11,3:7,p) = nan;
        end
        data_dis_ts(11,1,p) = squeeze(max_m); % Size of DCM
        data_dis_ts(11,2,p) = pres_dcm(squeeze(max_i)); % Depth of DCM
        data_dis_ts(11,8,p) = 100.*(nanmean(prof(DCM))-nanmean(prof(pres_dcm<tri_mld_p(p))))./nanmean(prof(pres_dcm<tri_mld_p(p))); 
        

     
    % Find C Chl Max
        prof = chlc;        
        [max_m, max_i] = max(prof);
        % Assign metrics if theres a DCM
        if pres_dcm(squeeze(max_i))>tri_mld_p(p)  && max_m> 1.1.*nanmean(prof(pres_dcm<tri_mld_p(p)))%  % IS DCM BELOW MLD
        dcm_thr =  nanmean(prof(pres_dcm<tri_mld_p(p))) + .5.* abs(nanmean(prof(pres_dcm<tri_mld_p(p)))- max_m);  % 50% of difference between MLD and max value 
        DCchM = find(prof>dcm_thr) ; 
        DCchM( pres_dcm(DCchM) < (tri_mld_p(p)-30)) = [] ;  % Dont condsider values 15 m above MLD
        data_dis_ts(12,3,p) = pres_dcm(DCchM(1));         % Top of DCM
        data_dis_ts(12,4,p) = pres_dcm(DCchM(end));       % Base of DCM
        data_dis_ts(12,5,p) = pres_dcm(DCchM(end)) - pres_dcm(DCchM(1));    % Width of DCM (m)
        data_dis_ts(12,6,p) = nanmean(prof(DCchM)).*(pres_dcm(DCchM(end)) - pres_dcm(DCchM(1))); % integrate of DCM (mg/m2)
        data_dis_ts(12,7,p) = 100.*(1./max_m - 1./nanmean(prof(pres_dcm<tri_mld_p(p))))./(1./nanmean(prof(pres_dcm<tri_mld_p(p)))); 
        else
        data_dis_ts(12,3:7,p) = nan;
        end
        data_dis_ts(12,1,p) = 1./squeeze(max_m);        % Size of DCM
        data_dis_ts(12,2,p) = pres_dcm(squeeze(max_i)); % Depth of DCM
        data_dis_ts(12,8,p) = 100.*(1./nanmean(prof(DCM))-1./nanmean(prof(pres_dcm<tri_mld_p(p))))./(1./nanmean(prof(pres_dcm<tri_mld_p(p)))); 
        
        
        % Sample plot
        if false
            figure(1)
            subplot(1,3,1) 
            plot( chl,pres_dcm);  
            if ~isnan(data_dis_ts(4,2,p))
                yline(data_dis_ts(4,2,p)); yline(data_dis_ts(4,3,p)); yline(data_dis_ts(4,4,p))
            end
            title( num2str(data_dis_ts(4,8,p)))
            xlabel('chl'); set(gca,'ydir','reverse')
               subplot(1,3,2) 
            plot( poc,pres_dcm);  
            if ~isnan(data_dis_ts(11,3,p))
                yline(data_dis_ts(11,2,p)); yline(data_dis_ts(11,3,p)); yline(data_dis_ts(11,4,p))
            end
            title( num2str(data_dis_ts(11,8,p)))
            xlabel('poc');set(gca,'ydir','reverse')
               subplot(1,3,3) 
            plot( chlc,pres_dcm);  
            if ~isnan(data_dis_ts(12,3,p))
                yline(data_dis_ts(12,2,p)); yline(data_dis_ts(12,3,p)); yline(data_dis_ts(12,4,p))
            end
            title( num2str(data_dis_ts(12,8,p)))
            xlabel('chl:c');set(gca,'ydir','reverse')
        end
        
        
        % Fill out dummies
        PRES(p,round(pres)) =pres;
        TEMP(p,round(pres)) = temp;
        SAL(p,round(pres)) = sal;
        
        % Remove bad values
          badi = find(tri_mld_p==0);
          tri_time_utc_p(badi)  = [];
          tri_time_loc_p(badi)  = [];
          tri_lat_p(badi)       = [];
          tri_lon_p(badi)       = [];
          tri_mld_p(badi)       = [];
          data_dis_ts(:,:,badi) = [];
          
        end
        end
      end
      
   % Transform lat lon to distance.
     tri_dis_p =[deg2km(distance('rh',tri_lat_p(1),tri_lon_p(1),avgStruct.latitude(1),avgStruct.longitude(1))), ...
                 deg2km(distance('rh',tri_lat_p(1),tri_lon_p(1),avgStruct.latitude(1),avgStruct.longitude(1)))+ ...
                 cumsum(deg2km(distance('rh',tri_lat_p(1:end-1),tri_lon_p(1:end-1),tri_lat_p(2:end),tri_lon_p(2:end))))];
          
     
      % save DCM and DBP profiles
      dum = file{i};
       FNAME = ['Tow_',dum(13),'_Leg_',dum(16),'.mat'] ;
       cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
        save(FNAME, 'data_dis_ts', 'tri_time_loc_p', 'tri_mld_p','tri_dis_p','-append')
        save(FNAME, 'tri_time_loc_p', 'tri_time_utc_p', 'tri_lat_p', 'tri_lon_p', 'tri_mld_p','-append')
        
      % Average Profile
           % Take average of all profiles in leg
           tri_time_loc(i)    = mean(tri_time_loc_p);
           tri_time_utc(i)    = mean(tri_time_utc_p);
           tri_lat(i)    = mean(tri_lat_p);
           tri_lon(i)    = mean(tri_lon_p);
       %
            % MLD 
            % tuncate to 100 for 31-33 legs
            pres = nanmean(PRES(:,1:200));  temp = nanmean(TEMP(:,1:200));  sal = nanmean(SAL(:,1:200)); 

            % Remove min dpeth
            pres(1:mind-1) = [];    temp(1:mind-1) = [];    sal(1:mind-1) = [];

            % remove nans
            nanind = any([isnan(pres);isnan(temp);isnan(sal)],1);
            pres(nanind) = [];      temp(nanind) = [];      sal(nanind) = [];

            % Find MLD
            findmld
            castmld_dg(i)  = gdmldp;          % density gradient method
            
            castmld_sp(i)  =  mixedsp ;        % from salinity algorithm
            castmld_tp(i)  =  mixedtp;         % from temperature algorithm
            castmld_dp(i)  =  mixeddp;         % from density algorithm
            
            castmld_tth(i) =  mldepthptmpp;   % temperature threshold method
            castmld_dth(i) =  mldepthdensp;   % density threshold method
            castmld_tg(i)  =  gtmldp;          % temperature gradient method
            
            
            tri_mld_all(i,:) = [gtmldp,mixedtp, mixedsp, mixeddp, mldepthptmpp, mldepthdensp, gdmldp];
    
       % plot
       if false
           for plott=1
        figure('Units','normalized','Position',[-.2 2 .5 1])
        clc
        set(gcf,'color','w');
        

        axPos = get(gca, 'Position');
        set(gca, 'Position', axPos + [0 0.2 0 -0.3]);
         line(temp,pres,'Color','k','linewidth',2)
         
         yline(castmld_dg(i),'m:','linewidth',2);     % D Grdient
         
         yline(castmld_sp(i),'r','linewidth',2);     % S Algorithm
         yline(castmld_tp(i),'k','linewidth',2);     % T Algorithm
         yline(castmld_dp(i),'m','linewidth',2);     % D Algorithm
         
         yline(castmld_tth(i),'k--','linewidth',2);     % T Th
         yline(castmld_tg(i),'k:','linewidth',2);     % T Grad
         
         yline(castmld_dth(i),'m--','linewidth',2);     % T Th

          set(gca,'ydir','rev');  
          ylim([10,200])
          ax1 = gca; % current axes
          ax1.XColor = 'k';
          ylabel('Depth (dbar)','interpreter','latex','fontsize',18)
          xlabel('Temperature (C)','interpreter','latex','fontsize',18)

         ax1_pos = ax1.Position;
         ax2 = axes('Position',ax1_pos,...
            'XAxisLocation','top',...
            'YAxisLocation','right',...
            'Color','none',...
            'xtick',[]);

        line(sal,pres,'Parent',ax2,'Color','r','linewidth',2)
         set(gca,'ydir','rev'); 
         ylim([10,200])

         ax3 = axes('position', (axPos .* [1 1 1 1e-3]) + [0 0.0 0 0], 'XColor', 'r', 'linewidth', 2);
        set(ax3, 'xlim', [ax3.XLim]);
        xlabel('Salinity (PSU)','interpreter','latex','fontsize',18)
        
        % Axis Five: Density
    ax1_pos = ax1.Position; 
    ax8 = axes('Position',ax1_pos,...
        'XAxisLocation','top',...
        'YAxisLocation','right',...
        'Color','none',...
        'xtick',[]);

        line(pden+1000,pres,'Parent',ax8,'Color','m','linewidth',2)
        set(gca,'ydir','rev'); 
        ylim([10,200])

        ax9 = axes('position', (axPos .* [1 1 1 1e-3]) - [0 0.06 0 0], 'XColor', 'm', 'linewidth', 2);
        set(ax9, 'xlim', [ax8.XLim]);
        xlabel('Density (kg/m3)','interpreter','latex','fontsize',18)

        title(ax1,['CTD Cast: ', num2str(i), '  | DG: ', num2str(castmld_dg(i)), ' (m), SA:' num2str(castmld_sp(i)),...
              ' (m), TA:' num2str(castmld_tp(i)),' (m), DA:' num2str(castmld_dp(i))],'interpreter','latex','fontsize',25)

           end
           end
%     
end

%% save
% cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
% save TRI_MLD tri_mld tri_min tri_max tri_mld_all unimld eyetest...
%            tri_lon tri_lat tri_time_loc tri_time_utc tri_time_loc_p
load TRI_MLD 

%% Load FiRE
 fprintf(['FIRe Files \n'])
for i = 9:14% 1:length(Fire_files) % 9:14 appear to be deployments 4 adn 5
    cd  '/Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus/Triaxus/processing/in2020_v08/FIRE'
    Fire_files= dir('*Triaxus*'); Fire_Files={Fire_files(:).name};
    
    % Pull data
    fire_time  = ncread(Fire_Files{i},'time') + datenum('30-Dec-1899 00:00:00') +11/24 ;  
          % 'days since 1899-12-30 00:00:00 UTC' converted to local time in days since Jan 01 2000
    fire_pres  = ncread(Fire_Files{i},'pressure')./10000;  % Pressure
          % Converted from Pascals to dBar
    fire_fv_fm = ncread(Fire_Files{i},'fv_fm');   % 'MaximumQuantum yield of photochemistry inPSII'
          % 'dimensionless'
    fire_Abs_sig = ncread(Fire_Files{i},'absolute_sigma');   % 'absolute functional absorption cross-section sigma'
          % 'Dimensionsless'
    fire_Rel_sig = ncread(Fire_Files{i},'relative_sigma');   % 'relative functional absorption cross-section sigma'
          % 'Dimensionsless'
    fprintf(['Start Time:', datestr(fire_time(1)), ' | End Time:', datestr(fire_time(end)), '\n'])
    
    % Save Fire Data With Ctd Data
    cd(TRIpath)
    File = dir('*g.mat'); file={File(:).name};
    fname = file{i+11};  % This moves through legs 4 adn 5 by hand...Not elegant
    save(fname, 'fire_time','fire_pres','fire_fv_fm','fire_Abs_sig','fire_Rel_sig','-append')
end
  

%  fprintf(['CTDs \n'])
% cd(TRIpath)
% File = dir('*g.mat'); file={File(:).name};
% for i = 20:25 %1:length(file)
%     load(file{i});
%     time = avgStruct.time;   % in utc times
%     fprintf(['Start Time:', datestr(time(1)), ' | End Time:', datestr(time(end)), '\n'])
% end
% 
%%
% cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
% FIRe = readtable('FIRe_4_5.csv');
% FIRe_day = squeeze(FIRe{:,3});
% FIRe_time = squeeze(FIRe{:,4});
% FIRe_FvFm = squeeze(FIRe{:,9});
% FIRe_Sigma = squeeze(FIRe{:,12});
% [X Y Z] = csvimport('FIRe_4_5.csv', 'columns', {'B', 'C', 'L'});
% 

%% Interpolator 
TRIpath = '/Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus/Triaxus/processing/in2020_v08/cap/cappro/avg';
cd(TRIpath)
File = dir('*g.mat'); file={File(:).name};

for TOWS =5:7
    fprintf(['\n Starting TOW: ',num2str(TOWS),'\n']) 
    
for siteDescription =1
% Sots 
% Tow 1
  if TOWS == 1
  filesN = [4,5,6];
  fname = 'Tow_1';
  
% Tow 2
  elseif TOWS == 2
  filesN = 7:10;
  fname = 'Tow_2';
  
% Tow 3
  elseif TOWS == 3
  filesN = 11:19;
  fname = 'Tow_3';
  
% Southern Site 1
% Tow 4
  elseif TOWS == 4
  filesN = [20,21,22];
  fname = 'Tow_4';
% Tow 5
  elseif TOWS == 5
  filesN = [23,24,25];
  fname = 'Tow_5'; 
% Southern Site 2
  % Tow 6
  elseif TOWS == 6
  filesN = [26,27];
  fname = 'Tow_6'; 
  % Tow 7
  elseif TOWS == 7
  filesN = [28,29,30];
  fname = 'Tow_7'; 
  end
end
 

for leg = 1:length(filesN)
    leg
clearvars data data_ts firedata firedata_raw
cd(TRIpath)
load(file{filesN(leg)});
FNAME = [fname,'_Leg_',num2str(leg)] ;   

% CTD Data
if true
for CTD_Vars = 1 

% variables
pres   = avgStruct.sensor.pressure.data;
time   = avgStruct.time + 11/24;           % in local times
lat    = avgStruct.latitude;
lon    = avgStruct.longitude;
chl    = avgStruct.sensor.chlorophyll.data;
dis_scl = 1;
dist_m =[0, cumsum(deg2km(distance('rh',lat(1:end-1),lon(1:end-1),lat(2:end),lon(2:end))).*dis_scl)];
if TOWS<4
chl    = chl_ctd_cf_sots(1) + chl_ctd_cf_sots(2).*chl;          % Correct based on ctd - shipboard correction for SOTS
else
chl    = chl_ctd_cf_south(1) + chl_ctd_cf_south(2).*chl;        % Correct based on ctd - shipboard correction for SOUTH 
end 
o2     = mean([avgStruct.sensor.oxygen.data;      avgStruct.sensor.oxygen_2.data]);
temp   = mean([avgStruct.sensor.temperature.data; avgStruct.sensor.temperature_2.data]);
trans  = avgStruct.sensor.transmissometer.data;
trans  = trans - trans_offset; % subtract offset with CTd Tranmisometer
% POC
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
if TOWS < 4 
    load('poc_corrections_sots.mat')
else
    load('poc_corrections_south.mat')
end

  % schallenbergMultiYearObservationsFluorescence2019
  c = -(1/.25) .* log(trans ./100);    % Calibrated for SAZ at SOTS using same Wetlabs CStar beam
                                       % transmissometer (25 cm path
                                       % length). Units of m^-1
  cp  = c - c_sw_ctd;                  % Contirbution from particulate equals total attenutation minus csw;  Units of m^-1
  bbp = 0.00345.*cp.^0.449;            % backscattering from cp;Units of m^-1
  POC_cs = 37601*bbp + 4.95;           % POC from Backscattering.  
                                       % Using schallenbergMultiYearObservationsFluorescence2019
                                       % relationship, calibrated to SOTS,
                                       % but with 2 transfomations
  POC_cs = POC_cs./12;                 % convert to uM (or umol/L) frpm mg/m3                       
                                       
  POC_wg = 33.5.*cp + 3.064;           % POC from Backscattering.  Units of mg/m3  
                                       % Using gardnerGlobalPOCSynthesis2003
                                       % relationship, calibrated for APFZ
                                       % but directly from cp. 
                                       
  POC_sol = cf_cp_poc_0.*cp;           % POC using SOLACe measurements
  % C_phyto = POC*??                   % Phytoplankton carbon estimated from bbp using Graff
                                       % et al. (2015) relationship, via
                                       % Haentjens et al. (2017) (from
                                       % Christina's paper)
                                      

% POC - Figure out regional SO regression and where to get beam attenuation
% coefficeint (c).
%poc    = 31.7.*trans + 0.785; % GLOBAL Regressing (No SO data) gardnerGlobalPOCSynthesis2003
                               % bishopTransmissometerMeasurementPOC1999
                               % schallenbergMultiYearObservationsFluorescence2019
sal    = mean([avgStruct.sensor.salinity.data;    avgStruct.sensor.salinity_2.data]);
PDen   = sw_pden(sal,temp,pres,0);
o2sol  = gsw_O2sol_SP_pt(sal,temp);       
o2sat  = 100*(o2-o2sol)./o2sol;


% query points
% xqs = linspace(time(1),time(end),100) ;% 100 timesteps but variable spacing
%  xqs = time(1):1/24/60/6:time(end);       % 10 second time steps 
xqs = 0:(.02.*dis_scl):dist_m(end);  
yqs = 1:150 ;

Time = xqs; 
Dist = xqs./dis_scl;
Pres = yqs;

% Include MLD
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
%load(FNAME)
%MLD =interp1(tri_time_loc_p,tri_mld_p,Dist); 
%MLDsm = movmean(MLD, 3*6,'omitnan');

% Salinty
% Depth resolved
var = sal;
[xq,yq] = meshgrid(xqs, yqs);
vq_sal = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(1,:,:) = vq_sal;
datanames{1} = 'Salinity ($PSU$)' ;


% Temperature 
% Depth resolved
var = temp;
[xq,yq] = meshgrid(xqs, yqs);
vq_temp = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(2,:,:) = vq_temp;
datanames{2} = 'Temperature ($^\circ C$)' ;
%   Time Series for Tmin --- not sure about this
% for t=1:length(xq);
%     prof = movmean(vq_temp(:,t),5,'omitnan');        % 5 Meter moving mean
%     if ~isnan(MLDsm(t))
%     %prof = prof(1:MLDsm(t));
%     [min_m, min_i] = min(prof);
%     data_ts(2,1,t) = squeeze(min_m); % Size of Tmin
%     data_ts(2,2,t) = squeeze(min_i); % Depth of Tmin
%     else
%     data_ts(2,:,t) = nan;
%     end  
% end
% data_ts(2,:,:) = movmean(squeeze(data_ts(2,:,:)),6*3,2,'omitnan');


%%%%% Potential density
var = PDen;
[xq,yq] = meshgrid(xqs, yqs);
vq_pden = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(3,:,:) = vq_pden;
datanames{3} = 'Potential Density ($kg/m3$)' ;


%%%%%% Chlorrophyl (mg/m3) 
%%%    Depth Resolved
var = chl;
[xq,yq] = meshgrid(xqs, yqs);
vq_chl = squeeze(griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq));
data(4,:,:) = vq_chl;
datanames{4} = {'Chlorophyll'; '($mg/m3$)'} ;
% Tim Series
for t=1:length(xqs);
    prof = movmean(vq_chl(:,t),5,'omitnan');        % 5 Meter moving mean 
    % vq_chl(:,t) = prof;                            % reasisng with smoothed
    dcm_thr = 1.0;%.*nanmean(prof(1:50));           % 30% higher than nan mean of top 50 m
    DCM = find(prof>dcm_thr) ;
    % Assign metrics if theres a DCM
    if~isempty(DCM) %&& DCM(1)>40;  % Assign metrics if theres a DCM
    [max_m, max_i] = max(prof);
    data_ts(4,1,t) = squeeze(max_m); % Size of DCM
    data_ts(4,2,t) = squeeze(max_i); % Depth of DCM
    data_ts(4,3,t) = DCM(1);         % Top of DCM
    data_ts(4,4,t) = DCM(end);       % Base of DCM
    data_ts(4,5,t) = length(DCM);    % Width of DCM (m)
    data_ts(4,6,t) = sum(prof(DCM)); % integrate of DCM (mg/m2)
    else
    data_ts(4,:,t) = nan;
    end
end
data_ts(4,:,:) = movmean(squeeze(data_ts(4,:,:)),6*3,2);

% Sample plot
% imagesc(Time,Pres,squeeze(vq_chl)); hold on;  colorbar; 
% plot(Time,MLDsm,'k','linewidth',2); hold on; 
% %plot(Time,squeeze(data_ts(2,2,:)),'r','linewidth',2); hold on; 
% plot(Time,squeeze(data_ts(4,2,:)),'k','linewidth',2); hold on
% plot(Time,squeeze(data_ts(4,3,:)),'k:','linewidth',2); hold on
% plot(Time,squeeze(data_ts(4,4,:)),'k:','linewidth',2); hold on

% Oxygen
% Depth Resolved
var = o2;
[xq,yq] = meshgrid(xqs, yqs);
vq_o2 = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(5,:,:) = vq_o2;
datanames{5} = 'Oxygen ($uM$)' ;

% Oxygen Saturation
% Depth Resolved
var = o2sat;
vq_o2sat = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(6,:,:) = vq_o2sat;
datanames{6} = 'Oxygen Saturation ($\%$)' ;

% Tranmisometer
% Depth Resolved
var = trans;
[xq,yq] = meshgrid(xqs, yqs);
vq_trans = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(7,:,:) = vq_trans;
datanames{7} = 'Transmisometer ($\%$)' ;
for t=1:length(xqs)
    prof = movmean(vq_trans(:,t),5,'omitnan');        % 5 Meter moving mean 
    %vq_trans(:,t) = prof;                            % reasisng with smoothed
    dcm_thr = 96.6;%.75.*nanmean(prof(1:50));           % 25% higher than nan mean of top 25 m
    DCM = find(prof<dcm_thr) ;
    % Assign metrics if theres a DCM
    if~isempty(DCM) % Assign metrics if theres a DCM
    [max_m, max_i] = min(prof);
    data_ts(7,1,t) = squeeze(max_m); % Size of DCM
    data_ts(7,2,t) = squeeze(max_i); % Depth of DCM
    data_ts(7,3,t) = DCM(1);         % Top of DCM
    data_ts(7,4,t) = DCM(end);       % Base of DCM
    data_ts(7,5,t) = length(DCM);    % Width of DCM (m)
    data_ts(7,6,t) = sum(prof(DCM)); % integrate of DCM (mg/m2)
    else
    data_ts(7,:,t) = nan;
    end
end
data_ts(7,:,:) = movmean(squeeze(data_ts(7,:,:)),6*3,2);

% Stratification
% Depth resolved
dens  = squeeze(data(3,:,:));
dp_dz = diff(dens) ; dp_dz = [ dp_dz(1,:);dp_dz];
p0    = 1026.7;  % mean of underwater surface density (ipden_uw)
N2    = (9.81./p0).* dp_dz;
data(8,:,:) = N2;   
datanames{8} = '$N^2$ ($s^{-2}$)'; % caxis([0,5e-5]);
% mahadevanEddyDrivenStratificationInitiates2012

% POC
var = POC_cs;
[xq,yq] = meshgrid(xqs, yqs);
vq_poc_cs = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(9,:,:) = vq_poc_cs;
datanames{9} = 'POC-CS (uM)' ;
for t=1:length(xqs)
    prof = movmean(vq_poc_cs(:,t),5,'omitnan');        % 5 Meter moving mean 
    %vq_trans(:,t) = prof;                            % reasisng with smoothed
    dcm_thr = 1.1.*nanmean(prof(1:50));           % 25% higher than nan mean of top 25 m
    DCM = find(prof>dcm_thr) ;
    % Assign metrics if theres a DCM
    if~isempty(DCM) % Assign metrics if theres a DCM
    [max_m, max_i] = max(prof);
    data_ts(9,1,t) = squeeze(max_m); % Size of DCM
    data_ts(9,2,t) = squeeze(max_i); % Depth of DCM
    data_ts(9,3,t) = DCM(1);         % Top of DCM
    data_ts(9,4,t) = DCM(end);       % Base of DCM
    data_ts(9,5,t) = length(DCM);    % Width of DCM (m)
    data_ts(9,6,t) = sum(prof(DCM)); % integrate of DCM (mg/m2)
    else
    data_ts(9,:,t) = nan;
    end
end
data_ts(9,:,:) = movmean(squeeze(data_ts(9,:,:)),6*3,2);

var = POC_wg;
[xq,yq] = meshgrid(xqs, yqs);
vq_poc_wg = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(10,:,:) = vq_poc_wg;
datanames{10} = 'POC-WG (uM)' ;
for t=1:length(xqs)
    prof = movmean(vq_poc_wg(:,t),5,'omitnan');        % 5 Meter moving mean 
    %vq_trans(:,t) = prof;                            % reasisng with smoothed
    dcm_thr = 1.25.*nanmean(prof(1:50));           % 25% higher than nan mean of top 25 m
    DCM = find(prof>dcm_thr) ;
    % Assign metrics if theres a DCM
    if~isempty(DCM) % Assign metrics if theres a DCM
    [max_m, max_i] = max(prof);
    data_ts(10,1,t) = squeeze(max_m); % Size of DCM
    data_ts(10,2,t) = squeeze(max_i); % Depth of DCM
    data_ts(10,3,t) = DCM(1);         % Top of DCM
    data_ts(10,4,t) = DCM(end);       % Base of DCM
    data_ts(10,5,t) = length(DCM);    % Width of DCM (m)
    data_ts(10,6,t) = sum(prof(DCM)); % integrate of DCM (mg/m2)
    else
    data_ts(10,:,t) = nan;
    end
end
data_ts(10,:,:) = movmean(squeeze(data_ts(10,:,:)),6*3,2);

var = POC_sol;
[xq,yq] = meshgrid(xqs, yqs);
vq_poc_sol = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(11,:,:) = vq_poc_sol;
datanames{11} = 'POC (uM)' ;
for t=1:length(xqs)
    prof = movmean(vq_poc_sol(:,t),5,'omitnan');        % 5 Meter moving mean 
    %vq_trans(:,t) = prof;                            % reasisng with smoothed
    dcm_thr = 1.25.*nanmean(prof(1:50));           % 25% higher than nan mean of top 25 m
    DCM = find(prof>dcm_thr) ;
    % Assign metrics if theres a DCM
    if~isempty(DCM) % Assign metrics if theres a DCM
    [max_m, max_i] = max(prof);
    data_ts(11,1,t) = squeeze(max_m); % Size of DCM
    data_ts(11,2,t) = squeeze(max_i); % Depth of DCM
    data_ts(11,3,t) = DCM(1);         % Top of DCM
    data_ts(11,4,t) = DCM(end);       % Base of DCM
    data_ts(11,5,t) = length(DCM);    % Width of DCM (m)
    data_ts(11,6,t) = sum(prof(DCM)); % integrate of DCM (mg/m2)
    else
    data_ts(11,:,t) = nan;
    end
end
data_ts(11,:,:) = movmean(squeeze(data_ts(11,:,:)),6*3,2);


% POC:Chl (SOL)
var = (POC_sol.*12)./chl;  % Convert POC to mg/m3
[xq,yq] = meshgrid(xqs, yqs);
vq_chl_c = griddata(squeeze(dist_m),squeeze(pres),squeeze(var),xq ,yq);
data(12,:,:) = vq_chl_c;
datanames{12} = 'C:Chl (g/g)' ;
for t=1:length(xqs)
    prof = movmean(vq_chl_c(:,t),5,'omitnan');       % 5 Meter moving mean 
    %vq_trans(:,t) = prof;                            % reasisng with smoothed
    dcm_thr = 1.25.*nanmean(prof(1:50));              % 25% higher than nan mean of top 25 m
    DCM = find(prof>dcm_thr) ;
    % Assign metrics if theres a DCM
    if~isempty(DCM) % Assign metrics if theres a DCM
    [max_m, max_i] = max(prof);
    data_ts(12,1,t) = squeeze(max_m); % Size of DCM
    data_ts(12,2,t) = squeeze(max_i); % Depth of DCM
    data_ts(12,3,t) = DCM(1);         % Top of DCM
    data_ts(12,4,t) = DCM(end);       % Base of DCM
    data_ts(12,5,t) = length(DCM);    % Width of DCM (m)
    data_ts(12,6,t) = sum(prof(DCM)); % integrate of DCM (mg/m2)
    else
    data_ts(12,:,t) = nan;
    end
end
data_ts(12,:,:) = movmean(squeeze(data_ts(12,:,:)),6*3,2);


% Time Series 

%  Save 
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
save(FNAME,'Time', 'Pres', 'datanames', 'data','data_ts', 'MLD','MLDsm','Dist','-append')

end
end

% FIRe Data
if false
for FIRe_Vars = 1 

% Fire Variable -  Which may not be correctly Saved (i.e. lined up)
% query points
% xqs = linspace(time(1),time(end),100) ;% 100 timesteps but variable spacing
xqs = fire_time(1):1/24/12:fire_time(end);       % 5 min time steps 
yqs = 1:2:200 ;

var = fire_fv_fm;
[xq,yq] = meshgrid(xqs, yqs);
vq_fv_fm = griddata(squeeze(fire_time),squeeze(fire_pres),squeeze(var),xq ,yq);
firedata(1,:,:) = vq_fv_fm;
firedata_raw(1,:,:) = fire_fv_fm;
firedatanames{1} = 'Fv\textbackslash Fm' ;

var = fire_Abs_sig;
[xq,yq] = meshgrid(xqs, yqs);
vq_abs_sig = griddata(squeeze(fire_time),squeeze(fire_pres),squeeze(var),xq ,yq);
firedata(2,:,:) = vq_abs_sig;
firedata_raw(2,:,:) = fire_Abs_sig;
firedatanames{2} = 'Absolute Sigma';

var = fire_Rel_sig;
[xq,yq] = meshgrid(xqs, yqs);
vq_rel_sig = griddata(squeeze(fire_time),squeeze(fire_pres),squeeze(var),xq ,yq);
firedata(3,:,:) = vq_rel_sig;
firedata_raw(3,:,:) = fire_Rel_sig;
firedatanames{3} = 'Relative Sigma';

fire_Time = xqs; 
fire_Pres = yqs;

%  Save 
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
save(FNAME,'fire_Time', 'fire_Pres', 'firedatanames', 'firedata','firedata_raw','fire_time','fire_pres', '-append')

end
end

end
end




%% CTD PLOTTER
% User Inputs
for CTDPlot = 1
% Variable options
% 1.  Salinity (PSU)      2. Temperature (C)          3. Potential Densisity (kg/m3) 
% 4.  Chlorophyll (mg/m3) 5. Oxygen (uM)              6. Oxygen Saturation (%)
% 7.  Transmisometer (%)  8. Stratification (1/s2)    9. POC (uM) form CS  
% 10. POC (uM) from WG    11. POC (uM) from Solace    12. Chl:C(from SOLACE)

% Pick 3 variables to plot from 3 legs.  
Tows = 7*[1,1,1]; % Presumably from the same tow, but not nessecarily
Legs = [1,2,3]; % Presumably sequential, but not nessecarily
Vars = [4,11,12]; % See options above

% Set directory to path where Triaxus Data is stored
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
 
% Create Figure
figure('Units','normalized','Position',[0 0 1.8 2])
clc
set(gcf,'color','w');


% Compile data to use consistent color bar across legs
clearvars alldata alldataV
for v = 1:3
   alldataV =[];
   for leg = 1:length(Legs)
   fname = ['Tow_',num2str(Tows(leg)), '_Leg_', num2str(Legs(leg)) ,'.mat'] ;   
   load(fname)
   vq =squeeze(data(Vars(v),:,:)); % Pull Variable from Data Matrix
   alldataV = [alldataV;vq(:)];
   end
   alldata(v,:) = alldataV;
end

% Loop through Legs
for leg = 1:length(Legs)
   fname = ['Tow_',num2str(Tows(leg)), '_Leg_', num2str(Legs(leg)) ,'.mat'] ;  
   tx = {['Triaxus Tow ', num2str(Tows(1))];datestr(Time(1),'mmm-DD') };
   % Load Main Data Set
   cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
   load(fname)
  
% Loop through Variables
for v = 1:3
vq =squeeze(data(Vars(v),:,:)); % Pull Variable from Data Matrix

% Plot
h=subplot(3,3,3*(v-1) + leg);

% Adjust Positions [left bottom width height]
posnew = get(h, 'Position'); 
posnew(1) = posnew(1) - 0.1; 
posnew(2) = posnew(2); 
posnew(3) = posnew(3)-.04;
posnew(4) = posnew(4)+.07; 

if leg==2; posnew(1)= posnew(1) - .1; end;
if leg==3; posnew(1)= posnew(1) - .2; end
% Bottom
if v==3; posnew(2)= posnew(2) - .06; end
if v==2; posnew(2)= posnew(2) - .04; end 
if v==1; posnew(2)= posnew(2) - .02; end

set(h, 'Position', posnew)
imagesc(Dist,Pres,vq); hold on
%plot(Time,MLDsm,'r','linewidth',2); hold on
plot(squeeze(tri_dis_p),squeeze(tri_mld_p),'r','linewidth',2); hold on


 if any(Vars(v)==[4,11,12])%,7,9,10,11,12])
 %plot(Time,squeeze(data_ts(Vars(v),2,:)),'k','linewidth',1); hold on; % Depth Max
 %plot(Time,squeeze(data_ts(Vars(v),3,:)),'k','linewidth',1); hold on; % Depth Max
 %plot(Time,squeeze(data_ts(Vars(v),4,:)),'k','linewidth',1); hold on; % Depth Max
     plot(squeeze(tri_dis_p),squeeze(data_dis_ts(Vars(v),2,:)),'k','linewidth',.5); hold on
     plot(squeeze(tri_dis_p),squeeze(data_dis_ts(Vars(v),3,:)),'k','linewidth',.5); hold on
     plot(squeeze(tri_dis_p),squeeze(data_dis_ts(Vars(v),4,:)),'k','linewidth',.5); hold on
     %MLD =interp1(tri_time_loc_p,tri_mld_p,Time); 
     %MLDsm = movmean(MLD, 3*6,'omitnan');
 end

% Pick Color map
if     Vars(v) == 1; cmap = [ 1 1 1; cmocean('matter')];
elseif Vars(v) == 2; cmap = [ 1 1 1; cmocean('thermal')]; 
elseif Vars(v) == 3; cmap = [ 1 1 1; cmocean('dense')]; 
elseif Vars(v) == 4; cmap = [ 1 1 1; cmocean('haline')]; 
elseif Vars(v) == 5; cmap = [1 1 1; brewermap(256,'ylorrd')]; 
elseif Vars(v) == 6; cmap = [ 1 1 1; (cmocean('balance'))]; %caxis([89,111]); 
elseif Vars(v) == 7; cmap = [ 1 1 1; cmocean('thermal')]; 
elseif Vars(v) == 8; cmap = cmocean('gray'); cmap =flipud(cmap);  
elseif Vars(v) == 9; cmap = [ 1 1 1; cmocean('thermal')];
elseif Vars(v) == 10; cmap = [ 1 1 1; cmocean('thermal')]; 
elseif Vars(v) == 11; cmap = [ 1 1 1; cmocean('thermal')]; 
elseif Vars(v) == 12; cmap = [ 1 1 1; flipud(cmocean('thermal'))]; 
end

% Colormap And Colorbar
colormap(gca,cmap)
if Vars(v)==6; rdum = prctile(squeeze(alldata(v,:)),[0,100]); range = [-max(abs(rdum)), max(abs(rdum))];
% elseif Vars(v)==8; range = [0,5e-5];
% elseif Vars(v)==10 ; range = [0,10];
 elseif Vars(v)==4 ; range = [0,2.4];      % South 1 [0,2] South 2[0,2.]
% elseif Vars(v)==9 ; range = [0,7.5];
 elseif Vars(v)==11 ; range = [0,9];        % South 1 [0,7] South 2[0.9]
 elseif Vars(v)==12 ; range = [30,100];    % South 1 [30,150] South [40,150]
else;    range = prctile(squeeze(alldata(v,:)),[0,100]);
end
caxis(range)
%cbr = colorbar; set(cbr, 'YTick', linspace(range(1), range(2), 5),'TickLabelInterpreter','latex');

% Axis
if leg ==2 || leg==3; set(gca,'ytick',[]); end
set(gca,'ydir','rev')
%datetick(gca,'x','HH:MM','keepticks') 
set(gca,'fontsize',20,'TickLabelInterpreter','latex')

% Title
if v ==1; title(['Leg ',num2str(leg)],'interpreter','latex','fontsize',20); end
%if leg==1; ylabel(datanames{Vars(v)},'interpreter','latex','fontsize',20); end

% Colorbar
if leg==3 || leg==2 && Tows(leg) ==6
cloc = get(gca,'Position');
h= colorbar(gca,'position', [cloc(1)+cloc(3)+.02, cloc(2), .015, cloc(4) ], ...
    'orientation','vertical','TickLabelInterpreter','latex','fontsize',25);
%set(get(h,'title'),'string', datanames{Vars(v)},'interpreter','latex','fontsize',10)
end

if leg==1 && v==1
% tx = {['Triaxus Tow ', num2str(Tows(1))];datestr(Time(1),'mmm-DD') };
%  text(.1,.1 , tx,'position', 'HorizontalAl','center',...
%      'interpreter','latex','fontsize',40,'Rotation',90)
end


end
end



% Title
% ax(2) = axes; 
% set(ax(2), 'Position',[.00   .00  .09  1]);
% set(gca,'box','off','ycolor','w')
% ht = text(.5,.53,tx,'Rotation',90,'interpreter','latex','FontSize',40,  'HorizontalAlignment', 'center');

end

%% FIRe PLOTTER Interpolated
% User Inputs
for Fire_Plot = 1

% Variable options
% 1. Fm/Fv               2. Absolute Sigma II        3. Relative Sigma II

% Pick 3 variables to plot from 3 legs.  
Tows = [4,4,4]; % Presumably from the same tow, but not nessecarily
Legs = [1,2,3]; % Presumably sequential, but not nessecarily
Vars = [1,2,3]; % See options above

% Set directory to path where Triaxus Data is stored
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
 
% Create Figure
figure('Units','normalized','Position',[0 0 2 2])
clc
set(gcf,'color','w');


% Compile data to use consistent color bar across legs
clearvars alldata alldataV
for v = 1:3
   alldataV =[];
   for leg = 1:3
   fname = ['Tow_',num2str(Tows(leg)), '_Leg_', num2str(Legs(leg)) ,'.mat'] ;   
   load(fname)
   vq =squeeze(firedata(Vars(v),:,:)); % Pull Variable from Data Matrix
   alldataV = [alldataV;vq(:)];
   end
   alldata(v,:) = alldataV;
end

% Loop through Legs
for leg = 1:3
   fname = ['Tow_',num2str(Tows(leg)), '_Leg_', num2str(Legs(leg)) ,'.mat'] ;  
   tx = {['Triaxus Tow ', num2str(Tows(1))];datestr(Time(1),'mmm-DD') };
   % Load Main Data Set
   cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
   load(fname)
  
% Loop through Variables
for v = 1:3
vq =squeeze(firedata(Vars(v),:,:)); % Pull Variable from Data Matrix

% Plot
h=subplot(3,3,3*(leg-1) + v);
posnew = get(h, 'Position'); posnew(2) = posnew(2) + 0.005; posnew(4) = posnew(4) + 0.02; 
set(h, 'Position', posnew)
plot(Time,MLD,'k','linewidth',2); hold on
imagesc(fire_Time,fire_Pres,vq); hold on
plot(Time,MLD,'k','linewidth',2); hold on


% Pick Color map
if     Vars(v) == 1; cmap = [ 1 1 1; cmocean('haline')];
elseif Vars(v) == 2; cmap = [ 1 1 1; cmocean('haline')]; 
elseif Vars(v) == 3; cmap = [ 1 1 1; cmocean('haline')]; 
% elseif Vars(v) == 4; cmap = [ 1 1 1; cmocean('haline')]; 
% elseif Vars(v) == 5; cmap = [1 1 1; brewermap(256,'ylorrd')]; 
% elseif Vars(v) == 6; cmap = [ 1 1 1; (cmocean('balance'))]; %caxis([89,111]); 
% elseif Vars(v) == 7; cmap = [ 1 1 1; cmocean('thermal')]; 
% elseif Vars(v) == 8; cmap = [ 1 1 1; cmocean('thermal')]; 
end

% Colormap And Colorbar
colormap(gca,cmap)
range = prctile(squeeze(alldata(v,:)),[5,95]);
caxis(range)
%cbr = colorbar; set(cbr, 'YTick', linspace(range(1), range(2), 5),'TickLabelInterpreter','latex');

% Axis
set(gca,'ydir','rev')
datetick(gca,'x','HH:MM','keepticks') 
set(gca,'fontsize',12,'TickLabelInterpreter','latex')

% Title
if leg ==1; title(firedatanames{Vars(v)},'interpreter','latex','fontsize',20); end
if v==1; ylabel(['Leg ',num2str(leg)],'interpreter','latex','fontsize',20); end

% Colorbar
if leg==3
cloc = get(gca,'Position');
h= colorbar(gca,'position', [cloc(1), cloc(2)-.085, cloc(3), .04], ...
    'orientation','horizontal','TickLabelInterpreter','latex');
%set(get(h,'title'),'string', {'[Chl.]';'($mg/m^3$)'},'interpreter','latex','fontsize',25)
end

if leg==1 && v==1
tx = {['Triaxus Tow ', num2str(Tows(1))];datestr(Time(1),'mmm-DD') };
text(.1,.1 , tx,'position', 'HorizontalAl','center',...
    'interpreter','latex','fontsize',40,'Rotation',90)
end

end
end

% Title
ax(2) = axes; 
set(ax(2), 'Position',[.00   .00  .1  1]);
set(gca,'box','off','ycolor','w')
ht = text(.5,.53,tx,'Rotation',90,'interpreter','latex','FontSize',50,  'HorizontalAlignment', 'center');

end



%% Regression plot for all 
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus

   met_names{1} = 'Max ';
   met_names{2} = 'Depth of Max ';
   met_names{3} = 'Top of ';
   met_names{4} = 'Base of ';
   met_names{5} = 'Thickness of ';
   met_names{6} = 'Inventory of ';
   met_names{7} = 'Deep Enhancement of ';
   met_names{8} = 'Deep Enhancement of ';
   
   
% var_names{4}  = 'Deep Chlorrophyl Maxmima ';
%    var_names{11} = 'Deep Biomass Maxmima ';
%    var_names{12} = 'C:Chl Minimum';
 
   var_names{4}  = 'Chlorrophyl ';
   var_names{11} = 'Max Biomass ';
   var_names{12} = 'C:Chl ';
   var_names{13} = 'MLD';
   
   for units=1
   met_units{4,1} = '$(mg/m^3)$'; met_units{11,1} = '$(uM)$';        met_units{12,1} = '$(-)$'; 
   met_units{4,2} = '$(m)$';      met_units{11,2} = '$(m)$';         met_units{12,2} = '$(m)$';   met_units{13,2} = '$(m)$'; 
   met_units{4,3} = '$(m)$';      met_units{11,3} = '$(m)$';         met_units{12,3} = '$(m)$'; 
   met_units{4,4} = '$(m)$';      met_units{11,4} = '$(m)$';         met_units{12,4} = '$(m)$';  
   met_units{4,5} = '$(m)$';      met_units{11,5} = '$(m)$';         met_units{12,5} = '$(m)$'; 
   met_units{4,6} = '$(mg/m^2)$'; met_units{11,6} = '$(umol/m^2)$';  met_units{12,6} = '$(g:g)$'; 
   met_units{4,7} = '$(\%)$';     met_units{11,7} = '$(\%)$';        met_units{12,7} = '$(\%)$'; 
   met_units{4,8} = '$(\%)$';     met_units{11,8} = '$(\%)$';        met_units{12,8} = '$(\%)$'; 
   end
   
   
   
% vars1 = [11, 11,  11, 12, 12, 12 ];
% mets1 = [2,  1,  8,   2,   1,  8];   
%    
% vars2 = [4, 4,  4, 4, 4, 4];
% mets2 = [2, 1,  8, 2, 1, 8 ];

% vars1 = [13, 13,  13, 13, 13, 13 ];
% mets1 = [2,  2,  2,   2,   2,  2];   
%    
% vars2 = [12, 12,  12, 12, 12,12];
% mets2 = [2, 1,  7, 5, 6, 7 ];

% DCM Fig 
% vars1 = [13, 11,  12, 13,  11, 12 ];
% mets1 = [2,  2,  2,   2,   8,  8];   
%    
% vars2 = [4, 4,  4, 4, 4, 4];
% mets2 = [2, 2,  2, 8, 8, 8 ];

vars1 = [11, 11,  11, 13,  11, 12 ];
mets1 = [2,  5,  8,   2,   8,  8];   
   
vars2 = [4, 4,  4, 4, 4, 4];
mets2 = [2, 5,  8, 8, 8, 8 ];

vars1 = [11, 11,  11, 11,  11, 11 ];
mets1 = [2,  5,  7,   1,   7,  7];   
   
vars2 = [4, 4,  4, 4, 4, 4];
mets2 = [2, 5,  8, 1, 8, 8 ];


%cmapSouth1 = brewermap(9, 'YlOrRd');
%cmapSouth2 = brewermap(9, 'PuBU');
msize =100;

% Make Figure
figure('Units','normalized','Position',[0 0 2 2])
clc
set(gcf,'color','w');

   for i=1%1:6
   h= subplot(2,3,i)
   
        pos = get(h, 'Position') ;
   if i ==2 || i==5
        posnew = pos; posnew(1) = posnew(1) + 0.015; set(h, 'Position', posnew)
   end
    if i ==3 || i==6
        posnew = pos; posnew(1) = posnew(1) + 0.02; set(h, 'Position', posnew)
   end

   % Both
   ALL_var1 =[]; ALL_var2 = [];
   
   % South 1
   ci=0; 
   for tow =5%:5
       var1=[];   var2=[];
     if tow ==6; lastleg =2; else lastleg = 3;end
      for leg = 1:lastleg 
       ci = ci+1;
        % Load and Assign
        fname = ['Tow_',num2str(tow), '_Leg_', num2str(leg) ,'.mat'] ;   
        load(fname)
        var1 = [var1; squeeze(data_dis_ts(vars1(i),mets1(i),:))]; 
        var2 = [var2; squeeze(data_dis_ts(vars2(i),mets2(i),:))]; 
        
        if vars1(i)==11 && mets1(i) == 5 % remove weird thickness anomalies
            var2(var1>40) =[];  var1(var1>40) =[];
        end
        
      end
        % Regression 
        mdl = fitlm(var1,var2); 
        coefs_reg(:,tow,i) = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
        xlims(:,tow-3) =  [min(var1(~isnan(var2))),max(var1(~isnan(var2)))]';
        ylims(:,tow-3) =  [min(var2),max(var2)]';

        % Scatter Plot 
        if tow==4
        scatter(var1,var2,msize, cmapSouth1(4,:), 'filled'); hold on
        else
        scatter(var1,var2,msize, cmapSouth1(8,:), 'filled'); hold on
        end
        
   end
   
   % South 2
%    ci=4;
%    for tow =6:7
%        var1=[];   var2=[];
%      if tow ==6; lastleg =2; else lastleg = 3;end
%       for leg = 1:lastleg 
%        ci = ci+1;
%         % Load and Assign
%         fname = ['Tow_',num2str(tow), '_Leg_', num2str(leg) ,'.mat'] ;   
%         load(fname)
%         var1 = [var1; squeeze(data_dis_ts(vars1(i),mets1(i),:))]; 
%         var2 = [var2; squeeze(data_dis_ts(vars2(i),mets2(i),:))]; 
%         if tow==7 && mets1(i) == 5 ; % remove weird thickness anomalies
%             var2(var1>18) =[];  var1(var1>18) =[];
%             var1(var2>18) =[];  var2(var2>18) =[];
%         end
%       end
%       
%         mdl = fitlm(var1,var2); 
%         coefs_reg(:,tow,i) = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
%         xlims(:,tow-3) =  [min(var1(~isnan(var2))),max(var1(~isnan(var2)))]';
%         ylims(:,tow-3) =  [min(var2),max(var2)]';
%        
%         % Scatter Plot
%         if tow ==6
%         scatter(var1,var2,msize, cmapSouth2(4,:), 'filled'); hold on
%         else
%         scatter(var1,var2,msize, cmapSouth2(8,:), 'filled'); hold on
%         end
%         set(gca,'fontsize',25,'TickLabelInterpreter','latex')
%    end
% 
%      % Plot regression
%      h(1) = plot(xlims(:,1), coefs_reg(1,4,i)+coefs_reg(2,4,i).*xlims(:,1),'color',cmapSouth1(4,:),'linewidth',3); hold on
%      h(2) = plot(xlims(:,2), coefs_reg(1,5,i)+coefs_reg(2,5,i).*xlims(:,2),'color',cmapSouth1(8,:),'linewidth',3); hold on
%      h(3) = plot(xlims(:,3), coefs_reg(1,6,i)+coefs_reg(2,6,i).*xlims(:,3),'color',cmapSouth2(4,:),'linewidth',3); hold on
%      h(4) = plot(xlims(:,4), coefs_reg(1,7,i)+coefs_reg(2,7,i).*xlims(:,4),'color',cmapSouth2(8,:),'linewidth',3); hold on
%      %if mets1(i)== 2 && mets2(i)== 2; plot(1:100,1:100,':k', 'linewidth',3); end 
%      if any(i==[1:3]); plot(1:100,1:100,':k', 'linewidth',3); end          
%      
%      ll= legend(h(1:4),[ '$r^2$: ',num2str(round(100*coefs_reg(3,4,i))/100),'  $\;$' ],...
%                    [ '$r^2$: ',num2str(round(100*coefs_reg(3,5,i))/100),'  $\;$' ],...
%                    [ '$r^2$: ',num2str(round(100*coefs_reg(3,6,i))/100),'  $\;$' ],...
%                    [ '$r^2$: ',num2str(round(100*coefs_reg(3,7,i))/100),'  $\;$' ],'fontsize',23,'interpreter','latex','location','northoutside','orientation','horizontal','box','off')
%        
        %pos = get(ll, 'Position') ;
        %posnew = pos; posnew(1) = posnew(1) + 0.03; set(ll, 'Position', posnew) 
        
%         legend(h(1:4),[ 'Dec 29: ',num2str(round(100*coefs_reg(3,4,i))/100) ],...
%                    [ 'Dec 30: ',num2str(round(100*coefs_reg(3,5,i))/100)],...
%                    [ 'Jan 02: ',num2str(round(100*coefs_reg(3,6,i))/100)],...
%                    [ 'Jan 11: ',num2str(round(100*coefs_reg(3,7,i))/100)],'fontsize',20,'interpreter','latex','location','northoutside','orientation','horizontal','box','off')
%      
%       
   % Labels and Lim 
     %xlabel([met_names{mets1(i)}, var_names{vars1(i)}, met_units{vars1(i), mets1(i)} ],'interpreter','latex','fontsize',35)
     %ylabel([met_names{mets2(i)}, var_names{vars2(i)}, met_units{vars2(i),mets2(i)} ],'interpreter','latex','fontsize',35)
     xlim([min(xlims(:)), max(xlims(:))])
     ylim([min(ylims(:)), max(ylims(:))])
     if vars2(i)==12 && mets1(i)~= 2; set(gca,'ydir','reverse'); end
     if vars1(i)==12 && mets1(i)~= 2; set(gca,'xdir','reverse'); end
   
   end
   
%% Regression plots
%     data_ts(7,1,t) = squeeze(max_m); % Size of DCM
%     data_ts(7,2,t) = squeeze(max_i); % Depth of DCM
%     data_ts(7,3,t) = DCM(1);         % Top of DCM
%     data_ts(7,4,t) = DCM(end);       % Base of DCM
%     data_ts(7,5,t) = length(DCM);    % Width of DCM (m)
%     data_ts(7,6,t) = sum(prof(DCM)); % integrate of DCM (mg/m2)
%     


figure('Units','normalized','Position',[0 0 2 2])
clc
set(gcf,'color','w');

cmap = brewermap([11], 'RdYlBu');
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus

 alldata_ts = []; allmld_ts = [];
 msize =100;
 % Version 1:
 i=0;
 for tow =4:7
    if tow ==6; lastleg =2; else lastleg = 3;end
   for leg = 1:lastleg 
       i = i+1;
   fname = ['Tow_',num2str(tow), '_Leg_', num2str(leg) ,'.mat'] ;   
   load(fname)
   
   subplot(2,3,1)
   scatter(squeeze(tri_mld_p),squeeze(data_dis_ts(4,1,:)),msize, cmap(i,:), 'filled'); hold on
   xlabel('MLD ($m$)','interpreter','latex','fontsize',20)
   ylabel('DCM ($mg/m^3$)','interpreter','latex','fontsize',20)
   
   subplot(2,3,2)
   scatter(squeeze(tri_mld_p),squeeze(data_dis_ts(4,6,:)),msize, cmap(i,:), 'filled'); hold on
   xlabel('MLD ($m$)','interpreter','latex','fontsize',20)
   ylabel('Integrated DCM ($mg/m^2$)','interpreter','latex','fontsize',20)
   
   subplot(2,3,3)
   scatter(squeeze(tri_mld_p),squeeze(data_dis_ts(4,2,:)),msize, cmap(i,:), 'filled'); hold on
   xlabel('MLD ($m$)','interpreter','latex','fontsize',20)
   ylabel('DCM Depth ($m$)','interpreter','latex','fontsize',20)
   
%    subplot(2,3,4)
%    scatter(squeeze(data_ts(4,1,:)), squeeze(data_ts(4,6,:)),10, cmap(i,:), 'filled'); hold on
%    xlabel('DCM ($mg/m^3$)','interpreter','latex','fontsize',20)
%    ylabel('Integrated DCM ($mg/m^2$)','interpreter','latex','fontsize',20)
%    
   subplot(2,3,4)
   scatter(squeeze(tri_mld_p),squeeze(data_dis_ts(12,1,:)),msize, cmap(i,:), 'filled'); hold on
   xlabel('MLD ($m$)','interpreter','latex','fontsize',20)
   ylabel('Chl:C max ','interpreter','latex','fontsize',20)
   
   
   subplot(2,3,5)
   scatter(squeeze(data_dis_ts(4,1,:)), squeeze(data_dis_ts(11,1,:)),msize, cmap(i,:), 'filled'); hold on
   xlabel('DCM ($mg/m^3$)','interpreter','latex','fontsize',20)
   ylabel('DBM Magnitude ($um/m^3$)','interpreter','latex','fontsize',20)
   %ylim([3.3,5])
  
   subplot(2,3,6)
   scatter(squeeze(data_dis_ts(4,2,:)), squeeze(data_dis_ts(11,2,:)),msize, cmap(i,:), 'filled'); hold on
   xlabel('DCM Depth ($m$)','interpreter','latex','fontsize',20)
   ylabel('$\sim$ DBM Depth ($m$)','interpreter','latex','fontsize',20)
   %xlim([60,100]);   ylim([60,100])
   
   alldata_ts = cat(3,alldata_ts,data_dis_ts);
   allmld_ts = cat(2,allmld_ts,tri_mld_p);
   end
 end

% add regression
 subplot(2,3,1)
 mdl = fitlm(squeeze(allmld_ts),squeeze(alldata_ts(4,1,:))); 
 coefs_reg = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
 xs = [min(squeeze(allmld_ts)), max(squeeze(allmld_ts))];
 plot(xs, coefs_reg(1)+coefs_reg(2).*xs,'k','linewidth',3)
 title(['$r^2$ = ',num2str(coefs_reg(3))],'fontsize',20,'interpreter','latex')
 
 subplot(2,3,2)
 mdl = fitlm(squeeze(allmld_ts),squeeze(alldata_ts(4,6,:))); 
 coefs_reg = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
 xs = [min(squeeze(allmld_ts)), max(squeeze(allmld_ts))];
 plot(xs, coefs_reg(1)+coefs_reg(2).*xs,'k','linewidth',3)
 title(['$r^2$ = ',num2str(coefs_reg(3))],'fontsize',20,'interpreter','latex')
 
 subplot(2,3,3)
 mdl = fitlm(squeeze(allmld_ts),squeeze(alldata_ts(4,2,:))); 
 coefs_reg = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
 xs = [min(squeeze(allmld_ts)), max(squeeze(allmld_ts))];
 plot(xs, coefs_reg(1)+coefs_reg(2).*xs,'k','linewidth',3)
 title(['$r^2$ = ',num2str(coefs_reg(3))],'fontsize',20,'interpreter','latex')
 
%  subplot(2,3,4)
%  mdl = fitlm(squeeze(alldata_ts(4,1,:)),squeeze(alldata_ts(4,6,:))); 
%  coefs1 = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
%  xs = [min(squeeze(alldata_ts(4,1,:))), max(squeeze(alldata_ts(4,1,:)))];
%  plot(xs, coefs1(1)+coefs1(2).*xs,'k','linewidth',3)
%  title(['$r^2$ = ',num2str(coefs1(3))],'fontsize',20)
 
 subplot(2,3,4)
 mdl = fitlm(squeeze(allmld_ts),squeeze(alldata_ts(12,1,:))); 
 coefs_reg = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
 xs = [min(squeeze(allmld_ts)), max(squeeze(allmld_ts))];
 plot(xs, coefs_reg(1)+coefs_reg(2).*xs,'k','linewidth',3)
 title(['$r^2$ = ',num2str(coefs_reg(3))],'fontsize',20,'interpreter','latex')

 subplot(2,3,5)
 mdl = fitlm(squeeze(alldata_ts(4,1,:)),squeeze(alldata_ts(11,1,:))); 
 coefs_reg = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
 xs = [min(squeeze(alldata_ts(4,1,:))), max(squeeze(alldata_ts(4,1,:)))];
 plot(xs, coefs_reg(1)+coefs_reg(2).*xs,'k','linewidth',3)
 %ylim([4.8,5.6])
 title(['$r^2$ = ',num2str(coefs_reg(3))],'fontsize',20,'interpreter','latex')

 subplot(2,3,6)
 mdl = fitlm(squeeze(alldata_ts(4,2,:)),squeeze(alldata_ts(11,2,:))); 
 coefs_reg = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
 xs = [min(squeeze(alldata_ts(4,2,:))), max(squeeze(alldata_ts(4,2,:)))];
 plot(xs, coefs_reg(1)+coefs_reg(2).*xs,'k','linewidth',3)
 title(['$r^2$ = ',num2str(coefs_reg(3))],'fontsize',20,'interpreter','latex')
 
 

%% FIRe PLOTTER Scatter
% User Inputs
for Fire_Plot = 1

% Variable options
% 1. Fm/Fv               2. Absolute Sigma II        3. Relative Sigma II

% Pick 3 variables to plot from 3 legs.  
Tows = [5,5,5]; % Presumably from the same tow, but not nessecarily
Legs = [1,2,3]; % Presumably sequential, but not nessecarily
Vars = [1,2,3]; % See options above

% Set directory to path where Triaxus Data is stored
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
 
% Create Figure
figure('Units','normalized','Position',[0 0 2 2])
clc
set(gcf,'color','w');


% Compile data to use consistent color bar across legs
clearvars alldata alldataV
for v = 1:3
   alldataV =[];
   for leg = 1:3
   fname = ['Tow_',num2str(Tows(leg)), '_Leg_', num2str(Legs(leg)) ,'.mat'] ;   
   load(fname)
   vq =squeeze(firedata(Vars(v),:,:)); % Pull Variable from Data Matrix
   alldataV = [alldataV;vq(:)];
   end
   alldata(v,:) = alldataV;
end

% Loop through Legs
for leg = 1:3
   fname = ['Tow_',num2str(Tows(leg)), '_Leg_', num2str(Legs(leg)) ,'.mat'] ;  
   tx = {['Triaxus Tow ', num2str(Tows(1))];datestr(Time(1),'mmm-DD') };
   % Load Main Data Set
   cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
   load(fname)
  
% Loop through Variables
for v = 1:3
vq =squeeze(firedata(Vars(v),:,:)); % Pull Variable from Data Matrix

% Plot
h=subplot(3,3,3*(leg-1) + v);
posnew = get(h, 'Position'); posnew(2) = posnew(2) + 0.005; posnew(4) = posnew(4) + 0.02; 
set(h, 'Position', posnew)
scatter(fire_time,fire_pres,100,squeeze(firedata_raw(v,:)),'filled');  hold on
plot(Time,MLD,'k','linewidth',4); hold on


% Pick Color map
if     Vars(v) == 1; cmap = [ 1 1 1; cmocean('haline')];
elseif Vars(v) == 2; cmap = [ 1 1 1; cmocean('haline')]; 
elseif Vars(v) == 3; cmap = [ 1 1 1; cmocean('haline')]; 
% elseif Vars(v) == 4; cmap = [ 1 1 1; cmocean('haline')]; 
% elseif Vars(v) == 5; cmap = [1 1 1; brewermap(256,'ylorrd')]; 
% elseif Vars(v) == 6; cmap = [ 1 1 1; (cmocean('balance'))]; %caxis([89,111]); 
% elseif Vars(v) == 7; cmap = [ 1 1 1; cmocean('thermal')]; 
% elseif Vars(v) == 8; cmap = [ 1 1 1; cmocean('thermal')]; 
end

% Colormap And Colorbar
colormap(gca,cmap)
range = prctile(squeeze(alldata(v,:)),[5,95]);
caxis(range)
%cbr = colorbar; set(cbr, 'YTick', linspace(range(1), range(2), 5),'TickLabelInterpreter','latex');

% Axis
set(gca,'ydir','rev')
datetick(gca,'x','HH:MM','keepticks') 
set(gca,'fontsize',12,'TickLabelInterpreter','latex')

% Title
if leg ==1; title(firedatanames{Vars(v)},'interpreter','latex','fontsize',20); end
if v==1; ylabel(['Leg ',num2str(leg)],'interpreter','latex','fontsize',20); end

% Colorbar
if leg==3
cloc = get(gca,'Position');
h= colorbar(gca,'position', [cloc(1), cloc(2)-.085, cloc(3), .04], ...
    'orientation','horizontal','TickLabelInterpreter','latex');
%set(get(h,'title'),'string', {'[Chl.]';'($mg/m^3$)'},'interpreter','latex','fontsize',25)
end

if leg==1 && v==1
tx = {['Triaxus Tow ', num2str(Tows(1))];datestr(Time(1),'mmm-DD') };
text(.1,.1 , tx,'position', 'HorizontalAl','center',...
    'interpreter','latex','fontsize',40,'Rotation',90)
end

end
end

% Title
ax(2) = axes; 
set(ax(2), 'Position',[.00   .00  .1  1]);
set(gca,'box','off','ycolor','w')
ht = text(.5,.53,tx,'Rotation',90,'interpreter','latex','FontSize',50,  'HorizontalAlignment', 'center');

end



