% this script is to plot the ctd transect
clear

pltadcp = 1; % switch to plot adcp velocity cross-shelf section

ctddir  = '~/research/cruise/2019_Brown_0512-0525/data/CTD_processed/';
figdir  = '~/research/cruise/2019_Brown_0512-0525/data/figs/ctd/';
fnpre   = 'RB1904';
fprint  = 1;

nfile  = [58:61 64:66];  % Transect 18, eastward; 
sname  = {'EW7','EW8','EW9','NS10','A12','EW10','NS4'};
time{1}= [2019 05 17];
time{2}= [2019 05 18];
tranum = [18 19];

mxdpth   = 150;
fsalt    = 34.5; % frontal density

dpthfile = '~/research/roms/bathymetry/ngdc2.mat';

% for add adcp data
adcpfile = '~/research/cruise/2019_Brown_0512-0525/data/adcp/os75bb_uv_blank8mbin4m_16to18.mat';
% read transect start time
addpath ~/research/cruise/2019_Brown_0512-0525/data/scripts/
rb1904_transect_startime


lineclr  = ['r','b'];
numstn   = length(nfile);
if length(time) == 1
    time{2} = time{1};
end
ffpre = [figdir 'ctd' num2str(nfile(1)) '-' num2str(nfile(numstn(1))) '_tran' num2str(tranum(end)) '_'];

for ifile = 1:length(nfile)
    if nfile(ifile) < 10
        filename = [ctddir fnpre '00' num2str(nfile(ifile)) '.mat'];
    elseif nfile(ifile) < 100
        filename = [ctddir fnpre '0' num2str(nfile(ifile)) '.mat'];
    else
        filename = [ctddir fnpre num2str(nfile(ifile)) '.mat'];
    end
    
    data = load(filename);
    
    if 0
        % plot the profiles
        figure(3)
        plot(data.iflor,data.idpth,'color',lineclr(ifile),'linewidth',2);hold on;
        if length(time) > 1
            httl = title([datestr(datenum(time{1}),'yyyy-mm-dd') ' - ' datestr(datenum(time{2}),'mm-dd')]);
        else
            httl = title(datestr(datenum(time{1}),'yyyy-mm-dd'));
        end
        ylabel('Depth (m)');
        xlabel('Chlorophyll (mg m^-^3)');
        grid on;
        
        
        figure(4)
        plot(data.iturb,data.idpth,'color',lineclr(ifile),'linewidth',2);hold on;
        if length(time) > 1
            httl = title([datestr(datenum(time{1}),'yyyy-mm-dd') ' - ' datestr(datenum(time{2}),'mm-dd')]);
        else
            httl = title(datestr(datenum(time{1}),'yyyy-mm-dd'));
        end
        ylabel('Depth(m)');
        xlabel('Turbidity (NTU)');
        grid on;
        
    end
    
    if ifile == 1
        idpth = data.idpth(:);
        lon   = data.ilon;
        lat   = data.ilat;
        wdpth = data.iwdpth;
        temp  = data.itemp(:);
        salt  = data.isalt(:);
        pden  = data.ipden(:);
        flor  = data.iflor(:);
        turb  = data.iturb(:);
        oxyg  = data.ioxyg(:);
        ntrt  = data.intrt(:);
        cdom  = data.icdom(:);
    else
        lon   = [lon data.ilon];
        lat   = [lat data.ilat];
        wdpth = [wdpth data.iwdpth];
        temp  = [temp data.itemp(:)];
        salt  = [salt data.isalt(:)];
        pden  = [pden data.ipden(:)];
        flor  = [flor data.iflor(:)];
        turb  = [turb data.iturb(:)];
        oxyg  = [oxyg data.ioxyg(:)];
        ntrt  = [ntrt data.intrt(:)];
        cdom  = [cdom data.icdom(:)];
    end
end

% re-order the stations according to the longitude
[lon,norder] = sort(lon);
lat   = lat(norder);
wdpth = wdpth(norder);
temp  = temp(:,norder);
salt  = salt(:,norder);
pden  = pden(:,norder);
flor  = flor(:,norder);
turb  = turb(:,norder);
oxyg  = oxyg(:,norder);
ntrt  = ntrt(:,norder);
cdom  = cdom(:,norder);
sname = sname(norder);

