%%%%%%%%%%%%%%%%%%%%%%
%% Compile MLD Data %%
%%%%%%%%%%%%%%%%%%%%%%
% Compile MLD data from all casts (CTD, XBT...TMR, Triaxus)
% MLD computing algorithm from Talle & holte 2008 
%     http://journals.ametsoc.org/doi/abs/10.1175/2009JTECHO543.1


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

% Load Interpolated Underway Data
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Underway
load iUnderway t_ncp ilat ilon time_uw
load SiteDeploymentDates

% Load CTD Data
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/CTD
load CTD_MLD 

% Load XBT Data
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/XBT
load XBT_MLD 

%     % remove first cast
%     xbt_mld = xbt_mld(2:end); xbt_time_loc = xbt_time_loc(2:end); 
%     xbt_lat = xbt_lat(2:end); xbt_lon = xbt_lon(2:end); 

% Load TMR Data
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles/TMR
load TMR_MLD 

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


%% %% %% Compile In Situ MLD %% %% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if  false
    for MakeInSitu = 1
% Concatinate
MLD      = [xbt_mld, ctd_mld, tmr_mld, tri_mld];  %   , tri_mld removed
MLD_min  = [xbt_min, ctd_min, tmr_min, tri_min];  %   , tri_min removed
MLD_max  = [xbt_max, ctd_max, tmr_max, tri_max];   %   , tri_max removed

time_mld_loc = [xbt_time_loc, casttime_loc, tmr_time_loc,tri_time_loc]; % tri_time_loc removed
time_mld_utc = [xbt_time_utc, casttime_utc, tmr_time_utc,tri_time_utc]; % tri_time_utc removed

% Sort by time
[time_mld_loc, time_ind] = sort(time_mld_loc);
MLD = MLD(time_ind);
MLD_min = MLD_min(time_ind);
MLD_max = MLD_max(time_ind);


% remove overlapping times  --- Not sure why they are though...(52 and 59; similiar mlds though) 
for j=1:2
for i = find(diff(time_mld_loc)==0)
time_mld_loc(i) = [];
MLD(i)      = [];
MLD_min(i)  = [];
MLD_max(i)  = [];
end
end

% interpolate in hourly data
%tinterp_mld = time_mld_loc(1):1/24:time_mld_loc(end);
iMLD          = interp1(time_mld_loc,MLD,t_ncp);
iMLDmin       = interp1(time_mld_loc,MLD_min,t_ncp);
iMLDmax       = interp1(time_mld_loc,MLD_max,t_ncp);
 
% smooth with 2 day window
siMLD = movmean(iMLD, 24*2,'omitnan');
siMLDmin = movmean(iMLDmin, 24*2,'omitnan');
siMLDmax = movmean(iMLDmax, 24*2,'omitnan');

  %% save
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles
save MLD tinterp_mld time_mld_loc time_mld_utc ...
        siMLD siMLDmin siMLDmax  ...
        xbt_time_utc casttime_utc tmr_time_utc tri_time_utc ...
        xbt_time_loc casttime_loc tmr_time_loc tri_time_loc ...
        xbt_mld ctd_mld tmr_mld tri_mld ...
        xbt_min ctd_min tmr_min tri_min ...
        xbt_max ctd_max tmr_max tri_max ...
        

    end
else 
cd /Users/roh029/Desktop/Data.nosync/SOLACE/Profiles
load MLD
end

%% %% %% Load Reanalysis MLD %% %% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if  false
    for MakeSOSE =1
%  SOSE Reference 

% mazloffEddyPermittingSouthernOcean2010a
% SOSE reference M. Mazloff, P. Heimbach, and C. Wunsch, 2010: "An Eddy-Permitting Southern Ocean State Estimate." J. Phys. Oceanogr., 40, 880–899. doi: 10.1175/2009JPO4236.1

% verdyDataAssimilatingModel2017
% B-SOSE reference A. Verdy and M. Mazloff, 2017: " A data assimilating model for estimating Southern Ocean biogeochemistry." J. Geophys. Res. Oceans.,  122, doi:10.1002/2016JC012650.

cd /Users/roh029/Desktop/Data.nosync/SOLACE/MLD/SOSE
NCsoseMLD = 'bsose_Iter135_2013to2019_5day_MLD.nc';

% Load File 
Lat_sose  = ncread(NCsoseMLD, 'YC');                                    % Degrees North                  (588x1)
Lon_sose  = ncread(NCsoseMLD, 'XC');                                    % Degrees East                   (2160x1)
time_sose = ncread(NCsoseMLD, 'time');                                  % Seconds since 2012-12-01       (511x1)
days_sose = double(time_sose)./60./60./24 -30;                           % Time to days since Jan 1 2000  (511x1)                         
time_sose = datenum('01-Dec-2012 00:00:00') + double(time_sose)./60./60./24;    % Time to days since Jan 1 2000  (511x1)                         
  
MLD_sose  = ncread(NCsoseMLD, 'BLGMLD'); % m                        (2160x588x511)

%% Grid
for y = 1:7
    y
    for d = 1:73
    MLD_CLIM(:,:,d,y) = MLD_sose(:,:,(y-1)*73+d); 
    end
end
    
