function [rlon,rlat,Etide]=ssh_tides(Fname,tide_start,time);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [rlon,rlat,Etide]=ssh_tides(Fname,tide_start,time)               %
%                                                                           %
% This function computes the sea surface height tidal contribution for      %
% the requested time. It uses the provided tidal constituents elevation     %
% amplitudes and phase from ROMS tidal forcing file to compute the total    %
% contribution.                                                             %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    Fname       ROMS tidal forcing NetCDF                                  %
%    tide_start  Reference time origin for tidal forcing (days)             %
%    time        Requested time (days) to compute                           %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    rlon        Grid longitude at rho-points.                              %
%    rlat        Grid latitude at rho-points.                               %
%    Etide       Sea surface height elevation due to tides (m).             %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Inquire about variables.

[vname,nvars]=nc_vname(Fname);

got_period=0;
got_Tamp=0;
got_Tphase=0;
got_rlon=0;
got_rlat=0;
got_rmask=0;

for n=1:nvars
  name=deblank(vname(n,:));
  switch name
    case 'tide_period'
      got_period=1;
    case 'tide_Eamp'
      got_Tamp;
    case 'tide_Ephase'
      got_Tphase=1;
    case 'lon_rho'
      got_rlon=1;
    case 'lat_rho'
      got_rlat=1;
    case 'mask_rho'
      got_rmask=1;
  end,
end,

day2sec=84600;
hour2sec=3600.0;
deg2rad=pi/180.0;

%----------------------------------------------------------------------------
% Read in tidal data.
%----------------------------------------------------------------------------

% Read in tidal period (hours) for each available constituent.

Tperiod=nc_read(Fname,'tide_period');
Tperiod=Tperiod.*hour2sec;

% Read in tidal elevation amplitude (m) and phase (degrees) for each
% constituent.

SSH_Tamp=nc_read(Fname,'tide_Eamp');
SSH_Tphase=nc_read(Fname,'tide_Ephase');
SSH_Tphase=SSH_Tphase.*deg2rad;

%----------------------------------------------------------------------------
% Read in grid data.
%----------------------------------------------------------------------------

rlon=nc_read(Fname,'lon_rho');
rlat=nc_read(Fname,'lat_rho');

if (got_rmask),
 rmask=nc_read(Fname,'mask_rho');
else,
 rmask=ones(size(rlon));
end,

%----------------------------------------------------------------------------
% Compute sea surface height elevation due to tidal forcing.
%----------------------------------------------------------------------------

[Im,Jm,NTC]=size(SSH_Tamp);

Etide=zeros(size(squeeze(SSH_Tamp(:,:,1))));

for itide=1:NTC,
  if (Tperiod(itide) > 0),
    omega=2.0*pi*(time-tide_start)*day2sec/Tperiod(itide);
    Etide(:,:)=Etide(:,:)+ ...
               SSH_Tamp(:,:,itide).*cos(omega-SSH_Tphase(:,:,itide)).*rmask;
  end,
end,

return
end
