function [f,z,t]=roms_zstation(fname,vname,istation);                                  
%function [f,z,t]=roms_zstation(fname,vname,istation);
%                         %
%                                                                           %
% This routine reads in SCRUM station NetCDF file and plots a time-series   %
% of the requested field profile.                                           %
%                                                                           %
%  On Input:                                                                %
%                                                                           %
%     fname    NetCDF filename without the extension suffix (string).       %
%     vname    NetCDF variable name to read (string).                       %
%     index    Station index to process (integer).                          %
%                                                                           %
%  On Output:                                                               %
%                                                                           %
%     t        Time of field (matrix).                                      %
%     z        Depths (m) of field (matrix).                                %
%     f        Requested field (matrix).                                    %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----------------------------------------------------------------------------
%  Initialize labels.
%----------------------------------------------------------------------------

slabel='  ';
stitle='  ';

%----------------------------------------------------------------------------
% Test input to see if it's in an acceptable form.
%----------------------------------------------------------------------------

if (nargin < 3),
  disp(' ');
  disp([setstr(7),'*** Error:  ZSTATION - too few arguments.',setstr(7)]);
  disp([setstr(7),'                       number of supplied arguments: ',...
       num2str(nargin),setstr(7)]);
  disp([setstr(7),'                       number of required arguments: 3',...
       setstr(7)]);
  disp(' ');
  return
end,


%----------------------------------------------------------------------------
% open file
%-------------------------
nc=netcdf(fname);


%----------------------------------------------------------------------------
%  Read in number of grid points.
%----------------------------------------------------------------------------
N=length(nc('s_rho'));
Np=N+1;

NS=length(nc('station'));
NT=length(nc('ocean_time'));

%----------------------------------------------------------------------------
%  Check if index of extraction is whitin bounds.
%----------------------------------------------------------------------------
if (istation < 1 | istation > NS),
  disp(' ');
  disp([setstr(7),'*** Error:  ZSTATION - illegal station index.',setstr(7)]);
  disp([setstr(7),'                       valid range:  1 <= index <= ',...
       num2str(Lp),setstr(7)]);
  disp(' ');
  return
end,

%----------------------------------------------------------------------------
%  Read in vertical S-coordinate parameters.
%----------------------------------------------------------------------------

hc=nc{'hc'}(:);
sc_r=nc{'s_rho'}(:);
Cs_r=nc{'Cs_r'}(:);
sc_w=nc{'s_w'}(:);
Cs_w=nc{'Cs_w'}(:);

%----------------------------------------------------------------------------
% Read in bathymetry and free-surface (meters) at RHO-points.
%----------------------------------------------------------------------------
 
h=nc{'h'}(istation);
zeta=nc{'zeta'}(:,istation);

%----------------------------------------------------------------------------
% Read in field.
%----------------------------------------------------------------------------
f=nc{vname}(:,istation,:);
f=squeeze(f);
t=nc{'ocean_time'}(:);

%----------------------------------------------------------------------------
%  Set zstation vertical coordinate.
%----------------------------------------------------------------------------

z=zeros(NT,N);
for k=1:N,
  z(:,k)=zeta(:).*(1.0+sc_r(k))+hc.*sc_r(k)+(h-hc)*Cs_r(k);
end,

%----------------------------------------------------------------------------
% Plot station.
%----------------------------------------------------------------------------
close(nc);
return