%------------------
% Purpose:
%------------------
% plot global M field
% plot topography
% plot ssh with following preprocessing:
% a) remove the running 15 week average
% b) optionally:	spatial filter in space
% make a movie of propagating ssh signals!

%------------------
% History:
% built from get_m_flex_filter, 	Nov. 16th/2007

%------------------
% flag_overlay 	= 1	overlay ssh with mdt
% flag_overlay  = 2	overlay ssh with bottom topography

flag_overlay = 2;

zlevels = [1:0.05:2.0];

y0 =25;
dy = 20;
y1 = y0+dy;

x0 = 315;
dx = 20/cosd((y0+y1)/2);
x1 = x0+dx;

load M_19940105_19991229_upd_highpass1_15weekanom
ireg = find(NbLongitudes>x0 & NbLongitudes<x1);
jreg = find(NbLatitudes>y0 & NbLatitudes<y1);
M1 = M(jreg,ireg);
clear M
load M_20000105_20051228_upd_highpass1_15weekanom
M2 = M(jreg,ireg);
g = find(~isnan(M1+M2));
tmp = corrcoef(M1(g),M2(g));
disp(['presistence of M in this region ' num2str(tmp(1,2))])

pause(5)
load ../../Bottom_Topography/topography.mat

fig = figure(1);
clf
hp = pcolor(NbLongitudes(ireg), NbLatitudes(jreg),M(jreg,ireg));
set(hp,'EdgeColor','none')
axis('square')

hold off 
[mdt lon_mdt lat_mdt] = load_mean_dynamic_topography(3);

imdt = find(lon_mdt>x0 & lon_mdt<x1);
jmdt = find(lat_mdt>y0 & lat_mdt<y1);
%contour(lon_mdt,lat_mdt,mdt,[-2:0.05:2],'-k')

%hp = pcolor(lon_mdt(imdt),lat_mdt(jmdt),mdt(jmdt,imdt));
%set(hp,'EdgeColor','none')



min_wks = 52;
 
flag_time_filter = 1;

%------------------
% spatial smoothing too?
%------------------
% flag_space_filter = 1	% use high pass filter


 flag_space_filter = 1;

%------------------
% currents:
%------------------
% Allow for switching between:

% flag_source = 1;	% AVISO ssha upd
% flag_source = 2;	% AVISO ssh upd

 flag_source = 2;
%------------------

% flag_date = 0;		% 2000-2005 inclusive
% flag_date = 1;		% 1994-1999 inclusive
% flag_date = 2;		% 2004-2006 inclusive

 flag_date = 0;

%------------------
% current data

switch flag_date
   case 0
%      filelist = 'ncfile_names_19991117_20060215.txt';
%      date_range = '20000105_20051228';   
      filelist = 'ncfile_names_2000_2005.txt';
      date_range = '20000105_20051228';   
   case 1
      filelist = 'ncfile_names_19931117_20000216.txt';
      date_range = '19940105_19991229';

   case 2
      filelist = 'ncfile_names_2004_2006.txt';
      date_range = '20040107_20061227';

   case 3
      filelist = 'ncfile_names_19990922_20060412.txt';
      date_range = '19990922_20060412';

end

switch flag_source

   case 1
      disp('Using upd 3 to 4-satellite AVISO ssha')
      filenamelen = 36;
%filelist = 'ncfile_names_2000_2005.txt';

      dir = '/disk/cg2/rscott/dt_upd_msla_ssha';
      file_hbar = strcat(dir,'/h_',date_range,'_ave_aviso_upd_highpass',int2str(flag_space_filter),'.mat');
      outfile_fin = strcat('movie_ssha_',date_range,'_highpass',int2str(flag_space_filter), ...
           '_15weekfilter',int2str(flag_time_filter),'_upd.mat')

   case 2
      disp('Using upd 3 to 4-satellite AVISO ssh')
      filenamelen = 36;
      dir = '/disk/staff/rscott/Aviso/dt_upd_madt_ssh';
      outfile_fin = strcat('movie_ssh_',date_range,'_highpass',int2str(flag_space_filter), ...
           '_15weekfilter',int2str(flag_time_filter),'_upd.mat')
   otherwise
      error('set flag_source')
end



pause(5)
%-------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% don't change below this line:

% initial arrays

