% Compare FVCOM 30 year ocean hindcast to current meter data
% found within specified time and bounding box contraints   
%
% Uses OpenSearch and OPeNDAP access routines from
% NCTOOLBOX: http://code.google.com/p/nctoolbox/

% Rich Signell (rsignell@usgs.gov) 
clear
tic
% Open a FVCOM model simulation dataset
url='http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3';
disp(sprintf('Reading grid information from %s ...',url))
nc=ncgeodataset(url);
% compare model to data in this time period:
jdstart=datenum([1985 8 1 0 0 0]);
jdstop=datenum([1985 11 1 0 0 0]);

jdm=nc.time('time');
lonc=nc.data('lonc');latc=nc.data('latc');
lonn=nc.data('lon');latn=nc.data('lat');
nv=nc.data('nv').';
% set up variable objects for depth and sigvar (no data read yet):
zvar=nc.variable('h');
sigvar=nc.variable('siglay');
%
% look for current meter data:  
u_std_name='eastward_sea_water_velocity';
v_std_name='northward_sea_water_velocity';

% look for data in this region:
%q.bbox=[min(lonn(:)) max(lonn(:)) min(latn(:)) max(latn(:))]; %whole domain
q.bbox=[-71.5 -65.0 40.0 45.0];  % Gulf of Maine
%q.bbox=[-70.9 -70.6  41.4 41.7];  % Buzzards Bay

% Use OpenSearch to find observational data that meets the contraints.
% The service endpoint in this case is a GI-CAT service that
% has harvested the NCISO metadata from three different THREDDS catalogs
% (WHOI, USGS and NOAA/NMFS time series).  

disp(sprintf('Finding data using OpenSearch...'))
%q.endpoint='http://geoport.whoi.edu/gi-cat/services/opensearch';
q.endpoint='http://geoport-dev.whoi.edu/gi-cat/services/opensearch';
q.string_text=u_std_name;
q.time_start=jdstart;
q.time_end=jdstop;

% 
[links,params]=opensearch(q);   % make the query
dap=links2dap(links); % find only the OPeNDAP links

%%
% eliminate any data with sampling rate less than 1 hour, since we capture
% both the "raw" data and "hourly averaged data" with this search. 
j=0;
for i=1:length(dap);
  ncd=ncgeodataset(dap{i});
  jdd=ncd.time('time');
  if length(jdd)>1,
    dthours=((jdd(end)-jdd(1))/(length(jdd)-1))*24;
    if dthours>=0.9,  % don't want "raw" data, only data with dt>=1hour
      j=j+1;
      b{j}=dap{i};
    end
  end
end
if j==0,return,end
dap=b;

%% loop through Time Series locations
for i=1:length(dap);
  disp(sprintf('Interpolating:(%i/%i):%s',i,length(dap),char(dap{i})))
  nci=ncgeodataset(dap{i});
  % find u and v variables.  We know they are in here because we
  % used 'eastward_seawater_velocity' for our opensearch query.
  % So loop through all variables looking for variables that match
  % the standard_names.
  vars=nci.variables;
  for j=1:length(vars);
    std_name=nci.attribute(vars{j},'standard_name');
    if strcmp(std_name,u_std_name)
      vart=nci.geovariable(vars{j});
      jds=vart.timewindowij(jdstart,jdstop);
      ui=vart.data(jds.index,1,1,1);
      % convert to m/s
      units=deblank(vart.attribute('units'));
      ui = ncunits(ui, units, 'm/s');
    elseif strcmp(std_name,v_std_name)
      vart=nci.geovariable(vars{j});
      jds=vart.timewindowij(jdstart,jdstop);
      vi=vart.data(jds.index,1,1,1);
      % convert to m/s
      units=deblank(vart.attribute('units'));
      vi = ncunits(vi, units, 'm/s');
    else
    end
  end
  %jdi=jds.time(1):1/24:jds.time(end);
  jdi=jds.time;
  loni=double(nci.data('lon'));
  lati=double(nci.data('lat'));
  wd=double(nci.attribute('WATER_DEPTH'));
  if isempty(wd),
    wdepth_obs(i)=nan;
  else
    wdepth_obs(i)=-wd;
  end
  dep_obs(i)=-double(nci.data('depth'));
  depi=nci.data('depth');
  
  % find FVCOM element center closest to moored lon,lat value
  indc=nearxy(lonc(:),latc(:),loni,lati);
  iele_mod(i)=indc;
  inodes=nv(indc,:);  % nodes surrounding this element
  for kk=1:3,
    zz(:,kk)=zvar.data(inodes(kk))*sigvar.data(:,inodes(kk));  % ignore zeta contrib to z 
    ww(kk)=-double(zvar.data(inodes(kk)));
  end
  wdepth_mod(i)=double(mean(ww(kk)));
  z=mean(zz,2);
  
  % if difference between dep_obs and wdepth_obs is less than 5 m, then
  % take observation from the height above model bottom
  % extract model data from closest layer to instrument depth
  
  if (dep_obs(i)-wdepth_obs(i)) < 5,
    k=near(z-wdepth_mod(i),dep_obs(i)-wdepth_obs(i));
  else
    k=near(z,dep_obs(i));
  end
  k_mod(i)=k;
  
  u_var_mod=nc.geovariable('u');
  v_var_mod=nc.geovariable('v');
  jdmod=u_var_mod.timewindowij(jdi(1),jdi(end));
  
  umod = u_var_mod.data(jdmod.index,k,indc); %extract data from model 
  vmod = v_var_mod.data(jdmod.index,k,indc);
  %remove duplicate time records from model results
  idup=find(diff(jdmod.time)==0);
  umod(idup)=[];  vmod(idup)=[];  jdmod.time(idup)=[];
  % interpolate onto observational time base
  umod=interp1(jdmod.time,umod,jdi);
  vmod=interp1(jdmod.time,vmod,jdi);
  
  % represent velocity as complex number (u + i*v)
  wmod{i}=complex(umod,vmod);
  wobs{i}=complex(ui,vi);
  dep_mod(i)=z(k);
  lon_obs(i)=loni;
  lat_obs(i)=lati;
  lon_mod(i)=lonc(indc);
  lat_mod(i)=latc(indc);
  jd{i}=jdi;
 
end
toc
save month.mat dap lon_mod lat_mod dep_mod lon_obs lat_obs dep_obs jd jdm wobs wmod k_mod iele_mod q jdstart jdstop