%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% "Triaxus" CTD  Data %%
%%%%%%%%%%%%%%%%%%%%%%
% Compile CTD Profiles from sites 1 and 2 for analysis with Triaxus Data

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


%% Determine indexes of CTD at each site:
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/CTD
load('CTD_MLD.mat')

cd /Users/roh029/Desktop/Data.nosync/SOLACE/CHL
load CHL cf_ctd CHL_ctd    % NPQ corrected profiles and correction factor to discrete flouresense
% This was made in CHL_Builder_plotter and regresses are discrete chl
% smaples against the sensor value if the sensor was inthe dark (ie below
% euphotic zone or at night) 
% EX. corrected values = cf_ctd(1) + cf_ctd(2).*CTD_Flourometer;

cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus


South1_CTDi = 41:63; 
South2_CTDi = 64:85;

% Triaxus tows
clearvars tri_i
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
load('Tow_1_Leg_1.mat', 'Time')
[dum,tri_i(1)]= min(abs(casttime_loc-Time(1)));
load('Tow_2_Leg_1.mat', 'Time')
[dum,tri_i(2)]= min(abs(casttime_loc-Time(1)));
load('Tow_3_Leg_1.mat', 'Time')
[dum,tri_i(3)]= min(abs(casttime_loc-Time(1)));
load('Tow_4_Leg_1.mat', 'Time')
[dum,tri_i(4)]= min(abs(casttime_loc-Time(1)));
load('Tow_5_Leg_1.mat', 'Time')
[dum,tri_i(5)]= min(abs(casttime_loc-Time(1)));
load('Tow_6_Leg_1.mat', 'Time')
[dum,tri_i(6)]= min(abs(casttime_loc-Time(1)));
load('Tow_7_Leg_1.mat', 'Time')
[dum,tri_i(7)]= min(abs(casttime_loc-Time(1)));


    
       
%% POC Correction Run Seperately for seperate sites!
for POC_conversion=1
% 1. Calibarate the two Transmisometer sensors by subtracting out 'c_sw' as
% measure at 186 m, which is the lowest the triaxus goes, 

% 2. Calibrate Site Specific CTD C_P against measure POC
Sots_CTD  = 3:39;  Sots_tri   = 4:19;   Sots_POC = 1:32;
South_CTD = 41:84; South_tri   = 20:30; South_POC = 33:68;
 
 fname = 'poc_corrections_sots' ;
 ctd_i = Sots_CTD; tri_i =Sots_tri ;  Station_i =  Sots_POC; cp_correct_range=100:150;

% fname = 'poc_corrections_south' ;
% ctd_i = South_CTD; tri_i =South_tri ;  Station_i =  South_POC;  cp_correct_range=100:120;

% Profile Data
all_trans_ctd = []; all_pres_ctd = [];
for i = ctd_i 
    all_trans_ctd = [all_trans_ctd,Profiles(i).trans]; 
    all_pres_ctd = [all_pres_ctd,Profiles(i).pres]; 
end

% Compute Beam Attenuation
all_c_ctd     = -(1/.25) .* log(all_trans_ctd ./100);
c_sw_ctd      = min(all_c_ctd);    % Minimum value recorded above lowest depth triaxus went
                                           % If Tri Trans gets fixed (so it
                                           % doesnt exceed 100 then use
                                           % all_c_ctd here)
cp_ctd    = all_c_ctd - c_sw_ctd;

% Compute POC
bbp = 0.00345.*cp_ctd.^0.449;            
POC_ctd_cs = 37601*bbp + 4.95;
POC_ctd_cs =POC_ctd_cs./12; % convert form gram to mol
POC_ctd_wg = 33.5.*cp_ctd + 3.064;

