%  upperEk;  an elementary upper ocean Ekman layer
%  Solve a 1-D diffusion flow numerically.  This
%  version has rotation, and thus makes an Ekman 
%  layer.  The IC is a state of rest.  The fluid is forced
%  by an imposed stress at the top of the column (z=0)
%  The lower boundary condition is free-slip to minimize
%  the effect of finite depth. 
%
%  There has been no effort to make this efficient or
%  the least bit sophisticated.  
%  Jim Price, Oct 99

clear
close all

%  display the text string str2
str2(1) =  {'upperEk:' };
str2(2) =  {'  '};
str2(3) =  {'An upper ocean Ekman layer driven by a  '};
str2(4) =  {'northward stress.  Solved by an elementary'};
str2(5) =  {'numerical method.  You will be asked to enter'};
str2(6) =  {'the latitude only.  To change other parameters,'}; 
str2(7) =  {'K, tau, dz, etc., see the script. '}; 
str2(8) =  {''}; 
str2(9) =  {''};
str2(10) = {''};
str2(11) = {'Jim Price, March, 2000.'};

hf3 = figure(1);
clf reset
% set(hf3,'Position',[50 200 550 550])
set(gca,'Visible','off')
text(0, 0.50, str2,'FontSize',13,'HorizontalAlignment','left')

dz = 2.;          %  the grid interval, m
L = 100.;         %  the grid extent (length), m  
z = 0:-dz:-L ;    %  the grid defined

stresBC = 1;
freeslip = 1;

ndays = 10.;     %  days to integrate
K = 1.e-2;       %  the diffusivity (MKS, appropriate for the upper ocean)
dt = 0.4*dz^2/K; %  the time step, set so that 
w = dt*K/dz^2    %  w is less then 0.5 for numerical stability
nstep = round(ndays*8.64e4/dt);

%  specify the latitude here
lat = input('Enter the latitude (degrees)  ');
f = 2*7.29e-5*sin(lat*pi/180);
IP = (2*pi/f)/8.64e4;

rho = 1025.;     %  nominal constant density of water
%  the (northward) stress
tauBC = 0.1/rho; %  imposed constant stress/density as upper BC
Ektrans = tauBC/f
Ekdepth = sqrt(2*K/f)

u = 0*z;         %  the dependent variable, a velocity
v = u;
[nn nz] = size(u)
uavg = zeros(1,nz);
vavg = uavg; 
navg = 0;

figure(1)
clf reset
title('upper ocean Ekman layer, \tau_y = 0.1 Pa')

set(0,'DefaultLineLineWidth',1.4)
set(0,'DefaultTextFontSize',12)
set(0,'DefaultAxesLineWidth',1.3)
set(0,'DefaultAxesFontSize',12)

%  begin time-stepping

 for j=1:nstep
       
 delsqu = [0 diff(u,2) 0];    
 delsqv = [0 diff(v,2) 0];

u = u + w*delsqu + dt*f*v;
v = v + w*delsqv - dt*f*u;

% apply the surface BC next
    u(1) = u(2);
    v(1) = v(2)  + dz*tauBC/K;
% the bottom BC (free-slip)
    u(nz) = u(nz-1);
    v(nz) = v(nz-1);
 
 %  compute and store some diagnostics

time(j) = (j-1)*dt;


%  sum the currents to make a time-average
if time(j)/8.64e4 >= 1.
   uavg = uavg + u;
   vavg = vavg + v;
   navg = navg + 1;   
end


%  time series of surface current and the transport
usurf(j) = u(1);
vsurf(j) = v(1);
mid = round(2*nz/3) -1;
umid(j) = u(mid);
vmid(j) = v(mid);
udeep(j) = u(nz-3);
vdeep(j) = v(nz-3);

transu(j) = (sum(u) - u(1))*dz;
transv(j) = (sum(v) - v(1))*dz;

topstress(j) = -K*(u(2) - u(1))/dz;
botstress(j) = -K*(u(nz) - u(nz-1))/dz;

nv1 = 2;
nv2 = nz-1;

topwork(j) = - (u(nv1)*K*(u(nv1+1) - u(nv1-1))/(2*dz) +...
  v(nv1)*K*(v(nv1+1) - v(nv1-1))/(2*dz))  ;
botwork(j) =  u(nv2)*K*(u(nv2+1) - u(nv2-1))/(2*dz) + ...
   v(nv2)*K*(v(nv2+1) - v(nv2-1))/(2*dz);

ds = 0.;
kes = 0.;
for m=nv1:nv2
   shru = (u(m+1) - u(m))/dz;
   shrv = (v(m+1) - v(m))/dz;
   ds = ds - dz*K*(shru^2 + shrv^2);
   kes = kes + 0.5*dz*(u(m)^2 + v(m)^2);
end
diss(j) = ds;
ke(j) = kes;

efold = 1 - exp(-1);
ssum = cumsum((sqrt(u.*u + v.*v)+1.e-4));
ze(j) = -interp1(ssum,z,efold*ssum(nz), 'linear');
zK(j) = sqrt(0.5*K*time(j));


%  vorticity budget over a volume from nv1 to nv2

