%load model mesh in 
%
% Saves out several .mat files with geometric information
% 
% mesh.mat contains
%   M the number of nodes, on which all scalers are stored (T,S,eta, H etc)
%   N the number of triangles, roughly twice M
%     also the number of U,V points, each on the centroid of the triangle
%   KBM1 the number of sigma levels on which U,V,T,S etc are stored.
%   
%   and the structure "mesh", whose feilds are
%   mesh.surround(N,3) the U,V (or centroid) points clockwise arround 
%                      a U,V point, called NBE in model and paper
%   mesh.trinodes(N,3) the t/s nodes surrounding a triangle, in clockwise order
%                      called NV in model and N(j) in paper
%   mesh.edges(M,6)  the up to 6 other nodes that each node is connected
%                    to. NaN if the edge does not exist.
%   mesh.nedges(M,1) the number of edges each node is connected to.
%   mesh.nodexy(M,2) the x and y position of the nodes, x is 1 coloumn
%   mesh.depth(M) the depth
%   mesh.uvnode(N,2) the x and y positions of the centroids of the triangles
%                    at which U and V are defined. (sorry for the poor name)
%   mesh.nodez(KBM1,M) the depth of each scaler, assuming free surface=0
%   mesh.zuv(KBM1,N)   the depth of each U,V point, assuming free surface=0
%   mesh.sigvec(KBM1) the sigma coordinate grid
%
% interpcoef.mat contains
% the interpolation coefficients to interpolate u,v and the scalers
% from where they are stored to any surrounding point. see the routines 
% window_sect_inf_make.f and window_sect_data_make.f for information, or
% secflux.m for how they are used.
%
% depthgrd.mat contains depth gridded onto a course rectangular grid
%    depthgrd.xvec is x dimension of grid
%    depthgrd.yvec is y dimension of grid
%    depthgrid.dpth is the depth
%
% depthgrd_fine.mat contains depth gridded onto a finer rectangular grid
%
%

clear all

%get M, the number of nodes, and N, the number of triangles
 filehist='../fvcom30oct/out/gom_history.cdf';
 nchist=netcdf(filehist,'nowrite');  
 t=nchist{'t1'};
 M=nchist{'m'}(1); %number of nodes
 N=nchist{'n'}(1); %number of triangles
 KBM1=size(t,2);
 close(nchist)

%get mesh data
 filemesh='../fvcom30oct/gom_grd.dat';
 fid=fopen(filemesh,'r');

%skip triangle data -- this will be read in below
 mesh.trinodes=zeros(N,3);
 for n=1:N
   line=fgetl(fid);
   jnk=sscanf(line,'%f %f %f %f %f');
 end

%read in node positions
 mesh.nodexy=zeros(M,2);
 for m=1:M
   line=fgetl(fid);
   jnk=sscanf(line,'%f %f %f %f');
   mesh.nodexy(m,:)=jnk(2:3);
 end
 fclose(fid);

%get connectivity data 
 filemesh='../fvcom01/outsms/mesh.inf';
 fid=fopen(filemesh,'r');
 fprintf('reading %s\n',fgetl(fid));

%read triangle data.  This data is used, and not the data in gom_grd, 
%because this data is in the clockwise order specified in the Chen et al. 
%paper
 mesh.trinodes=zeros(N,3);
 for n=1:N
   line=fgetl(fid);
   jnk=sscanf(line,'%f %f %f %f %f %f');
   mesh.trinodes(n,:)=jnk(4:6);
   mesh.surround(n,:)=jnk(1:3);
 end
 fclose(fid);

%get data on how each node is connected to three other nodes
 edges=zeros(M,6)*nan;
 nedges=zeros(M,1); %number of edges already found
 for n=1:N
   n;
   [edges,nedges]=insert_edge(mesh.trinodes(n,1),mesh.trinodes(n,2),...
       edges,nedges);
   [edges,nedges]=insert_edge(mesh.trinodes(n,2),mesh.trinodes(n,3),...
       edges,nedges);
   [edges,nedges]=insert_edge(mesh.trinodes(n,3),mesh.trinodes(n,1),...
       edges,nedges);
   [edges,nedges]=insert_edge(mesh.trinodes(n,2),mesh.trinodes(n,1),...
       edges,nedges);
   [edges,nedges]=insert_edge(mesh.trinodes(n,3),mesh.trinodes(n,2),...
       edges,nedges);
   [edges,nedges]=insert_edge(mesh.trinodes(n,1),mesh.trinodes(n,3),...
       edges,nedges);
 end
 mesh.edges=edges;
 mesh.nedges=nedges;

