%  A Matlab program to evaluate several different Ekman layer solutions.
%  by Jim Price, starting in 1996.   Not carefully tested in all respects,
%  so watch out!

 clear 

 kk = menu('Which case?','lotus3', 'EBC', '10N', '10Nsmall')

 nn = menu(' Choose your K', 'Specify K', 'Specify Q, tau', 'paramk(u*,f)')
  
 lat = input(' Enter the latitude (degrees)  ')
 f = 2.*7.292e-5*sin(lat*pi/180.)

 if nn == 1  %  just put k in

 k = input(' Enter the diffusivity (MKS units, e.g., 100e-4)  ')
 tau = input(' Enter the wind stress (Pa)')
 ustar = sqrt(tau/1028.)

 end

 if nn ==2   %  compute k by analogy with the SEL4 model

 Q = 600.;
 Q = input(' Enter the noon heat flux (W/m*m)')
 Qstar = 9.8*0.3*Q/(1025.*1025.*4100.)
 Pq = 10.*3.6e3;
 tau = input(' Enter the wind stress (Pa)  ')
 ustar = sqrt(tau/1028.)
 Ptau = (1/f)*sqrt(2 - 2*cos(f*Pq/2))

 dtrap = ustar*ustar*Ptau/sqrt(Qstar*Pq/2)
 
 d1 = 1.0;
 k = (f/2)*dtrap*dtrap*d1*d1

 end

if nn == 3   %  compute k by neutral, turbulent scaling 

 tau = input(' Enter the wind stress (Pa)  ')
 ustar = sqrt(tau/1028.)

 c2 = 0.015     %  this is the similarity coefficient that seems to work
 k = c2*ustar*ustar/f

end

 h = input(' Enter the thermcoline depth (m)  ')  

 dz = 1.;
 z = 0:dz:h;
 z = -z' - dz/2.;
 [n3, n31] = size(z);


 n = menu(' Choose your BC:','Infnite depth', 'Free slip at H', ...
   'Zero slip at H')


% The purest Ekman layer, the original model, infinite depth.

 if n == 1 
 
 de = sqrt(2.*k/f)
 uamp = (ustar^2)/(f*de)
 
 u = uamp*exp(z/de).*exp(i*z/de)*(1 + i);

 ur = real(u);
 vr = imag(u); 
 
 end
 
 
 if n == 2
 
% Free slip at z = H.
% In this case the shear is made zero at the depth h to give tau = 0.     
% Gonella's (1971) form follows (wind stress is made northward). 
 
 de = sqrt(2.*k/f);  
 deg = 2.*pi*sqrt(2.*k/f)    % This is not your usual Ekman layer depth!
 
 hn = h/deg;                % hn is the depth is made non-d by the Ekman depth.
 zabs = abs(z)/deg;         % z is positive down here.

 denom = cosh(4*pi*hn) - cos(4*pi*hn);         
 a = (ustar*ustar/sqrt(2.))*exp(i*pi/4.)*(exp(4.*pi*hn) - exp(i*4.*pi*hn))/denom;
 b = (ustar*ustar/sqrt(2.))*exp(i*pi/4.)*(exp(-i*4.*pi*hn) - exp(-4.*pi*hn))/denom;        
         
 u = a*exp(-2.*pi*(i + 1.)*zabs) + b*exp(2.*pi*(i + 1.)*zabs); 
 
 u = u/sqrt(2.*k*f);         % Dimensionalize the current, time scale = inertial period.
 ur = real(u);
 vr = imag(u);


%  JPs simplified form follows
                                                   
 deg = sqrt(2.*k/f)    % This is not your usual Ekman layer depth! 
 hn = h/deg;           % hn is the depth is made non-d by the Ekman
 zabs = abs(z)/deg;    % z is positive down here.

 denom = cosh(2.*hn) - cos(2.*hn);
         
 a = (ustar*ustar/sqrt(2.))*exp(i*pi/4.)*( exp(2.*hn) - ...
          exp(i*2.*hn) )/denom;
 b = (ustar*ustar/sqrt(2.))*exp(i*pi/4.)*( exp(-i*2.*hn) - ...
          exp(-2.*hn) )/denom;        
         
 u = a*exp(-(i + 1.)*zabs) + b*exp((i + 1.)*zabs); 
 
 u = u/sqrt(2.*k*f);   % Dimensionalize the current
                                                  
 ur = real(u);
 vr = imag(u);

