
% path to John's general roms slicer
addpath /home/jevans/roms/4dvar/gslice -end

% load the track coordinates previously generated with e.g.
% eac_create_cruisetrack
load('EACVOSXBT.mat')

% Get the control structure (Required to manage cruise track station 
% times that span multiple ROMS output files)

direc = '/moby/wilkin/roms/eac/out/grd4-1';
rexp = 'eac_his*';
ctl = jge_timectl ( direc, rexp );

base_date = parsetnc(nc_attget(ctl.files{1},'ocean_time','units'));
base_datenum = datenum(base_date); 
ctl.datenum = base_datenum + ctl.time/86400;

% modify the time properly
timtrk = timtrk*86400;  % convert to seconds

% Adjust timtrk to start around 19-Feb-2001 (totally arbitrary)
begin = min(find(ctl.datenum > datenum([2001 2 20])))
timtrk = ctl.time(begin) + timtrk;
datenumtrk = base_datenum + timtrk/86400;



% Drive the ship through ROMS to sample each variable to be included in 
% the twin experiment assimilation data set

f = 1; % counter for plots

for varlist = {'temp','salt','u','v'}
  
  varname = char(varlist);
  [slicevar, r, z, Ipos, Jpos] = ...
    gslice ( ctl, varname, lontrk(:), lattrk(:), timtrk(:) );
  
  eval([varname ' = slicevar;'])
  eval([varname '_zr = z;'])
  
  % switch Ipos,Jpos
  % STILL NEED TO CHECK THE FRACTIONAL VALUES AND WHETHER U,V ARE
  % STAGGERRED CORRECTLY
  eval([varname '_Ipos = Jpos;'])
  eval([varname '_Jpos = Ipos;'])
    
  figure(f); clf
  
  % plot versus time
  % pcolorjw(base_datenum+slicetime/86400,z,slicevar); datetick('x',1,'keeplimits')
  
  % plot versus index
  pcolorjw(repmat(1:length(lontrk),[size(slicevar,1) 1]),z,slicevar) 
  title([varname '   ' datestr(datenumtrk(1),1) ' to ' datestr(datenumtrk(end),1)])
  colorbar
  
  f=f+1;
    
end



