%
% Matlab script to process and plot CTD data from the Andaman Shelf.
%
% ckharris@vims.edu
% January 2018; slightly revised January 2020 for archiving.
%
% Data stored as *mat files within a directory for each station.
%

% These paths are local to CKHarris' computer; they allow access to
% some of the matlab functions used below.
addpath('~ckharris/Documents/MATLAB/TOOLBOX/m_cmg/trunk/RPSstuff');
addpath('~ckharris/Documents/MATLAB/TOOLBOX/CKH_STUFF');
addpath('~ckharris/Documents/MATLAB/TOOLBOX/seawater_ver3_3.1');

% change the mat_file depending on which station you want to process.
num_files=10;
fig_num=0;
all_locations=[];
all_surf_salt=[];
all_bot_salt=[];
all_mn_salt=[];
all_surf_temp=[];
all_bot_temp=[];
all_mn_temp=[];
all_surf_dens=[];
all_bot_dens=[];
all_mn_dens=[];
all_stations=[];
all_stations=[];
all_bot_turb=[];
all_max_turb=[];


for ifile=1:num_files
    switch ifile
      case 1
        mat_file='CTD_DATA/STATION1/ST001_018042_20171206_0024.mat';
      case 2
        mat_file='CTD_DATA/STATION2_3/STATION23_018042_20171207_0522.mat';
        %      case 3
        %        mat_file='CTD_DATA/STATION4_5/018042_20171205_2312.mat'; % data
                                                        % looks
                                                        % wonky
      case 3        
        mat_file='CTD_DATA/STATION6/018042_20171209_0039.mat'; 
        %mat_file='CTD_DATA/STATION8/018042_20171209_2339.mat'; % dataclear 
        
        %identical to station 6
      case 4
        mat_file='CTD_DATA/STATION11/018042_20171210_2310.mat'; 
      case 5
        mat_file='CTD_DATA/STATION22/018042_20171214_2057.mat';
      case 6
        mat_file='CTD_DATA/STATION27/018042_20171216_2150.mat';
      case 7
        mat_file='CTD_DATA/STATION28/018042_20171217_0151.mat'; 
      case 8
        mat_file='CTD_DATA/STATION29/018042_20171217_0459.mat'; 
      case 9
        mat_file='CTD_DATA/STATION30/018042_20171217_0641.mat'; 
      case 10
        mat_file='CTD_DATA/RiverTransect/018042_20171218_0410.mat'; 
    end % of switch to choose file name.
    

load(mat_file)

data=RBR.data;  % this should be a Nsamples X 10 data file.

% columns are: (1)conductivity, (2)temperature, (3)pressure,
% (4)voltage, (5)voltage, (6)depth, (7)salinity, (8)specific
% conductivity, (9)density anomaloy, (10)speed of sound.

dataname=RBR.datasetfilename;  % assume this is unique identifier
depth=data(:,6);         % units in meters
salinity=data(:,7);      % PSU
temp=data(:,2);          % degrees-C
conductivity=data(:,1);  % 
sp_cond=data(:,8);
pressure=data(:,3);
volt1=data(:,4);
volt2=data(:,5);

for jj=1:length(depth)
    density(jj)=sw_dens(salinity(jj),temp(jj),pressure(jj));
    %    dens_turb(jj)=density(jj)+volt2(jj);  % assumes 1 volt ~ 1g/L
end 


% choose the downcast and upcasts for each dataset manually.
switch str2num(dataname(end-7:end-4));
  case 0024
    num_casts=1;
    station=1;
    down=[111,591];
    lat_lon=[15.7717	96.0682];
  case 0522
    num_casts=2;  % stations 2 and 3
    station=[2 3];
    down=[545    1497
          150513 150923];
    lat_lon=[14.9604	96.21695;
             15.02445	96.31321667];
  case 0039
    num_casts=1; 
    station=6;
    down=[551 1028];
    lat_lon=[15.18853333	97.51615];
  case 2339
    num_casts=1; 
    station=8;
    down=[548 1028];
    lat_lon=[15.02513333	96.31085];
  case 2310
    num_casts=1; 
    station=11;
    down=[20804 22028]
    lat_lon=[15.21038333	95.57235];
  case 2057
    num_casts=1; 
    station=22;
    down=[1074 1804];
    lat_lon=[15.61483333	94.39885];
  case 2150
    num_casts=1;
    station=27;
    down=[5170,5638];
    lat_lon=[15.39673333	96.03738333];
  case 0151
    num_casts=1;
    station=28;
    down=[2750 3005];
    lat_lon=[15.65791667	96.46316667];
  case 0459
    num_casts=1;
    station=29;
    down=[7365 7650];
    lat_lon=[15.84901667	96.75523333];
  case 0641
    num_casts=1;
    station=[30];
    down=[2343 2790];
    lat_lon=[15.37292	96.85052833]
  case 0410
    num_casts=5;
    station=[9.1 9.2 9.3 9.4 9.5]; % river stations
    down=[ 4205        4411
        5187        5432
        6007        6230
        6561        6762
        7048        7223 ];
    lat_lon=[16.45368333	96.34128333
             16.53016667	96.25828333
             16.62795	        96.25681667
             16.71786667	96.225
             16.76451167	96.18057];
