

sstime(1) = ssttime;

readurl   = 1;  % switch to extract data from the UDel server; otherwise read local data

%url       = 'http://tds.maracoos.org/thredds/dodsC/AVHRR.nc';
url       = 'http://basin.ceoe.udel.edu/thredds/dodsC/avhrr_unfiltered_sst.nc';

% area for extracting data
lonlim  = [-76 -66];  % northeast
latlim  = [36 42];    % northeast
% area for plotting
plonlim = [-71.75 -69.75];
platlim = [39.25 41.1];    % northeast

clrlim  = [10 26];

if addadcp
    gps_tsg_ho_wgz2
    read_adcp
    % ADCP data
    tranum   = [transect];   % the transect number for adding adcp surface current
    quvscale = 1;   % quiver vector scale
    uvscale  = 0.2; %
    adcpfile = '/Users/hildeoliver/Documents/cruise_dashboard/adcp/os75bb_uv_blank8mbin8m.mat';
    % read transect start time
    addpath /Users/hildeoliver/Documents/cruise_dashboard/
    tn368_transect_startime
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];


if readurl  % read data from UDel server
    udtime = double(ncread(url,'time'))./86400 + datenum([1970 1 1]);
    udlon  = ncread(url,'lon');
    udlat  = ncread(url,'lat');
    % load('udel_sst_lonlat.mat');
    lonind = find(udlon>lonlim(1) & udlon < lonlim(2));
    latind = find(udlat>latlim(1) & udlat < latlim(2));
    lonlen = length(lonind);
    latlen = length(latind);
    lon    = udlon(lonind);
    lat    = udlat(latind);
else   % read local data
    matfile = dir(matdir);
    for iff = 3:length(matfile)
        matname = [matfile(iff).folder '/' matfile(iff).name];
        tmptime = load(matname,'time');
        udtime(iff-2) = tmptime.time;
    end
    tmplon  = load(matname,'lon');
    tmplat  = load(matname,'lat');
    lon     = tmplon.lon;
    lat     = tmplat.lat;
end

