% loadasciidrift.m
%
% load a kriged ASCII file (as provided by Mayra
% or Jessica) into Matlab.
% First, manually strip the header line, e.g.
%	sed -n '4,$p' filename > destinationname
% MAKE SURE it has the format:
%  ID TIME LAT LON TEMP(dgC) VE(cm/s) VN(cm/s) ...
%       SPD eLAT eLON eTEMP
%     - or -
%  ID MM DD YY LAT LON ...
% 
% Include drogue on/off status from the file
% specified by dirname.  If the header file has dates
% in the format MONTH- DAY- YEAR, use "nedit" to
% convert "- " to "  ".  (don't just remove "-", as
% there are negative numbers such as lats).
%   -> Don't forget to remove "Death codes" at EOF
%
% Death types:
%	1: ran aground
%	2: picked up by vessel
%	3. quit transmitting
%	4. unreliable/sporadic transmissions
%	5. bad battery voltage
%	6. inactive status
%	0. still active

clear all;
filename='24442.dat';
dirname='24442_dirfl.dat';

fprintf('Loading file ... ');
eval(['x=dlmread(''',filename,''');']);
fprintf('done.');


% replace bad placeholder 999.99 with NaN
q=find(x==999.999);
x(q)=NaN;

if (0)  % header: ID T LAT ...
    id=x(:,1);  t=x(:,2);  
    lat=x(:,3);  lon=x(:,4); sst=x(:,5); 
    u=x(:,6)/100; v=x(:,7)/100;
    elat=sqrt(x(:,9)); elon=sqrt(x(:,10)); 
    esst=sqrt(x(:,11));
elseif (1)    % header: ID MM DD YY LAT ...
    id=x(:,1);  
    t=juliand([x(:,4) x(:,2) x(:,3)]);
    lat=x(:,5);  lon=x(:,6);
    sst=x(:,7); u=x(:,8)/100; v=x(:,9)/100;
    elat=sqrt(x(:,11)); elon=sqrt(x(:,12)); 
    esst=sqrt(x(:,13));
end;
clear x;
% make sure lon is defined from -180 to 180.
q=find(lon>180);lon(q)=lon(q)-360;

plot(1,1,'.w');hold on;overlay;
plot(lon,lat,'.k','markersize',4);
hold off;
%keyboard;

% remove values with undefinied position
q=find(isnan(lat+lon));
id(q)=[];t(q)=[];lat(q)=[];lon(q)=[];
sst(q)=[];u(q)=[];v(q)=[];
elat(q)=[];elon(q)=[];esst(q)=[];
clear q;

% now load the direc file 
%   drglst_time=0 : buoy still has drogue.
%	    set drglst_time to 9e5 (Feb.11, 4443)
%   NOTE: drglst_time=1 means drogue status 
%	    uncertain from deployment (treat as undrogued)
eval(['x=textread(''',dirname,''');']);
clear dirname
if (0)  % header: time is JULIAN DRIFTER DAYS
        %   ID WMO EXP TYP Depl.TIME Depl.LAT Depl.LON 
        %       E.TIME DRO.OFF DEATH
    dir_id=x(:,1);
    dir_deploytime=x(:,5);
    dir_deploylat=x(:,6);
    dir_deploylon=x(:,7);
    dir_endtime=x(:,8);
    dir_drglst_time=x(:,9);
    dir_death=x(:,10);
    
    q=find(dir_drglst_time==0); % drogue never lost
    dir_drglst_time(q)=9e5;
    
else    % header: time is DATE (DD MM YY)
        %  ID WMO EXP DEP.DATE DEP.LAT DEP.LON 
        %       END.DATE DROGOFFDATE DEATH
    dir_id=x(:,1);
    dir_deploytime=juliand([x(:,6) x(:,4) x(:,5)]);
    dir_deploylat=x(:,7);
    dir_deploylon=x(:,8);
    dir_endtime=juliand([x(:,11) x(:,9) x(:,10)]);
    dir_drglst_time=juliand([x(:,14) x(:,12) x(:,13)]);
    dir_death=x(:,15);
    
    q=find(dir_drglst_time==juliand([0 0 0]));
    dir_drglst_time(q)=9e5;
end
clear x
q=find(dir_deploylon>180);
dir_deploylon(q)=dir_deploylon(q)-360;

% at this point, an individual drifter's trajectory
% including deployment, endtime and drogue lost time can
% be plotted with the script plotdirinfo.m

% flag data as unrecorded (-1), no drogue (0)
% or drogued (1)
disp('compiling with dir file; marking drogue status')
drogue=-1*ones(size(id));
for i=1:length(dir_id);
  fprintf(' %d of %d\n',i,length(dir_id));
  q=find(id==dir_id(i));
  drogue(q)=1;
  q2=find(t(q)>=dir_drglst_time(i));
  drogue(q(q2))=0;
end;
q=find(drogue<0);
if (length(q)) 
  disp('WARNING: some drifters not recorded in dir file!');
  disp('  Drogue status uncertain.')
  keyboard;
end;
clear dir_* q q2 i

if (strcmp(filename,'tropatl.dat'));
    q=find(lon<-60 & lat<0); truncatedrift;
    q=find(abs(lat)>25); truncatedrift;
    i=[-74.2 -76.1 -79.3 -81.4 -84.8 ...
            -85.8 -95.8 -135.9 -133.3];
    j=[-1.4 8.1 9.2 8.1 10.6 14.3 17.6 17.1 -2.1000];
    q=find(isinpoly(lon,lat,i,j));
    truncatedrift;
    clear i j;
end;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%disp('DROPPING ALL DRIFTERS AFTER 8 APRIL 2003');
%q=find(t>juliand([2003 4 8 0 0 0]));
%truncatedrift;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calculate wind speed and wind stress corresponding
% to each observation, using NCEP daily winds
year=gregoriand(t);
year=year(:,1);
uy=nanunique(year);
Wx=NaN*ones(size(lon));Wy=Wx;Tx=Wx;Ty=Wx;
for yr=min(uy):max(uy);
    q=find(year==yr);
    [tw,lonw,latw,wx,wy,taux,tauy]= ...
        loaddailywinds(yr);
    
     tt=t(q)+2443874;
     tt(find(tt<min(tw)))=min(tw);
     tt(find(tt>max(tw)))=max(tw);
     xx=lon(q);
     xx(find(xx<min(lonw)))=min(lonw);
     xx(find(xx>max(lonw)))=max(lonw);
     yy=lat(q);
     yy(find(yy<min(latw)))=min(latw);
     yy(find(yy>max(latw)))=max(latw);
     
     Wx(q)=interp3(lonw,latw,tw,wx, ...
         xx,yy,tt,'linear');
     Wy(q)=interp3(lonw,latw,tw,wy, ...
         xx,yy,tt,'linear');
     Tx(q)=interp3(lonw,latw,tw,taux, ...
         xx,yy,tt,'linear');
     Ty(q)=interp3(lonw,latw,tw,tauy, ...
         xx,yy,tt,'linear');
end;
clear year uy yr tt xx yy q;
clear tw lonw latw wx wy taux tauy;

% next, remove the slip. 
% Drogued: 7e-4 m/s per m/s of wind. (Niiler & Paduan, 1995)
% Undrogued: 7e-4+7.9e-3=8.6e-3 (", Pazan & Niiler 2000) 
q=find(drogue==1);
u(q)=u(q)-7e-4*Wx(q);
v(q)=v(q)-7e-4*Wy(q);
q=find(drogue==0);
u(q)=u(q)-8.6e-3*Wx(q);
v(q)=v(q)-8.6e-3*Wy(q);
clear q Wx Wy;

% run "lowpassdrifters.m".  This script
% runs a 5 day lowpass filter on the wind stresses
% and speeds, the latter to calculate a detided
% and de-inertialized (if that is that a word)
% u, v; it also NaNs the velocities and
% where elon or elat are greater than
% 0.03 degrees (approx. 10 cm/s max error for speed)
% or where esst>0.1 degC.
disp('Lowpassing the velocities.');
lowpassdrifters

% calculate the Ekman component of the
% velocity, using the Ralph & Niiler (1999)
% model with the updated Niiler WOCE book coefficient.
lateq=2;    % don't assume Ekman dynamics this close
            % to the equator.
A=0.081;    % empirical coefficient, units: sec^(-1/2)
            %  (updated value of beta=.065)
O=7.2921e-05;
tau=sqrt(Tx.^2+Ty.^2);  % magnitude of wind stress
f=2*O*sin(2*pi*lat/360);
q=find(abs(lat)<lateq);
f(q)=NaN;
load sfcdens;
rho=interp2(lon_coads,lat_coads',rho_coads,lon,lat);
ustar=sqrt(tau./rho);
uek=A*ustar./sqrt(abs(f));  % mag.of ekman velocity
Hstar=ustar./sqrt(abs(f))/A;  % ekman depth
zrat=15./Hstar;             % fractional ek.depth
% calculate angle of ek.velocity relative to wind
theta=sign(f)*45*(2*pi/360) .*(1+zrat/1.2);
% convert to angle relative to N/S by adding
% angle of wind with respect to due eastward
theta=atan2(Ty,Tx)-theta;
vek=uek.*sin(theta);
uek=uek.*cos(theta);
% set Ekman vel to 0 where it is NaN (very small wind
% stress, causing Ekman depth->0) over the ocean
q=find(isnan(uek) & isfinite(rho));
uek(q)=0;vek(q)=0;
q=find(abs(lat)<lateq);
uek(q)=NaN;vek(q)=NaN;

clear lateq A bt O tau f q lon_coads lat_coads rho_coads
clear ustar Hstar zrat theta ans rho Ty Tx

% Define u_g, v_g to be the lowpassed
% total speeds minus the Ekman speed
u_g=u-uek;
v_g=v-vek;
% clear uek vek - can be reconstructed from
% u_g and u.
clear uek vek

% calculate Lagrangian derivatives (D/Dt sst, u, v)
calcdeltas

% before calculating altimetry-derived geostrophic
% velocity anomalies, save everything and clear
% all variables but the necessary ones for the
% alt.calculations: lon, lat and t.
%disp(' ');
%disp('Saving data as "drift.mat"');
%save drift;
%clear all;
%load drift.mat 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
%load drift.mat;

% calculate Reynolds SST interpolated to drifter positions
calcreynoldssst

% save everything
clear filename
disp(' ');
disp('Saving data as "drift.mat"');
save drift
