% Interpolate a roms dataset from one bathymetry to another
%
% Assumes the horizontal grid doesn't alter
% just need to (a) do the vertical interpolation to new s-levels 
%              (b) handle possible extrapolation where hnew>hold


% the original his/avg etc file
grd_file_old = '/home/om/roms/latte/Data/roms_latte_grid_3b.nc';
infile = '/home/om/roms/latte/Data/espresso_2005mar_ini_alongcoast.nc';
infile = '/home/om/roms/latte/Data/espresso_2006feb_ini_alongcoast.nc';

% the new his/avg file:
%     created by (a) copying 'infile'
%                (b) altering h
%                (c) noting the new grd_file name
grd_file_new = '/home/om/roms/latte/Data/roms_latte_grid_3c.nc';
outfile = ...
    [ '/home/wilkin/roms/svn/omlab/wilkin/Projects/latte/Data/' ...
      'espresso_2005mar_ini_alongcoast_grd3c.nc'];
outfile = ...
    [ '/home/wilkin/roms/svn/omlab/wilkin/Projects/latte/Data/' ...
      'espresso_2006feb_ini_alongcoast_grd3c.nc'];


if exist(outfile)~=2
  reply = input([ outfile ' does not exist. Create it? (y/n): '],'s');
  if strcmp(reply,'n')
    return
  end  
  % create the new file by copying the old
  unixstr = [ 'cp ' infile ' ' outfile];
  disp([ 'Doing ' unixstr '...'])
  unix(unixstr)
  disp([ '... done'])
  disp([ 'Created ' outfile])
end

% note the new grid file name in the new history file
% nc.grd_file(:) = grd_file_new;
nc_attput(outfile,'global','grd_file',grd_file_new)



% copy new bathymetry to the new history file, overwriting the 
% old bathymetry inherited from the original file
%nc{'h'}(:) = grd.h;
%close(nc)
%disp([ 'Transferred new bathymetry to ' outfile])


% open the files for vertical interpolation
ncin = netcdf(infile);


% get grids
scoord = [5 0.4 50 30];
grd_old = roms_get_grid(grd_file_old,scoord,0,1);
grd_new = roms_get_grid(grd_file_new,scoord,0,1);

% Load inputs values into 'datain' and 'zin'
% interpolate to the new depths 'z' and store in 'data'
% [Don't need to interpolate the 2-D variables zeta, ubar, vbar]

for varlist = { 'temp','salt','u','v' }
  % for varlist = { 'temp','salt','u','v','NO3','NH4','phytoplankton','zooplankton','Ldetritus','Sdetritus','chlorophyll'}
  var = char(varlist);
  disp(['Doing ' var])
  tmp = ncin{var}(:);
  datain = reshape(tmp,[size(tmp,1) length(tmp(1,:))]);
  switch var
    case 'u'
      tmp = grd_old.z_u;
    case 'v'
      tmp = grd_old.z_v;
    otherwise
      tmp = grd_old.z_r;
  end
  zin = reshape(tmp,[size(tmp,1) length(tmp(1,:))]);
  
  % copy top and bottom rows of data so that if there are new z values out
  % of range of the inputs, the interpolation does not return NaN. What it
  % will give is an upward and downward copy of the shallowest and deepest
  % available data (it's rough, but we're only using this for an initial
  % condition)
  zin = [-10000*ones([1 length(zin)]); zin; ones([1 length(zin)])]; 
  datain = datain([1 1:end end],:);
  switch var
    case 'u'
      tmp = grd_old.z_u;
      water = find(grd_old.mask_u~=0);
    case 'v'
      tmp = grd_old.z_v;
      water = find(grd_old.mask_v~=0);
    otherwise
      tmp = grd_old.z_r;
      water = find(grd_old.mask_rho~=0);
  end
  z = reshape(tmp,[size(tmp,1) length(tmp(1,:))]);
  data = zeros(size(z));
  for i=water'
    data(:,i) = interp1(zin(:,i),datain(:,i),z(:,i));
    if rem(i,500)==0
      disp([ '  doing ' int2str(i) ' of ' int2str(length(data))])
    end  
  end
  
  % need to deal with NaNs
  
  % for near surface values, extrapolate upwards
  % find(isnan(data(30,:)==1))
  
  % [j,i]=ind2sub(size(grd.h),find(isnan(data(29,:))==1));

  % open the files for writing
  nc = netcdf(outfile,'write');
  nc{var}(:) = data;
  close(nc)

end


nc = netcdf(outfile,'write');
nc.details(:) = [nc.details(:) ' | Vertical interpolation ' ...
      'applied by ' which(mfilename) ' to convert to new ' ...
      grd_file_new ' bathymetry'];
close(nc)