for ifile = 1:length(sstime)
    
    % find the time index to read (note UDel time is generally several
    % minutes earlier than Rutgers
    ttind = find(udtime <= sstime(ifile));
    tind  = ttind(end);
    time  = udtime(tind);
    if abs(time - sstime(ifile)) < 30/(24*60)
        % only save the data when the UDel time is within half a hour
        % window around the Rutgers time
        if readurl
            mcsst = double(ncread(url,'mcsst',[lonind(1) latind(1) tind],[lonlen latlen 1]))';
        else
            matname = [matfile(tind+2).folder '/' matfile(tind+2).name];
            tmcsst  = load(matname,'mcsst');
            mcsst   = tmcsst.mcsst;
        end
        
        
        blonind1 = find(bdata.lon>plonlim(1) & bdata.lon < plonlim(2));
        blatind1 = find(bdata.lat>platlim(1) & bdata.lat < platlim(2));
        blon2    = bdata.lon(blonind1);
        blat2    = bdata.lat(blatind1);
        bdpth2   = bdata.h(blatind1(1):blatind1(end),blonind1(1):blonind1(end));
        
        plonind = find(lon > plonlim(1) & lon < plonlim(2));
        platind = find(lat > platlim(1) & lat < platlim(2));
        nlon   = lon(plonind);
        nlat   = lat(platind);
        nmcsst = mcsst(platind(1):platind(end),plonind(1):plonind(end));
        
        figure;clf;
        orient PORTRAIT
        colormap('jet');
        
        axes('pos',[0.05 0.1 0.85 0.8]);
        pcolorjw(nlon,nlat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon2,blat2,bdpth2,[40:10:100 200 500 1000 2000 3000]);
        set(hh,'color','k','linewidth',1);
        [cc,hh] = contour(blon2,blat2,bdpth2,[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.5:plonlim(2)],'ytick',[platlim(1):0.25:platlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15);
        %ylabel('Latitude');
        hclb = colorbar;
        set(hclb,'location','eastoutside');
        set(hclb,'pos',[0.75 0.2 0.02 0.6],'ytick',[clrlim(1):2:clrlim(2)],'tickdir','out','fontsize',15);
        text(plonlim(2)+0.1,41,'(^oC)','fontsize',15);
        %xlabel('Longitude');
        plot(mrlon,mrlat,'pk','markersize',6,'linewidth',2);
        plot(zdlon,zdlat,'dk','markersize',6,'linewidth',2);
        
        figname = [figdir 'sst_PA_' datestr(time,'yyyymmdd') '_' datestr(time,'HHMMSS') '.png'];
        tttxt = ['Sat SST on ' datestr(time)];
        htt = title(tttxt);
        set(htt,'fontsize',14,'fontname','Times');
        
        

        
        
                
        lonlim  = [-74 -66];  % northeast
        latlim  = [37 42];    % northeast
        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);
        bdpth   = bdata.h(blatind(1):blatind(end),blonind(1):blonind(end));
        
        plonind = find(lon > lonlim(1) & lon < lonlim(2));
        platind = find(lat > latlim(1) & lat < latlim(2));
        nlon   = lon(plonind);
        nlat   = lat(platind);
        nmcsst = mcsst(platind(1):platind(end),plonind(1):plonind(end));
        
        figure;clf;
        orient PORTRAIT
        colormap('jet');
        
        axes('pos',[0.05 0.1 0.85 0.8]);
        pcolorjw(nlon,nlat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon,blat,bdpth,[100 500 1000 2000 3000]);
        set(hh,'color','k','linewidth',1);
        
        pnc;
        pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
        xlim(lonlim);ylim(latlim);
        set(gca,'xtick',[lonlim(1):1:lonlim(2)],'ytick',[latlim(1):1:latlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15);
        %ylabel('Latitude');
        hclb = colorbar;
        set(hclb,'location','eastoutside');
        set(hclb,'ytick',[clrlim(1):2:clrlim(2)],'tickdir','out','fontsize',15);
        %xlabel('Longitude');
        plot(zdlon,zdlat,'dk','markersize',4,'linewidth',2);
        line([plonlim(1) plonlim(1)],[platlim(1) platlim(2)],'color',[0.5 0.5 0.5],'linewidth',2);
        line([plonlim(2) plonlim(2)],[platlim(1) platlim(2)],'color',[0.5 0.5 0.5],'linewidth',2);
        line([plonlim(1) plonlim(2)],[platlim(1) platlim(1)],'color',[0.5 0.5 0.5],'linewidth',2);
        line([plonlim(1) plonlim(2)],[platlim(2) platlim(2)],'color',[0.5 0.5 0.5],'linewidth',2);

        figname = [figdir 'sst_PA_' datestr(time,'yyyymmdd') '_' datestr(time,'HHMMSS') '.png'];
        tttxt = ['Sat SST on ' datestr(time)];
        htt = title(tttxt);
        set(htt,'fontsize',14,'fontname','Times');
        
        
        
                blonind1 = find(bdata.lon>plonlim(1) & bdata.lon < plonlim(2));
        blatind1 = find(bdata.lat>platlim(1) & bdata.lat < platlim(2));
        blon2    = bdata.lon(blonind1);
        blat2    = bdata.lat(blatind1);
        bdpth2   = bdata.h(blatind1(1):blatind1(end),blonind1(1):blonind1(end));
        
        plonind = find(lon > plonlim(1) & lon < plonlim(2));
        platind = find(lat > platlim(1) & lat < platlim(2));
        nlon   = lon(plonind);
        nlat   = lat(platind);
        nmcsst = mcsst(platind(1):platind(end),plonind(1):plonind(end));
        
        figure;clf;
        orient PORTRAIT
        colormap('jet');
        
        axes('pos',[0.05 0.1 0.85 0.8]);
        pcolorjw(nlon,nlat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon2,blat2,bdpth2,[40:10:100 200 500 1000 2000 3000]);
        set(hh,'color','k','linewidth',1);
        [cc,hh] = contour(blon2,blat2,bdpth2,[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.5:plonlim(2)],'ytick',[platlim(1):0.25:platlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15);
        %ylabel('Latitude');
        hclb = colorbar;
        set(hclb,'location','eastoutside');
        set(hclb,'pos',[0.75 0.2 0.02 0.6],'ytick',[clrlim(1):2:clrlim(2)],'tickdir','out','fontsize',15);
        text(plonlim(2)+0.1,41,'(^oC)','fontsize',15);
        %xlabel('Longitude');
        plot(mrlon,mrlat,'pk','markersize',6,'linewidth',2);
        plot(zdlon,zdlat,'dk','markersize',6,'linewidth',2);
        
        if addadcp
            gpsfile = ['/Users/hildeoliver/Documents/cruise_dashboard/gps/gps_',sprintf('%02d',day(ssttime)),'.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(-70.1,39.4,uvscale.*quvscale,0);
            set(hqvs,'autoscale','off','color','k','linewidth',2);
            hqvs2   = quiver(-70.1,39.4,uvscale.*quvscale,0);
            set(hqv2,'autoscale','off','color','w','linewidth',0.5);
            text(-70.1,39.3,[num2str(uvscale) ' m s^-^1'],'fontsize',15,'color','k');
            figname = [figdir 'sst_PA_' datestr(time,'yyyymmdd') '_' datestr(time,'HHMMSS') '_adcptrsct' num2str(tranum(1)) '-' num2str(tranum(end)) '.png'];
            tttxt = ['SST on ' datestr(time) ' + ADCP at transect ' num2str(tranum(1)) '-' num2str(tranum(end))];
        else
            figname = [figdir 'sst_PA_' datestr(time,'yyyymmdd') '_' datestr(time,'HHMMSS') '.png'];
            tttxt = ['SST on ' datestr(time)];
        end
        
        
    else
        disp(['Cannot find data for ' datestr(sstime(ifile)) ' and skipping...']);
    end
    
    
    
    
end