%get bathy data, and adjust as in model
 filemesh='../fvcom01/outsms/gom_dep.xy';
 fid=fopen(filemesh,'r');
 fprintf('reading \n%s\n',fgetl(fid));
 fprintf('reading %s\n',fgetl(fid));
 depth=zeros(M,1);
 for m=1:M
   line=fgetl(fid);
   jnk=sscanf(line,'%f %f %f');
   depth(m)=jnk(3);
 end
 mesh.depth=depth;
 fclose(fid);

%form u,v locations (the centroid of triangle)
 uvnodes=zeros(N,2);
 uvdepth=zeros(N,1);
 for n=1:N
   uvnode(n,1)=(mesh.nodexy(mesh.trinodes(n,1),1)+...
       mesh.nodexy(mesh.trinodes(n,2),1)+...
       mesh.nodexy(mesh.trinodes(n,3),1))/3;
   uvnode(n,2)=(mesh.nodexy(mesh.trinodes(n,1),2)+...
       mesh.nodexy(mesh.trinodes(n,2),2)+...
       mesh.nodexy(mesh.trinodes(n,3),2))/3;
   uvdepth(n)=(mesh.depth(mesh.trinodes(n,1))+...
       mesh.depth(mesh.trinodes(n,2))+...
       mesh.depth(mesh.trinodes(n,3)))/3;
 end
 mesh.uvnode=uvnode;
 mesh.uvdepth=uvdepth;


