%------------------
% Purpose:
%------------------
% Get global M field
% preprossing steps:	
% 1) remove 6-year mean
% 2) remove the running 15 week average
% 3) optionally:	spatial filter in space


%------------------
% History:
% built from get_m_flex_filtered.m	Feb. 27th, 2008
% built from get_m_flex.m, 	Nov. 16th/2007
% 
%------------------

 min_wks = 52;

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


 flag_filter = 0;

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

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

 flag = 1;
%------------------

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

   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

   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';

      outfile_fin = strcat('M_',date_range,'_upd_highpass',int2str(flag_filter),...
                           '_15weekanom_from_6yr_anom_ts.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';

      outfile_fin = strcat('M_',date_range,'_ref_highpass',int2str(flag_filter),...
                           '_15weekanom_from_6yr_anom.mat')

      load(strcat('../../Wind_power/uv_bar_',date_range,'_ref_smooth0_ts.mat'));

   otherwise
      error('set flag')
end

% get long term mean:

u_6yrave = u_bar;
v_6yrave = v_bar;
u_6yrave(uv_count<min_wks) = NaN;
v_6yrave(uv_count<min_wks) = NaN;
clear u_bar v_bar uv_count

load M_20000105_20051228_upd_highpass1_15weekanom_from_6yr_anom.mat

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

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


u_ts   = zeros(915,1080,15)*NaN;
v_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);

% initial arrays

ocean = find(ngood>52);

no = length(ocean);

usq   = zeros(no,nfiles);
vsq   = zeros(no,nfiles);
uv    = zeros(no,nfiles);
ngood = zeros(no,1);
unity = ones(no,1);




% 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);
     load(strcat(dir,'/',outfile))

% missing value replacement
     u = Grid_0001;
     u(Grid_0001>1e19) = NaN;
% convert to SI, orient so first index corresponds to latitude:
     u_ts(:,:,n+1) = u'/100 - u_6yrave;
% missing value replacement
     v = Grid_0002;
     v(Grid_0002>1e19) = NaN;
% convert to SI, orient so first index corresponds to latitude:
     v_ts(:,:,n+1) = v'/100 - v_6yrave;
  end


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

% loop over remaining files (Feb 23, 2000, is first)

  for n = 8:nfiles-7

% prepare the next set of 15 weeks
     
     u_ts(:,:,1:14) = u_ts(:,:,2:15);
     v_ts(:,:,1:14) = v_ts(:,:,2:15);

     infile = char(files(n+7));
     if (n/10 - fix(n/10)) == 0
        infile
     end
     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;
% convert to SI, orient so first index corresponds to latitude:
     u_ts(:,:,15) = u'/100 - u_6yrave;;
% missing value replacement
     v = Grid_0002;
     v(Grid_0002>1e19) = NaN;
% convert to SI, orient so first index corresponds to latitude:
     v_ts(:,:,15) = v'/100 - v_6yrave;

     u_bar = mean_notisnan(u_ts,3);
     v_bar = mean_notisnan(v_ts,3);

% extract the center week
% remove mean over surrounding 15 weeks
     u = squeeze(u_ts(:,:,7)) - u_bar;
     v = squeeze(v_ts(:,:,7)) - v_bar;
% high-pass filter in space
     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 

% just over the good ocean pts from 
% M_20000105_20051228_upd_highpass1_15weekanom_from_6yr_anom.mat

     u = u(ocean);
     v = v(ocean);
     usq(:,n) = u.^2;
     vsq(:,n) = v.^2;
     uv(:,n)  = u.*v;

end

g = find(ngood>=min_wks);

usq_bar  = mean_notisnan(usq,2);
vsq_bar  = mean_notisnan(vsq,2);
uv_bar   = mean_notisnan(uv,2);


M_bar = (usq_bar - vsq_bar) ./ (usq_bar + vsq_bar);


save(outfile_fin,'usq','vsq','uv','M_bar','ocean',...
                 'usq_bar','vsq_bar','uv_bar','NbLongitudes','NbLatitudes')