% this script is to plot the ctd transect
clear

pltadcp = 1; % switch to plot adcp velocity cross-shelf section
rmtide  = 1; % switch to remove tide in the adcp data

ctddir  = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/ctd_processed/';
figdir  = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/figs/ctd/';
fnpre   = 'tn368';
fprint  = 1;

% nfile   = [1:4 6:11 13:16];  % Transect 1, southward
% sname   = {'A5','A6','A7','A8','A9','A10','A11','A12','A13','A14','A15','A16','A17','A18'};
% time{1} = [2019 07 05];
% time{2} = [2019 07 06];
% tranum  = 1;
% caltrpt = 1; % switch to calculate transport
% trptind = [1:14]; % indices of the relatively straight section to compute the transport
%                      % the indices are obtained from the plot of lonl and latl 
% tnumsec = 100;       % number of intervals to separate the transport section into
% pltnh4  = 1; % switch to plot NH4 data
% nh4file = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/nh4/TN368_Transect1_Ammonium_raw_data.xlsx';

% nfile   = [17:29];  % Transect 3, southward CTD only
% sname   = {'A5','A6','A7','A8','A9','A10','A11','A12','A13','A14','A15','A16','A17'};
% time{1} = [2019 07 07];
% tranum  = 3;
% caltrpt = 0; % switch to calculate transport

% nfile   = [41 43:47 49:52];  % Transect 3, southward CTD only
% sname   = {'A5','A6','A7','A8','A9','A10','A11','A12','A13','A14'};
% time{1} = [2019 07 09];
% tranum  = 11;
% caltrpt = 0; % switch to calculate transport

% nfile   = [65:67 69:77 79 80];  % Transect 3, southward CTD only
% sname   = {'A18','A17','A16','A15','A14','A13','A12','A11','A10','A9','A8','A7','A6','A5'};
% time{1} = [2019 07 11];
% time{1} = [2019 07 12];
% tranum  = 26;
% caltrpt = 0; % switch to calculate transport
% pltnh4  = 1; % switch to plot NH4 data
% nh4file = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/nh4/TN368_Transect26_Ammonium_raw_data.xlsx';

% nfile   = [86:90 93:99];  % Transect 3, southward CTD only
% sname   = {'A17','A16','A15','A14','A13','A12','A11','A10','A9','A8','A7','A6'};
% time{1} = [2019 07 14];
% tranum  = 37;
% caltrpt = 0; % switch to calculate transport
% pltnh4  = 0; % switch to plot NH4 data
% nh4file = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/nh4/TN368_Transect26_Ammonium_raw_data.xlsx';

% nfile   = [103:112 122 121 120 119];  % Transect 3, southward CTD only
% sname   = {'A18','A17','A16','A15','A14','A13','A12','A11','A10','A9','A8','A7','A6','A5'};
% time{1} = [2019 07 16];
% time{2} = [2019 07 17];
% tranum  = [44 49];
% caltrpt = 0; % switch to calculate transport
% pltnh4  = 0; % switch to plot NH4 data
% nh4file = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/nh4/TN368_Transect26_Ammonium_raw_data.xlsx';

%nfile   = [125:134 135:140];  % Transect 3, southward CTD only
%sname   = {'A16','A15','A14','A13','A12','A11','A10','A9','A8','A7','A6','A5','L3','L4','L2','L1'};
nfile   = [125:136];  % Transect 3, southward CTD only
sname   = {'A16','A15','A14','A13','A12','A11','A10','A9','A8','A7','A6','A5'};
time{1} = [2019 07 17];
time{2} = [2019 07 18];
tranum  = [52 54];
caltrpt = 0; % switch to calculate transport
pltnh4  = 0; % switch to plot NH4 data
nh4file = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/nh4/TN368_Transect26_Ammonium_raw_data.xlsx';

mxdpth    = 150;
fsalt     = 34.5; % frontal density
fcoriolis = sw_f(40);
xrange    = [39.6 40.7]; % normal range covering A1-A18
% xrange    = [39.6 41.2];  % range covering LTER stations

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

no2offset = -3; % 

% for add adcp data
adcpfile = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/adcp/os75bb_uv_blank8mbin8m.mat';
% read transect start time
addpath ~/research/cruise/2019_Thompson_TN368_0703-0719/data/scripts/
tn368_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)) '_tran' num2str(tranum(1)) '-' 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

