% 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


nena_initial = 6


% the original his/avg etc file

infile = '../../in/nena_his_day458+.nc';
infile = '../../in/nena_natl_avg_001_oldh.nc';
infile = '../../in/nena_natl_init_day368.25+bio.nc';
% the new his/avg file:
%     created by (a) copying 'infile'
%                (b) altering h
%                (c) noting the new grd_file name
out_file = '../../in/nena_his_day458+_nb.nc';
out_file = '../../in/nena_natl_avg_001.nc';
out_file = '../../in/nena_natl_init_day368.25+bio-grd6.nc';
grd_file = [];

if exist('nena_initial')
  switch nena_initial
    
    case 2
      infile   = 'nena_natl_avg_001.5.nc';
      out_file = 'nena_natl_init_day368.25.nc';
      % the new grid file
      grd_file = 'roms_nena_grid_5.nc';
      
    case 6
      grd_file = 'roms_nena_grid_6.nc';
      800 convention place
    otherwise
      error('nena_initial option invalid')
  end
else
  warning(['Check infile and outfile etc option in ' mfilename])
  disp('paused...')
  pause
end

if exist(out_file)~=2
  reply = input([ out_file ' 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 ' ' out_file];
  disp([ 'Doing ' unixstr '...'])
  unix(unixstr)
  disp([ '... done'])
  disp([ 'Created ' out_file])
end

% open the new history file for writing
nc = netcdf(out_file,'write');

% open the new grid file
grd = roms_get_grid(grd_file);  

% note the new grid file name in the new history file
nc.grd_file(:) = grd_file;

% 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 ' out_file])

% open the files for vertical interpolation
ncin = netcdf(infile);
nc = netcdf(out_file,'write');

% get grids
grdin_file = ncin.grd_file(:);
if exist(grdin_file)~=2
  warning([ 'Unable to open ' grdin_file])
  grdin_file = grdin_file(1+max(findstr(grdin_file,'/')):end);
  disp([ 'Trying ' grdin_file])
  if exist(grdin_file)~=2
    warning([ 'Unable to open ' grdin_file])
  end  
end
grdin = roms_get_grid(grdin_file,infile,1,1);
grd = roms_get_grid(nc.grd_file(:),out_file,1,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 = grdin.z_u;
    case 'v'
      tmp = grdin.z_v;
    otherwise
      tmp = grdin.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.z_u;
      water = find(grd.mask_u~=0);
    case 'v'
      tmp = grd.z_v;
      water = find(grd.mask_v~=0);
    otherwise
      tmp = grd.z_r;
      water = find(grd.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));
  
  nc{var}(:) = data;
end

nc.details(:) = [nc.details(:) ' | Vertical interpolation ' ...
      'applied by ' which(mfilename) ' to convert to new ' ...
      grd_file ' bathymetry'];
  

close(nc)




