function oct_vert_profile(hisfile,gridfile,lon0,lat0,vname,tindex,coef,Yorig)
% OCT_VERT_PROFILE: oct_ prefixed copy of vert_profile using oct_ helpers

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
  tindex=1;
  disp(['Default time index: ',num2str(tindex)])
end
if nargin < 7
  coef=1;
  disp(['Default coef: ',num2str(coef)])
end
if nargin < 8
  Yorig=NaN;
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 using native oct_netcdf
ngid = oct_netcdf.open(gridfile,'NC_NOWRITE');
try
  h_full = oct_netcdf.getVar(ngid,oct_netcdf.inqVarID(ngid,'h'));
catch
  oct_netcdf.close(ngid);
  error(['h not found in ',gridfile])
end
oct_netcdf.close(ngid);
ncid = oct_netcdf.open(hisfile,'NC_NOWRITE');
try
  zeta_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'zeta'));
catch
  zeta_full = [];
end
if vname(1)=='u'
  zeta = mean(squeeze(zeta_full(tindex,J,I:I+1)),2);
  h = mean(squeeze(h_full(J,I:I+1)),2);
elseif vname(1)=='v'
  zeta = mean(squeeze(zeta_full(tindex,J:J+1,I)),2);
  h = mean(squeeze(h_full(J:J+1,I)),2);
else
  if ~isempty(zeta_full)
    zeta = squeeze(zeta_full(tindex,J,I));
  else
    zeta = 0;
  end
  h = squeeze(h_full(J,I));
end

% vertical coord parameters
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
  Nr = [];
end
s_coord = 1;
try
  VertCoordType = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'VertCoordType'));
catch
  VertCoordType = [];
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

% Read the variable
if vname(1)=='*'
  if strcmp(vname,'*Ke')
    u=mean(squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'u'))(tindex,:,J,I-1:I)),2);
    v=mean(squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'v'))(tindex,:,J-1:J,I)),2);
    var=coef.*0.5.*(u.^2+v.^2);
  elseif strcmp(vname,'*Speed')
    u=mean(squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'u'))(tindex,:,J,I-1:I)),2);
    v=mean(squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'v'))(tindex,:,J-1:J,I)),2);
    var=coef.*sqrt(u.^2+v.^2);
  elseif strcmp(vname,'*Rho')
    temp=squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'temp'))(tindex,:,J,I));
    salt=squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'salt'))(tindex,:,J,I));
    z=squeeze(zlevs(h,zeta,theta_s,theta_b,hc,Nr,'r',s_coord));
    var=coef*rho_eos(temp,salt,z);
  elseif strcmp(vname,'*Rho_pot')
    temp=squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'temp'))(tindex,:,J,I));
    salt=squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'salt'))(tindex,:,J,I));
    var=coef*rho_pot(temp,salt);
  elseif strcmp(vname,'*Chla')
    sphyto=squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'SPHYTO'))(tindex,:,J,I));
    lphyto=squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'LPHYTO'))(tindex,:,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
  try
    tmp = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,vname));
    var=coef*squeeze(tmp(tindex,:,J,I));
  catch
    oct_netcdf.close(ncid);
    error(['Variable ',vname,' not found in ',hisfile])
  end
end

N=length(var);
if N==Nr+1
  type='w';
else
  type='r';
end
Z=squeeze(zlevs(h,zeta,theta_s,theta_b,hc,Nr,type,s_coord));
oct_netcdf.close(ncid)

% Get the date
[day,month,year,imonth,thedate]=...
  oct_get_date(hisfile,tindex,Yorig);

plot(var,Z,'k')
hold on
plot(var,Z,'r.')
hold off
ylabel('Depth [m]')
xlabel([vname,' - ',thedate])
return