global float number timer temp pres lat lon u v wwww
[float, number,timer,temp,pres,lat, lon, u,v,wwww] = rfcread;
%plot(time,temp)
nn = length(lat);
[ n1, n2] = findgood( u, -99., 99. );
DuDt = zeros(nn,1);
DvDt = zeros(nn,1); 
uin  = zeros(nn,1); 
vin  = zeros(nn,1); 
n2 = n2 - 1;
j = 1;
while j < nn,
   Dt   = timer(j+1) - timer(j);
% --express the time in hours-----
   Dt = Dt * 24.;
   DuDt(j) = (u(j+1) - u(j))/ Dt;
   DvDt(j) = (v(j+1) - v(j))/ Dt;    
   uin(j)  = (u(j+1) + u(j))/ 2;
   vin(j)  = (v(j+1) + v(j))/ 2;
   f(j) = coriolis((lat(j+1)+lat(j))/2);
   f(j) = f(j) * 2 * pi;
   a1(j) = sqrt( DuDt(j)^2 + DvDt(j)^ 2);
   a2(j) = abs( f(j) * sqrt( uin(j)^2 + vin(j)^2)) ;
   tprime(j) = (timer(j+1) + timer(j))/2.;
   LHS(j)    = DuDt(j) - f(j)*vin(j) + (DvDt(j) + f(j)*uin(j)) * i;
   j = j + 1;
end
f(nn) = f(nn-1);
f = f';
n2 = n2 - 1;
i2 = figure ;
Ro = a1 ./ a2;
subplot(311);
plot(tprime(n1:n2), u(n1:n2), tprime(n1:n2), v(n1:n2),'--');
xlabel('time, days');
ylabel('Velocity, cm/s');
axis([ timer(1)  timer(nn)  -50.  50. ]);         
titlos = str2mat(['Float no ' float]);
title(titlos);
%
subplot(312);
  semilogy(tprime(n1:n2), Ro(n1:n2));
  axis([ timer(1)  timer(nn)  0.001  1 ]);
xlabel('time, days');
ylabel('Rossby number');
%
subplot(313);
plot(tprime, real(LHS),tprime, imag(LHS),'--' );
axis([ timer(1)  timer(nn)  -6  6 ]);             
xlabel('time, days');
ylabel('LHS, cm/s per hour');
pause
subplot(111);
% --- plot Squared Rossby number versus Squared LHS of momentum ---
% --- equations ---------------------------------------------------
fig2 = figure;
ampl = real(LHS) .^ 2 + imag(LHS) .^2;
loglog(ampl(n1:n2), Ro(n1:n2),'+');
xlabel('Squared LHS');
ylabel('Rossby number');
title(titlos); 
pause
% --- compute acceleration and velocity magnitudes to examine -----
% --- their correlations ------------------------------------------
fig3 = figure;
uampl = sqrt( u .^2 + v .^2 );
accel = sqrt( ampl );
fcoriol = abs(f .* uampl);
plot(uampl(n1:n2), fcoriol(n1:n2),uampl(n1:n2), accel(n1:n2),'+');
xlabel('Velocity amplitude (cm/s)');
ylabel('Total acceleration (cm/s per hour)');
title(titlos);  
pause
% --- compute and plot vector differences between Coriolis and ----
% --- relative accelerations --------------------------------------
fig4 = figure;
DX = DuDt + f.* v;
DY = DvDt - f.* u;
R  = sqrt( DX .^2 + DY .^2);
R  = abs( R ./ ( f .* uampl)) - 1.;
theta = atan2(DY, DX);
polar(theta(n1:n2), R(n1:n2), '+');
grid on;
title([ titlos, '  LHS/fU']);   
pause