%  MII

 deg = sqrt(2.*k/f)    % Ekman layer depth 
 hn = h/deg;           % hn is the depth is made non-d by the Ekman
 zabs = abs(z)/deg;    % z is positive down here.

 denom = cosh(2.*hn) - cos(2.*hn);
         
 a = exp(i*pi/4.)*( exp(2.*hn) - exp(i*2.*hn) )/denom;
 b = exp(i*pi/4.)*( exp(-i*2.*hn) - exp(-2.*hn) )/denom;        
         
 u =  (ustar^2/sqrt(4*k*f))*  ...
               (a*exp(-(i + 1.)*zabs) + b*exp((i + 1.)*zabs)); 
                                                   
 ur = real(u);
 vr = imag(u);

%  MIII

 deg = sqrt(2.*k/f)        % Ekman layer depth 
 alpha = h/deg;            % ratio h to deg
 hn = 1;                   % hn is the t-cline depth made non-d by itself!
 zabs = abs(z)/h;    % z is positive down here.

 denom = cosh(2.*alpha) - cos(2.*alpha);
         
 a = exp(i*pi/4.)*( exp(2.*alpha) - exp(i*2.*alpha) )/denom;
 b = exp(i*pi/4.)*( exp(-i*2.*alpha) - exp(-2.*alpha) )/denom;        
         
 u =  (ustar^2/(h*f*sqrt(2))*alpha*  ...
   ( a*exp(-(i + 1.)*zabs*alpha) + b*exp((i + 1.)*zabs*alpha)) ); 

 un = (ustar^2)/(h*f)
                                                   
 ur = real(u);
 vr = imag(u);

 end


 if n == 3 

% No slip at z = h. 
% The shallow water case in which the current is zero on a boundary at z = h.

 a = sqrt(f/(2*k))
 de = pi/a

 fact = ustar*ustar*de/(k*pi);
 denom = cosh(2*a*h) + cos(2*a*h) 
 alpha = fact*( cosh(a*h)*cos(a*h) + sinh(a*h)*sin(a*h) )/denom;
 beta = fact*(cosh(a*h)*cos(a*h) - sinh(a*h)*sin(a*h) )/denom; 
 
 hmz = h - abs(z);
 
 ur = alpha*sinh(a*hmz).*cos(a*hmz) - beta*cosh(a*hmz).*sin(a*hmz);
 vr = alpha*cosh(a*hmz).*sin(a*hmz) + beta*sinh(a*hmz).*cos(a*hmz); 
         
 end 
 

% Diagnostics and plotting.
       
% Check to see if the transport above z = h is Ekman transport.

 ektran = ustar*ustar/f
 eksum = sum(ur)*dz 
 ysum = sum(vr)*dz

 figure(1)
 clf reset
 
 subplot(2,2,1)
 plot(ur, z,'-', vr,z,'--')
 xlabel('Ekman layer currents')
 ylabel('depth, m')
 title('Ekman layer current profiles')
 grid
 
 subplot(2,2,2)
 z5 = 0:5:h;
 z5 = -(z5 + 0.5)';
 u5 = z5;
 v5 = z5;
 u5 = interp1(z,ur,z5);
 v5 = interp1(z,vr,z5);
 plot(ur, vr,'-', u5, v5,'*')
 axis('equal')
 axis('square')
 title('hodograph')
 grid

 subplot(2,2,3)
 x3 = ones(2,n3); 
 zero = 0.*x3;
 y3 = x3;   
 z3 = x3; 

 for m=1:n3;
 x3(1,m) = 0.;
 x3(2,m) = ur(m);
 y3(1,m) = 0.;
 y3(2,m) = vr(m);
 z3(1,m) = z(m);
 z3(2,m) = z(m); 
 end;

 plot3(x3,y3,z3,'g-'); 
 xlabel('east speed, m/s')
 ylabel('north speed, m/s')
 zlabel('depth, m') 

 clear x3 y3 z3;
 x3 = [0 0 0]; y3 = [0 0 0.03]; z3 = [-h 0 0]; 
 hold on
 plot3(x3, y3, z3, 'r+', x3, y3, z3, 'r-' )

 %  compute the Flatness, F

 zc = z; uc = ur; vc = vr;

 dz = 6. 
 zup = -46:2:-2;
 [m n] = size(zup);
 
 for j=2:n
 zj = zup(j) - dz/2.;
 vj = interp1(zc,vc,zj);
 uj = interp1(zc,uc,zj); 
 za = zj + dz/2.;
 va = interp1(zc,vc,za);
 ua = interp1(zc,uc,za);
 zb = zj - dz/2.;
 vb = interp1(zc,vc,zb);
 ub = interp1(zc,uc,zb);
 spd = sqrt(uj*uj + vj*vj);
 spda = sqrt(ua*ua + va*va);
 spdb = sqrt(ub*ub + vb*vb);
 dvdz = (spda - spdb)/dz;
 dtheta = acos(  (ua*ub + va*vb)/(spda*spdb) )/dz;
 F(j) = dvdz/(spd*dtheta);
 zF(j) = zj; 
 end

 F(1) = F(2); zF(1) = zF(2);

 subplot(2,2,4)
 axis([ 0 3 -60 0]); plot(F, zF); axis([ 0 3 -60 0])
 xlabel('Flatness')
 ylabel('depth, m')
  

 figure(2)    %  the three-d masterpiece used in the  paper
 clf reset 


 set(0,'DefaultTextFontSize',16)
 set(0,'DefaultAxesFontSize',16)
 set(0,'DefaultAxesLineWidth',1.4)
 set(0,'DefaultLineLineWidth',1.2)

 xlabel('0.01 m s^{-1}')
 ylabel('0.01 m s^{-1}')
 zlabel('depth, m')
  

