%% Script for Andrea
%
% The smoothed profiles are saved in the data structures as:
% data(1).denSmooth.value

clear all; clc;
filename = 'transect14_smooth';
filedir = '/Users/Alain/Desktop/OSCAHR_Andrea';
%filedir = '/home/andrea/MYMATLAB/???';

% Plots of N2, Richardson_phi, and Richardson_adcp

load([filedir filesep filename]);

nProf = length(data);
nDep = length(data(1).pres.value);

% Load necessary variables into matrices
pres = NaN(nDep,nProf); % matrix of pressure values (for plotting).
temp = NaN(nDep,nProf);
psal = NaN(nDep,nProf);
fluo = NaN(nDep,nProf);
dens = NaN(nDep,nProf);
dis  = NaN(1,nProf);
lon  = NaN(1,nProf);
lat  = NaN(1,nProf);

for j = 1:nProf
    pres(:,j) = data(j).pres.value;
    temp(:,j) = data(j).temp.value;
    psal(:,j) = data(j).psal.value;
    fluo(:,j) = data(j).analog.value;
    dens(:,j) = data(j).denSmooth.value+1000;
    lon(j)    = metadata(j).lon.value;
    lat(j)    = metadata(j).lat.value;
end

dis(2:end) = cumsum(m_idist(lon(1:end-1),lat(1:end-1),lon(2:end),lat(2:end)));
dis(1) = 0;
diskm = dis./1000;

N2 = NaN(nDep,nProf);
N2_window = 2;
RiB_window = 1;

g = 9.8; % m/s^2, for buoyancy
rho0 = 1024; % Reference density...
f = 2*7.2921e-5*sind(metadata(1).lat.value); % in degrees

for j = 1:nProf
    ii = find(~isnan(dens(:,j)));
    tempdens = dens(ii,j);
    dep = pres(ii,j);
    for k = 1+N2_window:length(ii)-N2_window
        dendiff = (tempdens(k+N2_window)-tempdens(k-N2_window));
        depdiff = (dep(k+2)-dep(k-2));
        N2(ii(k),j) = g*dendiff/(2*depdiff)/rho0;
    end
end

% plot:
figure(1)
pcolor(diskm,pres(:,1),N2); axis ij; shading interp;
set(gca,'fontsize',20); xlabel('Distance (km) ');
ylabel('Pressure (dBar) ');
axis([0 diskm(end) 0 100]);
h = colorbar;
ylabel(h,'Stratification');
title('Stratification'); n2_caxis = caxis;
hold on
contour(diskm,pres(:,1),dens,10,'w','linewidth',3);
caxis(n2_caxis);
hold off
drawnow
pause

% Richardson number derived from density.

% get horz buoyancy grad...
Bh = NaN(nDep,nProf); % Horizontal buoyancy gradient.
RiB = NaN(nDep,nProf); % Gradient Richardson number
dh = NaN(nDep,nProf); % Horizontal buoyancy difference
for j = 1+RiB_window:nDep-RiB_window
    for k = 2:nProf-1
        dh(j,k) = -g*((dens(j,k+1)-(dens(j,k-1))))/rho0;
        Bh(j,k) = (dh(j,k)/dis(k))^2;
        RiB(j,k) = N2(j,k)*f^2/Bh(j,k);
    end
end

% Remove RiB<=0
RiB(RiB<=0) = NaN;

figure(1)
pcolor(diskm,pres(:,1),log10(RiB)); axis ij; axis([0 diskm(end) 0 100]);
shading interp; set(gca,'fontsize',20); xlabel('Distance (km) ');
ylabel('Pressure (dBar) '); h = colorbar;
ylabel(h,'Richardson Number');
hold on
contour(diskm,pres(:,1),dens,10,'w','linewidth',3);
caxis([0 3]);
hold off
drawnow

% Note: may have to remove high values from RiB to get a good distribution...

%% The rest of this script is to match up ADCP profiles with MVP profiles,
%% and calculate the classical Richard #: N^2 / (du/dz)^2
% 
% % Get richardson thru ADCP!
% % It's N2 (have), and ADCP shear...
% 
% % CD to desktop/mio/rossby
% load ADCP_LDB_After.mat
% 
% % For each MVP cast, get nearest ADCP profile (lot more ADCP than MVP)
% equ_ind = zeros(1,96);
% 
% for j = 1:96
%     templon = MVP_lon(j); templat = MVP_lat(j);
%     tempdsit = m_idist(templon,templat,LDB_After_lon,LDB_After_lat);
%     equ_ind(j) = find(tempdsit==min(tempdsit),1,'first');
% end
% 
% RiADCP = NaN(42,96);
% 
% for j = 1:96
%     uprof = LDB_After_u(:,equ_ind(j));
%     vprof = LDB_After_v(:,equ_ind(j));
%     dep = LDB_After_depth;
%     ii = find(~isnan(uprof));
%     shearu = NaN(42,1);
%     shearv = NaN(42,1);
%     for k = 2:length(ii)-1
%         shearu(ii(k)) = (uprof(ii(k+1))-uprof(ii(k-1)))./(2*8);
%         shearv(ii(k)) = (vprof(ii(k+1))-vprof(ii(k-1)))./(2*8);
%     end
%     
%     shearU = (shearu.^2+shearv.^2);
%     
%     ii = find(~isnan(N2_again(:,j)));
%     N2interp = pchip(MVP_P(ii,1),N2_again(ii,j),-dep(1:42));
%     
%     RiADCP(1:42,j) = N2interp(1:42)./shearU(1:42);
% end
% 
% umat = NaN(42,96);
% vmat = NaN(42,96);
% 
% for j = 1:96
%     u(:,j) = LDB_After_u(1:42,equ_ind(j));
%     v(:,j) = LDB_After_v(1:42,equ_ind(j));
% end

    
    
    

        
        