dz   = abs(diff(data.idpth));
dz   = dz(1);
% !!! assuming the CTD data is on regular vertical grid!!!

% first two casts of SUNA has no value, remove the data
if tranum  == 1
    ntrt(:,1:2) = nan;
end

% re-order the stations according to the latitude
[lat,norder] = sort(lat);
lon   = lon(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
dlat  = (max(lat)-min(lat))/40;
ilat  = min(lat):dlat:max(lat);
ilon  = interp1(lat,lon,ilat);
itemp = griddata(lat,idpth,temp,ilat,idpth);
isalt = griddata(lat,idpth,salt,ilat,idpth);
ipden = griddata(lat,idpth,pden,ilat,idpth);
iflor = griddata(lat,idpth,flor,ilat,idpth);
iturb = griddata(lat,idpth,turb,ilat,idpth);
ioxyg = griddata(lat,idpth,oxyg,ilat,idpth);
intrt = griddata(lat,idpth,ntrt,ilat,idpth);
icdom = griddata(lat,idpth,cdom,ilat,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 lat0
for ii = 1:length(nfile)
    sdist(ii)  = sw_dist([lat(ii) lat0],[lon(ii) lon0],'km');
    if lat(ii) > lat0
        sdist(ii) = -sdist(ii);
    end
end

for ii = 1:length(ilat)
    idist(ii)  = sw_dist([ilat(ii) lat0],[ilon(ii) lon0],'km');
    if ilat(ii) > lat0
        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) > lat0
        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;
    trz2   = repmat(trz(:),   [1 length(trlon)]);
    if rmtide
        [utide,vtide] = otps_tidecurrent_function(trlon,trlat,trtime);
        tru = tru - repmat(utide(:)',[size(tru,1) 1]);
        trv = trv - repmat(vtide(:)',[size(trv,1) 1]);
    end
end

if pltnh4
    [nh4data,nh4txt,nh4raw] = xlsread(nh4file);
    nh4lon0   = nh4data(:,1);
    nh4lat0   = nh4data(:,2);
    nh4wdph0  = nh4data(:,3);
    nh4cast0  = nh4data(:,4);
    nh4nkno0  = nh4data(:,5);
    nh4depth0 = nh4data(:,6);
    nh4value0 = nh4data(:,7);
    nh4stn0   = nh4txt(2:end,1);
    
    [uqnh4stn,IA,IB] = unique(nh4stn0);
    clear nh4lat nh4lon nh4value
    for in = 1:length(uqnh4stn)
        nh4ind = find(strcmp(nh4stn0,uqnh4stn(in)));
        nh4lat(in) = nh4lat0(nh4ind(1));
        nh4lon(in) = nh4lon0(nh4ind(1));
        tnh4depth  = nh4depth0(nh4ind);
        tnh4value  = nh4value0(nh4ind);
        nh4value(:,in) = interp1(-tnh4depth,tnh4value,idpth);
    end
    [nh4lat,nh4ind2] = sort(nh4lat);
    nh4lon           = nh4lon(nh4ind2);
    nh4value         = nh4value(:,nh4ind2);
    inh4             = griddata(nh4lat,idpth,nh4value,ilat,idpth);
    nh4v             = griddata(nh4lat,idpth,nh4value,lat,idpth);
end

if caltrpt
    % pick out the section of the ship track for transport calculation
    tlonl   = lon(trptind);
    tlatl   = lat(trptind);
    tdpthc  = repmat(idpth(:),[1 length(trptind)]);
    ttemc   = temp(:,trptind);
    tsalc   = salt(:,trptind);
    tdenc   = pden(:,trptind);
    tflorc  = flor(:,trptind);
    tturbc  = turb(:,trptind);
    toxygc  = oxyg(:,trptind);
    tcdomc  = cdom(:,trptind);
    tntrtc  = ntrt(:,trptind);
    % !!! assume the nh4 has the same dimension as the ctd data !!! CHECK!
    tnh4v   = nh4value(:,trptind);  
        
    % to obtain coordinates of the regularly spaced points on
    % the straight transport section
    dtlon   = (tlonl(1)-tlonl(end))./tnumsec;
    dtlat   = (tlatl(1)-tlatl(end))./tnumsec;
    tlon    = tlonl(end):dtlon:tlonl(1);
    tlat    = tlatl(end):dtlat:tlatl(1);
    tlon2   = repmat(tlon(:)',[size(tdpthc,1) 1]);
    tlat2   = repmat(tlat(:)',[size(tdpthc,1) 1]);
    clear tdist
    for ip  = 1:length(tlon)
        % distance to the first point of the straight section
        tdist(ip) = sw_dist([tlat(ip) tlat(1)],[tlon(ip) tlon(1)],'km').*1000;
    end
    % 2D distance and depth arrays for interpolation purpose
    tdist2  = repmat(tdist(:)',[size(tdpthc,1) 1]);
    tz2     = repmat(idpth(:),[1 length(tlon)]);
    % angle of the straight transport line
    ssang   = atan((tlatl(1)-tlatl(end))./(tlonl(1)-tlonl(end)));
    % to map each point on the ship track onto the straight section by looking
    % for a point among a series of densely distributed ones on
    % the straight line that is closest to the orginal track point
    ttlon   = tlonl(end):dtlon./100:tlonl(1);
    ttlat   = tlatl(end):dtlat./100:tlatl(1);
    clear mtlonl mtlatl mtdist
    for ip  = 1:length(trptind)
        ttdist     = sqrt((tlonl(ip)-ttlon).^2 + (tlatl(ip)-ttlat).^2);
        [mindist,minind] = min(ttdist);
        mtlonl(ip) = ttlon(minind);
        mtlatl(ip) = ttlat(minind);
        mtdist(ip) = sw_dist([mtlatl(ip) tlat(1)],[mtlonl(ip) tlon(1)],'km').*1000;
    end
    % remove the points that have equal distance to the first point
    % of the straight line (for interpolation purpose)
    zoind             = find(diff(mtdist)==0);
    mtdist(zoind+1)   = [];
    tlonl(zoind+1)    = [];
    tlatl(zoind+1)    = [];
    tdpthc(:,zoind+1) = [];
    ttemc(:,zoind+1)  = [];
    tsalc(:,zoind+1)  = [];
    tdenc(:,zoind+1)  = [];
    tflorc(:,zoind+1) = [];
    tturbc(:,zoind+1) = [];
    toxygc(:,zoind+1) = [];
    tcdomc(:,zoind+1) = [];
    tntrtc(:,zoind+1) = [];

    % to map each adcp point onto the straight section by looking
    % for a point among a series of densely distributed ones on
    % the straight line that is closest to the orginal track point
    clear mvlonl mvlatl mvdist
    for ip  = 1:length(trind)
        ttdist     = sqrt((trlon(ip)-ttlon).^2 + (trlat(ip)-ttlat).^2);
        [mindist,minind] = min(ttdist);
        mvlonl(ip) = ttlon(minind);
        mvlatl(ip) = ttlat(minind);
        mvdist(ip) = sw_dist([mvlatl(ip) tlat(1)],[mvlonl(ip) tlon(1)],'km').*1000;
    end
    % remove the points that have equal distance to the first point
    % of the straight line (for interpolation purpose)
    zoind             = find(diff(mvdist)==0);
    mvdist(zoind+1)   = [];
    trlon(zoind+1)    = [];
    trlat(zoind+1)    = [];
    trz2(:,zoind+1) = [];
    tru(:,zoind+1)    = [];
    trv(:,zoind+1)    = [];

    
    mtdist2 = repmat(mtdist(:)',[size(tdpthc,1) 1]);
    mvdist2 = repmat(mvdist(:)',[size(trz2,1) 1]);
    % interpolate the data onto the regular distance grid
    Fu      = scatteredInterpolant(mvdist2(:),trz2(:),tru(:));
    trptu   = reshape(Fu(tdist2(:),tz2(:)),size(tdist2));
    Fv      = scatteredInterpolant(mvdist2(:),trz2(:),trv(:));
    trptv   = reshape(Fv(tdist2(:),tz2(:)),size(tdist2));
    Ft      = scatteredInterpolant(mtdist2(:),tdpthc(:),ttemc(:));
    trpttem = reshape(Ft(tdist2(:),tz2(:)),size(tdist2));
    Fs      = scatteredInterpolant(mtdist2(:),tdpthc(:),tsalc(:));
    trptsal = reshape(Fs(tdist2(:),tz2(:)),size(tdist2));
    Fd      = scatteredInterpolant(mtdist2(:),tdpthc(:),tdenc(:));
    trptden = reshape(Fd(tdist2(:),tz2(:)),size(tdist2));
    Ff      = scatteredInterpolant(mtdist2(:),tdpthc(:),tflorc(:));
    trptflo = reshape(Ff(tdist2(:),tz2(:)),size(tdist2));
    Ftb     = scatteredInterpolant(mtdist2(:),tdpthc(:),tturbc(:));
    trpttrb = reshape(Ftb(tdist2(:),tz2(:)),size(tdist2));
    Fo      = scatteredInterpolant(mtdist2(:),tdpthc(:),toxygc(:));
    trptoxy = reshape(Fo(tdist2(:),tz2(:)),size(tdist2));
    Fc      = scatteredInterpolant(mtdist2(:),tdpthc(:),tcdomc(:));
    trptcdm = reshape(Fc(tdist2(:),tz2(:)),size(tdist2));
    Fn      = scatteredInterpolant(mtdist2(:),tdpthc(:),tntrtc(:));
    trptntr = reshape(Fn(tdist2(:),tz2(:)),size(tdist2));
    
    % fill the gap on the surface for temp, salt, dens, fluo, turb and oxyg
    % and set the value below the bottom to 0
    for id = 1:length(tlon)
        if isnan(trpttem(1,id))
            nnnind = find(~isnan(trpttem(:,id)));
            trpttem(1:(nnnind(1)-1),id) = trpttem(nnnind(1),id);
        end
        if isnan(trptsal(1,id))
            nnnind = find(~isnan(trptsal(:,id)));
            trptsal(1:(nnnind(1)-1),id) = trptsal(nnnind(1),id);
        end
        if isnan(trptden(1,id))
            nnnind = find(~isnan(trptden(:,id)));
            trptden(1:(nnnind(1)-1),id) = trptden(nnnind(1),id);
        end
        if isnan(trptflo(1,id))
            nnnind = find(~isnan(trptflo(:,id)));
            trptflo(1:(nnnind(1)-1),id) = trptflo(nnnind(1),id);
        end
        if isnan(trpttrb(1,id))
            nnnind = find(~isnan(trpttrb(:,id)));
            trpttrb(1:(nnnind(1)-1),id) = trpttrb(nnnind(1),id);
        end
        if isnan(trptoxy(1,id))
            nnnind = find(~isnan(trptoxy(:,id)));
            trptoxy(1:(nnnind(1)-1),id) = trptoxy(nnnind(1),id);
        end
        if isnan(trptcdm(1,id))
            nnnind = find(~isnan(trptcdm(:,id)));
            trptcdm(1:(nnnind(1)-1),id) = trptcdm(nnnind(1),id);
        end
        if isnan(trptntr(1,id))
            nnnind = find(~isnan(trptntr(:,id)));
            if ~isempty(nnnind)
                trptntr(1:(nnnind(1)-1),id) = trptntr(nnnind(1),id);
            else
                trptntr(:,id) = 0;
            end
        end
    end
    % fill the gap on the bottom for temp, salt, dens, fluo, turb and oxyg
    trptdpth = griddata(bdata.lon',bdata.lat,bdata.h,tlon,tlat);
    for id = 1:length(tlon)
        if isnan(trpttem(end,id))
            nnnind = find(~isnan(trpttem(:,id)));
            trpttem((nnnind(end)+1:end),id) = trpttem(nnnind(end),id);
        end
        if isnan(trptsal(end,id))
            nnnind = find(~isnan(trptsal(:,id)));
            trptsal((nnnind(end)+1:end),id) = trptsal(nnnind(end),id);
        end
        if isnan(trptden(end,id))
            nnnind = find(~isnan(trptden(:,id)));
            trptden((nnnind(end)+1:end),id) = trptden(nnnind(end),id);
        end
        if isnan(trptflo(end,id))
            nnnind = find(~isnan(trptflo(:,id)));
            trptflo((nnnind(end)+1:end),id) = trptflo(nnnind(end),id);
        end
        if isnan(trpttrb(end,id))
            nnnind = find(~isnan(trpttrb(:,id)));
            trpttrb((nnnind(end)+1:end),id) = trpttrb(nnnind(end),id);
        end
        if isnan(trptoxy(end,id))
            nnnind = find(~isnan(trptoxy(:,id)));
            trptoxy((nnnind(end)+1:end),id) = trptoxy(nnnind(end),id);
        end
        if isnan(trptcdm(end,id))
            nnnind = find(~isnan(trptcdm(:,id)));
            trptcdm((nnnind(end)+1:end),id) = trptcdm(nnnind(end),id);
        end
        if isnan(trptntr(end,id))
            nnnind = find(~isnan(trptntr(:,id)));
            if ~isempty(nnnind)
                trptntr((nnnind(end)+1:end),id) = trptntr(nnnind(end),id);
            else
                trptntr(:,id) = 0;
            end
        end
        if isnan(trptu(end,id))
            nnnind = find(~isnan(trptu(:,id)));
            if ~isempty(nnnind)
                trptu((nnnind(end)+1:end),id) = trptu(nnnind(end),id);
            else
                trptu(:,id) = 0;
            end
        end
        if isnan(trptv(end,id))
            nnnind = find(~isnan(trptv(:,id)));
            if ~isempty(nnnind)
                trptv((nnnind(end)+1:end),id) = trptv(nnnind(end),id);
            else
                trptv(:,id) = 0;
            end
        end
        btmind = find(tz2(:,id) < - trptdpth(id));
        trpttem(btmind,id) = 0;
        trptsal(btmind,id) = 0;
        trptden(btmind,id) = 0;
        trptflo(btmind,id) = 0;
        trpttrb(btmind,id) = 0;
        trptoxy(btmind,id) = 0;
        trptcdm(btmind,id) = 0;
        trptntr(btmind,id) = 0;
        trptu(btmind,id)   = 0;
        trptv(btmind,id)   = 0;
        
    end
    
    % compute the horizontal density gradient and use the thermal
    % wind balance to fill the surface gap in the normal component
    % of the velocity
    drdx            = zeros(size(trptu));
    dtdist2         = diff(tdist2')';
    % smooth the density distribution
    trptden1        = shapiro2(trptden,2,2,50);
    drdx(:,2:end-1) = (trptden1(:,3:end)-trptden1(:,1:end-2)) ./ (dtdist2(:,1:end-1)+dtdist2(:,2:end));
    drdx(:,1)       = drdx(:,2);
    drdx(:,end)     = drdx(:,end-1);
    dnuvdz          = -9.81./fcoriolis./1025.*drdx;
    dnuvdzctr       = av2(dnuvdz);
    % velocity component of the section
    trptnuv         = -trptu.*sin(ssang) + trptv.*cos(ssang);
    ntrptnuv        = trptnuv;
    % fill the gap in the normal component of the velocity
    for id = 1:length(tlon)
        if isnan(ntrptnuv(1,id))
            nnnind = find(~isnan(trptnuv(:,id)));
            addnuv = flipud(cumsum(flipud(dnuvdzctr(1:(nnnind(1)-1),id).*dz)));
            ntrptnuv(1:(nnnind(1)-1),id) = ntrptnuv(nnnind(1),id)+addnuv;
        end
    end
    % for the nans in the middle of the water column, fill with 0
    uvnind = find(isnan(ntrptnuv));
    ntrptnuv(uvnind) = 0;
    
    % compute the volume transport
    shfwind              = find(trptsal <= 34.5);
    slpwind              = find(trptsal >  34.5 & tz2 < 100);
    dtdistctr            = zeros(size(ntrptnuv));
    dtdistctr(:,2:end-1) = (dtdist2(:,1:end-1)+dtdist2(:,2:end))./2;
    dtdistctr(:,1)       = dtdist2(:,1);
    dtdistctr(:,end)     = dtdist2(:,end);
    voltrpt              = sum(ntrptnuv(shfwind).*dz.*dtdistctr(shfwind));
    disp(['The downshelf volume transport is: ' num2str(voltrpt) ' m^3/s']);
    
end


%%%%%%%%%%%%%%%%%%%%%%%
figure(1);clf;
orient tall
colormap('jet');
subplot(311);
pcolorjw(ilat,idpth,itemp);
hold on
[cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
[cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
plot(ilat,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([xrange]);
grid on;
caxis([4 20]);
set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top','xdir','reverse');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lat,zeros(size(lat)),'vk','linewidth',2);
text(40.6,-mxdpth*0.9,'Temperature (^oC)','fontsize',15);
for ii = 1:numstn
    text(lat(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(ilat,idpth,isalt);
hold on
[cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
clabel(cc0,hh0);
[cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
plot(ilat,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([xrange]);
grid on;
caxis([31.5 36.5]);
%caxis([33.4 34.8]);
set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top','xdir','reverse');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lat,zeros(size(lat)),'vk','linewidth',2);
text(40.6,-mxdpth*0.9,'Salinity','fontsize',15);
for ii = 1:numstn
    text(lat(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end

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

% plot chl, turb, and oxy with latitude
figure(2); clf
orient tall
colormap('jet');
subplot(311);
pcolorjw(ilat,idpth,iflor);
hold on
[cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilat,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([xrange]);
grid on;
%caxis([0 1.5]);
caxis([0 3]);
set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top','xdir','reverse');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lat,zeros(size(lat)),'vk','linewidth',2);
text(40.6,-mxdpth*0.9,'Chlorophyll (mg m^-^3)','fontsize',15);
[cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lat(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(ilat,idpth,icdom);
hold on
[cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilat,-ibdpth,'k','linewidth',2);
caxis([0 3]);
ylim([-mxdpth 0]);
xlim([xrange]);
grid on;
set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top','xdir','reverse');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lat,zeros(size(lat)),'vk','linewidth',2);
text(40.6,-mxdpth*0.9,'CDOM (mg/m^3)','fontsize',15);
[cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lat(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end

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

% plot chl, turb, and oxy with latitude
figure(3); clf
orient tall
colormap('jet');
subplot(311);
pcolorjw(ilat,idpth,iflor);
hold on
[cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilat,-ibdpth,'k','linewidth',2);
ylim([-mxdpth 0]);
xlim([xrange]);
grid on;
%caxis([0 1.5]);
caxis([0 8]);
set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top','xdir','reverse');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lat,zeros(size(lat)),'vk','linewidth',2);
text(40.6,-mxdpth*0.9,'Chlorophyll (mg m^-^3)','fontsize',15);
[cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lat(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(ilat,idpth,iturb);
hold on
[cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
set(hh0,'color',[0.5 0.5 0.5]);
plot(ilat,-ibdpth,'k','linewidth',2);
caxis([0 .5]);
ylim([-mxdpth 0]);
xlim([xrange]);
grid on;
set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
set(gca,'xticklabel',[],'layer','top','xdir','reverse');
ylabel('Depth (m)');
pos = get(gca,'pos'); hclb = colorbar;
set(hclb,'location','eastoutside'); set(gca,'pos',pos);
set(hclb,'fontsize',13,'tickdir','out');
plot(lat,zeros(size(lat)),'vk','linewidth',2);
text(40.6,-mxdpth*0.9,'Turbidity (NTU)','fontsize',15);
[cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
set(hh,'color','k','linewidth',2);
for ii = 1:numstn
    text(lat(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
end

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

if pltadcp
    figure(4);clf;
    orient tall
    colormap('jet');
    subplot(311);
    pcolorjw(trlat,trz,tru);
    hold on
    [cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    [cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    plot(ilat,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
    xlim([xrange]);
    grid on;
    caxis([-0.5 0.5]);
    set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'xticklabel',[],'layer','top','xdir','reverse');
    ylabel('Depth (m)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    set(hclb,'fontsize',13,'tickdir','out');
    plot(lat,zeros(size(lat)),'vk','linewidth',2);
    text(40.6,-mxdpth*0.9,'Eastward velocity (m/s)','fontsize',15);
    for ii = 1:numstn
        text(lat(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(trlat,trz,trv);
    hold on
    [cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    clabel(cc0,hh0);
    [cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    plot(ilat,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
    xlim([xrange]);
    grid on;
    caxis([-0.2 0.2]);
    set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'xticklabel',[],'layer','top','xdir','reverse');
    ylabel('Depth (m)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    set(hclb,'fontsize',13,'tickdir','out');
    plot(lat,zeros(size(lat)),'vk','linewidth',2);
    text(40.6,-mxdpth*0.9,'Northward velocity (m/s)','fontsize',15);
    for ii = 1:numstn
        text(lat(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
    end
    
    subplot(313);
    %[cc0,hh0] = contourf(ilat,idpth,ipden-1000,[20:0.1:28]);
    pcolorjw(ilat,idpth,ipden-1000);
    set(hh0,'color',[0.5 0.5 0.5]);
    hold on
    [cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    [cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    caxis([25 27])
    plot(ilat,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
    xlim([xrange]);
    grid on;
    set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'layer','top','xdir','reverse');
    ylabel('Depth (m)');
    xlabel('Latitude (^oN)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    set(hclb,'fontsize',13,'tickdir','out');
    plot(lat,zeros(size(lat)),'vk','linewidth',2);
    text(40.6,-mxdpth*0.9,'\sigma_\theta(kg m^-^3)','fontsize',15);
    for ii = 1:numstn
        text(lat(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
    end
    drawnow;
    if fprint
        print('-dpng','-r300',[ffpre 'uvr_isohal_lat.png']);
    end
end

if pltnh4
    figure(5);clf;
    orient tall
    colormap('jet');
    subplot(311);
    pcolorjw(ilat,idpth,inh4);
    hold on
    [cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    [cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    plot(ilat,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
    xlim([xrange]);
    grid on;
    caxis([0 1000]);
    set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'xticklabel',[],'layer','top','xdir','reverse');
    ylabel('Depth (m)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    cpos = get(hclb,'pos');
    set(hclb,'fontsize',13,'tickdir','out','pos',[cpos(1:2) 0.01 cpos(4)]);
    plot(lat,zeros(size(lat)),'vk','linewidth',2);
    text(40.6,-mxdpth*0.9,'Ammonium (nmol/l)','fontsize',15);
    for ii = 1:numstn
        text(lat(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(ilat,idpth,inh4);
    hold on
    [cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    clabel(cc0,hh0);
    [cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    plot(ilat,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
    xlim([xrange]);
    grid on;
    caxis([0 50]);
    set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'xticklabel',[],'layer','top','xdir','reverse');
    ylabel('Depth (m)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    cpos = get(hclb,'pos');
    set(hclb,'fontsize',13,'tickdir','out','pos',[cpos(1:2) 0.01 cpos(4)]);
    plot(lat,zeros(size(lat)),'vk','linewidth',2);
    text(40.6,-mxdpth*0.9,'Ammonium (nmol/l)','fontsize',15);
    for ii = 1:numstn
        text(lat(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
    end
    
    subplot(313);
    pcolorjw(ilat,idpth,log10(inh4));
    set(hh0,'color',[0.5 0.5 0.5]);
    hold on
    [cc0,hh0] = contour(ilat,idpth,isalt,[33:0.5:36]);
    set(hh0,'color',[0.5 0.5 0.5]);
    [cc,hh] = contour(ilat,idpth,isalt,[fsalt fsalt+0.000001]);
    set(hh,'color','k','linewidth',2);
    caxis([0 3])
    plot(ilat,-ibdpth,'k','linewidth',2);
    ylim([-mxdpth 0]);
    xlim([xrange]);
    grid on;
    set(gca,'fontsize',15,'tickdir','out','xtick',[39.1:0.2:41],'ytick',[-mxdpth:mxdpth/5:0]);
    set(gca,'layer','top','xdir','reverse');
    ylabel('Depth (m)');
    xlabel('Latitude (^oN)');
    pos = get(gca,'pos'); hclb = colorbar;
    set(hclb,'location','eastoutside'); set(gca,'pos',pos);
    cpos = get(hclb,'pos');
    set(hclb,'fontsize',13,'tickdir','out','pos',[cpos(1:2) 0.01 cpos(4)]);
    plot(lat,zeros(size(lat)),'vk','linewidth',2);
    text(40.6,-mxdpth*0.9,'log_1_0(Ammonium (nmol/l))','fontsize',15);
    for ii = 1:numstn
        text(lat(ii)+0.01,mxdpth*0.05,sname{ii},'fontsize',15);
    end
    drawnow;
    if fprint
        print('-dpng','-r300',[ffpre 'nh4_isohal_lat.png']);
    end
    
    
end

