clear
prtflag   = 0; % flag
figdir    = '~/research/cruise/2019_Brown_0512-0525/data/figs/eddy/';
matdir    = '~/research/cruise/2019_Brown_0512-0525/data/figs/sst/mat/';
%files     = dir(matdir);
%files.name = ['udelsst_mab_20190515_2135.mat'];
%files.name = ['udelsst_mab_20190515_2348.mat'];
%files.name = ['udelsst_mab_20190516_0128.mat'];
%files.name = ['udelsst_mab_20190516_2124.mat'];
%files.name = ['udelsst_mab_20190516_2337.mat'];
%files.name = ['udelsst_mab_20190517_0947.mat'];
%files.name = ['udelsst_mab_20190517_1338.mat'];
%files.name = ['udelsst_mab_20190518_1327.mat'];
%files.name = ['udelsst_mab_20190519_0052.mat'];
files.name = ['udelsst_mab_20190520_0041.mat'];

% ADCP data
addadcp  = 1;
tranum   = 17:31;   % the transect number for adding adcp surface current
                    % *** NEED the start time of the (future) transect after tranum(end)
                    % which is provided in the last line (place holder) of rb1904_transect_startime.m 
quvscale = 0.5;   % quiver vector scale
uvscale  = 0.2;   %
adcpfile = '~/research/cruise/2019_Brown_0512-0525/data/adcp/os75bb_uv_blank8mbin4m_17to19.mat';

% read transect start time
addpath ~/research/cruise/2019_Brown_0512-0525/data/scripts/
rb1904_transect_startime

% underway data
addtsg  = 1;
tsgfile = '~/research/cruise/2019_Brown_0512-0525/data/gps/tsg_with_gps_17to19.mat';
addflr  = 1;
flrfile = '~/research/cruise/2019_Brown_0512-0525/data/gps/flr_with_gps_17to19.mat';

% area for plotting
plonlim = [-71.5 -70.5];
platlim = [39.5 40.5];    % northeast

clrlim = [4 19];

bdata   = load('~/research/roms/bathymetry/ngdc2.mat');
blonind = find(bdata.lon>plonlim(1) & bdata.lon < plonlim(2));
blatind = find(bdata.lat>platlim(1) & bdata.lat < platlim(2));
blon    = bdata.lon(blonind);
blat    = bdata.lat(blatind);
bdpth   = bdata.h(blatind(1):blatind(end),blonind(1):blonind(end));

% mooring locations
mrlon  = [-70.88333 -70.88333 -70.775   -70.8833 -70.88333 -70.775 -70.775];
mrlat  = [40.363333 40.226667 40.136667 40.09833 39.94     40.365  39.94];
mrname = ['IS','CI','CN','CO','OS','UI','UO'];
sname  = {'A1','A2','A3','A4','A5','A6','A7','A8','A9','A10','A11','A12','A13','A14','A15','A16','A17','A18'};
zdlon  = [-70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830];
zdlat  = [ 40.725  40.660  40.595  40.530  40.465  40.400  40.335  40.270  40.205  40.140  40.075  40.010  39.945  39.880  39.815  39.750  39.685  39.620];

% to read eddy chasing stations
eddypwf = '~/research/cruise/2019_Brown_0512-0525/data/eddy/eddy_waypoints_3.txt';
fid     = fopen(eddypwf);
wpind   = 0;
while ~feof(fid) %Check every line in the file
    wpind = wpind + 1;
    [llstr(wpind,:) count] = fscanf(fid,'%c%c%d    %d %f %c    %d %f %c\n',[1,9]);
