function [X,Z,VAR]=oct_get_section(fname,gname,lonsec,latsec,vname,tindex);
% OCT_GET_SECTION Procedural NetCDF variant of get_section
interp_type='linear';

if nargin < 1
  error('You must specify a file name')
end
if nargin < 2
  gname=fname;
  disp(['Default grid name: ',gname])
end
if nargin < 3
  lonsec=[12 18];
  disp(['Default longitude: ',num2str(lonsec)])
end
if nargin < 4
  latsec=-34;
  disp(['Default latitude: ',num2str(latsec)])
end
if nargin < 5
  vname='temp';
  disp(['Default variable to plot: ',vname])
end
if nargin < 6
  tindex=1;
  disp(['Default time index: ',num2str(tindex)])
end

% Find maximum grid angle size (dl)
[lat,lon,mask]=oct_read_latlonmask(gname,'r');
[M,L]=size(lon);
dl=1.5* max([max(max(abs(lon(2:M,:)-lon(1:M-1,:)))) ...
        max(max(abs(lon(:,2:L)-lon(:,1:L-1)))) ...
        max(max(abs(lat(2:M,:)-lat(1:M-1,:)))) ...
        max(max(abs(lat(:,2:L)-lat(:,1:L-1))))]);

% Read point positions
[type,vlevel]=get_type(fname,vname,10);
if (vlevel==0)
   error([vname,' is a 2D-H variable'])
   return
end
[lat,lon,mask]=oct_read_latlonmask(gname,type);
[M,L]=size(lon);

% Find minimal subgrids limits
minlon=min(lonsec)-dl; minlat=min(latsec)-dl;
maxlon=max(lonsec)+dl; maxlat=max(latsec)+dl;
sub=lon>minlon & lon<maxlon & lat>minlat & lat<maxlat;
if (sum(sum(sub))==0)
  error('Section out of the domain')
end
ival=sum(sub,1); jval=sum(sub,2);
imin=min(find(ival~=0)); imax=max(find(ival~=0));
jmin=min(find(jval~=0)); jmax=max(find(jval~=0));

% Get subgrids
lon=lon(jmin:jmax,imin:imax); lat=lat(jmin:jmax,imin:imax);
sub=sub(jmin:jmax,imin:imax); mask=mask(jmin:jmax,imin:imax);

% Put lat/lon of section in correct vector form
if (length(lonsec)==1)
  disp(['N-S section at longitude: ',num2str(lonsec)])
  if (length(latsec)==1)
    error('Need more points to do a section')
  elseif (length(latsec)==2)
    latsec=(latsec(1):dl:latsec(2));
  end
  lonsec=0.*latsec+lonsec;
elseif (length(latsec)==1)
  disp(['E-W section at latitude: ',num2str(latsec)])
  if (length(lonsec)==2)
    lonsec=(lonsec(1):dl:lonsec(2));
  end
  latsec=0.*lonsec+latsec;
elseif (length(lonsec)==2 & length(latsec)==2)
  Npts=ceil(max([abs(lonsec(2)-lonsec(1))/dl abs(latsec(2)-latsec(1))/dl]));
  if lonsec(1)==lonsec(2)
    lonsec=lonsec(1)+zeros(1,Npts+1);
  else
    lonsec=(lonsec(1):(lonsec(2)-lonsec(1))/Npts:lonsec(2));
  end
  if latsec(1)==latsec(2)
    latsec=latsec(1)+zeros(1,Npts+1);
  else
    latsec=(latsec(1):(latsec(2)-latsec(1))/Npts:latsec(2));
  end
elseif (length(lonsec)~= length(latsec))
  error('Section latitudes and longitudes are not of the same length')
end
Npts=length(lonsec);

% Get sub mask for each section point
sub=0*lon;
for i=1:Npts
  sub(lon>lonsec(i)-dl & lon<lonsec(i)+dl & lat>latsec(i)-dl & lat<latsec(i)+dl)=1;
end

