%  bottomEk:  solve for a bottom Ekman layer.
%  Solve a 1-D diffusion flow numerically.  This
%  version has rotation, and thus makes an Ekman 
%  layer.  The IC is a geostrophic current. The fluid 
%  is forced by the BC at the top of the column (z=0), 
%  the constant gesotrophic velocity (and the associated
%  pressure gradient) and by a no-slip condtion on the bottom. 
%
%  There has been no effort to make this efficient or
%  the least bit sophisticated. Public domain for all
%  personal, educational purposes. 
%  Jim Price, Oct 99.  modified Dec 2001.

clear
close all


%  display the text string str2
str2(1) =  {'bottomEk:' };
str2(2) =  {'  '};
str2(3) =  {'A bottom Ekman layer driven by an imposed  '};
str2(4) =  {'gesotrophic current, Uo, and 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) =  {'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')



Uo = 0.1;         %  the free stream (geostrophic) current
dz = 2.;          %  the grid interval, m
L = 100.;         %  the grid extent (length), m  
z = 0:-dz:-L ;    %  the grid defined
[nn nz] = size(z)


u = Uo*ones(1,nz);
v = 0*u;

stresBC = 0;    % this is for an imposed speed at the upper boundary
freeslip = 0;   % for a no-slip BC along the bottom

K = 0.5e-2;       %  the diffusivity (MKS)
  %  this K is appropriate for a turbulent, ocean BBL ocean
ndays = 10.;      %  the number of days to integrate
dt = 0.4*(dz^2)/K;       %  the time step, sec
nstep = ndays*8.64e4/dt;     %  number of time steps needed to take
w = dt*K/dz^2    %  must be less then 0.5 for numerical stability
rho = 1025.;     %  nominal constant density of water

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

PGy = Uo*f;      %  the pressure gradient needed to balance Uo

uavg = zeros(1,nz);
vavg = uavg; 
navg = 0;
jd = 0;

figure(1)
clf reset

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
    
%  compute the diffusion tendency
    
 for k=2:nz-1
    delsqu(k) = u(k-1) + u(k+1) - 2.*u(k);
    delsqv(k) = v(k-1) + v(k+1) - 2.*v(k);
 end
 delsqu(1) = 0.;
 delsqu(nz) = 0.;
 delsqv(1) = 0.;
 delsqv(nz) = 0.;
 
%  time step
 
u = u + w*delsqu + dt*f*v ;
v = v + w*delsqv + dt*PGy - dt*f*u;
   
% apply the upper boundary condition
    u(1) = Uo;
    v(1) = 0.;
 
 % apply the no-slip lower boundary condition
    u(nz) = 0.;
    v(nz) = 0.;
 
%  sum the current to make a time-average 
 
if dt*j/8.64e4 >= 1.
   uavg = uavg + u;
   vavg = vavg + v;
   navg = navg + 1;   
end 

% store some other diagnostics
 
 if mod(j,5) == 1
    
    jd = jd + 1;
    timed(jd) = (j-1)*dt/8.64e4;

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

%  the transport and top and bottom stresses

transu(jd) = (sum(u) - Uo*nz)*dz - 0.5*(u(nz)+u(1))*dz;
transv(jd) = sum(v)*dz - 0.5*(v(nz)+v(1))*dz;

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

ushr = (u(nz) - u(nz-1))/dz;
vshr = (v(nz) - v(nz-1))/dz;
bstrx(jd) = -K*ushr;
bstry(jd) = -K*vshr;

%  compute the vorticity and energy budgets over a portion of
%   the water column

nv1 = round(2*nz/3);
nv2 = nz - 1;

topwork(jd) = - (u(nv1)*K*(u(nv1+1) - u(nv1-1))/(2*dz) +...
  v(nv1)*K*(v(nv1+1) - v(nv1-1))/(2*dz))  ;
botwork(jd) =  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.;
ps = 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);
   ps = ps + dz*PGy*v(m);
end
diss(jd) = ds;
pgwork(jd) = ps;
ke(jd) = kes;


%  the vorticity budget over a volume from nv1 to nv2

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

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

%  plot profiles every nplot steps
nplot = 100;
if mod(jd,nplot) == nplot - 1
 hold on

disp('  time elapsed to this point is ')
 tdays = timed(jd)
    
 figure(1)  
 subplot(2,1,1)
 hold on
 plot(u, z)
    
 subplot(2,1,2)  
 hold on 
 plot(v, z)
   
