function oct_time_series(hisfile,gridfile,lon0,lat0,vname,vlevel,coef)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% oct_time_series: copy of time_series with oct_ prefix
%
% Defaults values
%
if nargin < 1
  error('You must specify a file name')
end
if nargin < 2
  gridfile=hisfile;
  disp(['Default grid name: ',gridfile])
end
if nargin < 3
  lon0=[];
end
if nargin < 4
  lat0=[];
end
if nargin < 5
  vname='temp';
  disp(['Default variable to plot: ',vname])
end
if nargin < 6
  vlevel=-10;
  disp(['Default vertical level: ',num2str(vlevel)])
end
if nargin < 7
  coef=1;
  disp(['Default coef: ',num2str(coef)])
end
%
% Get default values
%
if isempty(gridfile)
  gridfile=hisfile;
end
if vname(1)=='u'
  [lat,lon,mask]=oct_read_latlonmask(gridfile,'u');
elseif vname(1)=='v'
  [lat,lon,mask]=oct_read_latlonmask(gridfile,'v');
else
  [lat,lon,mask]=oct_read_latlonmask(gridfile,'r');
end
if isempty(lon0) | isempty(lat0)
  lat0=mean(mean(lat));
  lon0=mean(mean(lon));
end
%
% Find j,i indices for the profile
%
disp(['lon0 = ',num2str(lon0),' - lat0 = ',num2str(lat0)])
[J,I]=find((lat(1:end-1,1:end-1)<=lat0 & lat(2:end,2:end)>lat0 &...
            lon(2:end,1:end-1)<=lon0 & lon(1:end-1,2:end)>lon0)==1);
if isempty(I) |  isempty(J)
  disp('Warning: profile place not found')
  [M,L]=size(lon);
  I=round(L/2);
  J=round(M/2);
end
disp(['I = ',int2str(I),' J = ',int2str(J)])
lon1=lon(J,I);
lat1=lat(J,I);
disp(['lon1 = ',num2str(lon1),' - lat1 = ',num2str(lat1)])
%
% get the vertical levels
%
% Open file using native oct_netcdf interface
ncid = oct_netcdf.open(hisfile,'NC_NOWRITE');

surface_vars = {'zeta','ubar','vbar','sustr','svstr','shflux',... 
                'swflux','shflx_rsw','swrad','shflx_rlw','shflx_sen',... 
                'shflx_lat','sst_skin'};
if any(strcmp(vname,surface_vars))
  try
    tmp = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,vname));
  catch
    error(['Variable ',vname,' not found in ',hisfile])
  end
  var = coef*squeeze(tmp(:,J,I));
elseif vlevel>0
%
% Read the variable
%
  if vname(1)=='*'
    if strcmp(vname,'*Ke') || strcmp(vname,'*Speed')
      u_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'u'));
      v_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'v'));
      u = mean(squeeze(u_full(:,vlevel,J,I-1:I)),2);
      v = mean(squeeze(v_full(:,vlevel,J-1:J,I)),2);
      if strcmp(vname,'*Ke')
        var = coef.*0.5.*(u.^2+v.^2);
      else
        var = coef.*sqrt(u.^2+v.^2);
      end
    elseif strcmp(vname,'*Rho_pot')
      temp_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'temp'));
      salt_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'salt'));
      temp = squeeze(temp_full(:,vlevel,J,I));
      salt = squeeze(salt_full(:,vlevel,J,I));
      var = coef*rho_pot(temp,salt);
    elseif strcmp(vname,'*Chla')
      sphyto_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'SPHYTO'));
      lphyto_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'LPHYTO'));
      sphyto = squeeze(sphyto_full(:,vlevel,J,I));
      lphyto = squeeze(lphyto_full(:,vlevel,J,I));
      theta_m  = 0.020;
      CN_Phyt  = 6.625;
      var = coef*theta_m*(sphyto+lphyto)*CN_Phyt*12.;
      var(var<=0)=NaN;
    else
      disp('Sorry not implemented yet')
      oct_netcdf.close(ncid);
      return
    end
  else
    tmp = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,vname));
    var = coef*squeeze(tmp(:,vlevel,J,I));
  end