% get coefficients of objective analysis
londata=lon(sub==1); latdata=lat(sub==1);
coef=oacoef(londata,latdata,lonsec,latsec,100e3);

% Get the mask
mask=mask(sub==1);
m1=griddata(londata,latdata,mask,lonsec,latsec,'nearest');
londata=londata(mask==1); latdata=latdata(mask==1);

% Get vertical levels
ngid = netcdf.open(gname,'NC_NOWRITE');
try
  h = netcdf.getVar(ngid,netcdf.inqVarID(ngid,'h'));
catch
  h = [];
end
netcdf.close(ngid);
hmin=min(min(h));
h=h(jmin:jmax,imin:imax);
h=h(sub==1); h=h(mask==1);
h=griddata(londata,latdata,h,lonsec,latsec,interp_type);

ncid = netcdf.open(fname,'NC_NOWRITE');
% hmorph
try
  hm_tmp = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'hmorph'));
  if ndims(hm_tmp)==3, hm = squeeze(hm_tmp(tindex,:,:)); else hm = hm_tmp; end
catch
  hm = [];
end
if ~isempty(hm)
 h=hm; hmin=min(min(h)); h=h(jmin:jmax,imin:imax); h=h(sub==1); h=h(mask==1);
 h=griddata(londata,latdata,h,lonsec,latsec,interp_type);
end

% zeta
try
  zeta_tmp = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'zeta'));
  if ndims(zeta_tmp)==3
    zeta = squeeze(zeta_tmp(tindex,jmin:jmax,imin:imax));
  else
    zeta = zeta_tmp;
  end
catch
  zeta = [];
end
if isempty(zeta)
  zeta=0.*h;
else
  zeta = zeta(sub==1); zeta = zeta(mask==1);
  zeta = griddata(londata,latdata,zeta,lonsec,latsec,interp_type);
end

% vertical coords
try, theta_s = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_s')); catch, theta_s = []; end
if isempty(theta_s)
  try, theta_s = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_s')); catch, theta_s = []; end
  try, theta_b = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_b')); catch, theta_b = []; end
  try, Tcline = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'Tcline')); catch, Tcline = []; end
else
  try, theta_b = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_b')); catch, theta_b=[]; end
  try, Tcline = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'Tcline')); catch, Tcline=[]; end
end
if isempty(Tcline)
  try, hc = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'hc')); catch, hc = []; end
else
  hmin=min(min(h)); hc=min(hmin,Tcline);
end

% N and s_coord
try
  [~,N] = netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'s_rho'));
catch
  N = 0;
end
s_coord=1;
try
  VertCoordType = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'VertCoordType');
catch
  VertCoordType = '';
end
if isempty(VertCoordType)
  try vtrans = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'Vtransform')); if ~isempty(vtrans), s_coord=vtrans; end; catch, end
elseif strcmp(VertCoordType,'NEW'), s_coord=2; end
if s_coord==2, hc=Tcline; end

Z=squeeze(zlevs(h,zeta,theta_s,theta_b,hc,N,type,s_coord));
[N,Nsec]=size(Z);

% Loop vertical levels
VAR=0.*Z;
for k=1:N
  try
    var_full = netcdf.getVar(ncid,netcdf.inqVarID(ncid,vname));
    % var_full expected dims [time z y x] or [z y x]
    if ndims(var_full)==4
      var_slice = squeeze(var_full(tindex,k,jmin:jmax,imin:imax));
    elseif ndims(var_full)==3
      var_slice = squeeze(var_full(k,jmin:jmax,imin:imax));
    else
      var_slice = var_full;
    end
  catch
    var_slice = [];
  end
  var_slice = var_slice(sub==1); var_slice = var_slice(mask==1);
  var_grid = griddata(londata,latdata,var_slice,lonsec,latsec,interp_type);
  VAR(k,:) = m1 .* var_grid;
end
netcdf.close(ncid);

% Distances
dist=spheric_dist(latsec(1),latsec,lonsec(1),lonsec)/1e3;
X=squeeze(tridim(dist,N));
return
