% Reads in the WaveWatchIII (WW3) data and get Hs and Fp for Jim Edson's 
% buoy ocated at 38.318N, -65.  It was deployed from November 1, 2005
% to January 31, 2007.

% Global files are 30 arc-minute data, or .5deg spatial resolution (~50km)
% ATTNW files (North West Atlantic) are 10 min spatial resolution (~16km)  <-- use those files here

% time:units = "days since 1990-01-01 00:00:00" ;
% time:conventions = "relative julian days with decimal part (as parts of the day )" ;
%    hs       
%            Size:       720x317x248 (global) or 343x211x248(attnw)
%            Dimensions: longitude,latitude,time
  
% 09/2019 - LB

%% Input variables
DataPath='/Volumes/PISTON/WW3_4Jim_C4.0/ATNW/'; %path of where the WW3 netcdf files reside
% year_in=2005;month_in=11;  %start of deployment
% year_out=2007;month_out=01; %end of deployment
%buoy lat/lon ref
lonref=-65;
latref=38.318;

%% reads in file and extract appropriate time series at buoy location
fnames_hs=dir(fullfile(DataPath,'*_hs.nc'));
fnames_fp=dir(fullfile(DataPath,'*_fp.nc'));

for ii=1:length(fnames)
    [out_hs] = read_netcdfFile(fullfile(DataPath,fnames_hs(ii).name)); 
    [out_fp] = read_netcdfFile(fullfile(DataPath,fnames_fp(ii).name));
    year=str2num(fnames_hs(ii).name(14:17));
    WW3_doy=out_hs.time+datenum(1990,1,1)-datenum(year,1,0);

    latind=find(abs(out_hs.latitude-latref)==min(abs(out_hs.latitude-latref))) ; %lat indice in WW3 closest to latref
    lonind=find(abs(out_hs.longitude-lonref)==min(abs(out_hs.longitude-lonref))) ; %lon indice in WW3 closest to lonref
    
    WW3_hs=squeeze(out_hs.hs(lonind,latind,:));
    WW3_fp=squeeze(out_fp.fp(lonind,latind,:));
    
    WW3_all=[WW3_doy.*0+year WW3_doy WW3_hs WW3_fp];
    
    % save copy in text format
    save(fullfile(DataPath,'HsFp_38N65W_11012005_01312007.txt'),'WW3_all','-ascii','-tabs','-append');
    
    clear WW3_all
end    

    
%% do plots
%map    
lonmin = -75; lonmax = -55; latmin = 33; latmax = 45;
[lonm,latm] = meshgrid(out.longitude,out.latitude);
f = figure;
h = worldmap([latmin latmax],[lonmin lonmax]);
s4 = geoshow(latm,lonm,out.hs(:,:,ii)','displaytype','texturemap');  % <-USE GEOSHOW
% caxis([-250 250])
h=colorbar;
set(get(h,'ylabel'),'String','Hs (m)');
title(['WW3 Hs @ 38.318N/65W ' datestr(WW3_doy(end)+datenum(year,1,0))], 'fontsize',14,'fontweight','bold');
geoshow(latref, lonref, 'Marker','+','MarkerEdgeColor','red')

%Full time series
WW3_all=load(fullfile(DataPath,'HsFp_38N65W_11012005_01312007.txt'));
WW3_time=WW3_all(:,2)+datenum(WW3_all(:,1),1,0);
figure;
subplot(2,1,1);
    plot(WW3_time,WW3_all(:,3));
    xlim([WW3_time(1) WW3_time(end)]);
    grid;
    datetick('x','yyyy/mm/dd','keeplimits')
    xlabel('Time (yyyy/mm/dd)');ylabel ('Hs (m)'); 
    title('WW3 Hs and Fp @ 38.318N/65W', 'fontsize',14,'fontweight','bold');
subplot(2,1,2);
    plot(WW3_time,WW3_all(:,4));
    xlim([WW3_time(1) WW3_time(end)]);
    grid;
    datetick('x','yyyy/mm/dd','keeplimits')
    xlabel('Time (yyyy/mm/dd)');ylabel ('Fp (s^{-1})'); 
    
