%
% load bank150 grid
%
[in,x,y,z,bnd]=loadgrid('bank150');
[y,x]=xy2ll(x,y);
size(in);
ne=ans(1);
size(x);
nn=ans(1);
%
% plot bank150 boundary and isobaths
%
normalfigure;
bndo=plotbnd(x,y,bnd);
set(bndo,'Color','black')
cval=[60 100 150];
h=lcontour2(in,x,y,z,cval);
set(h,'Color','black')
hold on
%
% load sea soar data
% 4 m bins with centers at 13-17-21-...-201-205-209
%
load oc299_uv.mat
%
% plot (lon,lat) locations
%
plot(txy(:,2),txy(:,3),'bo');
%
% ensure u and v contain identical NaN entries
%
unan=isnan(u);
vnan=isnan(v);
max(max(unan-vnan))
%
% vertically integrate non-nan velocity observations
%
nobs=0;
for i=1:length(txy)
   for j=1:length(depth)
      if unan(i,j) == 0
         jtop=j;
         ztop=depth(jtop)-2.0;
         break
      end
   end
   if j == length(depth);
   else;
      nobs=nobs+1;
      aobs(nobs,1)=txy(i,1);
      aobs(nobs,2)=txy(i,2);
      aobs(nobs,3)=txy(i,3);
      aobs(nobs,6)=ztop;
      jbot=length(depth);
      aobs(nobs,7)=depth(jbot)+2.0;
      for j=jtop+1:length(depth)
         if unan(i,j) == 1
            jbot=j-1;
            aobs(nobs,7)=depth(jbot)+2.0;
            break
         end
      end
      usum=0.0;vsum=0.0;
      for j=jtop:jbot
         usum=usum+u(i,j);
         vsum=vsum+v(i,j);
      end
      aobs(nobs,4)=usum/(jbot-jtop+1);
      aobs(nobs,5)=vsum/(jbot-jtop+1);
   end;
end
%
% Plot depth range of observations
%
normalfigure;
plot(aobs(:,1),aobs(:,7),'b-.')
plot(aobs(:,1),aobs(:,6),'r-')
legend('Maximum Depth','Minimum Depth')
xlabel('Yearday')
%
% Plot vector field
%
normalfigure;
plot(aobs(:,2),aobs(:,3),'ro');
hold on
bndo=plotbnd(x,y,bnd);
set(bndo,'Color','black')
cval=[60 100 150];
h=lcontour2(in,x,y,z,cval);
set(h,'Color','black')
hv=vecplot2(aobs(:,2),aobs(:,3),aobs(:,4),aobs(:,5),1.0,'m/s');
set(hv,'Color','r');
%
% Horizontally filter to remove observations outside bank150
%
load bank150.elev.ncd
bndnodes=bank150(:,1);
isin=isinpoly(aobs(:,2),aobs(:,3),x(bndnodes),y(bndnodes));
iobs=0;
for i=1:nobs
   if isin(i) == 1
      iobs=iobs+1;
      aiobs(iobs,:)=aobs(i,:);
   end
end
%
% Set bathymetric depth at observational locations to model bathymetry
%
[smin,nclosest]=closestnode(x,y,aiobs(:,2),aiobs(:,3));
aiobs(:,8)=z(nclosest);
scrange(aiobs(:,8)-aiobs(:,7))
%
% Plot depth range of observations and model bathymetry
%
landscapefigure;
plot(aiobs(:,1),aiobs(:,8),'g--')
plot(aiobs(:,1),aiobs(:,7),'b-.')
plot(aiobs(:,1),aiobs(:,6),'r-')
legend('Model Bathymetry','Maximum Depth','Minimum Depth')
xlabel('Yearday')
%
% Plot vector field
%
landscapefigure;
plot(aiobs(:,2),aiobs(:,3),'ro');
hold on
bndo=plotbnd(x,y,bnd);
set(bndo,'Color','black')
cval=[60 100 150];
h=lcontour2(in,x,y,z,cval);
set(h,'Color','black')
hv=vecplot2(aiobs(:,2),aiobs(:,3),aiobs(:,4),aiobs(:,5),1.0,'m/s');
set(hv,'Color','r');
%
% Convert observational coordinates to (x,y) adjust vertical coordinates
%
[aiobs(:,2),aiobs(:,3)]=ll2xy(aiobs(:,3),aiobs(:,2));
aiobs(:,6)=-aiobs(:,6);
aiobs(:,7)=-aiobs(:,7);
aiobs(:,8)=-aiobs(:,8);
%
% Output .m3d file
%
fname='oc299_sea_soar.m3d';
fmt=' %f %f %f %f %f %f %f %f\n';
data=aiobs;
gridname='bank150';
year=1997;
ncol=8;
gmt='gmt';
write_obs(fname,fmt,data,gridname,year,ncol,gmt);