% interpolate the data onto fine grid
dlon  = (max(lon)-min(lon))/40;
ilon  = max(lon):-dlon:min(lon);
ilat  = interp1(lon,lat,ilon);
itemp = griddata(lon,idpth,temp,ilon,idpth);
isalt = griddata(lon,idpth,salt,ilon,idpth);
ipden = griddata(lon,idpth,pden,ilon,idpth);
iflor = griddata(lon,idpth,flor,ilon,idpth);
iturb = griddata(lon,idpth,turb,ilon,idpth);
ioxyg = griddata(lon,idpth,oxyg,ilon,idpth);
intrt = griddata(lon,idpth,ntrt,ilon,idpth);
icdom = griddata(lon,idpth,cdom,ilon,idpth);

bdata  = load(dpthfile);
ibdpth = griddata(bdata.lon',bdata.lat,bdata.h,ilon(:),ilat(:));
%if ilat(1)<ilat(end)
%    % if shoreward transect, flip the depth and use the onshore end as the
%    % start
%    %ibdpth = flipud(ibdpth(:));
%    lat0  = ilat(end);
%    lon0  = ilon(end);
%else
%    lat0  = ilat(1);
%    lon0  = ilon(1);
%end
% A1 station
lat0 = 40.725;
lon0 = -70.830;
% distance from the lon0
for ii = 1:length(nfile)
    sdist(ii)  = sw_dist([lat(ii) lat0],[lon(ii) lon0],'km');
    if lon(ii) > lon0
        sdist(ii) = -sdist(ii);
    end
end

for ii = 1:length(ilon)
    idist(ii)  = sw_dist([ilon(ii) lat0],[ilon(ii) lon0],'km');
    if ilon(ii) > lon0
        idist(ii) = -idist(ii);
    end
end
for ii = 1:length(lat)
    dist(ii)  = sw_dist([lat(ii) lat0],[lon(ii) lon0],'km');
    if lat(ii) > lon0
        dist(ii) = -dist(ii);
    end
    
end

if pltadcp
    adcp   = load(adcpfile);
    trind  = find(adcp.adcptime > tstime(tranum(1)) & adcp.adcptime < tstime(tranum(end)+1));
    tru    = adcp.u(:,trind);
    trv    = adcp.v(:,trind);
    trlon  = adcp.adcplon(trind);
    trlat  = adcp.adcplat(trind);
    trtime = adcp.adcptime(trind);
    trz    = -adcp.zc;
end

%%%%%%%%%%%%%%%%%%%%%%%
figure(1);clf;
orient tall
colormap('jet');
subplot(311);
pcolorjw(ilon,idpth,itemp);
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
plot(ilon,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
caxis([4 20]);
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'Temperature (^oC)','fontsize',15);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end
httl = title(['Transect ' num2str(tranum(1)) '-' num2str(tranum(end)) ' -- ' datestr(datenum(time{1}),'yyyy-mm-dd') ' - ' datestr(datenum(time{2}),'mm-dd')]);
set(httl,'pos',[40.16 15]);

subplot(312);
pcolorjw(ilon,idpth,isalt);
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
clabel(cc0,hh0);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
plot(ilon,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
caxis([31.5 36.5]);
%caxis([33.4 34.8]);
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'Salinity','fontsize',15);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end

subplot(313);
%[cc0,hh0] = contourf(ilon,idpth,ipden-1000,[20:0.1:28]);
pcolorjw(ilon,idpth,ipden-1000);
set(hh0,'color',[0.5 0.5 0.5]);
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
caxis([25 27])
plot(ilon,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'layer','top');
ylabel('Depth (m)');
xlabel('Longitude (^oW)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'\sigma_\theta(kg m^-^3)','fontsize',15);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end
drawnow;
if fprint
    print('-dpng','-r300',[ffpre 'tsr_isohal_lon.png']);
end

% plot chl, turb, and oxy with lonitude
figure(2); clf
orient tall
colormap('jet');
subplot(311);
pcolorjw(ilon,idpth,iflor);
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilon,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
%caxis([0 1.5]);
caxis([0 3]);
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'Chlorophyll (mg m^-^3)','fontsize',15);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end
httl = title(['Transect ' num2str(tranum(1)) '-' num2str(tranum(end)) ' -- ' datestr(datenum(time{1}),'yyyy-mm-dd') ' - ' datestr(datenum(time{2}),'mm-dd')]);
set(httl,'pos',[40.16 15]);

subplot(312);
pcolorjw(ilon,idpth,icdom);
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilon,-ibdpth,'k','linewidth',2);
caxis([0 3]);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'CDOM (mg/m^3)','fontsize',15);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end

