% script to interpolate a roms dataset from one S coordinate to another
% so that a new run with different S coordinate parameters
% [theta_s theta_b Tcline N] can be initialized
%
% 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

% NOTE: before this m-file is executed, take these steps at the system prompt:
% 1. create a CDL header of the old file
%   ncdump -h old_his.nc > new_ini.cdl
% 2. edit the new_ini.cdl and change s_rho and s_w to set the new number of layers
% 3. generate a new init.nc file
%   ncgen -o new_ini.nc new_ini.cdl
%   ncks -A old_his.nc new_ini.nc
%
%  specify new scoord parameters:
% specify new parameters here.  N must match that set in file_out.
%%
%theta_s=4;
%theta_b=0.9;

theta_s=5;
theta_b=0.1;
Tcline=4;
N=30;
tstep=21;
[sc_r,Cs_r,sc_w,Cs_w]=scoord0(theta_s,theta_b,Tcline,N)
%%

%file_in = 'sed017_his_0037_16.nc';
file_in='http://stellwagen.er.usgs.gov/cgi-bin/nph-dods/models/adria/adria03_sed017_his_0037.nc';
file_out = 'new_ini.nc';

%copyfile(file_in,file_out);

% get the old grid info
grd_in = roms_get_grid(file_in,file_in,tstep,1);  

scoord_in=[grd_in.theta_s grd_in.theta_b grd_in.Tcline grd_in.N]



scoord_out=[theta_s theta_b Tcline N];
%grd_out=roms_get_grid(file_in,scoord_out,0,1);

% open the files for vertical interpolation
nc_in = netcdf(file_in);
ocean_time=nc_in{'ocean_time'}(:);
jd=ocean_time/(24*3600)+2440000;
nc_out = netcdf(file_out,'write');
% it would be nice if dimension resizing worked, but it does not seem to
%s_rho=nc_out('s_rho');
%s_rho(:)=N;    % resizing dimension
%s_w=nc_out('s_w');
%s_w(:)=N+1;    % resizing dimension
nc_out{'theta_s'}(:)=theta_s;
nc_out{'theta_b'}(:)=theta_b;
nc_out{'Tcline'}(:)=Tcline;
nc_out{'s_rho'}(:)=sc_r;
nc_out{'s_w'}(:)=sc_w;
nc_out{'Cs_r'}(:)=Cs_r;
nc_out{'Cs_w'}(:)=Cs_w;
nc_out{'ocean_time'}(1)=ocean_time(tstep);
close(nc_out);

grd_out=roms_get_grid(file_out,file_out,1,1);

nc_out = netcdf(file_out,'write');


% 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 = nc_in{var}(tstep,:,:,:);
  data_in = reshape(tmp,[size(tmp,1) length(tmp(1,:))]);
  switch var
    case 'u'
      tmp = grd_in.z_u;
    case 'v'
      tmp = grd_in.z_v;
    otherwise
      tmp = grd_in.z_r;
  end
  z_in = 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)
  z_in = [-10000*ones([1 length(z_in)]); z_in; ones([1 length(z_in)])]; 
  data_in = data_in([1 1:end end],:);
  switch var
    case 'u'
      tmp = grd_out.z_u;
      water = find(grd_out.mask_u~=0);
    case 'v'
      tmp = grd_out.z_v;
      water = find(grd_out.mask_v~=0);
    otherwise
      tmp = grd_out.z_r;
      water = find(grd_out.mask_rho~=0);
  end
  z_out = reshape(tmp,[size(tmp,1) length(tmp(1,:))]);
  data_out = zeros(size(z_out));
  for i=water'
    data_out(:,i) = interp1(z_in(:,i),data_in(:,i),z_out(:,i));
    if rem(i,500)==0
      disp([ '  doing ' int2str(i) ' of ' int2str(length(data_out))])
    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_out{var}(1,:,:,:) = data_out;
end
nc_out{'zeta'}(1,:,:)=nc_in{'zeta'}(tstep,:,:);
nc_out{'ubar'}(1,:,:)=nc_in{'ubar'}(tstep,:,:);
nc_out{'vbar'}(1,:,:)=nc_in{'vbar'}(tstep,:,:);

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

close(nc_in);
close(nc_out);