if kk == 1
 zc = [  -5 -10 -15 -25]
 zebc = [-0.1 -5 -10 -15 -20 -25 -30 -35 -40 -45 -50]
 zflag = [0 1 1 1 0 1 0 0 0 0 0 0 0]
 zn = length(zebc); 
 dz = 5
 text(-1, 6, 9, 'a','fontsize',20);
 zb = -50 
end

if kk == 2
 zc = -[ 8 12 16 20 24 28 32 36 40 44 48 50];   %  ebc depths
 zebc = -[0.1 4 8 12 16 20 24 28 32 36 40 44 48 50];
 zflag = [0 0 1 1 1 1 1 1 1 1 1 1 1 1]
 zn = length(zebc);
 dz = 4
 zb = -50.
 text(-1, 6, 9, 'b','fontsize',20); end   %  ebc text
end

if kk >= 3
 zc = [ -20 -30 -40 -50 -60  -70 -80]
 zebc = [ -0.1 -10 -20 -30 -40 -50 -60  -70 -80]
 zflag = [0 0 1 1 1 1 1 1 1]
 zn = length(zebc);
 text(-1, 6, 9, 'c','fontsize',20)
 dz = 10
 zb = -80
end

  axis([ -1 6 -1 6 zb 0])  
for m=1:zn-1;
 hold on
 m4 = (m-1)*4 + 1;
 arrow53( [0 0 zebc(m)], [100.*ur(m4) 100.*vr(m4) zebc(m)], ...
    0.04, 2.0, zflag(m))

 end
% add wind stress
 arrow53([0 0 3], [0 2.5 3], 0.04, 2.0)
% add a central stem
 arrow53( [0 0 0], [0 0 zb], 0.0, 2.0 ) 
view([-4 -5 2.])
 xyl = 6


 xx = [zb xyl xyl]; yy = [xyl xyl -1]; zz = [zb zb zb]; % lower rear
 hold on; plot3(xx, yy, zz, ':'); clear xx yy zz
 xx = [xyl xyl]; yy = [xyl xyl]; zz = [zb 0]; % back corner
 hold on; plot3(xx, yy, zz, ':'); clear xx yy zz
 xx = [xyl xyl]; yy = [-1 -1]; zz = [zb 0];  %  rt front edge
 hold on; plot3(xx, yy, zz, 'LineWidth', 1.0); clear xx yy zz
 xx = [zb xyl xyl]; yy = [xyl xyl zb]; zz = [0 0 0]; % lower rear
 hold on; plot3(xx, yy, zz, 'LineWidth', 1.0); clear xx yy zz  %  top rear
 grid

 
 figure(3)    % yet another three-d masterpiece, all normed-up

 clf reset 

 un = (ustar^2)/(h*f)

 xlabel('U/U_n')
 ylabel('V/U_n')
 zlabel('depth/H')
 axis([ -1 5 -1 5 -50/h 0]) 
 for m=1:3:n3;
 hold on
 arrow53( [0 0 z(m)/h], [ur(m)/un vr(m)/un z(m)/h], 0.04, 2.0)
 end
% add wind stress
 arrow53([0 0 3/h], [0 2.5 3/h], 0.04, 2.0)
% add a central stem
 arrow53( [0 0 0], [0 0 -50/h], 0., 1.7) 
 view([-4 -5 2.])
 grid

 xx = [-1 5 5]; yy = [5 5 -1]; zz = [-1 -1 -1]; % lower rear
 hold on; plot3(xx, yy, zz, ':'); clear xx yy zz
 xx = [5 5]; yy = [5 5]; zz = [-1 0]; % back corner
 hold on; plot3(xx, yy, zz, ':'); clear xx yy zz
 xx = [5 5]; yy = [-1 -1]; zz = [-1 0];  %  rt front edge
 hold on; plot3(xx, yy, zz, 'LineWidth', 1.0); clear xx yy zz
 xx = [-1 5 5]; yy = [5 5 -1]; zz = [0 0 0]; % lower rear
 hold on; plot3(xx, yy, zz, 'LineWidth', 1.0); clear xx yy zz  %  top rear