end

% plot the downcast
for i_cast=1:num_casts
    fig_num=fig_num+1
    figure(fig_num)
    dn_cast=down(i_cast,1):down(i_cast,2);
clf
sethelv
subplot(141)
plot(salinity(dn_cast),-depth(dn_cast),'linewidth',2);
title(['Station: ',num2str(station(i_cast)),' File: ',dataname(1:6)])
ylabel('Depth m')
xlabel('Salinity')

subplot(142)
plot(temp(dn_cast),-depth(dn_cast),'linewidth',2);
xlabel('Temperature C')

subplot(143)
plot(volt2(dn_cast),-depth(dn_cast),'linewidth',2);
hold on; 
%plot(volt1(dn_cast),-depth(dn_cast),'linewidth',2,'color','r');
xlabel('Voltage')

subplot(144)
plot(density(dn_cast),-depth(dn_cast),'linewidth',2);
%hold on;
%plot(dens_turb(dn_cast),-depth(dn_cast),'linewidth',2,'color','r');
xlabel('Density kg/m^3')

wysiwyg

all_surf_salt=[all_surf_salt; salinity(dn_cast(1))];
all_bot_salt=[all_bot_salt; salinity(dn_cast(end))];
all_mn_salt=[all_mn_salt; mean(salinity(dn_cast))];

all_surf_temp=[all_surf_temp; temp(dn_cast(1))];
all_bot_temp=[all_bot_temp; temp(dn_cast(end))];
all_mn_temp=[all_mn_temp; mean(temp(dn_cast))];

all_surf_dens=[all_surf_dens; density(dn_cast(1))];
all_bot_dens=[all_bot_dens; density(dn_cast(end))];
all_mn_dens=[all_mn_dens; mean(density(dn_cast))];

all_bot_turb=[all_bot_turb; volt2(dn_cast(end))];
all_max_turb=[all_max_turb; max(volt2(dn_cast))];

end 

