function oct_getpot(clmname,grdname)
% OCT_GETPOT compute potential temperature using octave-netcdf
ncid_g = netcdf.open(grdname,'NC_NOWRITE');
h = netcdf.getVar(ncid_g,netcdf.inqVarID(ncid_g,'h'));
netcdf.close(ncid_g);

ncid = netcdf.open(clmname,'NC_WRITE');
theta_s = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_s'));
theta_b = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_b'));
hc  = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'hc'));
[~,N] = netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'s_rho'));
% Determine tlen if present
try
  [~,tlen] = netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'tclm_time'));
catch
  tlen = 1;
end

try
  vtransform = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'Vtransform'));
catch
  vtransform = 1;
end

P = -1e-4*1025*9.81*zlevs(h,0.*h,theta_s,theta_b,hc,N,'r',vtransform);

for l=1:tlen
  disp(['   getpot: Time index: ',num2str(l),' of total: ',num2str(tlen)])
  % read temp and salt for time index l
  varid_temp = netcdf.inqVarID(ncid,'temp');
  varid_salt = netcdf.inqVarID(ncid,'salt');
  T = netcdf.getVar(ncid,varid_temp, [0 0 0 l-1], [Inf Inf Inf 1]);
  S = netcdf.getVar(ncid,varid_salt, [0 0 0 l-1], [Inf Inf Inf 1]);
  T = squeeze(T);
  S = squeeze(S);
  % compute potential temperature
  Tout = theta(S,T,P);
  % write back
  netcdf.putVar(ncid,varid_temp, [0 0 0 l-1], size(Tout), Tout);
end

netcdf.close(ncid);
return