else
  % Open grid file with native oct_netcdf
  ngid = oct_netcdf.open(gridfile,'NC_NOWRITE');
  if vname(1)=='u'
    zeta_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'zeta'));
    zeta = mean(squeeze(zeta_full(:,J,I:I+1)),2);
    h_full = oct_netcdf.getVar(ngid,oct_netcdf.inqVarID(ngid,'h'));
    h = mean(squeeze(h_full(J,I:I+1)),2);
  elseif vname(1)=='v'
    zeta_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'zeta'));
    zeta = mean(squeeze(zeta_full(:,J:J+1,I)),2);
    h_full = oct_netcdf.getVar(ngid,oct_netcdf.inqVarID(ngid,'h'));
    h = mean(squeeze(h_full(J:J+1,I)),1);
  else
    zeta_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'zeta'));
    zeta = squeeze(zeta_full(:,J,I));
    h_full = oct_netcdf.getVar(ngid,oct_netcdf.inqVarID(ngid,'h'));
    h = squeeze(h_full(J,I));
  end
  oct_netcdf.close(ngid)
  % Read vertical coordinate parameters (robustly)
  try
    theta_s = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'theta_s'));
  catch
    theta_s = [];
  end
  if isempty(theta_s)
    try
      theta_s = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'theta_s'));
      theta_b = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'theta_b'));
      Tcline  = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'Tcline'));
    catch
      theta_s = [];
    end
  else
    try
      theta_b = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'theta_b'));
      Tcline  = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'Tcline'));
    catch
      Tcline = [];
    end
  end
  if isempty(Tcline)
    try
      hc = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'hc'));
    catch
      hc = [];
    end
  else
    hmin = min(min(h));
    hc = min(hmin,double(Tcline));
  end
  % number of s_rho
  try
    s_rho = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'s_rho'));
    Nr = length(s_rho);
  catch
    % fallback: inquire dimension 's_rho'
    try
      [~, Ndims, ~, ~] = oct_netcdf.inqVar(ncid, oct_netcdf.inqVarID(ncid,'temp'));
      Nr = Ndims;
    catch
      Nr = [];
    end
  end
  s_coord = 1;
  VertCoordType = [];
  try
    VertCoordType = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'VertCoordType'));
  catch
  end
  if isempty(VertCoordType)
    try
      vtrans = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'Vtransform'));
      if ~isempty(vtrans), s_coord = vtrans; end
    catch
    end
  elseif strcmp(char(VertCoordType(:)'),'NEW')
    s_coord = 2;
  end
  if s_coord==2
    hc = Tcline;
  end
%
  if vname(1)=='*'
    if strcmp(vname,'*Ke') || strcmp(vname,'*Speed')
      u_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'u'));
      v_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'v'));
      u = mean(squeeze(u_full(:,:,J,I-1:I)),3);
      v = mean(squeeze(v_full(:,:,J-1:J,I)),3);
      if strcmp(vname,'*Ke')
        var2 = coef.*0.5.*(u.^2+v.^2);
      else
        var2 = coef.*sqrt(u.^2+v.^2);
      end
    elseif strcmp(vname,'*Rho_pot')
      temp_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'temp'));
      salt_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'salt'));
      temp = squeeze(temp_full(:,:,J,I));
      salt = squeeze(salt_full(:,:,J,I));
      var2 = coef*rho_pot(temp,salt);
    elseif strcmp(vname,'*Chla')
      sphyto_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'SPHYTO'));
      lphyto_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'LPHYTO'));
      sphyto = squeeze(sphyto_full(:,:,J,I));
      lphyto = squeeze(lphyto_full(:,:,J,I));
      theta_m  = 0.020;
      CN_Phyt  = 6.625;
      var2 = coef*theta_m*(sphyto+lphyto)*CN_Phyt*12.;
      var2(var2<=0)=NaN;
    else
      disp('Sorry not implemented yet')
      oct_netcdf.close(ncid);
      return
    end
  else
    tmp2 = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,vname));
    var2 = coef*squeeze(tmp2(:,:,J,I));
  end
%
  [T,N]=size(var2);
  if N==Nr+1
    type='w';
  else
    type='r';
  end
%
  var=0*(1:T);
%
  for l=1:T
    Z=squeeze(zlevs(h,squeeze(zeta(l,:,:)),theta_s,theta_b,hc,Nr,type,s_coord));
    var(l)=interp1(Z,var2(l,:),vlevel);
  end
end
%
% Get the time
% Attempt to read scrum_time with native oct_netcdf
try
  time_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'scrum_time'));
  time = double(time_full)/(24*3600);
catch
  time = [];
end
oct_netcdf.close(ncid)
%
plot(time,var,'k')
hold on
plot(time,var,'r.')
hold off
xlabel('Time [days]')
if strcmp(vname,'zeta') | strcmp(vname,'ubar') | strcmp(vname,'vbar') 
  ylabel([vname])
elseif vlevel>0
  ylabel([vname,' - level = ',num2str(vlevel)])
else 
  ylabel([vname,' - z = ',num2str(vlevel)])
end
return