
% calceddyterms.m
%
% calculate the eddy components up, vp, sstp
% from u, v, sst and Ubar, Vbar, Tbar.
% also calculate the Lagrangian time scales
% Tlu, Tdu, Tlv, Tdv, Tlt, Tdt.
%
% assumes the mean is stored in a file called
% "meandrift_*"

clear all;
%filename='trop_atlantic'; meanfilename='tropatl';
filename='atlantic'; meanfilename='atl';
%filename='pacific_93-03'; meanfilename='wpac';
%filename='trop_indian';meanfilename='tropind';

fprintf('loading file %s.mat\n',filename);
eval(['load ',filename,' t sst u* v* id lon lat drogue']);
dt=median(diff(t(1:100)));

% trajectory length (days) for time scale calculations
lengthtraj=90;  
lengthtraj=round(lengthtraj)/dt;   % convert to data points

fprintf('loading file meandrift_%s.mat\n',meanfilename);
eval(['load meandrift_',meanfilename, ...
        ' Lon Lat U*barI V*barI', ...
	    ' TbarI NumObs* *I_ann* *I_san*']);
% truncate outside bounds of Lon, Lat
q=find(lon<min(Lon)-median(abs(diff(Lon))) | ...
    lon>max(Lon)+median(abs(diff(Lon))));
truncatedrift;
q=find(lat<min(Lat)-median(abs(diff(Lat))) | ...
    lat>max(Lat)+median(abs(diff(Lat))));
truncatedrift

% define primed vectors
up=NaN*ones(size(u));
vp=up;sstp=up;ubar=up;vbar=up;sstbar=up;
ugbar=up;vgbar=up;
u_gp=up; v_gp=up;
u_gp_tot=up; v_gp_tot=vp;
% initialize Lagrangian scales
Tlu=up; Tlv=up; Tdu=up; Tdv=up;
%Tlt=up; Tdt=up;

qid=find(diff(id));
qid=[0;qid;length(id)];

