% Optimally interpolate a set of CTD station data on standard levels
%
% Assumes:
%    station data have been loaded and interpolated to standard levels
%       and are present as the variables:
%    zstd, Lon, Lat, Time (jd convention), 
%    Temp, Salt dimensioned [depth,station] 
% and
%    X,Y is the target lon,lat grid
%    clim_day is the target "climatological day" given as a Gregorian 
%      date vector with the year being irrelevant, 
%      e.g. clim_day = [NaN 7 1 0 0 0] for July-1.

load MJJA_NODC_stdlev

% preallocate the output mapped Data field (and expected Error)
T = NaN*ones([length(zstd) size(X)]);
% S = NaN*ones([length(zstd) size(X)]);
% E = T;

% prepare data for the OI

% convert times to numbers of days from climatological July 1 (this has to
% be consistent with the selection criteria applied in
% nmfs_clim_getdata/nmfs_subset when extacting the subset of data to be
% mapped

% Create a time variable Dday which is number of days from the nominal 
% climatological analysis date set by clim_day
tmp = gregorian(Time(:)); 
year_of_obs = tmp(:,1);
num = length(year_of_obs);
Dday = Time(:)-julian(year_of_obs,clim_day(2)*ones([num 1]),...
    clim_day(3)*ones([num 1])); 

% step through standard levels
% for stdlev = fliplr(1:length(zstd)) % flip to do small datasets first
for stdlev = 3 % just do one level for the purpose of speed testing
  
  disp([ 'Standard level z = ' num2str(zstd(stdlev))])

  % temperature ----------------------------------------------------
  
  disp([ ' Temperature'])
  % make a temporary position-data matrix
  tmp = [Lat(:) Lon(:) Dday(:) Temp(stdlev,:)'];

  % only retain data within +/- this many days
  too_old = 45.0;
  clip = find(abs(Dday) >= too_old);
  tmp(clip,:) = [];
            
  % then remove entries with missing data
  bad = find(any(isnan(tmp'))==1);
  tmp(bad,:) = [];
  
  % separate the position and data
  posdata = tmp(:,1:3);
  data = tmp(:,4);
  
  % OI
  posgrid = [Y(:) X(:) zeros(size(X(:)))]; % time=zero at analysis date 

  % DO TWO OIS WITH DECREASING LENGTH/TIME SCALES
  %
  % first OI at larger length/time covariance scales for 'mean' or 
  % background field
  lonscl = 6;
  latscl = 6;
  timscl = 60;
  e = 0.1; 
  s = 1;
  detailstr1 = [ ...
      ' | OI parameters: lonscl = ' num2str(lonscl) ...
      '; latscl = ' num2str(latscl) '; timscl = ' num2str(timscl) ...
      '; e = ' num2str(e) '; s = ' num2str(s)];
  disp(detailstr1)
  ax_range = [min(X(:)) max(X(:)) min(Y(:)) max(Y(:))];
  subregions = gensubregions(ax_range,0,[1 1]);
  vout = oi_3d_opt(data,posgrid,posdata,subregions,...
      lonscl,latscl,timscl,e,s);
  T(stdlev,:,:) = reshape(vout.vgrid,size(X));
  % E(stdlev,:,:) = reshape(vout.vgriderror,size(X));
  
end

if ~exist('details')
   details = [];
end
oidetails = detailstr1;

% save

lon = X;
lat = Y;
save testfile lon lat T oidetails clim_day zstd