%create a full xyz list of node locations
%
%for node, their are KBM1 depths, so we must include all these depths in ...
%the list of locations
 zall=repmat((1:KBM1)',1,M);
 zuv=repmat((1:KBM1)',1,N);

%make z in depth coordinates, not the discrete equivalent of sigman coords.
 dsig=1/(KBM1+1);
 sigvec=(0:(KBM1-1))'*dsig+0.5*dsig;
 for m=1:M
   zall(:,m)=-sigvec*mesh.depth(m);
 end
 for n=1:N
   zuv(:,n)=-sigvec*mesh.uvdepth(n);
 end
 mesh.nodez=zall;
 mesh.zuv=zuv;
 mesh.sigvec=sigvec;
 
 save mesh mesh M N KBM1

%load the interpolation coefficients to interpolate u,v and the scalers
%from where they are stored to any surrounding point. see the routines 
%window_sect_inf_make.f and window_sect_data_make.f for information
 filemesh='../fvcom01/outsms/shape.inf';
 fid=fopen(filemesh,'r');

 interpcoef.au=zeros(N,4);
 interpcoef.av=zeros(N,4);
 fprintf('reading %s\n',fgetl(fid));
 for n=1:N
   line=fgetl(fid);
   jnk=sscanf(line,'%f %f %f %f %f %f %f %f');
   interpcoef.au(n,:)=jnk(1:4);
   interpcoef.av(n,:)=jnk(5:8);
 end
 fclose(fid);

 filemesh='../fvcom01/outsms/awxcof.inf';
 fid=fopen(filemesh,'r');

 interpcoef.aw0=zeros(N,3);
 interpcoef.awx=zeros(N,3);
 interpcoef.awy=zeros(N,3);
 fprintf('reading %s\n',fgetl(fid));
 for n=1:N
   line=fgetl(fid);
   jnk=sscanf(line,'%f %f %f %f %f %f %f %f %f');

   interpcoef.aw0(n,1)=jnk(1);
   interpcoef.awx(n,1)=jnk(2);
   interpcoef.awy(n,1)=jnk(3);
   
   interpcoef.aw0(n,2)=jnk(4);
   interpcoef.awx(n,2)=jnk(5);
   interpcoef.awy(n,2)=jnk(6);
   
   interpcoef.aw0(n,3)=jnk(7);
   interpcoef.awx(n,3)=jnk(8);
   interpcoef.awy(n,3)=jnk(9);
   
 end
 fclose(fid);

 save interpcoef interpcoef

%the following computes the depth gradient on u/v points, and the
%bottom drag coefficient at each u/v point.  NOTE THAT THE LATER DEPENDS
%ON THE CODE, AND IS NOT EXACTLY CORRECT BELOW.

%compute depth gradient
hgrad.x=zeros(N,1);
hgrad.y=zeros(N,1);
for n=1:N
  %compute gradient, by representing the depth surrounding 
  %the u,v point as a plane a*x+b*y+c=h
  x1=mesh.nodexy(mesh.trinodes(n,1),1);
  x2=mesh.nodexy(mesh.trinodes(n,2),1);
  x3=mesh.nodexy(mesh.trinodes(n,3),1);
  y1=mesh.nodexy(mesh.trinodes(n,1),2);
  y2=mesh.nodexy(mesh.trinodes(n,2),2);
  y3=mesh.nodexy(mesh.trinodes(n,3),2);
  h1=mesh.depth(mesh.trinodes(n,1));
  h2=mesh.depth(mesh.trinodes(n,2));
  h3=mesh.depth(mesh.trinodes(n,3));
 
  %maple rocks!
  a = -(-y3*h2+y1*h2-y1*h3-y2*h1+y3*h1+y2*h3)/...
      (y2*x1-x3*y2+x3*y1+x2*y3-y1*x2-x1*y3);
  c = (-x3*y2*h1+x3*y1*h2+x2*y3*h1+y2*x1*h3-y1*x2*h3-x1*y3*h2)...
      /(y2*x1-x3*y2+x3*y1+x2*y3-y1*x2-x1*y3);
  b = (x3*h1+x1*h2-x2*h1-x1*h3-x3*h2+x2*h3)/(y2*x1-x3*y2+x3*y1+x2*y3-y1*x2-x1*y3);

  hgrad.x(n)=a;
  hgrad.y(n)=b;
end

%compute Cd, the bottom drag coefficient, on u,v points
%note that the calculation is not exact, because in the model the 
%free surface elevation is included in the depth, and here it is not.
fprintf('\nNOTE THAT BOTTOM DRAG COEF CALCULATION\n');
fprintf('MUST MATCH THAT IN BOTTOM_ROUGHNESS IN FVCOM.F\n\n');
Cd=zeros(N,1);
for i=1:N
  d1=mesh.uvdepth(i); %not exactly correct - should be h+el
  if(d1<=40.) 
    z0=3.e-3;
  elseif (d1>40.0) & (d1<=70.0) 
    z0=3.e-3*exp(-(d1-40.)/8.8204);
  elseif (d1>70.0) & (d1<=100.0) 
    z0=1.e-4*exp(-(d1-70.)/13.0288);
  else
    z0=1.e-5;
  end 

  ztemp=(0.5/KBM1)*d1/z0;   %ONLY CORRECT FOR EVEN SIGMA SPACEING!
  cbc = .16/(log(ztemp))^2;

  Cd(i)=cbc;
end  

save hgrad hgrad
save Cd Cd
    
%the following produces depth on a 2D regular grid, and the arrays 
%needed to interpolate T/S and other scalers onto a 3D regular grid,
%for ease of visualization.

%define the regular grid
 dx=3000;
 xvec=507948:dx:1774300;
 yvec=-656500:dx:352646;
 zvec=-300:10:0;
 
%make the 2D depth array, and save
 dpth=griddata(mesh.nodexy(:,1),mesh.nodexy(:,2),mesh.depth,xvec,yvec');
 depthgrd.xvec=xvec;
 depthgrd.yvec=yvec;
 depthgrd.dpth=dpth;
 save depthgrd depthgrd

%define the regular grid
 dx=1000;
 xvec=507948:dx:1774300;
 yvec=-656500:dx:352646;
 zvec=-300:10:0;
 
%make the 2D depth array, and save
 dpth=griddata(mesh.nodexy(:,1),mesh.nodexy(:,2),mesh.depth,xvec,yvec');
 depthgrd.xvec=xvec;
 depthgrd.yvec=yvec;
 depthgrd.dpth=dpth;
 save depthgrd_fine depthgrd


 
 
%plot nodes
plot(mesh.nodexy(:,1)/1e3,mesh.nodexy(:,2)/1e3,'b.',...
     mesh.uvnode(:,1)/1e3,mesh.uvnode(:,2)/1e3,'r.');

hold on
for m=1:M
  for ne=1:nedges(m)
    plot([mesh.nodexy(m,1) mesh.nodexy(edges(m,ne),1)],...
	[mesh.nodexy(m,2) mesh.nodexy(edges(m,ne),2)],'k-')
  end
end
hold off