% loop through, one drifter at a time
for ind=2:length(qid);

  fprintf('%d of %d\n',ind,length(qid));
  a=qid(ind-1)+1:qid(ind);

  uCa=interp2(Lon,Lat,UI_ann(:,:,1,1),lon(a),lat(a),'linear');
  uSa=interp2(Lon,Lat,UI_ann(:,:,2,1),lon(a),lat(a),'linear');
  uCsa=interp2(Lon,Lat,UI_san(:,:,1,1),lon(a),lat(a),'linear');
  uSsa=interp2(Lon,Lat,UI_san(:,:,2,1),lon(a),lat(a),'linear');
  vCa=interp2(Lon,Lat,VI_ann(:,:,1,1),lon(a),lat(a),'linear');
  vSa=interp2(Lon,Lat,VI_ann(:,:,2,1),lon(a),lat(a),'linear');
  vCsa=interp2(Lon,Lat,VI_san(:,:,1,1),lon(a),lat(a),'linear');
  vSsa=interp2(Lon,Lat,VI_san(:,:,2,1),lon(a),lat(a),'linear');
  tCa=interp2(Lon,Lat,TI_ann(:,:,1,1),lon(a),lat(a),'linear');
  tSa=interp2(Lon,Lat,TI_ann(:,:,2,1),lon(a),lat(a),'linear');
  tCsa=interp2(Lon,Lat,TI_san(:,:,1,1),lon(a),lat(a),'linear');
  tSsa=interp2(Lon,Lat,TI_san(:,:,2,1),lon(a),lat(a),'linear');
  ugCa=interp2(Lon,Lat,UGI_ann(:,:,1,1),lon(a),lat(a),'linear');
  ugSa=interp2(Lon,Lat,UGI_ann(:,:,2,1),lon(a),lat(a),'linear');
  ugCsa=interp2(Lon,Lat,UGI_san(:,:,1,1),lon(a),lat(a),'linear');
  ugSsa=interp2(Lon,Lat,UGI_san(:,:,2,1),lon(a),lat(a),'linear');
  vgCa=interp2(Lon,Lat,VGI_ann(:,:,1,1),lon(a),lat(a),'linear');
  vgSa=interp2(Lon,Lat,VGI_ann(:,:,2,1),lon(a),lat(a),'linear');
  vgCsa=interp2(Lon,Lat,VGI_san(:,:,1,1),lon(a),lat(a),'linear');
  vgSsa=interp2(Lon,Lat,VGI_san(:,:,2,1),lon(a),lat(a),'linear');
  
  cosTa=cos(2*pi*t(a)/365.24); 
  sinTa=sin(2*pi*t(a)/365.24);
  cosTsa=cos(4*pi*t(a)/365.24);
  sinTsa=sin(4*pi*t(a)/365.24);
  
  uS=(uCa.*cosTa+uSa.*sinTa)+(uCsa.*cosTsa+uSsa.*sinTsa);
  vS=(vCa.*cosTa+vSa.*sinTa)+(vCsa.*cosTsa+vSsa.*sinTsa);
  tS=(tCa.*cosTa+tSa.*sinTa)+(tCsa.*cosTsa+tSsa.*sinTsa);
  ugS=(ugCa.*cosTa+ugSa.*sinTa)+(ugCsa.*cosTsa+ugSsa.*sinTsa);
  vgS=(vgCa.*cosTa+vgSa.*sinTa)+(vgCsa.*cosTsa+vgSsa.*sinTsa);
  clear uSa uCa uSsa uCsa vSa vCa vSsa vCsa tSa tCa tSsa tCsa;
  clear ugSa ugCa ugSsa ugCsa vgSa vgCa vgSsa vgCsa;

  ubar(a)=interp2(Lon,Lat,UbarI,lon(a),lat(a),'linear');
  vbar(a)=interp2(Lon,Lat,VbarI,lon(a),lat(a),'linear');
  sstbar(a)=interp2(Lon,Lat,TbarI,lon(a),lat(a),'linear');
  ugbar(a)=interp2(Lon,Lat,UGbarI,lon(a),lat(a),'linear');
  vgbar(a)=interp2(Lon,Lat,VGbarI,lon(a),lat(a),'linear');
  
  up(a)=u(a)-uS-ubar(a);
  vp(a)=v(a)-vS-vbar(a);
  sstp(a)=sst(a)-tS-sstbar(a);
  u_gp(a)=u(a)-ugS-ugbar(a);
  v_gp(a)=v(a)-vgS-vgbar(a);
  u_gp_tot(a)=u(a)-ugbar(a);
  v_gp_tot(a)=v(a)-vgbar(a);
  
  if (length(isfinite(up(a)))>=lengthtraj);
    uu=up(a(2:length(a)-1));
    vv=vp(a(2:length(a)-1));
    sstt=sstp(a(2:length(a)-1));
    tt=t(a(2:length(a)-1))-t(a(2));
    tlu=NaN*ones(size(tt));tdu=tlu;
    tlv=tlu;tdv=tdu; tlt=tlu; tdt=tlu;
    % calculate time scales every 10 days
    for j=round(lengthtraj/2):round(10/dt): ...
          length(uu)-floor(lengthtraj/2)-1;
      jj=j-round(lengthtraj/2)+1: ...
        j+round(lengthtraj/2)+1;
      if (isfinite(mean(uu(jj)+vv(jj))) ...
              & max(diff(tt(jj)))==t(2)-t(1))
        [r,tau]=normcorr(uu(jj),tt(jj));
        [tdu(j),tlu(j)]=lagrangetime(r,tau);
        [r,tau]=normcorr(vv(jj),tt(jj));
        [tdv(j),tlv(j)]=lagrangetime(r,tau);
      end;
      %if (isfinite(mean(sstt(jj))) ...
      %        & max(diff(tt(jj)))==t(2)-t(1))
      %  [r,tau]=normcorr(sstt(jj),tt(jj));
      %  [tdt(j),tlt(j)]=lagrangetime(r,tau);
      %end;    % no NaNs in uu, vv
    end; % loop every 10 days
    Tdu(a(2:length(a)-1))=tdu;
    Tlu(a(2:length(a)-1))=tlu;
    Tdv(a(2:length(a)-1))=tdv;
    Tlv(a(2:length(a)-1))=tlv;
    %Tdt(a(2:length(a)-1))=tdt;
    %Tlt(a(2:length(a)-1))=tlt;

  end; % length>lengthtraj
end;  % loop over drifters

% convert lengthtraj back to days
lengthtraj=round(lengthtraj*dt);

clear qid ind a q k l uS vS tS sinT cosT ugS vgS
clear Lon Lat UbarI VbarI TbarI *_ann *_san
clear NumObs* thresh ans dt uu vv tt sstt
clear tlu tlv tdu tdv tlt tdt j jj ans r tau

fprintf('Saving eddyterms_%s.mat\n',meanfilename);
eval(['save eddyterms_',meanfilename, ...
        ' up ubar vp vbar sstp sstbar ', ...
        ' u_gp* v_gp* ugbar vgbar ', ...
        ' Tl* Td* lengthtraj lon lat t id drogue']);

% before calculating altimetry-derived geostrophic
% velocity anomalies, clear 
% all variables but the necessary ones for the
% alt.calculations: lon, lat and t.
clear u* v* sst* T* drogue id
eval(['load eddyterms_',meanfilename,' lon lat t']);
% calculate AVISO velocities interpolated to
% drifter positions; this routine also creates
% weekly files of altimetry for the region
% spanned by the drifter observations
calcavisovels
% re-load the other variables
eval(['load eddyterms_',meanfilename]);
% save the lot
fprintf('Saving eddyterms_%s.mat\n',meanfilename);
eval(['save eddyterms_',meanfilename, ...
        ' up ubar vp vbar sstp sstbar ', ...
        ' u_gp* v_gp* ugbar vgbar u_av v_av ', ...
        ' Tl* Td* lengthtraj lon lat t id drogue']);

