% Simple script for plotting some features of the ROMS RISSAGA simulations

% simple control of default behaviour for 
if exist('pause_cmd')~=1
  pause_cmd = 'pause'; % set to 'drawnow' at command to animate sequence
end
if exist('ax')~=1
  ax = 'tight'; % can override with zoom area axis range
end


% Name output file to plot; Load grid 
f = 'his_rissaga_021.nc';
gf = nc_attget(f,'global','grd_file');
g = roms_get_grid(gf,f);
try
  g.lon_coast = nc_varget(gf,'lon_coast');
  g.lat_coast = nc_varget(gf,'lat_coast');
catch
  g.lon_coast = NaN;
  g.lat_coast = NaN;
  warning(['Trouble loading coast data from ' gf ])
end

% Select the grid point at which to plot time series
if exist('Ipos')~=1 | exist('Jpos')~=1
  figure
  pcolorjw(g.lon_rho,g.lat_rho,g.h.*g.mask_rho_nan)
  [Jpos,Ipos] = closest(g.lon_rho,g.lat_rho,ginput(1));
  hold on
  plot(g.lon_rho(Jpos,Ipos),g.lat_rho(Jpos,Ipos),'ks')
  plot(g.lon_coast,g.lat_coast)
  hold off
end

% Load data for plotting ---------------------------------------

% Load model time and free surface solution
zeta=nc_varget(f,'zeta');
ot=nc_varget(f,'ocean_time');

% Reconstruct the air pressure forcing
clear pair
for t=1:length(ot)
  p=ana_pair(ot(t),g,[7 20*pi/180 30 18e3],1);
  pair(t,:,:)=p;
end

% std deviation of sea level variability
zeta_std = squeeze(std(zeta,1,1));

figure
% time series at a single grid point
plot(ot/60,zeta(:,Jpos,Ipos),ot/60,(pair(:,Jpos,Ipos)-1020)/70)
set(gcf,'pos',[18   392   560   420]);

figure
% wysiwyg
set(gcf,'pos',[610 238 768 576])
% step through time showing model response, with air pressure
% and sea level time series
for t=1:length(ot)

  hax(1) = subplot(211);  
  data = squeeze(zeta(t,:,:))    ;  %./zeta_std;
  pcolorjw(g.lon_rho,g.lat_rho,g.mask_rho_nan.*data)
  hold on
  han = plot(g.lon_rho(Jpos,Ipos),g.lat_rho(Jpos,Ipos),'kx');
  set(han,'markersize',5,'linew',2)
  ylim = get(gca,'ylim');
  set(gca,'DataAspectRatio',[1 cos(mean(ylim)*pi/180) 1]);
  axis(ax)
  % axis tight
  plot(g.lon_coast,g.lat_coast)
  hold off
  caxis(0.05*[-1 1])
  colorbar
  
  hax(2) = subplot(212);
  T = 1:t;
  % rescale pair simply to have it plot on the zeta axis range
  plot(ot(T)/60,zeta(T,Jpos,Ipos),ot(T)/60,(pair(T,Jpos,Ipos)-1020)/70)
  set(gca,'xlim',[0 max(ot)/60],'ylim',[-0.2 0.2])
  xlabel('minutes')
  grid on
  
  set(hax(1),'pos',[0.1 0.3 0.775 0.6]);
  set(hax(2),'pos',[0.1 0.08 0.775 0.15])
  eval(pause_cmd) % could be 'drawnow' or 'pause' to control how frames advance 
  
end