avgvort(j) = u(nv1) - u(nv2);
topvflux(j) = -K*delsqu(nv1)/dz^2;
botvflux(j) = -K*delsqu(nv2)/dz^2;
netvflux(j) = botvflux(j) - topvflux(j);
tippingx(j) = f*(u(nv1) - u(nv2));

avgvorty(j) = v(nv1) - v(nv2);
topvfluxy(j) = -K*delsqv(nv1)/dz^2;
botvfluxy(j) = -K*delsqv(nv2)/dz^2;
netvfluxy(j) = botvfluxy(j) - topvfluxy(j);
tippingy(j) = f*(v(nv1) - v(nv2)); 

%  plot profiles every nplot steps
nplot = 200;
if mod(j,nplot) == nplot - 1
   hold on
    
 figure(1)  
 subplot(2,1,1)
 hold on
 plot(u, z)
 ylabel('depth, m')
 xlabel('east current, m/sec')
 title('Upper ocean Ekman layer, \tau_y = 0.1 Pa')
    
 subplot(2,1,2)  
 hold on 
 plot(v, z)
 ylabel('depth, m')
 xlabel('north current, m/sec')
 drawnow

end   % if on whether to plot this time step

end   %  loop on j for time step


%  compute and plot diagnostics

%  time-integrate the top and bottom stress to compare with the
%   transport

totstr = cumsum(topstress-botstress)*dt;
%botstr = cumsum(botstress)*dt;

figure(4)
clf reset
subplot(2,1,1)
plot(time, ze, time, zK,':')
xlabel('time, s')
ylabel('bl thickness, m')
legend('e-folding', 'sqrt(0.5*K*t)')



%  energy budget
figure(6)
clf reset

%subplot(2,1,1)
dkedt = diff(ke)/dt;
tdkedt = cumsum(diff(time));
plot(time, topwork,':', time, botwork,'-.', ...
    time, diss, '--', tdkedt, dkedt)
legend('top work rate', 'bot work rate', 'dissipation rate','dKE/dt')
ylabel('Work rate, m*m/s*s*s')
xlabel('time, s')
title('energy budget, 20-30 m')
% axis([0 4e6 -5e-5 5e-5])

%  vorticity budget
figure(7)
clf reset

subplot(2,1,1)
dvortdt = diff(avgvort)/dt;
timedv = cumsum(diff(time));


plot(time, -netvflux, '-', timedv, dvortdt, '-.', time, tippingx, '--')
%  note, the sign on net flux is negative ????

ylabel('vorticity fluxes,x')
legend('-div(flux)', 'dvort/dt','tipping')
title('Overall vorticity budget')
axis([0 6e5 -1.e-5 1.e-5])

subplot(2,1,2)
dvortdty = diff(avgvorty)/dt;
timedv = cumsum(diff(time));
plot(time, netvfluxy, '-',  ...
   timedv, dvortdty, '-.', time, tippingy, '--')
xlabel('time, s')
ylabel('vorticity fluxes, y')

legend('-div(flux)','dvort/dt', 'tipping')
axis([0 6e5 -1.e-5 1.e-5])


figure(15)
clf reset
subplot(2,1,1)
plot(time, usurf, time, umid, time, udeep)
ylabel('east, m/s')
grid
title('east and north currents, various depths')
subplot(2,1,2)
plot(time, vsurf, time, vmid, time, vdeep)
ylabel('north, m/s')
xlabel('time, s')
grid


figure(11)
clf reset

plot(transu, transv)
hold on
plot([0 Ektrans], [0 0],'r+')
Tx = mean(transu);
Ty = mean(transv);
plot([0 Tx], [0 Ty], 'g')
xlabel('U trans ')
ylabel('V trans')
% axis([-0.1 0.15 -0.25 0.])
axis('equal')
title('Transport with time, and average')

figure(12)
clf reset

uavg = uavg/navg;
vavg = vavg/navg;

plot(uavg, z, vavg, z)
xlabel('east and north currents, m/sec')
ylabel('depth, m')
title('Mean current profiles')

figure(13)
clf reset

 x3 = ones(2,nz); 
 zero = 0.*x3;
 y3 = x3;   
 z3 = x3; 

 for m=1:nz;
 x3(1,m) = 0.;
 x3(2,m) = uavg(m);
 y3(1,m) = 0.;
 y3(2,m) = vavg(m);
 z3(1,m) = z(m);
 z3(2,m) = z(m); 
 end;

hp =  plot3(x3,y3,z3,'b-'); 
set(hp,'LineWidth',1.3)
 xlabel('east speed, m/s')
 ylabel('north speed, m/s')
 zlabel('depth, m') 
 
 %  show the wind stress direction
 clear x3 y3 z3;
 x3 = [0 0 0]; y3 = [0 0 0.1]; z3 = [-L 0 0]; 
 hold on
 hp =  plot3(x3, y3, z3, 'r-' );
 set(hp,'LineWidth',1.5)
 grid
 

figure(14)
clf reset

% plot(uavg, vavg, '*')
hold on
plot([0 0], [0 00.06], 'r')

for j=1:nz
   uh = [0 uavg(j)];
   vh = [0 vavg(j)];
   hold on 
   plot(uh, vh, 'LineWidth', 1.3)
end

xlabel('east speed, m/s')
ylabel('north speed, m/s')
title('Hodograph, dz = 2.5 m')
