clear

fdir   = '~/research/cruise/2019_Brown_0512-0525/data/vpr/';
%flist  = {'05141745.y19'};
%trlist = {5};
%tmlist = [2019 05 14];
%flist  = {'05161836.y19'};
%trlist = {13};
%tmlist = [2019 05 16];
flist  = {'05191512.y19'};
trlist = {26:30};
tmlist = [2019 05 16];

fnpre     = 'rb1904';
algno2    = 1; % flag to align O2 data with the rest of the variables
minoffset = 30;  % the maximum offset for the O2 alignment for the downcast
maxoffset = 30;  % the maximum offset for the O2 alignment for the downcast

figdir = ['~/research/cruise/2019_Brown_0512-0525/data/figs/vpr/aligned_' num2str(minoffset) '-' num2str(maxoffset) 's/'];
%figdir = ['~/research/cruise/2019_Brown_0512-0525/data/figs/vpr/aligned_test/'];

mxdpth  = 150; % maximum depth
fdens   = 1026.5; % frontal density
fdens2  = 1026.2; % frontal density
fsalt   = 34.5;
Nd      = 200; % number of points used for depth resolution

sname = {'A2','A3','A4','A5','A6','A7','A8','A9','A10','A11','A12','A13','A14','A15','A16','A17'};
slon  = [-70.830 -70.830 -70.830 -70.830 -70.829 -70.830 -70.829 -70.828 -70.829 -70.828 -70.826 -70.830 -70.830 -70.830 -70.830 -70.830];
slat  = [40.664 40.598 40.532 40.466  40.400  40.336  40.268  40.203  40.140  40.073  40.010  39.943  39.880  39.815  39.750  39.686];