% Exmaple Plot
% imagesc(Lon_sose,Lat_sose,squeeze(MLD_CLIM(:,:,30,1))'); colorbar; caxis([0,200]); set(gca, 'ydir','normal')

% Create Climatology
MLD_sose_clim = squeeze(mean(MLD_CLIM,4));
doy = (1:73)*5;
% reorganize to start June 24 !(doy 175)
MLD_sose_clim_aust = cat(3,MLD_sose_clim(:,:,35:end),MLD_sose_clim(:,:,1:34));
time_sose_clim_aus =  datenum('19-June-2020 00:00:00',date_frmt) + doy;


% Interpolate along boat track
 MLD_sose  =  interp3(Lat_sose,Lon_sose,time_sose_clim_aus,MLD_sose_clim_aust,ilat ,ilon,t_ncp);
 %squeeze(mean(mean(interp3(Lat_sose,Lon_sose,time_sose_clim_aus,MLD_sose_clim_aust, xq , yq, tq))));
 
% Smooth
siMLD_sose = movmean(MLD_sose, 24,'omitnan');

    % Plot single day distribution
    %imagesc(Lon_sose,(Lat_sose),(squeeze(MLD_sose_clim(:,:,70))')); colorbar; set(gca,'ydir','normal')
    % Plot single place TS
    % plot(squeeze(MLD_sose_clim(850,427,:))) @ -50Nand 139E
    


%% Save 
cd /Users/roh029/Desktop/Data.nosync/SOLACE/MLD/SOSE
save MLD_SOSE_clim MLD_SOSE_clim Lat_sose Lon_sose doy MLD_sose_clim MLD_sose_clim_aust time_sose_clim_aus 
save MLD_SOSE siMLD_sose MLD_sose t_ncp

    end
else   
    cd /Users/roh029/Desktop/Data.nosync/SOLACE/MLD/SOSE
    load MLD_SOSE siMLD_sose MLD_sose t_ncp
end


%% PLOT
    figure('Units','normalized','Position',[0 1 .9 .8])
    set(gcf,'color','w');

    subplot(2,1,1)
    yyaxis right
    plot(t_ncp,ilon,'linewidth',2); hold on
    ylabel('Longitude ($^\circ S$)','interpreter','latex')
    yyaxis left
    plot(t_ncp,ilat,'linewidth',2); hold on
        xline(de(3),'k','linewidth',2); hold on;    xline(de(4),'k','linewidth',2); hold on    % SOTS
        xline(de(7),'k:','linewidth',2); hold on;   xline(de(8),'k:','linewidth',2); hold on   % South 1
        xline(de(9),'k-.','linewidth',2); hold on;  xline(de(10),'k-.','linewidth',2); hold on % South 2
    xlim([t_ncp(1),t_ncp(end)])
    datetick(gca,'x','mmm-dd','keeplimits')   
    set(gca,'fontsize',15,'TickLabelInterpreter','latex')
    ylabel('Latitude ($^\circ S$)','interpreter','latex')


    
 %
 subplot(2,1,2)
    %plot(tinterp_mld,iMLD,'k','linewidth',2); hold on
    scatter(casttime_loc,ctd_mld,'k','linewidth',1); hold on
    scatter(xbt_time_loc,xbt_mld,'kx','linewidth',1); hold on
    scatter(tmr_time_loc,tmr_mld,'k*','linewidth',1); hold on
    scatter(tri_time_loc,tri_mld,'k^','linewidth',1); hold on
    % MLD Reanalysis averaged onto box at each station
     %  plot(t_sots,sots_mld_hy,'b','linewidth',1); hold on;       plot(t_sots,sots_mld_so,'b:','linewidth',1); hold on; 
     %  plot(t_south1,south1_mld_hy,'b','linewidth',1); hold on;   plot(t_south1,south1_mld_so,'b:','linewidth',1); hold on;
     %  plot(t_south2,south2_mld_hy,'b','linewidth',1); hold on;   plot(t_south2,south2_mld_so,'b:','linewidth',1); hold on;
    plot(t_ncp,siMLD,'r','linewidth',2); hold on
    plot(t_ncp,siMLDmin,'k','linewidth',1); hold on
    plot(t_ncp,siMLDmax,'k','linewidth',1); hold on
   
    
    
    %errorbar(casttime_loc,ctd_mld,ctd_min, ctd_max,'o') ; hold on
    
    
        xline(de(3),'k','linewidth',2); hold on;    xline(de(4),'k','linewidth',2); hold on    % SOTS
        xline(de(7),'k:','linewidth',2); hold on;   xline(de(8),'k:','linewidth',2); hold on   % South 1
        xline(de(9),'k-.','linewidth',2); hold on;  xline(de(10),'k-.','linewidth',2); hold on % South 2

    set(gca,'ydir','rev')
    ylim([0,175])
    xlim([t_ncp(1),t_ncp(end)])
    datetick(gca,'x','mmm-dd','keeplimits')  
    set(gca,'fontsize',15,'TickLabelInterpreter','latex')
    ylabel('Mixed Layer Depth (m)','interpreter','latex')
    legend('CTD','XBT','TMR','Triaxus','fontsize',18,'interpreter','latex','location','southeast')


    
    
%%     
%%%% Appendix A
% JAMESTEC might be old, use the SIO one i got processing scripts from?

% Clim MLD
% Monthly Climatologies used by liSatelliteEstimatesNet2016 
% Hosoda, S., T. Ohira, K. Sato, and T. Suga (2010), Improved description of global mixed‐layer 
% depth using Argo profiling floats, J. Oceanogr., 66(6), 773–787.
% hosodaImprovedDescriptionGlobal2010
% Download: http://www.jamstec.go.jp/ARGO/argo_web/ancient/MILAGPV/index_e.html


%fprintf('\n MLD Data has been loaded and processed \n')

    
    