end
fclose(fid);
for iw = 1:wpind
    ecsname{iw} = [char(llstr(iw,1)) char(llstr(iw,2)) num2str(llstr(iw,3))];
    latdgr(iw)  = llstr(iw,4);
    latmnt(iw)  = llstr(iw,5);
    latnsl(iw)  = llstr(iw,6);
    londgr(iw)  = llstr(iw,7);
    lonmnt(iw)  = llstr(iw,8);
    lonewl(iw)  = llstr(iw,9);
    ecslat(iw)  = latdgr(iw) + latmnt(iw)./60;
    ecslon(iw)  = londgr(iw) + lonmnt(iw)./60;
    if latnsl(iw) == 83
        ecslat(iw) = -ecslat(iw);
    elseif latnsl(iw) ~= 78
        errmes = ['Latitude letter error in file: ' eddypwf];
        error(errmes);
    end
    if lonewl(iw) == 87
        ecslon(iw) = -ecslon(iw);
    elseif lonewl(iw) ~= 69
        errmes = ['Latitude letter error in file: ' eddypwf];
        error(errmes);
    end
end

for ifile = 1:length(files)
    
    % find the time index to read (note UDel time is generally several
    % minutes earlier than Rutgers
    fname = files(ifile).name;
    if length(fname) > 10 && strcmp(fname(end-2:end),'mat')
        mcsst  = load([matdir fname]);
        lon    = mcsst.lon;
        lat    = mcsst.lat;
        nmcsst = mcsst.mcsst;
        time   = mcsst.time;
                    
        figure(1);clf;
        orient landscape
        colormap('jet');
        
        axes('pos',[0.09 0.2 0.4 0.75]);
        pcolorjw(lon,lat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon,blat,bdpth,[40:10:100 200 500 1000 2000 3000]);
        set(hh,'color','k','linewidth',1);
        [cc,hh] = contour(blon,blat,bdpth,[100 200]);
        set(hh,'color','k','linewidth',2);
        clabel(cc,hh,'color',[0.5 0.5 0.5],'fontsize',15,'labelspacing',350);
        pnc;
        pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
        xlim(plonlim);ylim(platlim);
        hclb = colorbar;
        set(hclb,'location','southoutside');
        set(hclb,'pos',[0.1 0.13 0.38 0.02],'ytick',[clrlim(1):2:clrlim(2)],'tickdir','out','fontsize',15);
        xlabel(hclb,'Temperature','fontsize',15);
        set(gca,'xtick',[plonlim(1):0.25:plonlim(2)],'ytick',[platlim(1):0.25:platlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15);
        plot(mrlon,mrlat,'pk','markersize',6,'linewidth',2);
        plot(zdlon,zdlat,'dk','markersize',6,'linewidth',2);
        plot(ecslon,ecslat,'ok','markersize',6,'linewidth',2);
        tttxt   = ['SST on ' datestr(time)];
        
        if addadcp
            gpsfile = '~/research/cruise/2019_Brown_0512-0525/data/gps/gps_17to19.mat';
            gps     = load(gpsfile);
            gpslon2 = gps.lon(1:60:end); % decimated gps for computing average speed
            gpslat2 = gps.lat(1:60:end);
            gpstim2 = datenum(gps.gpsdate(1:60:end));
            dglon2  = (gpslon2(2:end)-gpslon2(1:end-1))./180.*pi;
            dglat2  = (gpslat2(2:end)-gpslat2(1:end-1))./180.*pi;
            dgtime2 = (gpstim2(2:end)-gpstim2(1:end-1)).*86400;
            dgdist2 = sqrt(dglon2.^2 + dglat2.^2) .* rearth;
            gspeed2 = dgdist2 ./ dgtime2;
            gpstimc = (gpstim2(1:end-1) + gpstim2(2:end))./2;            
            
            adcp   = load(adcpfile);
            trind  = find(adcp.adcptime > tstime(tranum(1)) & adcp.adcptime < tstime(tranum(end)+1));
            tru    = adcp.u(:,trind);
            trv    = adcp.v(:,trind);
            trlon0 = adcp.adcplon(trind);
            trlat0 = adcp.adcplat(trind);
            trtime = adcp.adcptime(trind);
            surfu0 = tru(1,:);
            surfv0 = trv(1,:);
            
            trsped = interp1(gpstimc,gspeed2,trtime);
            smvind = find(trsped > 1);
            surfu  = surfu0(smvind);
            surfv  = surfv0(smvind);
            trlon  = trlon0(smvind);
            trlat  = trlat0(smvind);
            
            hqv    = quiver(trlon(:),trlat(:),surfu(:).*quvscale,surfv(:).*quvscale);
            set(hqv,'autoscale','off','color','k','linewidth',2);
            hqv2   = quiver(trlon(:),trlat(:),surfu(:).*quvscale,surfv(:).*quvscale);
            set(hqv2,'autoscale','off','color','w','linewidth',0.5);
            hqvs   = quiver(-71.45,39.6,uvscale.*quvscale,0);
            set(hqvs,'autoscale','off','color','k','linewidth',2);
            hqvs2   = quiver(-71.45,39.6,uvscale.*quvscale,0);
            set(hqv2,'autoscale','off','color','w','linewidth',0.5);
            text(-71.49,39.55,[num2str(uvscale) ' m s^-^1'],'fontsize',15,'color','k');
        end

        axes('pos',[0.51 0.2 0.4 0.75]);
        pcolorjw(lon,lat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon,blat,bdpth,[40:10:100 200 500 1000 2000 3000]);
        set(hh,'color','k','linewidth',1);
        [cc,hh] = contour(blon,blat,bdpth,[100 200]);
        set(hh,'color','k','linewidth',2);
        clabel(cc,hh,'color',[0.5 0.5 0.5],'fontsize',15,'labelspacing',350);
        pnc;
        pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
        hclb = colorbar;
        set(hclb,'location','southoutside');
        set(hclb,'pos',[0.52 0.13 0.38 0.02],'ytick',[clrlim(1):2:clrlim(2)],'tickdir','out','fontsize',15);
        xlabel(hclb,'Temperature','fontsize',15);
        %text(plonlim(2)+0.02,platlim(1)-0.2,'(^oC)','fontsize',15);
        xlim(plonlim);ylim(platlim);
        set(gca,'xtick',[plonlim(1):0.25:plonlim(2)],'ytick',[platlim(1):0.25:platlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15,'yticklabel',[]);
        plot(mrlon,mrlat,'pk','markersize',6,'linewidth',2);
        plot(zdlon,zdlat,'dk','markersize',6,'linewidth',2);
        plot(ecslon,ecslat,'ok','markersize',6,'linewidth',2);
        tttxt   = ['SST on ' datestr(time)];
        if addtsg
            tsg     = load(tsgfile);
            tsgind  = find(datenum(tsg.date_tsg) > tstime(tranum(1)) & datenum(tsg.date_tsg) < tstime(tranum(end)+1));
            tsgtemp = tsg.t1(tsgind);
            tsgsalt = tsg.s(tsgind);
            nnind   = find(tsgsalt < 30); % remove the salinity data less than 30
            tsgsalt(nnind) = [];
            tsgtemp(nnind) = [];
            tsgind(nnind)  = [];
            tsglon  = tsg.lon_tsg(tsgind);
            tsglat  = tsg.lat_tsg(tsgind);
            tsgtime = datenum(tsg.date_tsg(tsgind));
            hc = colourplot(tsglon(1:10:end),tsglat(1:10:end),tsgtemp(1:10:end));            
            tttxt   = [tttxt ' + ADCP + TSG (T) at transect ' num2str(tranum(1)) '-' num2str(tranum(end))];
            caxis(clrlim);
        end
        htt = title(tttxt);
        set(htt,'fontsize',14,'fontname','Times','pos',[plonlim(1) platlim(2)+0.03]);
        if prtflag
            figname = [figdir 'ssteddy_PA_' datestr(time,'yyyymmdd') '_' datestr(time,'HHMMSS') '_adcp_tsgtemp.png'];
            print('-dpng','-r300',figname);
            disp(['Printed figure file: ' figname]);
        end
        
        figure(2);clf;
        orient landscape
        colormap('jet');
        
        hax1 = axes('pos',[0.09 0.2 0.4 0.75]);
        pcolorjw(lon,lat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon,blat,bdpth,[40:10:100 200 500 1000 2000 3000]);
        set(hh,'color','k','linewidth',1);
        [cc,hh] = contour(blon,blat,bdpth,[100 200]);
        set(hh,'color','k','linewidth',2);
        clabel(cc,hh,'color',[0.5 0.5 0.5],'fontsize',15,'labelspacing',350);
        pnc;
        pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
        xlim(plonlim);ylim(platlim);
        set(gca,'xtick',[plonlim(1):0.25:plonlim(2)],'ytick',[platlim(1):0.25:platlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15);
        plot(mrlon,mrlat,'pk','markersize',6,'linewidth',2);
        plot(zdlon,zdlat,'dk','markersize',6,'linewidth',2);
        plot(ecslon,ecslat,'ok','markersize',6,'linewidth',2);
        if addtsg
            hax2 = axes('pos',[0.09 0.2 0.4 0.75]);
            caxis(hax2,[32.5 36]);
            hc = colourplot(tsglon(1:10:end),tsglat(1:10:end),tsgsalt(1:10:end));
            pnc;
            pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
            xlim(plonlim);ylim(platlim);
            hclb = colorbar('peer',hax2);
            set(hclb,'location','southoutside');
            set(hclb,'pos',[0.1 0.13 0.38 0.02],'ytick',[32:1:37],'tickdir','out','fontsize',15);
            xlabel(hclb,'Salinity','fontsize',15);
            set(hax2,'visible','off');
        end
        
        hax3 = axes('pos',[0.51 0.2 0.4 0.75]);
        pcolorjw(lon,lat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon,blat,bdpth,[40:10:100 200 500 1000 2000 3000]);
        set(hh,'color','k','linewidth',1);
        [cc,hh] = contour(blon,blat,bdpth,[100 200]);
        set(hh,'color','k','linewidth',2);
        clabel(cc,hh,'color',[0.5 0.5 0.5],'fontsize',15,'labelspacing',350);
        pnc;
        pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
        xlim(plonlim);ylim(platlim);
        set(gca,'xtick',[plonlim(1):0.25:plonlim(2)],'ytick',[platlim(1):0.25:platlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15,'yticklabel',[]);
        plot(mrlon,mrlat,'pk','markersize',6,'linewidth',2);
        plot(zdlon,zdlat,'dk','markersize',6,'linewidth',2);
        plot(ecslon,ecslat,'ok','markersize',6,'linewidth',2);
        tttxt   = ['SST on ' datestr(time)];
        if addflr
            flr     = load(flrfile);
            flrind  = find(datenum(flr.date_flr) > tstime(tranum(1)) & datenum(flr.date_flr) < tstime(tranum(end)+1));
            flrflr  = flr.flr(tsgind);
            flrlon  = flr.lon_flr(flrind);
            flrlat  = flr.lat_flr(flrind);
            flrtime = datenum(flr.date_flr(flrind));
            
            hax4 = axes('pos',[0.51 0.2 0.4 0.75]);
            caxis([0 3]);
            hc = colourplot(flrlon(1:10:end),flrlat(1:10:end),flrflr(1:10:end));
            pnc;
            pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
            xlim(plonlim);ylim(platlim);
            set(hax4,'visible','off');
            tttxt   = [tttxt ' + TSG (S & Chl) at transect ' num2str(tranum(1)) '-' num2str(tranum(end))];
            caxis(hax4,[0 3]);
            hclb = colorbar('peer',hax4);
            set(hclb,'location','southoutside');
            set(hclb,'pos',[0.52 0.13 0.38 0.02],'ytick',[0:1:3],'tickdir','out','fontsize',15);
            xlabel(hclb,'Fluorescence','fontsize',15);
        end
        htt2 = title(hax3,tttxt);
        set(htt2,'fontsize',14,'fontname','Times','pos',[plonlim(1) platlim(2)+0.03]);
        if prtflag            
            figname = [figdir 'ssteddy_PA_' datestr(time,'yyyymmdd') '_' datestr(time,'HHMMSS') '_tsgsalt_flr.png'];
            print('-dpng','-r300',figname);
            disp(['Printed figure file: ' figname]);
        end
    end
    
end