subplot(313);
pcolorjw(ilon,idpth,ioxyg); % change the unit from ml/l to mu mol/l
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilon,-ibdpth,'k','linewidth',2);
caxis([200 350]);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'layer','top');
ylabel('Depth (m)');
xlabel('Longitude (^oW)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'Oxygen (\mu mol/l)','fontsize',15);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end
drawnow;
if fprint
    print('-dpng','-r300',[ffpre 'cco_isohal_lon.png']);
end

% plot chl, turb, and oxy with lonitude
figure(3); clf
orient tall
colormap('jet');
subplot(311);
pcolorjw(ilon,idpth,iflor);
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilon,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
%caxis([0 1.5]);
caxis([0 3]);
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'Chlorophyll (mg m^-^3)','fontsize',15);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end
httl = title(['Transect ' num2str(tranum(1)) '-' num2str(tranum(end)) ' -- ' datestr(datenum(time{1}),'yyyy-mm-dd') ' - ' datestr(datenum(time{2}),'mm-dd')]);
set(httl,'pos',[40.16 15]);

subplot(312);
pcolorjw(ilon,idpth,iturb);
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilon,-ibdpth,'k','linewidth',2);
caxis([0 .5]);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'Turbidity (NTU)','fontsize',15);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end

subplot(313);
pcolorjw(ilon,idpth,intrt);
hold on
[cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilon,-ibdpth,'k','linewidth',2);
%caxis([0 5]);
caxis([-2 5]);
ylim([-mxdpth 0]);
xlim([min(ilon)-.05 max(ilon)+.05]); 
grid on;
set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'layer','top');
ylabel('Depth (m)');
xlabel('Longitude (^oW)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lon,zeros(size(lon)),'vk','linewidth',2);
text(-71.32,-mxdpth*0.9,'Nitrate (\mu mol/l)','fontsize',15);
[cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end
drawnow;
if fprint
    print('-dpng','-r300',[ffpre 'ctn_isohal_lon_2.png']);
end

if pltadcp
    figure(4);clf;
    orient tall
    colormap('jet');
    subplot(311);
    pcolorjw(trlon,trz,tru);
    hold on
    [cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    [cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    plot(ilon,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
     xlim([min(ilon)-.05 max(ilon)+.05]); 
    grid on;
    caxis([-0.5 0.5]);
    set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'xticklabel',[],'layer','top');
    ylabel('Depth (m)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    set(hclb,'fontsize',13,'tickdir','out');
    plot(lon,zeros(size(lon)),'vk','linewidth',2);
    text(-71.32,-mxdpth*0.9,'Eastward velocity (m/s)','fontsize',15);
    for ii = 1:numstn
        text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
    end
    httl = title(['Transect ' num2str(tranum(1)) '-' num2str(tranum(end)) ' -- ' datestr(datenum(time{1}),'yyyy-mm-dd') ' - ' datestr(datenum(time{2}),'mm-dd')]);
    set(httl,'pos',[40.16 15]);
    
    subplot(312);
    pcolorjw(trlon,trz,trv);
    hold on
    [cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    clabel(cc0,hh0);
    [cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    plot(ilon,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
     xlim([min(ilon)-.05 max(ilon)+.05]); 
    grid on;
    caxis([-0.2 0.2]);
    set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'xticklabel',[],'layer','top');
    ylabel('Depth (m)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    set(hclb,'fontsize',13,'tickdir','out');
    plot(lon,zeros(size(lon)),'vk','linewidth',2);
    text(-71.32,-mxdpth*0.9,'Northward velocity (m/s)','fontsize',15);
    for ii = 1:numstn
        text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
    end
    
    subplot(313);
    %[cc0,hh0] = contourf(ilon,idpth,ipden-1000,[20:0.1:28]);
    pcolorjw(ilon,idpth,ipden-1000);
    set(hh0,'color',[0.5 0.5 0.5]);
    hold on
    [cc0,hh0] = contour(ilon,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    [cc,hh] = contour(ilon,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    caxis([25 27])
    plot(ilon,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
     xlim([min(ilon)-.05 max(ilon)+.05]); 
    grid on;
    set(gca,'fontsize',15,'tickdir','out','xtick',[-71.34:0.05:-71.08],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'layer','top');
    ylabel('Depth (m)');
    xlabel('Longitude (^oW)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    set(hclb,'fontsize',13,'tickdir','out');
    plot(lon,zeros(size(lon)),'vk','linewidth',2);
    text(-71.32,-mxdpth*0.9,'\sigma_\theta(kg m^-^3)','fontsize',15);
    for ii = 1:numstn
        text(lon(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
    end
    drawnow;
    if fprint
        print('-dpng','-r300',[ffpre 'uvr_isohal_lon.png']);
    end
    
end

