% Example script for defining a cruisetrack through a ROMS domain for twin
% experiments on assimilation of subsurface profile data from, e.g.
% WOCE/TOGA Volunteer Observing Ship (VOS) program XBT/XCTD profiles, 
% underway ADCP, or AUVs. 

% East Australia Current application

% Coordinates for way points of Tasman Sea lines in high resolution 
% repeat VOS XBT/XCTD sections

% Sydney
xp(1) = 151+15/60;
yp(1) = -33-52/60;
% Wellington
xp(2) = 172+45/60;
yp(2) = -39-45/60;
% North Cape
xp(3) = 172+15/60;
yp(3) = -34-15/60;
% Noumea
xp(4) = 167+45/60;
yp(4) = -23-15/60;
% Brisbane
xp(5) = 153;
yp(5) = -27;
% WOCE PR11 is along 30S
xp(6) = 153;
yp(6) = -30;
xp(7) = 163;
yp(7) = -30;
% WOCE PR13 is along 43S
xp(8) = 148;
yp(8) = -43;
xp(9) = 163;
yp(9) = -43;


% load the roms grid structure (required by roms_cruisetrack to define the
% perimeter of the psi points grid)
if exist('g')~=1
  g = grd('eac4');
end

% plot the routes over the ROMS bathymetry
clf
hold on
pcolorjw(g.lon_rho,g.lat_rho,g.mask_rho_nan.*g.h)
amerc
axis([145.6029  179.6727  -47.0782  -20.5115])

plot(xp([1 2]),yp([1 2]),'-') % PX34
% plot(xp([1 3]),yp([1 3]),'-') % Sydney to North Cape not a standard route
plot(xp([1 4]),yp([1 4]),'-') %
plot(xp([5 4]),yp([5 4]),'-') % PX30/31
plot(xp([6 7]),yp([6 7]),'-') % WOCE PR11
plot(xp([8 9]),yp([8 9]),'-') % WOCE PR13

% draw the coastline
% gebco_eez(0,'k')

% label the way points
clear han
k=1;
han(k) = text(xp(1),yp(1),'Sydney ','horizontalalignment','right');k=k+1;
%han(k) = text(xp(2),yp(2),' Wellington');k=k+1;
%han(k) = text(xp(3),yp(3),' North Cape');k=k+1;
%han(k) = text(xp(4),yp(4),' Noumea');k=k+1;
han(k) = text(xp(5),yp(5),'Brisbane ','horizontalalignment','right');k=k+1;
han(k) = text(xp(6),yp(6),'PR11 ','horizontalalignment','right');k=k+1;
han(k) = text(xp(8),yp(8),'PR13 ','horizontalalignment','right');k=k+1;
set(han,'color','r','fontsize',10)


% ------------------------------------------------------------------------


% The cruise path is defined by vectors lonpath and latpath (which could be
% specified as waypoints, or created interactively with ginput)

% Suppose a single vessel followed the route:
% Hobart (PR13) - Wellington - Sydney - Noumea - Brisbane - then allong PR11
% lonpath = xp([8 9 2 1 4 5 6 7]);
% latpath = yp([8 9 2 1 4 5 6 7]);

opt = 'EACVOSXBT';
%opt = 'EACWOCE';
%opt = 'interactive';

switch opt
  
  case 'EACVOSXBT'
    opt
    % Suppose a single vessel followed the route of the high-resolution VOS
    % XBT/XCTD lines from
    % Wellington - Sydney - Noumea - Brisbane - then along PR11
    lonpath = xp([2 1 4 5 6 7]);
    latpath = yp([2 1 4 5 6 7]);

  case 'interactive'
    opt
    % To generate a path interactively
    [lonpath,latpath] = ginput;

  case 'EACWOCE'  
    opt
    % the WOCE Tasman Sea seasonal repeat hydrography
    lonpath = [152 163 xp(3) xp(1) 149 155 165 165 xp(8)];
    latpath = [-28 -28 yp(3) yp(3) -40.5 -38 -38 yp(8) yp(8)];
       
end

% compute the modeled cruisetrack following the way points in
% lonpath,latpath
dkm = 30; % nominal station spacing
speed = 15; % knots

% addpth /home/wilkin/matlab -end % path to my roms tools
% grd, roms_get_grid and roms_cruisetrack
[lontrk,lattrk,timtrk] = roms_cruisetrack(lonpath,latpath,g,dkm,speed);

% Output timtrk is the nominal time in days at which the stations would be
% sampled at the given ship speed (relative to timtrk(1)=0 at the start of
% the cruise)

% Plot the sampling locations
hold on
plot(lontrk,lattrk,'m*')
hold off

notes = ['Idealized EAC cruise track for twin experiments created by ' ...
  which(mfilename) ' on ' datestr(now)];

disp(['Saving track to ' opt '.mat'])
save(opt,'lontrk','lattrk','timtrk','notes')





