% Make MLD timeseries for OUTPACE cruise, merging the CTD's and MVP
% datasets. The definition used is 0.03 sigma-theta more than 10 m depth.
% Since temperature is more reliable than density for the MVP, here we use
% the temperature change equivalent for the range of T-S space we're in.

% Values of practical salinity and temperature at 10 m
psal10 = [];
temp10 = [];

% Loaded up CTD data w/ the script, now cycle through the casts.
for j = 1:length(CTD)
    tempt = CTD(j).temp;
    temps = CTD(j).psal;
    tempp = CTD(j).pres;
    ii = find(tempp==10,1,'first');
    psal10 = [psal10 temps(ii)];
    temp10 = [temp10 tempt(ii)];
end

clearvars -except psal10 temp10

% Since MVP salinity is suspect, we won't use the data to find the range of
% salinity values to consider.

% % Now load up the MVP data...
% mvp_dir = '/Users/Alain/Desktop/MIO/MVP_Dry2/L2_OUTPACE/';
% lda = 'AdeV_uCTD_slow_LDA_ALL/';
% ldb = 'AdeV_uSVP_slow_LDB_ALL/';
% ldc = 'AdeV_uSVP_slow_LDC_ALL/';
% 
% cd([mvp_dir lda]);
% load Transect_01_10_qc_2.mat
% for j = 1:length(data)
%     tempt = data(j).temp.value;
%     temps = data(j).psal.value;
%     tempp = data(j).pres.value;
%     ii = find(tempp==10,1,'first');
%     psal10 = [psal10 temps(ii)];
%     temp10 = [temp10 tempt(ii)];
% end
%     
% cd([mvp_dir ldb]);
% load Transect_19_30_qc_2.mat
% for j = 1:length(data)
%     tempt = data(j).temp.value;
%     temps = data(j).psal.value;
%     tempp = data(j).pres.value;
%     ii = find(tempp==10,1,'first');
%     psal10 = [psal10 temps(ii)];
%     temp10 = [temp10 tempt(ii)];
% end
% 
% cd([mvp_dir ldc]);
% load Transect_31_33_qc_2.mat
% for j = 1:length(data)
%     tempt = data(j).temp.value;
%     temps = data(j).psal.value;
%     tempp = data(j).pres.value;
%     ii = find(tempp==10,1,'first');
%     psal10 = [psal10 temps(ii)];
%     temp10 = [temp10 tempt(ii)];
% end

% Find the 'normal' values of salinity and temperature
Salval = nanmedian(psal10);
Temval = nanmedian(temp10);

rho_val = gsw_rho(Salval,Temval,0);

% get the temperature change for this rho to increase by 0.03...
rho_03 = rho_val+0.03;
t_trial = linspace(Temval,Temval-3,150);
r_trial = gsw_rho(Salval.*ones(1,length(t_trial)),t_trial,0);
ii = find(r_trial<=rho_03,1,'last');
temdiff = Temval - t_trial(ii);
clearvars -except temdiff

%% Now go back to the datasets, and make a list of Time, Lon, Lat, and MLD

% Do CTD's first...
Time = [];
Lon = [];
Lat = [];
MLD = [];

for j = 1:length(CTD)
    Time = [Time CTD(j).Time];
    Lon = [Lon CTD(j).Lon];
    Lat = [Lat CTD(j).Lat];
    temp = CTD(j).pres;
    temt = CTD(j).temp;
    ii = find(temp<=10,1,'last'); % Find the 10 dBar position
    if isempty(ii)
        ii = 1; % Or use the first observation...
    end
% Check to see if temperature is NaN... if so, go to first non-NaN value
    if isnan(temt(ii))
        posval = find(~isnan(temt));
        ii = posval(1);
        tval = temt(ii);
    else
        tval = temt(ii);
    end
    tval2 = tval-temdiff;
    iii = find(temt>=tval2,1,'last'); % find where the temperature changes
    MLD = [MLD temp(iii+1)];
    
end

% Now for the MVP Transects...
mvp_dir = '/Users/Alain/Desktop/MIO/MVP_Dry2/L2_OUTPACE/';
lda = 'AdeV_uCTD_slow_LDA_ALL/';
ldb = 'AdeV_uSVP_slow_LDB_ALL/';
ldc = 'AdeV_uSVP_slow_LDC_ALL/';

cd([mvp_dir lda]);
load Transect_01_10_qc_2.mat
for j = 1:length(data)
    Time = [Time metadata(j).date.value];
    Lon = [Lon metadata(j).lon.value];
    Lat = [Lat metadata(j).lat.value];
    temp = data(j).pres.value;
    temt = data(j).temp.value;
    ii = find(temp<=10,1,'last');
    if isempty(ii)
        ii = 1;
    end
% Check to see if temperature is NaN... if so, go to first non-NaN value
    if isnan(temt(ii))
        posval = find(~isnan(temt));
        ii = posval(1);
        tval = temt(ii);
    else
        tval = temt(ii);
    end
    tval2 = tval-temdiff;
    iii = find(temt>=tval2,1,'last'); % find where the temperature changes
    MLD = [MLD temp(iii+1)];
    
end

cd([mvp_dir ldb]);
load Transect_19_30_qc_2.mat
for j = 1:length(data)
    Time = [Time metadata(j).date.value];
    Lon = [Lon metadata(j).lon.value];
    Lat = [Lat metadata(j).lat.value];
    temp = data(j).pres.value;
    temt = data(j).temp.value;
    ii = find(temp<=10,1,'last');
    if isempty(ii)
        ii = 1;
    end
% Check to see if temperature is NaN... if so, go to first non-NaN value
    if isnan(temt(ii))
        posval = find(~isnan(temt));
        ii = posval(1);
        tval = temt(ii);
    else
        tval = temt(ii);
    end
    tval2 = tval-temdiff;
    iii = find(temt>=tval2,1,'last'); % find where the temperature changes
    MLD = [MLD temp(iii+1)];


end

cd([mvp_dir ldc]);
load Transect_31_33_qc_2.mat
for j = 1:length(data)
    Time = [Time metadata(j).date.value];
    Lon = [Lon metadata(j).lon.value];
    Lat = [Lat metadata(j).lat.value];
    temp = data(j).pres.value;
    temt = data(j).temp.value;
    ii = find(temp<=10,1,'last');
    if isempty(ii)
        ii = 1;
    end
% Check to see if temperature is NaN... if so, go to first non-NaN value
    if isnan(temt(ii))
        posval = find(~isnan(temt));
        ii = posval(1);
        tval = temt(ii);
    else
        tval = temt(ii);
    end
    tval2 = tval-temdiff;
    iii = find(temt>=tval2,1,'last'); % find where the temperature changes
    MLD = [MLD temp(iii+1)];


end

% Sort list into chronological order...
[Time ii] = sort(Time);
Lon = Lon(ii);
Lat = Lat(ii);
MLD = MLD(ii);
T = table(Time',Lon',Lat',MLD','VariableNames',{'Date','Longitude','Latitude','MLD'});
T.Date = datetime(T.Date,'ConvertFrom','datenum');
cd /Users/Alain/Desktop/MIO/Outpace/MLD/
writetable(T,'OUTPACE_MLD.txt');