drawnow

end   % if on whether to plot this time step

end  %  if on whether to store diagnostics this step

end   %  loop on j for time step

%  finish figure 1
subplot(2,1,1)
xlabel('U, m/s')
ylabel('Z, m')
delt = nplot*dt/(L^2/K);
tit = strcat('1-D rotating flow, numerically, dt = ',num2str(delt), ...
   ' (non-d)');
title(tit)

subplot(2,1,2)
xlabel('V, m/s')
ylabel('Z, m')

%  plot diagnostics

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

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

figure(2)
clf reset

subplot(2,1,1)
plot(timed, 1000.*bstrx, timed, 1000.*bstry)
ylabel('bottom stress, m^2 s^{-2}')
legend('bot str x','bot str y')
axis([0 6 -0.0 0.1])

subplot(2,1,2)
plot(timed, transu, timed, transv)
ylabel('transport (-Uo)')
xlabel('time, days')
legend('east','north')

figure(3)
clf reset

plot(transu, transv)
xlabel('U trans')
ylabel('V trans')
axis([-1.5 0 0 1.5])
if lat == 10
   axis([-3 0 0 3])
end
axis('square')
grid

%  the ageo-transports in the x-direction often match up rather poorly

disp(' the time-mean ageostrophic transports are') 
trax = mean(transu)
tray = mean(transv)

disp(' the ageo transports from the final currents are')
txua = sum(u)*dz - L*Uo
tyva = sum(v)*dz

disp('  the time-avgeraged bottom stresses are ')
bxa = 1030*mean(bstrx)
bya = 1030*mean(bstry)

disp(' the Ekman transports computed from the bottom stress are')
trxek = -bya/(f*1030)
tryek = bxa/(f*1030)

disp(' the geostrophic drag coefficient is')
Cd = sqrt(bstrx(jd)^2 + bstry(jd)^2)/(Uo^2)

%  energy budget
figure(6)
clf reset

dkedt = diff(ke)/(8.64e4*(timed(2)- timed(1)));
tdkedt = cumsum(diff(timed));
plot(timed, topwork,':', timed, botwork,'-.', ...
    timed, pgwork,':', timed, diss, '--', tdkedt, dkedt)
legend('top work rate', 'bot work rate', 'P grad work', ...
  'dissipation rate','dKE/dt')
ylabel('Work rate, m*m/s*s*s')
xlabel('time, days')
title('Energy budget, lower third')
axis([0 6 -1e-5  1e-5])


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


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

figure(12)
clf reset

plot(uavg, z, 'b', vavg, z, 'b')
grid
xlabel('east and north current, m/s')
ylabel('depth, m')
title('numerical (blue) and analytical (red) current profiles')

De = sqrt(K/(f/2));
ut = Uo*(1 - exp(-(z+L)/De).*cos((z+L)/De));
vt = Uo*exp(-(z+L)/De).*sin((z+L)/De);
hold on
plot(ut, z,'r+', vt, z,'r+')



%  integrate these average currents to find the bl transport
trxa = 0.;
trya = 0.;
clear za uta vta 
Lbl = 200.;
dza = 0.2;
uta = zeros(1,1000);
vta = uta;
za = uta; zn = uta;

for m=1:1000
   za(m) = -Lbl + (dza*m - dza);
   zn(m) = (Lbl + za(m))/De;  
   uta(m) = Uo*(1. - exp(-zn(m))*cos(zn(m)));
   vta(m) = Uo*exp(-zn(m))*sin(zn(m));   
trxa = trxa + uta(m)*dza;
trya = trya + vta(m)*dza;
end

disp(' ageostrophic transports from an analytic solution are')
trxa = trxa - Lbl*Uo
trya


figure(13)
clf reset

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

 for m=1:2: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 geo current direction
 clear x3 y3 z3;
 x3 = [0 0 Uo]; y3 = [0 0 0.0]; z3 = [-L 0 0]; 
 hold on
 hp =  plot3(x3, y3, z3, 'r-' );
 set(hp,'LineWidth',1.5)
 grid
 axis([-0.05 0.15 -0.05 0.15 -L 0])
 view(-20., 30.)
 
figure(14)
clf reset

hold on 

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

xlabel('east speed, m/s')
ylabel('north speed, m/s')
title('Hodograph')
axis('equal')