%rdSAMOS_nc_ENRR
% Plot SAMOS data captured off SAMOS web site

% D Wolfe 3/4/2016
clear
clf
close all

pltflg = 0; %flag to plot the data: 1=yes, 0=no
pngflg = 0; %flag to make a png of the plot: 1=yes, 0=no
cruise = 'ENRR_2016';


inpath = 'c:/Data/ElNino/data/SAMOS/netcdf/WTEC_';
outpath  = 'c:/Data/ElNino/data/SAMOS/images/';

% Change for date/time you want to process
by = 2016;
bm = 3;
bd = 30;
ey = 2016;
em = 3;
ed = 31;


if 0  %( nargin == 0 )
    bys = input('Input the beginning year to process (2007) ','s');
    bms = input('Input the beginning the month to process (1-12) ','s');
    bds = input('Input the beginning day to process (1-31) ','s');
    eys = input('Input the endding year to process (2007) ','s');
    ems = input('Input the endding the month to process (1-12) ','s');
    eds = input('Input the endding day to process (1-31) ','s');
    by = str2double( bys );
    bm = str2double( bms );
    bd = str2double( bds );
    ey = str2double( eys );
    em = str2double( ems );
    ed = str2double( eds );
end

bdn = datenum([by bm bd 0 0 0]);
edn = datenum([ey em ed 0 0 0]);
j = 1;
for idn = bdn:edn
    dv = datevec(idn);
    %     yday = datenum( dv(1), dv(2), dv(3) );
    jul_day = datenum(dv(1),dv(2),dv(3),0,0,0)-datenum(dv(1)-1,12,31,00,00,00);
    jd = sprintf('%03d', jul_day);
    if idn == bdn
        xmin = jul_day;
    end
    
    %     WTEC_20160217v10001
    fname = sprintf('%s%04d%02d%02dv10001.nc', inpath, dv(1), dv(2), dv(3) )
    fname2 = sprintf('%s%04d%02d%02dv10002.nc', inpath, dv(1), dv(2), dv(3) )
    
    %     fname = 'c:/Data/ElNino/data/SAMOS/netcdf/WTEC_20160217v10001.nc'
    %      fprintf(1,'*** Processing SAMOS daily %d-%02d-%02d***\n',...
    %          dv(1),dv(2),dv(3));
    
    fid = fopen( fname );
    if( fid > 0 )
        fclose( fid );
    else
        fprintf(1,'No SAMOS1 file: %s\n',fname);
        continue;
    end
    
    ncid = netcdf.open(fname,'NC_NOWRITE');
    % From Jeremy
    %     nc = netcdf(fname, 'nowrite');                        % Open NetCDF file.
    %      description = ncid.description(:);
    %ncdisp(fname)
    % Global attribute.
    %      variables = var(ncid);    % Get variable data.
    
    [tname,tlen]=netcdf.inqDim(ncid,3);
    tid = netcdf.inqVarID(ncid,'time_of_day');
    Ta = netcdf.inqVarID(ncid,'T');
    Ts = netcdf.inqVarID(ncid,'TS');
    dat = netcdf.inqVarID(ncid,'date');
    pres = netcdf.inqVarID(ncid,'P');
    Rh = netcdf.inqVarID(ncid,'RH');
    
    time = netcdf.getVar(ncid,tid,'double'); %hhmmss
    airt = netcdf.getVar(ncid,Ta,'double');
    sst = netcdf.getVar(ncid,Ts,'double');
    DT = netcdf.getVar(ncid,dat,'double'); %YYYYMMDD
    Ps = netcdf.getVar(ncid,pres,'double'); %pressure
    rh = netcdf.getVar(ncid,Rh,'double'); %relative humidity
    
    
    ln = netcdf.inqVarID(ncid,'lon');
    lt = netcdf.inqVarID(ncid,'lat');
    lon = netcdf.getVar(ncid,ln,'double');
    lat = netcdf.getVar(ncid,lt,'double');
    
    TWs = netcdf.inqVarID(ncid,'SPD');  %FWD earth relative wind direction
    TWd = netcdf.inqVarID(ncid,'DIR');  %FWD earth relative wind direction
    tws = netcdf.getVar(ncid,TWs,'double');
    twd = netcdf.getVar(ncid,TWd,'double');
    
    RWs = netcdf.inqVarID(ncid,'PL_WSPD');  %FWD earth relative wind direction
    RWd = netcdf.inqVarID(ncid,'PL_WDIR');  %FWD earth relative wind direction
    rws = netcdf.getVar(ncid,RWs,'double');
    rwd = netcdf.getVar(ncid,RWd,'double');
    
    Crs = netcdf.inqVarID(ncid,'PL_CRS');  %FWD earth relative wind direction
    Hd = netcdf.inqVarID(ncid,'PL_HD');  %FWD earth relative wind direction
    hd = netcdf.getVar(ncid,Hd,'double'); %Platform course over ground
    crs = netcdf.getVar(ncid,Crs,'double'); % Platform heading over ground
    Pspd = netcdf.inqVarID(ncid,'PL_SPD');  %Platform  speed over ground
    pspd = netcdf.getVar(ncid,Pspd,'double'); %Platform speed over ground
    
    Rsw = netcdf.inqVarID(ncid,'RAD_SW');  %Incoming short-wave
    rsw = netcdf.getVar(ncid,Rsw,'double'); %Incoming short-wave
    
    
    %     calculate Year day
    dtt = datevec(datenum(idn)); % date vector
    hh = fix(time/10000);
    t1 = (time-hh*10000);
    mm = fix(t1/100);
    
    jdd = datenum(dtt(1),dtt(2),dtt(3),hh,mm,0)-datenum(dtt(1)-1,12,31,00,00,00);
    
    netcdf.close(ncid);
    
    %     tit1 = sprintf('%s%s%04d%02d%02d', cruise(1:4), ' ', dv(1), dv(2), dv(3) );
    tit1 = 'ENRR 2016';
    
    
    figure(1)
    plot(jdd,airt,'r')
    if idn == bdn
        hold
    end
    plot(jdd,sst,'b')
    title(tit1);
    ylabel('Temp (C)');
    ylim([15 30 ]);
    legend('Air T','SST','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_TaSST_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    
    figure(2)
    plot(jdd,tws,'r')
    if idn == bdn
        hold
    end
    plot(jdd,rws,'b')
    plot(jdd,pspd,'g')
    title(tit1);
    ylabel('Wind Spd (ms^-^1)');
    ylim([ 0 30]);
    legend('True','Rel','SOG','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_Spd_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(6)
    plot(jdd,rws,'b')
    if idn == bdn
        hold
    end
    plot(jdd,tws,'r')
    title(tit1);
    ylabel('Wind Spd (ms^-^1)');
    ylim([ 0 30]);
    legend('Rel','True','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_SpdTR_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(3)
    plot(jdd,twd,'r')
    if idn == bdn
        hold
    end
    plot(jdd,rwd,'b')
    plot(jdd,crs,'k')
    plot(jdd,hd,'g')
    title(tit1);
    ylabel('Dir (deg)');
    ylim([0 360]);
    legend(' W True',' W Rel','Course','Heading','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_Dir_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(7)
    plot(jdd,twd,'r.')
    if idn == bdn
        hold
    end
    plot(jdd,rwd,'b.')
    title(tit1);
    ylabel('Dir (deg)');
    ylim([0 360]);
    legend(' W True',' W Rel','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_DirTR_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(4)
    plot(jdd,rh,'k')
    title(tit1);
    ylabel('RH (%)')
    ylim([0 100]);
    if idn == bdn
        hold
    end
    if pltflg
        print_buffer = [outpath, cruise, '_RH_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(5)
    plot(jdd,Ps,'k')
    if idn == bdn
        hold
    end
    title(tit1);
    ylabel('Pres (mb)')
    ylim([1000 1020]);
    if pltflg
        print_buffer = [outpath, cruise, '_P_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(8)
    subplot(2,1,1),plot(jdd,tws,'k')
    ylim([0 25]);
    title(tit1);
    ylabel('Wind Spd T (ms^-^1)');
    if idn == bdn
        hold
    end
    subplot(2,1,2),plot(jdd,twd,'k')
    ylabel('Wind Dir T');
    ylim([ 0 360]);
    if idn == bdn
        hold
    end
    if pltflg
        print_buffer = [outpath, cruise, '_WndSpdDirT_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(9)
    subplot(2,1,1),plot(jdd,rws,'k')
    ylim([ 0 25]);
    title(tit1);
    ylabel('Wind Spd Rel (ms^-^1)');
    if idn == bdn
        hold
    end
    subplot(2,1,2),plot(jdd,rwd,'k')
    ylabel('Wind Dir Rel');
    ylim([ 0 360]);
    if idn == bdn
        hold
    end
    if pltflg
        print_buffer = [outpath, cruise, '_WndSpdDirR_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(10)
    plot(jdd,rsw,'k')
    title(tit1);
    ylabel('SW radiation (Wm^-^2)')
    ylim([0 1300]);
    if idn == bdn
        hold
    end
    if pltflg
        print_buffer = [outpath, cruise, '_SW_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    
    %     figure
    %     plot(lon-360.,lat)
    %         title(tit1);
    %     xlabel('Longitude')
    %     ylabel('Latitude');
    %      print_buffer = [outpath, cruise, '_LatLon' jd, '.png'];
    %     print('-dpng ', print_buffer);
    
    %clear all
    %retval = 0;
    
    %% check for a 2nd file
    fid2 = fopen(fname2);
    
    if( fid2 > 0 )
        fclose( fid2 );
    else
        fprintf(1,'No SAMOS2 file: %s\n',fname2);
        continue;
    end
    ncid2 = netcdf.open(fname2,'NC_NOWRITE');
    % From Jeremy
    %     nc = netcdf(fname, 'nowrite');                        % Open NetCDF file.
    %      description = ncid.description(:);
    %ncdisp(fname)
    % Global attribute.
    %      variables = var(ncid);    % Get variable data.
    
    [tname,tlen]=netcdf.inqDim(ncid2,3);
    tid = netcdf.inqVarID(ncid2,'time_of_day');
    Ta = netcdf.inqVarID(ncid2,'T');
    Ts = netcdf.inqVarID(ncid2,'TS');
    dat = netcdf.inqVarID(ncid2,'date');
    pres = netcdf.inqVarID(ncid2,'P');
    Rh = netcdf.inqVarID(ncid2,'RH');
    
    time = netcdf.getVar(ncid2,tid,'double'); %hhmmss
    airt = netcdf.getVar(ncid2,Ta,'double');
    sst = netcdf.getVar(ncid2,Ts,'double');
    DT = netcdf.getVar(ncid2,dat,'double'); %YYYYMMDD
    Ps = netcdf.getVar(ncid2,pres,'double'); %pressure
    rh = netcdf.getVar(ncid2,Rh,'double'); %relative humidity
    
    
    ln = netcdf.inqVarID(ncid2,'lon');
    lt = netcdf.inqVarID(ncid2,'lat');
    lon = netcdf.getVar(ncid2,ln,'double');
    lat = netcdf.getVar(ncid2,lt,'double');
    
    TWs = netcdf.inqVarID(ncid2,'SPD');  %FWD earth relative wind direction
    TWd = netcdf.inqVarID(ncid2,'DIR');  %FWD earth relative wind direction
    tws = netcdf.getVar(ncid2,TWs,'double');
    twd = netcdf.getVar(ncid2,TWd,'double');
    
    RWs = netcdf.inqVarID(ncid2,'PL_WSPD');  %FWD earth relative wind direction
    RWd = netcdf.inqVarID(ncid2,'PL_WDIR');  %FWD earth relative wind direction
    rws = netcdf.getVar(ncid2,RWs,'double');
    rwd = netcdf.getVar(ncid2,RWd,'double');
    
    Crs = netcdf.inqVarID(ncid2,'PL_CRS');  %FWD earth relative wind direction
    Hd = netcdf.inqVarID(ncid2,'PL_HD');  %FWD earth relative wind direction
    hd = netcdf.getVar(ncid2,Hd,'double'); %Platform course over ground
    crs = netcdf.getVar(ncid2,Crs,'double'); % Platform heading over ground
    Pspd = netcdf.inqVarID(ncid2,'PL_SPD');  %Platform  speed over ground
    pspd = netcdf.getVar(ncid2,Pspd,'double'); %Platform speed over ground
    
    Rsw = netcdf.inqVarID(ncid,'RAD_SW');  %Incoming short-wave
    rsw = netcdf.getVar(ncid,Rsw,'double'); %Incoming short-wave
    
    
    %     calculate Year day
    dtt = datevec(datenum(idn)); % date vector
    hh = fix(time/10000);
    t1 = (time-hh*10000);
    mm = fix(t1/100);
    
    jdd = datenum(dtt(1),dtt(2),dtt(3),hh,mm,0)-datenum(dtt(1)-1,12,31,00,00,00);
    
    netcdf.close(ncid2);
    
    %     tit1 = sprintf('%s%s%04d%02d%02d', cruise(1:4), ' ', dv(1), dv(2), dv(3) );
    tit1 = 'ENRR 2016';
    
    
    figure(1)
    plot(jdd,airt,'r')
    if idn == bdn
        hold
    end
    plot(jdd,sst,'b')
    title(tit1);
    ylabel('Temp (C)');
    ylim([15 30 ]);
    legend('Air T','SST','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_TaSST_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(2)
    plot(jdd,tws,'r')
    if idn == bdn
        hold
    end
    plot(jdd,rws,'b')
    plot(jdd,pspd,'g')
    title(tit1);
    ylabel('Wind Spd (ms^-^1)');
    ylim([ 0 30]);
    legend('True','Rel','SOG','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_Spd_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(6)
    plot(jdd,tws,'r')
    if idn == bdn
        hold
    end
    plot(jdd,rws,'b')
    title(tit1);
    ylabel('Wind Spd (ms^-^1)');
    ylim([ 0 30]);
    legend('True','Rel','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_SpdTR_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(3)
    plot(jdd,twd,'r')
    if idn == bdn
        hold
    end
    plot(jdd,rwd,'b')
    plot(jdd,crs,'k')
    plot(jdd,hd,'g')
    title(tit1);
    ylabel('Dir (deg)');
    ylim([0 360]);
    legend(' W True',' W Rel','Course','Heading','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_Dir_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(7)
    plot(jdd,twd,'r.')
    if idn == bdn
        hold
    end
    plot(jdd,rwd,'b.')
    title(tit1);
    ylabel('Dir (deg)');
    ylim([0 360]);
    legend(' W True',' W Rel','Location','Best');
    if pltflg
        print_buffer = [outpath, cruise, '_DirTR_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(4)
    plot(jdd,rh,'k')
    title(tit1);
    ylabel('RH (%)')
    ylim([0 100]);
    if idn == bdn
        hold
    end
    if pltflg
        print_buffer = [outpath, cruise, '_RH_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(5)
    plot(jdd,Ps,'k')
    if idn == bdn
        hold
    end
    title(tit1);
    ylabel('Pres (mb)')
    ylim([1000 1020]);
    if pltflg
        print_buffer = [outpath, cruise, '_P_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(8)
    subplot(2,1,1),plot(jdd,tws,'k')
    ylim([ 0 30]);
    title(tit1);
    ylabel('Wind Spd T(m^-^1)');
    if idn == bdn
        hold
    end
    subplot(2,1,2),plot(jdd,twd,'k.')
    ylabel('Wind Dir T');
    ylim([ 0 360]);
    if idn == bdn
        hold
    end
    if pltflg
        print_buffer = [outpath, cruise, '_WndSpdDirt_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(9)
    subplot(2,1,1),plot(jdd,rws,'k')
    title(tit1);
    ylabel('Wind Spd Rel (m^-^1)');
    ylim([ 0 30]);
    if idn == bdn
        hold
    end
    subplot(2,1,2),plot(jdd,rwd,'k.')
    ylabel('Wind Dir Rel');
    ylim([ 0 360]);
    if idn == bdn
        hold
    end
    if pltflg
        print_buffer = [outpath, cruise, '_WndSpdDirR_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(10)
    plot(jdd,rsw,'k')
    title(tit1);
    ylabel('SW radiation (Wm^-^2)')
    ylim([0 1300]);
    if idn == bdn
        hold
    end
    if pltflg
        print_buffer = [outpath, cruise, '_SW_' jd, '.png'];
        print('-dpng ', print_buffer);
    end
    %
    
    
end  % loop through all files
xmax = (fix(max(jdd)))+1;

figure(1)
xlabel('Year Day (UTC)');
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_TaSST.png'];
print('-dpng ', print_buffer);

figure(2)
xlabel('Year Day (UTC)');
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_Spd.png'];
print('-dpng ', print_buffer);

figure(3)
xlabel('Year Day (UTC)');
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_Dir.png'];
print('-dpng ', print_buffer);

figure(4)
xlabel('Year Day (UTC)');
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_RH.png'];
print('-dpng ', print_buffer);

figure(5)
xlabel('Year Day (UTC)');
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_P.png'];
print('-dpng ', print_buffer);

figure(6)
xlabel('Year Day (UTC)');
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_SpdTR.png'];
print('-dpng ', print_buffer);

figure(7)
xlabel('Year Day (UTC)');
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_DirTR.png'];
print('-dpng ', print_buffer);

figure(8)
xlabel('Year Day (UTC)');
subplot(2,1,1)
xlim([xmin xmax])
ylim([0 25]);
grid
subplot(2,1,2)
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_WndSpdDirR.png'];
print('-dpng ', print_buffer);

figure(9)
xlabel('Year Day (UTC)');
subplot(2,1,1)
xlim([xmin xmax])
ylim([ 0 25]);
grid
subplot(2,1,2)
xlim([xmin xmax])
grid
print_buffer = [outpath, cruise, '_WndSpdDirT.png'];
print('-dpng ', print_buffer);

figure(10)
xlabel('Year Day (UTC)');
xlim([xmin xmax]);
grid
print_buffer = [outpath, cruise, '_SW.png'];
print('-dpng ', print_buffer);