function [] = oct_add_tidal_data(tidename,grdname,frcname,Ntides,tidalrank,...
                             Yorig,year,month,coastfileplot, ...
                             makeplot,pot_tides, ...
                             sal_tides,salname)
% oct_add_tidal_data: procedural netcdf replacement for add_tidal_data

interp_method = 'linear';
Roa           = 0.;
makeplot      = 0;

if nargin < 12
  sal_tides=0;
end

% Get start time
date_mjd=mjd(year,month,1);
[pf,pu,t0,phase_mkB]=egbert_correc(date_mjd,0.,0.,0.);
deg=180.0/pi;
rad=pi/180.0;

% phase correction
t0=t0-24.*mjd(Yorig,1,1);

% Read CROCO grid using octave-netcdf
ng = netcdf.open(grdname,'NC_NOWRITE');
lonr = netcdf.getVar(ng, netcdf.inqVarID(ng,'lon_rho'));
latr = netcdf.getVar(ng, netcdf.inqVarID(ng,'lat_rho'));
lonu = netcdf.getVar(ng, netcdf.inqVarID(ng,'lon_u'));
latu = netcdf.getVar(ng, netcdf.inqVarID(ng,'lat_u'));
lonv = netcdf.getVar(ng, netcdf.inqVarID(ng,'lon_v'));
latv = netcdf.getVar(ng, netcdf.inqVarID(ng,'lat_v'));
rangle = netcdf.getVar(ng, netcdf.inqVarID(ng,'angle'));
h = netcdf.getVar(ng, netcdf.inqVarID(ng,'h'));
rmask = netcdf.getVar(ng, netcdf.inqVarID(ng,'mask_rho'));
netcdf.close(ng);

% Read TPXO file
nctides = netcdf.open(tidename,'NC_NOWRITE');
periods = netcdf.getVar(nctides, netcdf.inqVarID(nctides,'periods'));
% components may be a char array stored as string variable
try
  cmpt = netcdf.getVar(nctides, netcdf.inqVarID(nctides,'components'));
catch
  cmpt = '';
end
netcdf.close(nctides);
Nmax = length(periods);
Ntides = min([Nmax Ntides]);

% Prepare forcing file
components = '';
for i=1:Ntides
  components = [components, cmpt(3*tidalrank(i)-2:3*tidalrank(i)-1), ' '];
end

oct_nc_add_tides(frcname,Ntides,date_mjd,components);

% Open forcing file for write
ncfrc = netcdf.open(frcname,'NC_WRITE');

% Loop on tidal components
for itide=1:Ntides
  it=tidalrank(itide);
  disp(['Processing tide : ',num2str(itide),' of ',num2str(Ntides)])
  % write tide_period
  varid_tp = netcdf.inqVarID(ncfrc,'tide_period');
  netcdf.putVar(ncfrc,varid_tp,itide-1,periods(it));

  correc_amp=pf(it);
  correc_phase=-phase_mkB(it)-pu(it)+360.*t0./periods(it);

  % SSH
  disp('  ssh...')
  ur = oct_ext_data_tpxo(tidename,'ssh_r',it,lonr,latr,'r',Roa);
  ui = oct_ext_data_tpxo(tidename,'ssh_i',it,lonr,latr,'r',Roa);
  ei = complex(ur,ui);
  varid_Ephase = netcdf.inqVarID(ncfrc,'tide_Ephase');
  varid_Eamp = netcdf.inqVarID(ncfrc,'tide_Eamp');
  count = [1 size(lonr,1) size(lonr,2)];
  start = [itide-1 0 0];
  netcdf.putVar(ncfrc,varid_Ephase,start,count,mod(-deg*angle(ei)+correc_phase,360.0));
  netcdf.putVar(ncfrc,varid_Eamp,start,count,abs(ei)*correc_amp);

  % U
  disp('  u...')
  ur = oct_ext_data_tpxo(tidename,'u_r',it,lonr,latr,'u',Roa);
  ui = oct_ext_data_tpxo(tidename,'u_i',it,lonr,latr,'u',Roa);
  ei = complex(ur,ui);
  upha = mod(-deg*angle(ei)+correc_phase,360.0);
  uamp = abs(ei)*correc_amp;

  % V
  disp('  v...')
  ur = oct_ext_data_tpxo(tidename,'v_r',it,lonr,latr,'v',Roa);
  ui = oct_ext_data_tpxo(tidename,'v_i',it,lonr,latr,'v',Roa);
  ei = complex(ur,ui);
  vpha = mod(-deg*angle(ei)+correc_phase,360.0);
  vamp = abs(ei)*correc_amp;

  % Ellipse
  disp('Convert to tidal ellipse parameters...')
  [major,eccentricity,inclination,phase] = ap2ep(uamp,upha,vamp,vpha);
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Cmin'), start, count, major.*eccentricity);
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Cmax'), start, count, major);
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Cangle'), start, count, inclination);
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Cphase'), start, count, phase);

  if pot_tides
    disp('Process equilibrium tidal potential...')
    if periods(it)<13.0
      Pamp = correc_amp*A(it)*B(it)*cos(rad*latr).^2;
      Ppha = mod(-2.*lonr+correc_phase,360.0);
    elseif periods(it)<26.0
      Pamp = correc_amp*A(it)*B(it)*sin(2.*rad*latr);
      Ppha = mod(-lonr+correc_phase,360.0);
    else
      Pamp = correc_amp*A(it)*B(it)*(1-1.5*cos(rad*latr).^2);
      Ppha = mod(correc_phase,360.0);
    end
    if sal_tides & it<9
      disp('Process tidal loading and self-attraction potential...')
      [SALamp,SALpha] = oct_ext_data_sal(grdname,salname,'tide_SALamp','tide_SALpha',it);
      SALamp = SALamp*correc_amp;
      SALpha = mod(SALpha+correc_phase,360.0);
      Ptot = Pamp.*exp(1i*Ppha*rad) + SALamp.*exp(1i*SALpha*rad);
      Pamp = abs(Ptot);
      Ppha = deg*angle(Ptot);
      Ppha(Pamp<0.0001)=0.;
      Ppha = mod(Ppha,360.0);
    end
    netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Pamp'), start, count, Pamp);
    netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Pphase'), start, count, Ppha);
  end
end

netcdf.close(ncfrc);

% Make plots (if requested)
if makeplot==1
  warning off
  if Ntides>=5, tiderg=[1 5]; else tiderg=1; end
  for itide=tiderg
    figure
    oct_plot_tide(grdname,frcname,itide,0.5,2,coastfileplot)
  end
  warning on
  if pot_tides
    ncfrc = netcdf.open(frcname,'NC_NOWRITE');
    domaxis=[min(min(lonr)) max(max(lonr)) min(min(latr)) max(max(latr))];
    Pamp = squeeze(netcdf.getVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Pamp'), 0, [1 size(lonr,1) size(lonr,2)]) );
    figure
    m_proj('mercator','lon',[domaxis(1) domaxis(2)],'lat',[domaxis(3) domaxis(4)]);
    m_pcolor(lonr,latr,Pamp);
    shading flat; colorbar
    m_gshhs_i('color','r'); m_gshhs_i('patch','k'); m_grid('box','fancy','xtick',5,'ytick',5,'fontsize',7);
    title('Potential Tides: M2 Amplitude [m]','fontsize',16)
    netcdf.close(ncfrc);
  end
end
end
