% Converts wind data in NetCDF format to SWAN ASCII wind forcing files on a ROMS grid
%
% These files will be created in a sub-directory called "uvswan"

% Define netcdf forcing file for wind
filename='/d2/rps/meteo/coamps/coamps_wind.nc';

% Output SWAN forcing files for this period:
start=[2002 9 10 0 0 0];  
stop=[2003 6 18 0 0 0];

% read NCOM/SWAN grid information
%load /d3/rps/meteo/lami/whole/ncom_grid.mat lon lat depth
%load ncom_grid.mat lon lat depth
nc=netcdf('/d2/rps/coarse_input/adria02_grid4a.nc');
lon=nc{'lon_rho'}(:);
lat=nc{'lat_rho'}(:);
mask=nc{'mask_rho'}(:);
depth=nc{'h'}(:);
ang=nc{'angle'}(:);  % radians, even though grid file says "degrees"
close(nc)
rlon=lon;
rlat=lat;
d=depth;
iland=find(mask==0);
iwater=find(mask==1);

% make the SWAN "input" directory (which we will fill with wind files)
!mkdir uvswan 

% Read the wind model output
nc=netcdf(filename);
jd=nc{'time'}(:);
jd=jd+2440000;
alon=nc{'lon'}(:);
alat=nc{'lat'}(:);
g=gregorian(jd);

% Find time indices within the specified range
nind=tind(jd,start,stop);

% set restart to 1 to load in data calculated previously,
% set restart to 0 to recalculate the points inside (if either
% the MET or ROMS grids change in any way)

%restart=1;
restart=0;

% Set up to interpolate on ij, rather than lon/lat
[nx,ny]=size(alon);
[ii,jj]=meshgrid(1:ny,1:nx);

if (restart)
  load forcing_cache_lowres.mat ib jb 
else
  ib=griddata(alon,alat,ii,rlon,rlat);
  jb=griddata(alon,alat,jj,rlon,rlat);
  save forcing_cache_lowres.mat ib jb
end

%for n=1:length(jd);
for n=1:length(nind);
  nn=nind(n);
  u=nc{'U10',1}(nn,:,:);
  v=nc{'V10',1}(nn,:,:);
  w=u+sqrt(-1)*v;
  wi=zeros(size(d));
  wi(iwater)=interp2(ii,jj,w,ib(iwater),jb(iwater));
  % rotate the wind into XI, ETA coordinate directions
  wir=wi.*exp(-sqrt(-1)*ang);   %rotated into model grid coordinates
  %arot=90;
  %wir=wir.*exp(sqrt(-1)*arot*pi/180);
  wir=wir.';
  uv=[real(wir(:)) imag(wir(:))];
  fil=sprintf('uvswan2/uvswan_%4.4d%2.2d%2.2d%2.2d',g(nn,1:4))
  fid=fopen(fil,'w');
  form='%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n';
  fprintf(fid,form,uv);
  fclose(fid);
end
close(nc);