%  Reads WAVEWATCH SPC output files
function [lat,lon,depth,freq,dir,dates,Efth,U10,Udir,UC,Cdir,names,gridname,nk,nth,np,nvar,varnames,ndateso]=readWWSPC(filename,ndates)
% 
fid=fopen(filename,'r');
% 
% Reads header and size of arrays
%
tline=fgetl(fid);
data=textscan(tline,'%23c%d%d%d %32c',1); %(A23,X,3I6,X,A32)
%        
nk=data{2};
nth=data{3};
np=data{4};
gridname=data{5};
%
nvar=1;
varnames=['E(,)    '];
%
%if (gridname(1)~=char(34) & gridname(1)~=char(39)) ;
if (size(data{5},1)==0) ; 
    % In this case source terms are also expected
    data=textscan(tline,'%23c%d%d%d %1c %1c %1c %1c %1c %1c',1); %(A23,X,3I6,X,A32)
    nvar=6;
    varna=['E(,)    ', 'Sin(,)  ', 'Snl(,)  ' , 'Sds(,)  ', 'Sbt(,)  ','Stot(,) '];
    varnames=reshape(varna,8,6)';
    invar=zeros(6,1);
    for i=1:6
       if (data{4+i}=='T') 
           invar(i)=1;
       end
    end
    nvar=sum(invar);
    varn='';
    for i=1:6
       if (invar(i)==1) 
           varn=[varn varnames(i,:)];
       end
    end
    varnames=reshape(varn,8,nvar)';
end
freq=zeros(nk,1);
for i=1:nk
  data=textscan(fid,'%f',1);
  freq(i)=data{1};
end
%data=textscan(fid,'%f',nth);
%dir=data{1};
dir=zeros(nth,1);
for i=1:nth
  data=textscan(fid,'%f',1);
  dir(i)=data{1};
end

%
%
names=[''];
lat=zeros(np,1);
lon=zeros(np,1);
depth1=zeros(ndates,np,1);
U101=zeros(ndates,np,1);
Udir1=zeros(ndates,np,1);
UC1=zeros(ndates,np,1);
Cdir1=zeros(ndates,np,1);
Efth1=zeros(ndates,np,nvar,nk,nth);
dates1=zeros(ndates,1);
%
% Loop on times and points
%
error=1;
itime=0;
while (error > 0 & itime < ndates)
   data=textscan(fid,'%4d%2d%2d%2d%2d%2d',1);
   error=size(data{1});
   if (error > 0) 
      itime=itime+1;
      iy=double(data{1});
      im=double(data{2});
      id=double(data{3});
      ih=double(data{4});
      imin=double(data{5});
      is=double(data{6});
      dates1(itime)=datenum(iy,im,id,ih,imin,is);
      for ip=1:np
      data=textscan(fid,'%1c%10c%1c%f%f%f%f%f%f%f',1);
      if (itime==1)
          names=[names data{2}];
          lat(ip)=data{4};
          lon(ip)=data{5};
      end
      depth1(itime,ip)=data{6};
      U101(itime,ip)=data{7};
      Udir1(itime,ip)=data{8};
      UC1(itime,ip)=data{9};
      Cdir1(itime,ip)=data{10};
      for ivar=1:nvar
          data=textscan(fid,'%f',nk*nth);
          Efth1(itime,ip,ivar,:,:)=reshape(data{1},nk,nth);
      end
   end
   end
end
ndateso=itime;
depth=zeros(ndateso,np,1);
U10=zeros(ndateso,np,1);
Udir=zeros(ndateso,np,1);
UC=zeros(ndateso,np,1);
Cdir=zeros(ndateso,np,1);
Efth=zeros(ndateso,np,nvar,nk,nth);
dates=zeros(ndateso,1);
%
% Reduces the time dimension to ndateso
%
depth(1:ndateso,:,:)=depth1(1:ndateso,:,:);
U10(1:ndateso,:,:)=U101(1:ndateso,:,:);
Udir(1:ndateso,:,:)=Udir1(1:ndateso,:,:);
UC(1:ndateso,:,:)=UC1(1:ndateso,:,:);
Cdir(1:ndateso,:,:)=Cdir1(1:ndateso,:,:);
Efth(1:ndateso,:,:,:,:)=Efth1(1:ndateso,:,:,:,:);
dates(1:ndateso)=dates1(1:ndateso);
if (ndates==ndateso) 
    fprintf('%s %d %s \n','WARNING: file may contain more than ',ndates,' time steps.');
end
fclose(fid);