all_stations=[all_stations; station'];
all_locations=[all_locations; lat_lon];
% need to clear out variables.
if (length(ifile)>1)
    keep ifile num_files fig_num all_locations all_surf_salt all_bot_salt ...
    all_stations all_mn_salt all_surf_temp all_bot_temp all_mn_temp ...
    all_surf_dens all_bot_dens all_mn_dens all_bot_turb all_max_turb
end

end % of loop through the files

% make a location map showing stratification
figure
clf
sethelv    % simple function to set fonts in plot
wysiwyg    % function to set plotting window size
load ayerwaddy_bathy.mat   % bathymetry file developed based on
                           % ETOPO data

%make land nan
zbathy_orig=zbathy;
land=find(zbathy<0);
zbathy(land)=NaN;

norm_strat=(all_bot_salt>all_surf_salt);
norm_strat(1)=0;  % this has a mid-water-column peak

volt_thresh=1;
high_turb=find(all_bot_turb>volt_thresh);

pcolorjw(xbathy,ybathy,log10(zbathy+eps));  
set(gca,'ylim',[14.5,18])
set(gca,'xlim',[93.5,98]);
dasp(16);
caxis(log10([1 1000]));
hold on
hturb=plot(all_locations(high_turb,2),all_locations(high_turb,1),'ko')
set(hturb,'markersize',20)
set(hturb,'linewidth',2.5);
hnorm=plot(all_locations(find(norm_strat),2),all_locations(find(norm_strat),1),'ks');
set(hnorm,'markersize',15);
set(hnorm,'markerfacecolor','b')
hinv=plot(all_locations(find(~norm_strat),2),all_locations(find(~norm_strat),1),'kh');
set(hinv,'markersize',15);
set(hinv,'markerfacecolor','r')
[cc,hh]=contour(xbathy,ybathy,zbathy_orig,[-1:1:1],'k');
set(hh,'linewidth',2);
[cc,hh]=contour(xbathy,ybathy,zbathy_orig,[10:10:100 1000],'k');
set(hh,'linewidth',1)
title('CTD Locations')
text(94.57,17.75,'Color: bathymetry')
legend([hnorm hinv hturb],'Normal stratification','Fresher at depth?',...
                    'High turbidty')
htxt=text(all_locations(1:end-5,2)+0.1,all_locations(1:end-5,1),num2str(all_stations(1:end-5)))


hclb=colorbar;
set(hclb,'Yticklabel',['1m   ';'3m   ';'10m  ';'30m  ';'100m ';...
                    '300m ';'1000m']);



% make a location map showing density stratification
figure
clf
sethelv
wysiwyg

norm_strat=(all_bot_dens>all_surf_dens);
dens_diff=all_bot_dens - all_surf_dens;
% the density difference for CTD Cast 1 is actually indicating less
% dense water under more dense water because there is a mid-water
% column density maximum that does not show up from bottom -
% surface. 
dens_diff(1)=-.7;
%norm_strat(1)=0;  % this has a mid-water-column peak

pcolorjw(xbathy,ybathy,log10(zbathy+eps));  
set(gca,'ylim',[14.5,18])
set(gca,'xlim',[93.5,98]);
dasp(16);
caxis(log10([1 1000]));

hold on
[cc,hh]=contour(xbathy,ybathy,zbathy_orig,[-1 0 1],'k');
set(hh,'linewidth',2)
[cc,hh]=contour(xbathy,ybathy,zbathy_orig,[10:10:100 1000],'k');
set(hh,'linewidth',1)

if (0)
hnorm=plot(all_locations(find(norm_strat),2),all_locations(find(norm_strat),1),'ko');
set(hnorm,'markersize',10);
set(hnorm,'markerfacecolor','b')
hinv=plot(all_locations(find(~norm_strat),2),all_locations(find(~norm_strat),1),'ks');
set(hinv,'markersize',15);
set(hinv,'markerfacecolor','r')
else
    thresh=0.4;
    inv_strat=find(dens_diff<0);
    well_mixed=find((dens_diff>=0)&(dens_diff<=thresh))
    norm_strat=find(dens_diff>thresh);
    %    [hdt]=scatter(all_locations(:,2),all_locations(:,1),...
    %          200,dens_diff,'filled');
    hturb=plot(all_locations(high_turb,2),all_locations(high_turb,1),'ko')
    set(hturb,'markersize',20)
    set(hturb,'linewidth',2.5);
    hinv=plot(all_locations(inv_strat,2),all_locations(inv_strat,1),'kv')
    set(hinv,'markersize',15);
    set(hinv,'markerfacecolor','r')
    hmix=plot(all_locations(well_mixed,2),all_locations(well_mixed,1),'ko')
    set(hmix,'markersize',15);
    set(hmix,'markerfacecolor','m')
    hnorm=plot(all_locations(norm_strat,2),all_locations(norm_strat,1),'^')
    set(hnorm,'markersize',15);
    set(hnorm,'markerfacecolor','b')
    legend([hinv hmix hnorm hturb],'Less dense at depth?','Well mixed?',...
                        'Denser at depth','High Turbidity')
end

title('CTD Locations')
text(94.57,17.75,'Color: bathymetry')
htxt=text(all_locations(1:end-5,2)+0.1,all_locations(1:end-5,1),num2str(all_stations(1:end-5)))


hclb=colorbar;
set(hclb,'Yticklabel',['1m   ';'3m   ';'10m  ';'30m  ';'100m ';...
                    '300m ';'1000m']);



% make a location map just showing location
figure
clf
sethelv
wysiwyg

pcolorjw(xbathy,ybathy,log10(zbathy+eps));  
set(gca,'ylim',[14.5,18])
set(gca,'xlim',[93.5,98]);
dasp(16);
caxis(log10([1 1000]));

hold on
[cc,hh]=contour(xbathy,ybathy,zbathy_orig,[-1 0 1],'k');
set(hh,'linewidth',2)
[cc,hh]=contour(xbathy,ybathy,zbathy_orig,[10:10:100 1000],'k');
set(hh,'linewidth',1)

hsta=plot(all_locations(:,2),all_locations(:,1),'ko');
set(hsta,'markersize',10);
set(hsta,'markerfacecolor','y')


title('CTD Locations')
text(94.57,17.75,'Color: bathymetry')
htxt=text(all_locations(1:end-5,2)+0.1,all_locations(1:end-5,1),num2str(all_stations(1:end-5)))


hclb=colorbar;
set(hclb,'Yticklabel',['1m   ';'3m   ';'10m  ';'30m  ';'100m ';...
                    '300m ';'1000m']);