% 1. Find and calibrate c_p from south 1 and 2 for ctd and triaxus
TRIpath = '/Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus/Triaxus/processing/in2020_v08/cap/cappro/avg';
cd(TRIpath)
all_trans_tri = []; all_pres_tri = [];
File = dir('*g.mat'); file={File(:).name};
for i =tri_i; load(File(i).name); 
    all_trans_tri = [all_trans_tri,avgStruct.sensor.transmissometer.data]; 
    all_pres_tri = [all_pres_tri,avgStruct.sensor.pressure.data];
end

% Find Correction factor for Triaxus and CTD
clearvars trans_ctd_mean trans_tri_mean CTD_cp
  for i=cp_correct_range
      trans_tri_mean(i-99) = nanmean(all_trans_tri(find(all_pres_tri<i & all_pres_tri>=i-1)));
      trans_ctd_mean(i-99) = nanmean(all_trans_ctd(find(all_pres_ctd<i & all_pres_ctd>=i-1)));
  end
  
  mdl = fitlm(squeeze(trans_tri_mean),squeeze(trans_ctd_mean)); 
  coefs1 = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted] ;
  trans_offset =  mean(trans_tri_mean) - mean(trans_ctd_mean);
  if true
  for samplefig=1
      figure('Units','normalized','Position',[0 0 2 2])
      set(gcf,'color','w');
      plot(mdl)
      set(gca,'fontsize',25,'TickLabelInterpreter','latex')
      xlabel('Tranmittance (Triaxus)','interpreter','latex','fontsize',35)
      ylabel('Tranmittance (CTD)','interpreter','latex','fontsize',35)
      title(['$r^{2} = $ ', num2str(coefs1(3)),' |  slope = ', num2str(coefs1(2)) ,' | mean offset = ', num2str(trans_offset),],'interpreter','latex','fontsize',35)
  end
  end
  
  % Add correction
  %all_trans_tri = + all_trans_tri - trans_offset;
  
  % Compute beam attenuation (and POC using Christinas realtionship
    all_c_tri = -(1/.25) .* log(all_trans_tri ./100); 
    cp_tri    = all_c_tri - c_sw_ctd;
    
  % Compute POC
    bbp = 0.00345.*cp_tri.^0.449;            
    POC_tri_cs = 37601*bbp + 4.95;
    POC_tri_cs =POC_tri_cs./12; % convert form gram to mol
    POC_tri_wg = 33.5.*cp_tri + 3.064;

% 2. c_p To POC relationship from staitions 1 and 2. 
% See Email from cathryn: FW: SOLACE HPLC data.  This is where
% POC_pig_CTD_data.mat files comes from
% LOAD C:CHL using HPLC Pigment Data
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Carbon
POC=readtable('POC_Summary.xlsx');


   POC_flag=POC{Station_i, 13};
   POC_Cast = POC{Station_i, 1}; 
   POC_Depth = POC{Station_i, 2}; 
   POC_dis_kg = POC{Station_i, 6};   % umol/kg         -- Discrete Bottle Fitlrations
       POC_temp=POC{Station_i, 4};
       POC_sal=POC{Station_i, 3};    % NOTE some salinity vlaues were missing in cathryns spreadsheet so i use value above in same profile.  marked in red inPOC_Summary sheet
       POC_dens =sw_pden(POC_sal,POC_temp,POC_sal.*0, 0);
       POC_dis =  POC_dis_kg.*POC_dens./1000;   % umol/L OR mmol/m3
%
for i=1:length(POC_Cast)
    ctd_trans = Profiles(POC_Cast(i)).trans;
    ctd_cp = -(1/.25) .* log(ctd_trans ./100) - c_sw_ctd ; 
    ctd_pres = Profiles(POC_Cast(i)).pres;
    depth_i = find( abs(ctd_pres-POC_Depth(i)) < 5 ) ;   % Take average sensor reading with in 2 m of bottle sample
    if isempty(depth_i);
         depth_i = find( abs(ctd_pres-POC_Depth(i)) < 7.5 ) ;   % Take average sensor reading with in 2 m of bottle sample
    end
    % Average vales with 5m of bottle sample on sensors
    CTD_cp(i)  = nanmean(ctd_cp(depth_i));
end

% Remove Flagged data
POC_dis_redflag = POC_dis(POC_flag==1);
CTD_cp_redflag  = CTD_cp(POC_flag==1);
POC_dis_qc      = POC_dis(POC_flag==0);
CTD_cp_qc       = CTD_cp(POC_flag==0);
%
% Regressions
cf_cp_poc = [ones(length(CTD_cp_qc),1) CTD_cp_qc']\POC_dis_qc;
cf_cp_poc_0 =  CTD_cp_qc'\POC_dis_qc;  % Force through 0 using only qc'd data
cf_cp_poc_0_all =  CTD_cp'\POC_dis;    % Force through 0 using all data
mdl_qc = fitlm(CTD_cp_qc',POC_dis_qc,'Intercept',false); 
mdl_all = fitlm(CTD_cp',POC_dis,'Intercept',false); 

% Transform
POC_tri_solace = cf_cp_poc_0.*cp_tri;
POC_ctd_solace = cf_cp_poc_0.*cp_ctd;

% Plot Regression
if true
    for PlotReg = 1
% figure('Units','normalized','Position',[0 0 2 2])
% clc
% set(gcf,'color','w');
figure(1)      
subplot(1,2,1)
scatter(CTD_cp_qc,POC_dis_qc,150,'b','filled'); hold on
scatter(CTD_cp_redflag,POC_dis_redflag,150,'b'); hold on
%h(1:3)=plotregression(CTD_cp,POC_dis); hold on
% r=0.898
% r2 =0.8064;
xs = [min(CTD_cp),max(CTD_cp)];
%plot(xs,cf_cp_poc(1) + cf_cp_poc(2).*xs,'r')
h(1)=plot(xs,cf_cp_poc_0.*xs,'b','linewidth',2);
h(2)=plot(xs,cf_cp_poc_0_all.*xs,'b:','linewidth',2);
xlim([0,max(CTD_cp)])
set(gca,'fontsize',20,'TickLabelInterpreter','latex')
xlabel('$c_p$ $(m^{-1})$','interpreter','latex','fontsize',25)
ylabel('POC $(\mu M)$','interpreter','latex','fontsize',25)
%title(['$r^{2} = $', num2str(mdl_qc.Rsquared.Adjusted),' | slope = ', num2str(mdl_qc.Coefficients.Estimate)],'interpreter','latex','fontsize',35)
%legend(h(1:2),'QCd Data','All Data','interpreter','latex','fontsize',35) 
%
subplot(1,2,2)
scatter(cp_ctd,POC_ctd_solace,'b'); hold on
scatter(cp_ctd,POC_ctd_cs,'m'); hold on
scatter(cp_ctd,POC_ctd_wg,'g'); hold on
set(gca,'fontsize',20,'TickLabelInterpreter','latex')
%legend('Solace Site Samples', 'SOTS (Schallenberg et al)', ' APFZ (Gardner et al)')
xlabel('$c_p$ (CTD)','interpreter','latex','fontsize',25)
ylabel('POC (uM)','interpreter','latex','fontsize',25)
xlim([0,.4])
ylim([0,15])



    end
end

end
% Save Corrections
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus
save(fname,'cf_ctd', 'cf_cp_poc_0','cf_cp_poc_0_all', 'c_sw_ctd', 'trans_offset')

%% Sample Plot
for splots =1
figure('Units','normalized','Position',[0 0 2 2])
clc
set(gcf,'color','w');

subplot(2,2,1)
scatter(all_trans_ctd,all_pres_ctd,'r'); hold on
 set(gca,'fontsize',25,'TickLabelInterpreter','latex')
xlim([90,100]); 
ylim([0,200])
set(gca,'ydir','rev'); 
title('CTD','interpreter','latex','fontsize',35)
xlabel('Transmitttance','interpreter','latex','fontsize',25)
ylabel('Depth (m)','interpreter','latex','fontsize',25)

subplot(2,2,3)
scatter(all_trans_tri,all_pres_tri,'b');
 set(gca,'fontsize',25,'TickLabelInterpreter','latex')
%xlim([90,100]); 
xlim([90,100]); 
ylim([0,200])
set(gca,'ydir','rev'); 
title('Triaxus','interpreter','latex','fontsize',35)
xlabel('Transmitttance','interpreter','latex','fontsize',25)
ylabel('Depth (m)','interpreter','latex','fontsize',25)

% subplot(2,2,1)
% scatter(cp_ctd,all_pres_ctd,'r'); hold on
% xlim([0,0.4]); 
% ylim([0,200])
% set(gca,'ydir','rev'); 
% title('CTD','interpreter','latex','fontsize',15)
% xlabel('$c_p$  $(m^{-1})$','interpreter','latex','fontsize',15)
% ylabel('Depth (m)','interpreter','latex','fontsize',15)
% 
% subplot(2,2,3)
% scatter(cp_tri,all_pres_tri,'b');
% %xlim([90,100]); 
% xlim([0,0.4]); 
% ylim([0,200])
% set(gca,'ydir','rev'); 
% title('CTD','interpreter','latex','fontsize',15)
% xlabel('$c_p$  $(m^{-1})$','interpreter','latex','fontsize',15)
% ylabel('Depth (m)','interpreter','latex','fontsize',15)

subplot(2,2,2)
%scatter((cf_cp_poc(1) + cf_cp_poc(2).*cp_ctd),all_pres_ctd,'r'); hold on
scatter((cf_cp_poc_0.*cp_ctd),all_pres_ctd,'r'); hold on
xlim([0,8]); 
ylim([0,200])
set(gca,'ydir','rev'); 
title('Triaxus','interpreter','latex','fontsize',15)
xlabel('POC (uM)','interpreter','latex','fontsize',15)
ylabel('Depth (m)','interpreter','latex','fontsize',15)

subplot(2,2,4)
%scatter(cf_cp_poc(1) + cf_cp_poc(2).*cp_tri,all_pres_tri,'b');
scatter(cf_cp_poc_0.*cp_tri,all_pres_tri,'b');
%xlim([90,100]); 
xlim([0,8]); 
ylim([0,200])
set(gca,'ydir','rev'); 
title('Triaxus','interpreter','latex','fontsize',15)
xlabel('POC (uM)','interpreter','latex','fontsize',15)
ylabel('Depth (m)','interpreter','latex','fontsize',15)

end

%% Site Profiles
% Create Figure
figure('Units','normalized','Position',[0 0 2 2])
clc
set(gcf,'color','w');
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/Trixaus

% SOTS 
% Stationi = 3:39;
% tx = 'SOTS';
% load poc_corrections_sots.mat

% % SOUTH 1 
 South1_CTDi = 41:63; 
 Stationi = 41:63;
 tx = 'SOUTH 1';
 load poc_corrections_south.mat
% 
% % SOUTH 2
% South2_CTDi = 64:85;
% Stationi = 64:85;
% tx = 'SOUTH 2';
% load poc_corrections_south.mat


maxdepth = 500;

cmap = brewermap(140,'RdYlBu');
cday = casttime_loc(Stationi); cday = cday-cday(1); cday = 1+round(cday*10);


for i = Stationi
subplot(2,3,1)
if any(tri_i-i==0)
line(Profiles(i).temp,Profiles(i).pres,'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',3)
else
line(Profiles(i).temp,Profiles(i).pres,'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',0.5)
end
set(gca,'ydir','rev');  
ylim([0,maxdepth]);
if maxdepth >1050; xlim([0,4.5]); end 
set(gca,'fontsize',15,'TickLabelInterpreter','latex')
ylabel('Depth (dbar)','interpreter','latex','fontsize',18)
xlabel('Temperature (C)','interpreter','latex','fontsize',18)

subplot(2,3,2)
if any(tri_i-i==0)
line(Profiles(i).sal,Profiles(i).pres,'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',3)
else
line(Profiles(i).sal,Profiles(i).pres,'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',0.5)
end
set(gca,'ydir','rev');
ylim([0,maxdepth])
if maxdepth >1050; xlim([33.79,34.45]); end 
set(gca,'fontsize',15,'TickLabelInterpreter','latex')
ylabel('Depth (dbar)','interpreter','latex','fontsize',18)
xlabel('Salinty (psu)','interpreter','latex','fontsize',18)

subplot(2,3,3)
if any(tri_i-i==0)
line(Profiles(i).par,Profiles(i).pres,'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',3)
else
line(Profiles(i).par,Profiles(i).pres,'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',0.5)
end
set(gca,'ydir','rev');
if maxdepth >1050; xlim([0,1000]); end 
set(gca,'xscale','log')
ylim([0,maxdepth])
set(gca,'fontsize',15,'TickLabelInterpreter','latex')
ylabel('Depth (dbar)','interpreter','latex','fontsize',18)
xlabel('PAR (uE/m2/s)','interpreter','latex','fontsize',18)

subplot(2,3,4)
if any(tri_i-i==0)
line(CHL_ctd(i,:),Profiles(i).pres(1:500),'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',3)
else
line(CHL_ctd(i,:),Profiles(i).pres(1:500),'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',0.5)
end
set(gca,'ydir','rev');
ylim([0,maxdepth])
if maxdepth >1050; xlim([0,1.1]); end 
set(gca,'fontsize',15,'TickLabelInterpreter','latex')
ylabel('Depth (dbar)','interpreter','latex','fontsize',18)
xlabel('Chl Flo (mg/m3)','interpreter','latex','fontsize',18)

subplot(2,3,5)
ctd_trans = Profiles(i).trans; 
ctd_cp = -(1/.25) .* log(ctd_trans ./100) - c_sw_ctd;
if any(tri_i-i==0)
line(cf_cp_poc_0.*ctd_cp,Profiles(i).pres,'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',3)
else
line(cf_cp_poc_0.*ctd_cp,Profiles(i).pres,'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',0.5)
end
set(gca,'ydir','rev');
ylim([0,maxdepth]); 
if maxdepth >1050; xlim([0,9]); end 
xlim([0,15])
set(gca,'fontsize',15,'TickLabelInterpreter','latex')
ylabel('Depth (dbar)','interpreter','latex','fontsize',18)
xlabel('POC (uM)','interpreter','latex','fontsize',18)

subplot(2,3,6)
ctd_trans = Profiles(i).trans(1:500); 
ctd_cp = -(1/.25) .* log(ctd_trans ./100) - c_sw_ctd;
if any(tri_i-i==0)
line(CHL_ctd(i,:)./(12.*cf_cp_poc_0.*ctd_cp),...
      Profiles(i).pres(1:500),'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',3)
else
line(CHL_ctd(i,:)./(12.*cf_cp_poc_0.*ctd_cp),...
      Profiles(i).pres(1:500),'Color',cmap(cday(i-Stationi(1)+1),:),'linewidth',0.5)
end
set(gca,'ydir','rev');  
ylim([0,maxdepth])
if maxdepth >1050; xlim([0,.05]); end 
set(gca,'fontsize',15,'TickLabelInterpreter','latex')
ylabel('Depth (dbar)','interpreter','latex','fontsize',18)
xlabel('Chl:C','interpreter','latex','fontsize',18)
end

colormap(cmap); caxis([0,10])
cloc = get(gca,'Position');
h= colorbar(gca,'position', [cloc(1)+cloc(3)+.02, cloc(2), .015, cloc(4).*2.45 ], ...
    'orientation','vertical','TickLabelInterpreter','latex');
set(get(h,'ylabel'),'string', 'Days Since First Cast at Station','interpreter','latex','fontsize',20)


% 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');




%%
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles
load('MLD.mat', 'ctd_mld','casttime_loc')

South1_CTDi = 41:63; 
South2_CTDi = 64:85;

Stationi = South1_CTDi;
fname = 'DCM_ctd_mets_South_1';
clearvars  ctd_met_ts


for i = Stationi

ctd_trans = Profiles(i).trans; 
ctd_cp    = -(1/.25) .* log(ctd_trans ./100) - c_sw_ctd;
ctd_poc   = cf_cp_poc_0.*ctd_cp;   ctd_poc=movmean(ctd_poc,20,'omitnan');
ctd_chl   = CHL_ctd(i,:);          % Already smoothed with 20m moving mean in Chl_Build_plotter
ctd_pres  = Profiles(i).pres; 
mld       = max(ctd_mld(i),50); 
time      = casttime_loc(i);
ctd_chlc      = ctd_chl./(12.*ctd_poc(1:500));  

p = (i+1)-Stationi(1);

ctd_met_ts(13,2,p) = mld; 
ctd_met_ts(14,1,p) = time; 

% DCM
prof = ctd_chl(1:150);   
        [max_m, max_i] = max(prof);
        % Assign metrics if theres a DCM
        if pres(squeeze(max_i))>mld %&& max_m> 1.05.*nanmean(prof(pres<mld)) %  % IS DCM BELOW MLD
        dcm_thr =  nanmean(prof(pres<mld)) + .5.* abs(max_m-nanmean(prof(pres<mld)));  % 50% of difference between MLD and max value 
        DCM = find(prof>dcm_thr) ; 
        DCM(pres(DCM) < (mld-30)) = [] ;  % Dont condsider values 15 m above MLD
        ctd_met_ts(4,1,p) = squeeze(max_m); % Size of DCM
        ctd_met_ts(4,2,p) = pres(squeeze(max_i)); % Depth of DCM
        ctd_met_ts(4,3,p) = pres(DCM(1));         % Top of DCM
        ctd_met_ts(4,4,p) = pres(DCM(end));       % Base of DCM
        ctd_met_ts(4,5,p) = pres(DCM(end)) - pres(DCM(1));    % Width of DCM (m)
        ctd_met_ts(4,6,p) = nanmean(prof(DCM)).*(pres(DCM(end)) - pres(DCM(1))); % integrate of DCM (mg/m2)
        ctd_met_ts(4,7,p) = 100.*abs(nanmean(prof(pres<mld))- max_m)./nanmean(prof(pres<mld)); 
        else
        ctd_met_ts(4,1:8,p) = nan;
        end
    
     % Find DBM 
     prof = ctd_poc(1:150);        
        [max_m, max_i] = max(prof);
        % Assign metrics if theres a DCM
         if pres(squeeze(max_i))>mld %%&& max_m> 1.1.*nanmean(prof(pres<mld)) %  % IS DCM BELOW MLD
        dcm_thr =  nanmean(prof(pres<mld)) + .5.* abs(max_m-nanmean(prof(pres<mld)));  % 50% of difference between MLD and max value 
        DCM = find(prof>dcm_thr) ; 
        DCM(pres(DCM) < (mld-30)) = [] ;  % Dont condsider values 15 m above MLD
        ctd_met_ts(11,1,p) = squeeze(max_m); % Size of DCM
        ctd_met_ts(11,2,p) = pres(squeeze(max_i)); % Depth of DCM
        ctd_met_ts(11,3,p) = pres(DCM(1));         % Top of DCM
        ctd_met_ts(11,4,p) = pres(DCM(end));       % Base of DCM
        ctd_met_ts(11,5,p) = pres(DCM(end)) - pres(DCM(1));    % Width of DCM (m)
        ctd_met_ts(11,6,p) = nanmean(prof(DCM)).*(pres(DCM(end)) - pres(DCM(1))); % integrate of DCM (mg/m2)
        ctd_met_ts(11,7,p) = 100.*abs(nanmean(prof(pres<mld))- max_m)./nanmean(prof(pres<mld)); 
        else
        ctd_met_ts(11,1:8,p) = nan;
        end

     
    % Find Chl C Max
    prof = ctd_chlc(1:150);        
        [max_m, max_i] = max(prof);
         % Assign metrics if theres a DCM
        if pres(squeeze(max_i))>mld %%&& max_m> 1.1.*nanmean(prof(pres<mld)) %  % IS DCM BELOW MLD
        dcm_thr =  nanmean(prof(pres<mld)) + .5.* abs(max_m-nanmean(prof(pres<mld)));  % 50% of difference between MLD and max value 
        DCM = find(prof>dcm_thr) ; 
        DCM(pres(DCM) < (mld-30)) = [] ;  % Dont condsider values 15 m above MLD
        ctd_met_ts(12,1,p) = squeeze(max_m); % Size of DCM
        ctd_met_ts(12,2,p) = pres(squeeze(max_i)); % Depth of DCM
        ctd_met_ts(12,3,p) = pres(DCM(1));         % Top of DCM
        ctd_met_ts(12,4,p) = pres(DCM(end));       % Base of DCM
        ctd_met_ts(12,5,p) = pres(DCM(end)) - pres(DCM(1));    % Width of DCM (m)
        ctd_met_ts(12,6,p) = nanmean(prof(DCM)).*(pres(DCM(end)) - pres(DCM(1))); % integrate of DCM (mg/m2)
        ctd_met_ts(12,7,p) = 100.*abs(nanmean(prof(pres<mld))- max_m)./nanmean(prof(pres<mld)); 
        else
        ctd_met_ts(12,1:8,p) = nan;
        end


end

cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/CTD
save(fname, 'ctd_met_ts')

%% Regression plot for all 
   met_names{1} = 'Magnitude of ';
   met_names{2} = 'Depth of ';
   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 ';
   
   var_names{4}  = 'Deep Chlorrophyl Maxmima ';
   var_names{11} = 'Deep Biomass Maxmima ';
   var_names{12} = 'Chl:C Maxima';
   var_names{13} = 'MLD';
   var_names{14} = 'Time'; 
   
   for units=1
   met_units{4,1} = '$(mg/m^3)$'; met_units{11,1} = '$(uM)$';        met_units{12,1} = '$(-)$';   met_units{14,1} = 'Day'; 
   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} = '$(-)$'; 
   met_units{4,7} = '$(\%)$';     met_units{11,7} = '$(\%)$';        met_units{12,7} = '$(\%)$'; 
   end
   
   
vars1 = [11, 11,  11, 12, 12, 12 ];
mets1 = [2,  1,  7,   2,   1,  7];   
   
vars2 = [4, 4,  4, 4, 4, 4];
mets2 = [2, 1,  7, 2, 1, 7 ];

% 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 ];


% Vs Time
vars1 = [14, 14,  14, 13, 13, 13 ];
mets1 = [1,  1,  1,   2,   2,  2];   
   
vars2 = 12*[1, 1,  1, 1, 1,1];
mets2 = [1, 6,  7,  1,  6, 7 ];


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

cmap = brewermap(100,'RdYlBu');
cday = casttime_loc(Stationi); cday = cday-cday(1); cday = 1+round(cday*10);

msize =100;

% Make Figure
figure('Units','normalized','Position',[0 0 2 2])
clc
set(gcf,'color','w');
   for i=1:6
   subplot(2,3,i)
           
   % Both
   ALL_var1 =[]; ALL_var2 = [];
   
   % South 1 
        load('DCM_ctd_mets_South_1.mat');
        site_prof1_time = ctd_met_ts(14,1,1);
        var1 = squeeze(ctd_met_ts(vars1(i),mets1(i),:)); 
        var2 = squeeze(ctd_met_ts(vars2(i),mets2(i),:)); 
        South1_var1 = [var1]; ALL_var1 =[ALL_var1;var1];
        South1_var2 = [var2]; ALL_var2 =[ALL_var2;var2];

        % Scatter Plot 
        scatter(var1,var2,msize, squeeze(ctd_met_ts(14,1,:)- site_prof1_time), 'filled'); hold on
       % scatter(var1,var2,msize, 'r', 'filled'); hold on
   
   % South 2
        load('DCM_ctd_mets_South_2.mat');
         site_prof2_time = ctd_met_ts(14,1,1);
        var1 = squeeze(ctd_met_ts(vars1(i),mets1(i),:)); 
        var2 = squeeze(ctd_met_ts(vars2(i),mets2(i),:)); 
        South2_var1 = [var1]; ALL_var1 =[ALL_var1;var1];
        South2_var2 = [var2]; ALL_var2 =[ALL_var2;var2];

        % Scatter Plot
        scatter(var1,var2,msize, squeeze(ctd_met_ts(14,1,:)-site_prof1_time), 'filled'); hold on
       % scatter(var1,var2,msize, 'b', 'filled'); hold on
   
   % Add Regressions 
     % South 1
     mdl = fitlm(South1_var1,South1_var2); 
     coefs_reg(:,1,i) = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
     % South 1
     mdl = fitlm(South2_var1,South2_var2); 
     coefs_reg(:,2,i) = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
     % Both
     mdl = fitlm(ALL_var1,ALL_var2); 
     coefs_reg(:,3,i) = [mdl.Coefficients.Estimate; mdl.Rsquared.Adjusted]; 
    
     % Plot regression
     xs = [min(squeeze(ALL_var1)), max(squeeze(ALL_var1))];
     h(1) = plot(xs, coefs_reg(1,1,i)+coefs_reg(2,1,i).*xs,'color',[1,0,0,.4],'linewidth',2); hold on
     h(2) = plot(xs, coefs_reg(1,2,i)+coefs_reg(2,2,i).*xs,'color',[0,0,1,.4],'linewidth',2); hold on
     %h(3) = plot(xs, coefs_reg(1,3,i)+coefs_reg(2,3,i).*xs,'color',[0,0,0,.7],'linewidth',2); hold on
     legend(h(1:2),['South 1: $r^2$ = ',num2str(round(100*coefs_reg(3,1,i))/100) ],...
                   [ 'South 2: $r^2$ = ',num2str(round(100*coefs_reg(3,2,i))/100)],...
                   [ 'Both: $r^2$ = ',num2str(round(100*coefs_reg(3,3,i))/100)],'fontsize',15,'interpreter','latex','location','best')

   % Labels and Lim 
     if vars1(i)==14
     xlabel([met_units{vars1(i), mets1(i)} ],'interpreter','latex','fontsize',15)
     xline(site_prof2_time)
     datetick('x', 'dd', 'keeplimits')
     else
     xlabel([met_names{mets1(i)}, var_names{vars1(i)}, met_units{vars1(i), mets1(i)} ],'interpreter','latex','fontsize',15)
     end
     ylabel([met_names{mets2(i)}, var_names{vars2(i)}, met_units{vars2(i),mets2(i)} ],'interpreter','latex','fontsize',15)
     xlim([min(squeeze(ALL_var1(~isnan(ALL_var2)))), max(squeeze(ALL_var1(~isnan(ALL_var2))))])
     ylim([min(squeeze(ALL_var2(~isnan(ALL_var1)))), max(squeeze(ALL_var2(~isnan(ALL_var1))))])

   % Colorbar
     colormap(cmap)
     caxis([0,17])
   end

   % colormap(cmap); caxis([0,10])
     cloc = get(gca,'Position');
     h= colorbar(gca,'position', [cloc(1)+cloc(3)+.02, cloc(2), .015, cloc(4).*2.35 ], ...
    'orientation','vertical','TickLabelInterpreter','latex');
     set(get(h,'ylabel'),'string', 'Days Since First Cast at South Station 1','interpreter','latex','fontsize',20)


   