usq   = zeros(915,1080);
vsq   = zeros(915,1080);
uv    = zeros(915,1080);
ngood = zeros(915,1080);
unity = ones(915,1080);

ssh_ts   = zeros(915,1080,15)*NaN;

% geophysical constants:
 omega = 2*pi/(3600*23+56*60);
 g=9.8;
 a = 6371e3*2*pi/360;

% all file names:
  files = textread(strcat(dir,'/',filelist),'%s');
  nfiles = size(files,1);
%? my_movie = NaN*ones(
% loop over first 14 files (fill weeks 2 through 15)

  disp('the first 14 files')

  for n = 1:14
     infile = char(files(n))
     outfile = strcat(infile(1:filenamelen),'.mat');
% convert to mat file
%     getnetcdffile(infile,outfile);
     if exist(strcat(dir,'/',outfile),'file')
        load(strcat(dir,'/',outfile))
     else
        Grid_0001=getnc(strcat(dir,'/',infile),'Grid_0001'); 
     end
% missing value replacement
     ssh = Grid_0001;
     ssh(Grid_0001>1e19) = NaN;
% convert to SI, orient so first index corresponds to latitude:
     ssh_ts(:,:,n+1) = ssh'/100;

  end



 
  disp('reading week 8 of 2000')
  pause(5)

% loop over remaining files (Feb 23, 2000, is first)
pp = 0;
  for n = 8:nfiles-7
     pp = pp+1;
% prepare the next set of 15 weeks
     
     ssh_ts(:,:,1:14) = ssh_ts(:,:,2:15);

     infile = char(files(n+7))
     outfile = strcat(infile(1:filenamelen),'.mat');
     if exist(strcat(dir,'/',outfile),'file')
        load(strcat(dir,'/',outfile))
     else
        Grid_0001=getnc(strcat(dir,'/',infile),'Grid_0001'); 
     end
% missing value replacement
     ssh = Grid_0001;
     ssh(Grid_0001>1e19) = NaN;
% convert to SI, orient so first index corresponds to latitude:
     ssh_ts(:,:,15) = ssh'/100 ;
%
     clear ssh Grid_0001
% extract the center week
% remove mean over surrounding 15 weeks
     if flag_time_filter == 1
        ssh = squeeze(ssh_ts(:,:,7)) - mean(ssh_ts,3);
     else
        ssh = squeeze(ssh_ts(:,:,7));
     end
% high-pass filter in space
     if flag_space_filter  == 1
        ssh_sm = gaussian_smooth(ssh,19,6);
	ssh = ssh - ssh_sm;
     end
     fig = figure(1);
     clf

     contourf(NbLongitudes(ireg), NbLatitudes(jreg), ssh(jreg,ireg), zlevels) %,'-k')
%
%     contourf(NbLongitudes(ireg), NbLatitudes(jreg), ssh(jreg,ireg), [-0.5:0.02:0.5]) %,'-k')

%     hp = pcolor(NbLongitudes(ireg), NbLatitudes(jreg),M(jreg,ireg));
%     hp = pcolor(NbLongitudes(ireg), NbLatitudes(jreg),ssh(jreg,ireg));
%    hp = pcolor(lon_mdt(imdt),lat_mdt(jmdt),mdt(jmdt,imdt));
%     set(hp,'EdgeColor','none')
%     caxis([-0.25 0.25])
%     caxis([0 2.5])



     colorbar
     hold on

     switch flag_overlay
        case 1
           contour(lon_mdt(imdt),lat_mdt(jmdt),mdt(jmdt,imdt),[-2:0.05:2],'-k');
        case 2
           contour(NbLongitudes(ireg), NbLatitudes(jreg), ...
	      TopoInterp(jreg,ireg),[-8000:250:0],'-w','LineWidth',2)
	   
     end
     caxis([min(zlevels) max(zlevels)])
     axis('square')
     hold on
     basemap(y0,dy,x0,dx,0,'k')
     hold on

     my_movie(pp) = getframe(gcf);

end

movie2avi(my_movie,'ssh_over_topog_v1.avi','fps',5)

save(outfile_fin,'my_movie','NbLongitudes','NbLatitudes', 'x0','x1','y0','y1',...
        'dx','dy','ireg','jreg','TopoInterp','mdt','lon_mdt','lat_mdt', ...
	'flag_date','flag_time_filter','flag_space_filter', ...
		 'flag_source','flag_overlay', 'flag_source')