%% Program to plot Ekman spiral
%% Input section
clear all; close all;
i=sqrt(-1) % To make sure i is the complex i
h=120; % Water depth
nut=.02; % (constant) turbulence viscosity
phi=51; % Latitude (deg.)
rho=1000; % Water density
tauwx=0; % Wind shear stress x-direction
tauwy=1; % Wind shear stress y-direction
tau=tauwx+tauwy*i; % Shear stress vector
f=2*7.27e-5*sin(phi*pi/180); % Coriolis coefficient
dek=sqrt(2*nut/f) % Ekman depth
dz=5; % Vertical step size
zek=[-h:dz:0]; % Vertical coordinate
%% Define some additional matrices for plotting purposes
xek=zeros(size(zek));
yek=zeros(size(zek));
taux=zeros(size(zek));
tauy=zeros(size(zek));
taux(end)=tauwx;
tauy(end)=tauwy;
%% Solution of Ekman spiral
s=tau*dek/rho/(1+i)/nut*exp((1+i)/dek*zek);
uek=real(s);
vek=imag(s);
%% Plot results
figure(1);
quiver3(xek,yek,zek,uek,vek,0,'linewidth',2);
hold on
quiver3(xek,yek,zek,taux,tauy,0,'r','linewidth',2);
title('Ekman spiral')
xlabel('u');ylabel('v');
zlabel('z')
print('-dpng','ekman.png')
print('-depsc','ekman.eps')
