%------------------
% smoothing?
%------------------
% flag_filter = 1	% use high pass filter


 flag_filter = 0;

 min_wks = 52;
%------------------
% currents:
%------------------
% Allow for switching between:

% flag = 1;	% AVISO absolute currents upd
% flag = 2;	% AVISO absolute currents ref

 flag = 2;
%------------------

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

 flag_date = 1;

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

switch flag_date
   case 0
      filelist = 'ncfile_names_2000_2005.txt';
      date_range = '20000105_20051228';   
   case 1
      filelist = 'ncfile_names_1994_1999.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

   case 1
      disp('Using updated AVISO anomalies u,v')
      filenamelen = 37;
%      dir = '/disk/staff/rscott/Aviso/dt_upd_madt_UV';
      dir = '/disk/cg2/rscott/dt_upd_madt_UV';
%      file_uvbar = strcat(dir,'/uv_',date_range,'_ave_aviso_upd.mat');
      outfile_fin = strcat('M_',date_range,'_upd_highpass',int2str(flag_filter),'.mat')
      load(strcat('../../Wind_power/uv_bar_',date_range,'_upd_smooth0.mat'));
   case 2
      disp('Using ref 2-satellite AVISO absolute u,v')
      filenamelen = 37;
      dir = '/disk/staff/rscott/Aviso/dt_ref_madt_UV';
%     file_uvbar = strcat(dir,'/uv_',date_range,'_ave_aviso_ref_highpass',int2str(flag_filter),'.mat');
      file_uvbar = strcat(dir,'/uv_',date_range,'_ave_aviso_ref.mat');
%/disk/staff/rscott/Aviso/dt_ref_madt_UV/uv_19940105_19991229_ave_aviso_ref.mat
      outfile_fin = strcat('M_',date_range,'_ref_highpass',int2str(flag_filter),'.mat')
      load(file_uvbar)
   otherwise
      error('set flag')
end


pause(15)
%-------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);

% 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);

% loop over files
  for n = 1:nfiles
     infile = char(files(n))
     outfile = strcat(infile(1:filenamelen),'.mat');
% convert to mat file
%     getnetcdffile(infile,outfile);
     load(strcat(dir,'/',outfile))

% missing value replacement
     u = Grid_0001;
     u(Grid_0001>1e19) = NaN;
% missing value replacement
     v = Grid_0002;
     v(Grid_0002>1e19) = NaN;
% convert to SI, orient so first index corresponds to latitude:
% remove mean
     u = u'/100 - u_bar;
     v = v'/100 - v_bar;
% high-pass filter
     if flag_filter  == 1
        u_sm = gaussian_smooth(u,19,6);
        v_sm = gaussian_smooth(v,19,6);
	
	u = u - u_sm;
	v = v - v_sm;
     end 
     good = find(~isnan(u+v));
     usq(good) = usq(good) + u(good).^2;
     vsq(good) = vsq(good) + v(good).^2;
     uv(good)  = uv(good)  + u(good).*v(good);
     ngood(good) = ngood(good) + unity(good);

% compute std errors


end


g = find(ngood>=min_wks);

usq(g)  = usq(g) ./ngood(g);
vsq(g)  = vsq(g) ./ngood(g);
uv(g)   = uv(g)  ./ngood(g);

b = find(ngood<min_wks);
usq(b) = NaN;
vsq(b) = NaN;
uv(b)  = NaN;

M = (usq - vsq) ./ (usq + vsq);

source_info = 'generated with get_m_flex.m, includes possiblity of gaussian spatial filter set with flag_filter';
save(outfile_fin,'usq','vsq','uv','M','ngood','NbLongitudes','NbLatitudes','flag_filter','source_info','min_wks','flag_date')