dpthfile = '~/research/roms/bathymetry/ngdc2.mat';
bdata    = load(dpthfile);
lonlim   = [-71.1 -70.8];
latlim   = [39.6 40.6];
blonind  = find(bdata.lon >= lonlim(1) & bdata.lon <= lonlim(2));
blatind  = find(bdata.lat >= latlim(1) & bdata.lat <= latlim(2));
blon     = bdata.lon(blonind);
blat     = bdata.lat(blatind);
bathy    = -bdata.h(blatind(1):blatind(end),blonind(1):blonind(end));
blon2    = repmat(blon(:)',[length(blat) 1]);
blat2    = repmat(blat(:), [1 length(blon)]);

% A1 station
lat0 = 40.725;
lon0 = -70.830;

% distance of the stations from the lat0
for ii = 1:length(slon)
    stdist(ii)  = sw_dist([slat(ii) lat0],[slon(ii) lon0],'km');
    if slat(ii) > lat0
        stdist(ii) = -stdist(ii);
    end
end

for iv = 1:length(flist)
    dfile  = [fdir flist{iv}];
    data   = load(dfile);
    ttime  = tmlist(iv,:);
    tranum = trlist{iv};
    ffpre  = [figdir 'vpr' num2str(iv) '_tran' num2str(tranum(1)) '-' num2str(tranum(end)) '_'];
    
    % start to read VPR data and sort the data in time ascending order.
    vprtimestep = data(:,1);
    [indi,indj] = sort(vprtimestep,'ascend');
    gmt         = data(indj,2);
    date        = data(indj,3);
    latitude    = data(indj,4); % format is ddmm.mmmm
    longitude   = data(indj,5); % format is ddmm.mmmm
    ddlat       =  fix(latitude/100);
    latitude    = ddlat + (latitude-100*ddlat)/60;
    ddlon       =  fix(longitude/100);
    longitude   = ddlon + (longitude-100*ddlon)/60;
    maxdepth    = data(indj,6);
    turb        = data(indj,11);% unit is ug/L, identical to mg/m^3
    fluo        = data(indj,10); % NOT SURE if 10 and 11 are flipped, yes, should be 10!
    par         = data(indj,12);
    oxyg        = data(indj,13);
    salt        = data(indj,14);
    temp        = data(indj,15);
    depth       = data(indj,16);
    
    % transfer the time into seconds from April 16, 2018, 00:00 am
    %hr = fix(gmt/10000);
    %mm = fix((gmt - 10000*hr)/100);
    %ss = gmt-10000*hr - 100*mm;
    %time = floor(date/10000-16)*3600*24 + hr*3600 + mm*60 + ss;
    
    hr   = floor(gmt./10000);
    mm   = floor((gmt - 10000.*hr)./100);
    ss   = gmt - 10000.*hr - 100.*mm;
    tday = floor(date./10000);
    tmmm = floor((date - 10000.*tday)./100);
    yyyy = 2000 + date - 10000.*tday - 100.*tmmm;
    time = datenum([yyyy(:) tmmm(:) tday(:) hr(:) mm(:) ss(:)]);
    
    fluo        = data(indj,10); % NOT SURE if 10 and 11 are flipped, yes, should be 10!
    par         = data(indj,12);
    oxyg        = data(indj,13);
    salt        = data(indj,14);
    temp        = data(indj,15);
    
    % interpolate the data onto regular time grid
    itime = min(time):1/86400:max(time);
    idepth= interp1(time,depth,itime);
    itemp = interp1(time,temp,itime);
    isalt = interp1(time,salt,itime);
    ioxyg = interp1(time,oxyg,itime);
    ipar  = interp1(time,par, itime);
    ifluo = interp1(time,fluo,itime);
    iturb = interp1(time,turb,itime);
    ilat  = interp1(time,latitude,itime);
    ilon  = interp1(time,longitude,itime);
    
    time  = itime;
    depth = idepth;
    temp  = itemp;
    salt  = isalt;
    oxyg  = ioxyg;
    par   = ipar;
    fluo  = ifluo;
    turb  = iturb;
    latitude  = ilat;
    longitude = ilon;
    dens = sw_pden(salt,temp,abs(depth),1);
    
    if 0
        % scatter plot to check
        scatter(time,-depth,[],dens-1000)
        caxis([26 27])
        colormap(jet)
        colorbar
    end
    
    % transfer lat lon to earth sytem x, y and z
    altitude   = zeros(size(latitude));
    [xs,ys,zs] = lla2ecef(latitude,longitude,altitude);
    % transfer x,y with respective to cross shelf and along shelf
    phi        = pi/6;
    % along-shelf velocity (northeastward poistive)
    % x         = xs.*sin(phi) + ys.* cos(phi);
    x          = xs;
    % cross-shelf velocity (onshore/northwestward positive)
    % y         = -xs.* cos(phi) + ys.* sin(phi);
    y          = ys;
    % Now y is cross shelf direction!!!
    
    % Get the local maxima
    % downcast is clean
    [indid,indjd] = findpeaks(depth,'MinPeakProminence',20);
    % upcast at the top is noisy
    [indiu,indju] = findpeaks(-depth,'MinPeakProminence',5);
    
    % Now the data has been splitted into downcast and upcast
    % for each maxdepth, find the nearest top point
    timeup = time(indju);
    clear indi indj
    upup = []; updn = [];
    for i = 1:length(indjd)
        tmark   = time(indjd(i));
        updn(i) = indjd(i);
        % now find the nearest time at surface
        % WGZ: WHY DO YOU DO THIS??? Why can't you just use the next indju?
        [indi,indj] = min(abs(timeup-tmark));
        time(indju(indj));
        % index = indju(indj) in original time index
        if (time(indju(indj))<tmark && indj<length(timeup)) % up is back, need to advance
            upup = [upup indju(indj+1)];
        elseif(time(indju(indj))<tmark && indj>=length(timeup)) % there is no up cast then
            upup = [upup 0];
        else
            upup = [upup indju(indj)];
        end
    end
    if(upup(end)==0)
        upup = upup(1:end-1);
        updn = updn(1:end-1);
    end
    if upup(end) < length(depth)-200
        % if the top of the last upcast is substantially away for the last
        % sample, the recovery upcast is likely omitted. So add the
        % recovery upcast to the end (remove the last 60 points);
        upup = [upup length(depth)-60];
        updn = [updn indjd(end)];
    end
    
    % What about the downcast
    clear indi indj
    downup = []; downdn = [];
    for i = 1:length(indjd)
        tmark = time(indjd(i));
        downdn(i) = indjd(i);
        % now find the nearest time at surface
        [indi,indj] = min(abs(timeup-tmark));
        % index = indju(indj) in original time index
        if (time(indju(indj))>tmark && indj>1) % up is back, need to advance
            downup = [downup indju(indj-1)];
        elseif(time(indju(indj))>tmark && indj<=1) % there is no up cast then
            downup = [downup 0];
        else
            downup = [downup indju(indj)];
        end
    end
    if(downup(1)==0)
        downup = downup(2:end);
        downdn = downdn(2:end);
    end
    %     if downup(end) > 100
    %         % if the top of the first downcast is substantially away for the
    %         % first sample, the launch downcast is likely omitted. So add the
    %         % launch downcast to the end
    %         downup = [1 downup];
    %         downdn = [indjd(1) downdn];
    %     end
    
    
    if 0
        % check the pick of maximum depth
        figure('Position',[100 100 650 350]);
        subplot('Position',[0.1 0.1 0.85 0.85]);
        scatter(time, -depth,[],temp);
        hold on
        plot(time(indju),-depth(indju),'b-^','markerSize',7);
        hold on
        plot(time(indjd),-depth(indjd),'g-^','markerSize',7);
        % Now the downcast is downup(i) to downdn(i)
        % Now the upcast is updn(i) to upup(i)
        hold on
        for i = 1:length(updn)
            hold on
            plot(time(updn(i):upup(i)),-depth(updn(i):upup(i)),'k-','lineWidth',2);
        end
        for i = 1:length(downdn)
            hold on
            plot(time(downup(i):downdn(i)),-depth(downup(i):downdn(i)),'r-','lineWidth',2);
        end
        hold off
        ylabel('$t$(s)','interpreter','latex','fontsize',15);
        ylabel('$h$(m)','interpreter','latex','fontsize',15);
        title('Data following the time and trajectory of RUV');
    end
    
    ym=[];xm=[];lonm=[];latm=[];timup=[];
    % Let's first try to obtain from the upcast
    yup    = zeros(Nd,length(updn));
    dpthup = yup;
    salup  = yup;
    temup  = yup;
    denup  = yup;
    turbup = yup;
    parup  = yup;
    florup = yup;
    oxygup = yup;
    uoffset= [];
    for i = 1:length(updn)
        ind = updn(i):upup(i);
        dd  = depth(ind);
        [indx,indy] = unique(dd);
        dd  = indx;
        
        tt   = time(ind(indy));
        latt = latitude(ind(indy));
        lonn = longitude(ind(indy));
        xx   = x(ind(indy));
        yy   = y(ind(indy));
        ssa   = salt(ind(indy));
        tem   = temp(ind(indy));
        den   = dens(ind(indy));
        turbw = turb(ind(indy));
        parr  = par(ind(indy));
        flo   = fluo(ind(indy));
        %oxy   = oxyg(ind(indy));
        
        % notice that time is not equal dist, so need interpolation
        ym(i)    = interp1(tt,yy,(tt(1)+tt(end))/2,'linear');
        xm(i)    = interp1(tt,xx,(tt(1)+tt(end))/2,'linear');
        lonm(i)  = interp1(tt,lonn,(tt(1)+tt(end))/2,'linear');
        latm(i)  = interp1(tt,latt,(tt(1)+tt(end))/2,'linear');
        % Should use the interpolated depth here
        timup(i) = (tt(1)+tt(end))/2;  % time of the cast; added by Weifeng Zhang
        %dpm = interp1(time(indjd),depth(indjd),(tt(1)+tt(end))/2,'linear');
        %dpthup(:,i) = linspace(0,2.*dpm,Nd);
        dpm = depth(indjd(i));    % Corrected by Weifeng Zhang
        dpthup(:,i) = linspace(0,dpm,Nd);
        salup(:,i)  = interp1(dd,ssa,dpthup(:,i),'linear');
        temup(:,i)  = interp1(dd,tem,dpthup(:,i),'linear');
        denup(:,i)  = interp1(dd,den,dpthup(:,i),'linear');
        turbup(:,i) = interp1(dd,turbw,dpthup(:,i),'linear');
        parup(:,i)  = interp1(dd,parr,dpthup(:,i),'linear');
        florup(:,i) = interp1(dd,flo,dpthup(:,i),'linear');
        %oxygup(:,i) = interp1(dd,oxy,dpthup(:,i),'linear');
        
        % align the oxygen time series with the other variables
        % to minimize the mismatch between this cast and the two neighboring downcast;
        if algno2 == 1 && length(downup) > i
            tind1 = downup(i):downdn(i);
            tdd1  = depth(tind1);
            [tindx1,tindy1] = unique(tdd1);
            tdd1  = tindx1;
            ttt1  = time(tind1(tindy1));
            
            tind2 = downup(i+1):downdn(i+1);
            tdd2  = depth(tind2);
            [tindx2,tindy2] = unique(tdd2);
            tdd2  = tindx2;
            ttt2  = time(tind2(tindy2));
            clear epsilon;
            for io = minoffset:maxoffset
                uoxy     = oxyg(ind(indy)+io);
                uoxygup  = interp1(dd,uoxy,dpthup(:,i),'linear');
                doxy1    = oxyg(tind1(tindy1)+io);
                doxy2    = oxyg(tind2(tindy2)+io);
                toxygdn1 = interp1(tdd1,doxy1,dpthup(:,i),'linear');
                toxygdn2 = interp1(tdd2,doxy2,dpthup(:,i),'linear');
                epsilon(io-minoffset+1) = nanmean((2.*uoxygup-toxygdn1-toxygdn2).^2);
            end
            [mineps, tuoffset] = min(epsilon);
            uoffset(i) = tuoffset+minoffset-1;
            
        else
            uoffset(i) = 0;
        end
        oxy    = oxyg(ind(indy)+uoffset(i));
        oxygup(:,i) = interp1(dd,oxy,dpthup(:,i),'linear');
    end
    uoffset
    lonup = repmat(lonm,Nd,1);
    latup = repmat(latm,Nd,1);
    yup   = repmat(ym,Nd,1);
    xup   = repmat(xm,Nd,1);
    
    % for the downcast
    lonm=[]; latm=[]; xm=[]; ym=[]; timdn=[];
    ydn    = zeros(Nd,length(downdn));
    dpthdn = ydn;
    saldn  = ydn;
    temdn  = ydn;
    dendn  = ydn;
    turbdn = ydn;
    pardn  = ydn;
    flordn = ydn;
    oxygdn = ydn;
    doffset= [];
    for i = 1:length(downdn)
        
        ind = downup(i):downdn(i);
        
        dd  = depth(ind);
        [indx,indy] = unique(dd);
        dd  = indx;
        tt   = time(ind(indy));
        latt = latitude(ind(indy));
        lonn = longitude(ind(indy));
        xx   = x(ind(indy));
        yy   = y(ind(indy));
        
        ssa   = salt(ind(indy));
        tem   = temp(ind(indy));
        den   = dens(ind(indy));
        turbw = turb(ind(indy));
        parr  = par(ind(indy));
        %oxy   = oxyg(ind(indy));
        flo   = fluo(ind(indy));
        % notice that time is not equal dist, so need interpolation
        ym(i) = interp1(tt,yy,(tt(1)+tt(end))/2,'linear');
        xm(i) = interp1(tt,xx,(tt(1)+tt(end))/2,'linear');
        lonm(i)     = interp1(tt,lonn,(tt(1)+tt(end))/2,'linear');
        latm(i)     = interp1(tt,latt,(tt(1)+tt(end))/2,'linear');
        % Should use the interpolated depth here
        timdn(i)    = (tt(1)+tt(end))/2;  % time of the cast; added by Weifeng Zhang
        %dpm = interp1(time(indjd),depth(indjd),(tt(1)+tt(end))/2,'linear');
        %dpthdn(:,i) = linspace(0,2.*dpm,Nd);
        dpm         = depth(indjd(i));   % corrected by Weifeng Zhang
        dpthdn(:,i) = linspace(0,dpm,Nd);
        saldn(:,i)  = interp1(dd,ssa,dpthdn(:,i),'linear');
        temdn(:,i)  = interp1(dd,tem,dpthdn(:,i),'linear');
        dendn(:,i)  = interp1(dd,den,dpthdn(:,i),'linear');
        turbdn(:,i) = interp1(dd,turbw,dpthdn(:,i),'linear');
        pardn(:,i)  = interp1(dd,parr,dpthdn(:,i),'linear');
        flordn(:,i) = interp1(dd,flo,dpthdn(:,i),'linear');
        %oxygdn(:,i) = interp1(dd,oxy,dpthdn(:,i),'linear');
        
        
        % align the oxygen time series with the other variables
        % to minimize the mismatch between this cast and the next downcast;
        if algno2 == 1 && length(updn) >= i && i > 1
            tind1 = updn(i-1):upup(i-1);
            tdd1  = depth(tind1);
            [tindx1,tindy1] = unique(tdd1);
            tdd1  = tindx1;
            ttt1  = time(tind1(tindy1));
            
            tind2 = updn(i):upup(i);
            tdd2  = depth(tind2);
            [tindx2,tindy2] = unique(tdd2);
            tdd2  = tindx2;
            ttt2  = time(tind2(tindy2));
            clear epsilon;
            for io = minoffset:maxoffset
                doxy     = oxyg(ind(indy)+io);
                doxygdn  = interp1(dd,doxy,dpthdn(:,i),'linear');
                uoxy1    = oxyg(tind1(tindy1)+io);
                uoxy2    = oxyg(tind2(tindy2)+io);
                toxygdn1 = interp1(tdd1,uoxy1,dpthdn(:,i),'linear');
                toxygdn2 = interp1(tdd2,uoxy2,dpthdn(:,i),'linear');
                epsilon(io-minoffset+1)  = nanmean((2.*doxygdn-toxygdn1-toxygdn2).^2);
            end
            [mineps, tdoffset] = min(epsilon);
            doffset(i) = tdoffset + minoffset - 1;
        else
            doffset(i) = 0;
        end
        oxy    = oxyg(ind(indy)+doffset(i));
        oxygdn(:,i) = interp1(dd,oxy,dpthdn(:,i),'linear');
        
    end
    doffset
    londn = repmat(lonm,Nd,1);
    latdn = repmat(latm,Nd,1);
    ydn   = repmat(ym,Nd,1);
    xdn   = repmat(xm,Nd,1);
    
    xc    = [xup,xdn];
    yc    = [yup,ydn];
    timc  = [timup(:)',timdn(:)'];   % added by Weifeng Zhang
    latc  = [latup,latdn];
    lonc  = [lonup,londn];
    dpthc = [dpthup,dpthdn];
    salc  = [salup,saldn];
    temc  = [temup,temdn];
    denc  = [denup,dendn];
    turbc = [turbup,turbdn];
    parc  = [parup,pardn];
    oxygc = [oxygup,oxygdn];
    florc = [florup,flordn];
    
    % the following two line only helps to figure out the sort
    %indc = cat(2,fix((upup+updn)/2),fix((downup+downdn)/2));
    % sort by latc, checked, the same as above, so use above
    %indc = latc(1,:);
    % Using time variable to sort the cast (by Weifeng Zhang);
    indc = timc;
    
    [indcx,indcy] = sort(indc);
    xc    = xc(:,indcy);
    yc    = yc(:,indcy);
    latc  = latc(:,indcy);
    lonc  = lonc(:,indcy);
    dpthc = dpthc(:,indcy);
    salc  = salc(:,indcy);
    temc  = temc(:,indcy);
    denc  = denc(:,indcy);
    turbc = turbc(:,indcy);
    parc  = parc(:,indcy);
    oxygc = oxygc(:,indcy);
    florc = florc(:,indcy);
    
    latl  = latc(1,:);
    lonl  = lonc(1,:);
    sdist = [];
    for ii = 1:length(latl)
        sdist(ii) = sw_dist([latl(ii) lat0],[lonl(ii) lon0],'km'); % in km
    end
    
    distc = repmat(sdist,Nd,1);
    
    % for downcast only
    latld  = latdn(1,:);
    lonld  = londn(1,:);
    sdistd = [];
    for ii = 1:length(latld)
        sdistd(ii) = sw_dist([latld(ii) lat0],[lonld(ii) lon0],'km'); % in km
    end
    distcd = repmat(sdistd,Nd,1);
    
    
    % Now map the measured data to bathmetry
    dpl    = griddata(bdata.lon',bdata.lat,bdata.h,lonl,latl);
        
    % compute the cross-shelf gradients
    ddistc = diff(distc,1,2);
    drdy   = -diff(denc,1,2)./ddistc./1e3;
    dtdy   = -diff(temc,1,2)./ddistc./1e3;
    dsdy   = -diff(salc,1,2)./ddistc./1e3;
    ddpthc = dpthc(:,1:end-1);
    ddist  = av2(distc')';
    
    
    % plot 3d surface for along-shelf and cross-shelf mixed tows   
    if 1
        % plot the transect with x-axis as latitude
        figure(3);clf;
        colormap('jet');
        orient tall
        
        subplot(311);
        %pcolorjw(latc,-dpthc,temc);
        surface(lonc,latc,-dpthc,temc);
        caxis([4 20]);
        view(55,46);
        hold on; grid on;
        %[cc0,hh0] = contour(latc,-dpthc,salc,[33:0.5:36]);
        %set(hh0,'color',[0.5 0.5 0.5]);
        %[cc,hh] = contour(latc,-dpthc,denc,[fdens fdens+0.000001]);
        %set(hh,'color','k','linewidth',2);
        %[cc,hh] = contour(latc,-dpthc,salc,[fsalt fsalt+0.000001]);
        %set(hh,'color','k','linewidth',2,'linestyle','-');
        zlim([-mxdpth 0]);
        ylim([39.6 40.7]);
        hh = colorbar;
        set(gca,'fontsize',15);
        shading interp
        hmsh = mesh(blon2(1:1:end,1:1:end),blat2(1:1:end,1:1:end),bathy(1:1:end,1:1:end),'edgecolor','k');
        plot3(lonl,latl,-dpl,'k','lineWidth',2);
        %camlight; lighting phong
        plot3(slon,slat,zeros(size(slat)),'vk','linewidth',2,'markersize',4);
        for ii = 1:length(slon)
            %text(slon(ii),slat(ii)+0.01,7.5,sname{ii},'fontsize',13);
        end
        zlabel('Depth (m)','fontsize',15);
        text(40.52,-125,'Temperature ($^\circ{C}$)','interpreter','latex','fontsize',15);
        httl = title(['Transect ' num2str(tranum(1)) '-' num2str(tranum(end)) ' -- towed VPR on ' datestr(datenum(ttime))]);
        set(httl,'pos',[40.1 15],'fontsize',13);
        
        subplot(312);
        %pcolorjw(latc,-dpthc,salc);
        surface(lonc,latc,-dpthc,salc);
        view(55,46);
        caxis([31.5 36.5]);
        hold on; grid on
        %[cc0,hh0] = contour(latc,-dpthc,salc,[33:0.5:36]);
        %set(hh0,'color',[0.5 0.5 0.5]);
        %clabel(cc0,hh0);
        %[cc,hh] = contour(latc,-dpthc,denc,[fdens fdens+0.000001]);
        %set(hh,'color','k','linewidth',2);
        %[cc,hh] = contour(latc,-dpthc,salc,[fsalt fsalt+0.000001]);
        %set(hh,'color','k','linewidth',2,'linestyle','-');
        zlim([-mxdpth 0]);
        ylim([39.6 40.7]);
        hh = colorbar;
        set(gca,'fontsize',15);
        shading interp
        hmsh = mesh(blon2(1:1:end,1:1:end),blat2(1:1:end,1:1:end),bathy(1:1:end,1:1:end),'edgecolor','k');
        plot3(lonl,latl,-dpl,'k','lineWidth',2);
        plot3(slon,slat,zeros(size(slat)),'vk','linewidth',2,'markersize',4);
        for ii = 1:length(slon)
            %text(slon(ii),slat(ii)+0.01,7.5,sname{ii},'fontsize',13);
        end
        zlabel('Depth (m)','fontsize',15);
        text(40.52,-125,'Salinity','interpreter','latex','fontsize',15);
        
        subplot(313);
        %[cc0,hh0] = contourf(latc,-dpthc,denc-1000,[20:0.1:28]);
        surface(lonc,latc,-dpthc,denc-1000);
        view(55,46)
        caxis([24 28]);
        hold on; grid on;
        %set(hh0,'color',[0.5 0.5 0.5]);
        %clabel(cc0,hh0);
        %[cc0,hh0] = contour(latc,-dpthc,salc,[33:0.5:36]);
        %set(hh0,'color',[0.5 0.5 0.5]);
        %[cc,hh] = contour(latc,-dpthc,denc,[fdens fdens+0.000001]);
        %set(hh,'color','k','linewidth',2);
        %[cc,hh] = contour(latc,-dpthc,salc,[fsalt fsalt+0.000001]);
        %set(hh,'color','k','linewidth',2,'linestyle','-');
        zlim([-mxdpth 0]);
        ylim([39.6 40.7]);
        hh = colorbar;
        set(gca,'fontsize',15);
        shading interp
        hmsh = mesh(blon2(1:1:end,1:1:end),blat2(1:1:end,1:1:end),bathy(1:1:end,1:1:end),'edgecolor','k');
        plot3(lonl,latl,-dpl,'k','lineWidth',2);
        plot3(slon,slat,zeros(size(slat)),'vk','linewidth',2,'markersize',4);
        for ii = 1:length(slon)
            %text(slon(ii),slat(ii)+0.01,7.5,sname{ii},'fontsize',13);
        end
        xlabel('Longitude (^oC)','fontsize',15);
        ylabel('Latitude (^oC)','fontsize',15);
        zlabel('Depth (m)','fontsize',15);
        text(40.52,-125,'$$\sigma_\theta \; (kg\; m^{-3})$$','interpreter','latex','fontsize',15);
        print('-dpng','-r300',[ffpre 'tsr_isohal_lat.png']);
        
        figure(4);clf;
        orient tall
        colormap('jet');
        
        subplot(311);
        %pcolorjw(latc,-dpthc,florc/250);
        %pcolorjw(latdn,-dpthdn,flordn/250); % plot downcast only
        surface(londn,latdn,-dpthdn,flordn/250); % plot downcast only
        view(55,46);
        caxis([0 3]);
        hold on; grid on;
        %[cc0,hh0] = contour(latc,-dpthc,salc,[33:0.5:36]);
        %set(hh0,'color',[0.5 0.5 0.5]);
        %[cc,hh] = contour(latc,-dpthc,denc,[fdens fdens+0.000001]);
        %set(hh,'color','k','linewidth',2);
        %[cc,hh] = contour(latc,-dpthc,salc,[fsalt fsalt+0.000001]);
        %set(hh,'color','k','linewidth',2,'linestyle','-');
        zlim([-mxdpth 0]);
        ylim([39.6 40.7]);
        hh = colorbar;
        set(gca,'fontsize',15);
        shading interp
        hmsh = mesh(blon2(1:1:end,1:1:end),blat2(1:1:end,1:1:end),bathy(1:1:end,1:1:end),'edgecolor','k');
        plot3(lonl,latl,-dpl,'k','lineWidth',2);
        plot3(slon,slat,zeros(size(slat)),'vk','linewidth',2,'markersize',4);
        for ii = 1:length(slon)
            %text(slon(ii),slat(ii)+0.01,7.5,sname{ii},'fontsize',13);
        end
        zlabel('Depth (m)','fontsize',15);
        text(40.52,-120,'Chlorophyl ($$\mu g/l$$)','interpreter','latex','fontsize',15);
        text(40.52,-135,'(Florescence column / 250)','fontsize',15);
        httl = title(['Transect ' num2str(tranum(1)) '-' num2str(tranum(end)) ' -- towed VPR on ' datestr(datenum(ttime))]);
        set(httl,'pos',[40.1 15],'fontsize',13);
        
        subplot(312);
        %pcolorjw(latc,-dpthc,turbc);
        surface(lonc,latc,-dpthc,turbc);
        view(55,46);
        caxis([0 250]);
        hold on; grid on;
        %[cc0,hh0] = contour(latc,-dpthc,salc,[33:0.5:36]);
        %set(hh0,'color',[0.5 0.5 0.5]);
        %[cc,hh] = contour(latc,-dpthc,denc,[fdens fdens+0.000001]);
        %set(hh,'color','k','linewidth',2);
        %[cc,hh] = contour(latc,-dpthc,salc,[fsalt fsalt+0.000001]);
        %set(hh,'color','k','linewidth',2,'linestyle','-');
        zlim([-mxdpth 0]);
        ylim([39.6 40.7]);
        hh = colorbar;
        set(gca,'fontsize',15);
        shading interp
        hmsh = mesh(blon2(1:1:end,1:1:end),blat2(1:1:end,1:1:end),bathy(1:1:end,1:1:end),'edgecolor','k');
        plot3(lonl,latl,-dpl,'k','lineWidth',2)
        plot3(slon,slat,zeros(size(slat)),'vk','linewidth',2,'markersize',4);
        for ii = 1:length(slon)
            %text(slon(ii),slat(ii)+0.01,7.5,sname{ii},'fontsize',13);
        end
        zlabel('Depth (m)','fontsize',15);
        text(40.52,-125,'Turbidity','interpreter','latex','fontsize',15);
        
        subplot(313);
        %pcolorjw(latc,-dpthc,oxygc./10);
        %pcolorjw(latdn,-dpthdn,oxygdn./10);   % plot downcast oxygen only
        surface(londn,latdn,-dpthdn,oxygdn./10);   % plot downcast oxygen only
        view(55,46);
        %pcolorjw(latup,-dpthup,oxygup./10);   % plot downcast oxygen only
        caxis([200 350]);
        hold on; grid on
        %[cc0,hh0] = contour(latc,-dpthc,salc,[33:0.5:36]);
        %set(hh0,'color',[0.5 0.5 0.5]);
        %[cc,hh] = contour(latc,-dpthc,denc,[fdens fdens+0.000001]);
        %set(hh,'color','k','linewidth',2);
        %[cc,hh] = contour(latc,-dpthc,salc,[fsalt fsalt+0.000001]);
        %set(hh,'color','k','linewidth',2,'linestyle','-');
        zlim([-mxdpth 0]);
        ylim([39.6 40.7]);
        hh = colorbar;
        set(gca,'fontsize',15);
        shading interp
        hmsh = mesh(blon2(1:1:end,1:1:end),blat2(1:1:end,1:1:end),bathy(1:1:end,1:1:end),'edgecolor','k');
        plot3(lonl,latl,-dpl,'k','lineWidth',2);
        plot3(slon,slat,zeros(size(slat)),'vk','linewidth',2,'markersize',4);
        for ii = 1:length(slon)
            %text(slon(ii),slat(ii)+0.01,7.5,sname{ii},'fontsize',13);
        end
        xlabel('Longitude (^oC)','fontsize',15);
        ylabel('Latitude (^oC)','fontsize',15);
        zlabel('Depth (m)','fontsize',15);
        text(40.52,-125,'Oxygen (ml/l)','interpreter','latex','fontsize',15);
        print('-dpng','-r300',[ffpre 'cto_isohal_lat.png']);
        drawnow;
    end    
end


