function oct_make_tides(tidename,grdname,frcname,Ntides,tidalrank,tidalranklist)
% oct_make_tides: procedural netcdf variant of make_tides.m
% This wrapper calls oct_ helpers where appropriate.

% For simplicity this function mirrors the original make_tides behavior
% but uses the oct_ prefixed helpers created alongside it.

% We'll assume `crocotools_param` is already in the path and sets variables.

crocotools_param

pot_tides=1;
if ~exist('sal_tides','var')
  sal_tides=0;
end

% Get start time
date_mjd=mjd(Ymin,Mmin,Dmin);
[pf,pu,t0,phase_mkB]=egbert_correc(date_mjd,Hmin,Min_min,Smin);
deg=180.0/pi; rad=pi/180.0;

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

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

% Read TPXO
nctides = netcdf.open(tidename,'NC_NOWRITE');
periods = netcdf.getVar(nctides, netcdf.inqVarID(nctides,'periods'));
cmpt = netcdf.getVar(nctides, netcdf.inqVarID(nctides,'components'));
netcdf.close(nctides);
Nmax = length(periods);
Ntides = min([Nmax Ntides]);

% prepare forcing file using oct_nc_add_tides
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');

for itide=1:Ntides
  it=tidalrank(itide);
  disp(['Processing tide : ',num2str(itide),' of ',num2str(Ntides)])
  % set tide_period variable
  varid = netcdf.inqVarID(ncfrc,'tide_period');
  netcdf.putVar(ncfrc,varid,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);
  % write variables
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Ephase'), [itide-1 0 0], [1 size(lonr,1) size(lonr,2)], mod(-deg*angle(ei)+correc_phase,360.0));
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Eamp'),   [itide-1 0 0], [1 size(lonr,1) size(lonr,2)], 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;

  % Convert to tidal ellipses
  disp('Convert to tidal ellipse parameters...')
  [major,eccentricity,inclination,phase] = ap2ep(uamp,upha,vamp,vpha);
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Cmin'),   [itide-1 0 0], [1 size(lonr,1) size(lonr,2)], major.*eccentricity);
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Cmax'),   [itide-1 0 0], [1 size(lonr,1) size(lonr,2)], major);
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Cangle'), [itide-1 0 0], [1 size(lonr,1) size(lonr,2)], inclination);
  netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Cphase'), [itide-1 0 0], [1 size(lonr,1) size(lonr,2)], 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'), [itide-1 0 0], [1 size(lonr,1) size(lonr,2)], Pamp);
    netcdf.putVar(ncfrc, netcdf.inqVarID(ncfrc,'tide_Pphase'), [itide-1 0 0], [1 size(lonr,1) size(lonr,2)], Ppha);
  end
end

netcdf.close(ncfrc);
end
