%  Solve a 1-D Couette flow numerically.  
%  The IC is a state of rest.  The fluid is forced
%  by the BC at the top of the column (z=0), either
%    a constant velocity, or,
%    a constant stress, determined by stresBC.
%  The lower boundary condition is either 
%    u = 0 (as if a fixed wall and no-slip), 
%    or dudz = 0 (as if free slip).
%
%  There has been no effort to make this efficient or
%  the least bit sophisticated. The plotting is mostly
%  done in the dimensional variables, mainly because
%  the BC options make it hard to use non-d variables.  
%
%  Jim Price, Oct 99;  Modified in Dec 01.  
%
%  This is public domain.
clear
close all

set(0,'DefaultLineLineWidth',1.2)
set(0,'DefaultTextFontSize',13)
set(0,'DefaultAxesLineWidth',1.1)
set(0,'DefaultAxesFontSize',13)

%  the text can go to about here                       > X
str2(1) =  {'Couette:' };
str2(2) =  {'  '};
str2(3) =  {'Solve for diffusive flow in a channel, sometimes  '};
str2(4) =  {'called Couette flow, by numerical methods. '};
str2(5) =  {'The only input is required to specify the upper'};
str2(6) =  {'and lower boundary conditions.  The upper BC '}; 
str2(7) =  {'drives the flow by either an imposed speed or'}; 
str2(8) =  {'an imposed stress;  the lower BC is either '};
str2(9) =  {'no-slip or free-slip.  Solutions are compared to'};
str2(10) = {'a similarity solution that assumes an infinite'};
str2(11) = {'domain, and various energy and vorticity '};
str2(12) = {'budgets are evaluated.  The dimensional values'};
str2(13) = {'are consistent with a turbulent upper ocean.'};
str2(14) = {''};
str2(15) = {'This code is public domain.'};
str2(16) = {''};
str2(17) = {'Jim Price, March, 2000.'};

hf3 = figure(10);
clf
set(hf3,'Position',[100 100 600 600])
set(gca,'Visible','off')
text(0, 0.50, str2,'FontSize',13,'HorizontalAlignment','left')

%  set some parameters
%   here we use dimensional variables, with scales
%   like those of a 'turbulent' upper ocean

dz = 1.;          %  the grid interval, m
L = 50.           %  the grid extent (length), m  
z = 0:-dz:-L ;    %  the grid defined
s = 0*z;          %  the dependent variable, think of it as velocity
[nn ns] = size(s)
delsq = s;

nstep = 3000;    %  the number of time steps
dt = 25.;        %  the time step, sec
K = 1.e-2;       %  the diffusivity (MKS, appropriate for the upper ocean)
w = dt*K/dz^2    %  must be less then 0.5 for numerical stability
rho = 1025.;     %  nominal constant density of water


%  define the type of BC:
Uo = 0.1;           %  an imposed speed upper BC
tauBC = 0.1/1025;   %  imposed constant stress/density as upper BC

stresBC = menu('Choose the upper BC:', 'Imposed speed', ...
'Imposed stress')
stresBC = stresBC - 1;

freeslip = menu('Choose the lower BC:', 'No slip', 'Free slip')
freeslip = freeslip - 1;


figure(1)
clf reset

km = 1:ns-2;
kp = km + 2;
kmid = 2:ns-1;

%  begin time-stepping;  the 'scheme' used here is forward in time
%   and centered in space - very crude, and it works pretty well!

for j=1:nstep
   
%  estimate the Laplacian
 delsq(kmid) = s(km) + s(kp) - 2.*s(kmid);
   
 delsq(1) = 0.;
 delsq(ns) = 0.;
 
 %  step forward
 
 s = s + w*delsq;
   
 % now set the upper boundary condition
 if stresBC == 0
    s(1) = Uo;
 else  
    s(1) = s(2) + dz*tauBC/K;
 end
 
 % the lower boundary condition
 if freeslip == 0
    s(ns) = 0.;
 else
    s(ns) = s(ns-1);
 end

%  compute and store some diagnostics

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

topstress(j) = -K*(s(2) - s(1))/dz;
botstress(j) = -K*(s(ns) - s(ns-1))/dz;
trans(j) = (sum(s) - s(1))*dz;

%  do a volume integrals over these limits 
nv1 = round(20/dz);
nv2 = round(30/dz);

topwork(j) = -s(nv1)*K*(s(nv1+1) - s(nv1-1))/(2*dz);
botwork(j) =  s(nv2)*K*(s(nv2+1) - s(nv2-1))/(2*dz);
ds = 0.;
kes = 0.;
for m=nv1:nv2
   ds = ds - dz*K*(((s(m+1)-s(m))/dz)^2);
   kes = kes + 0.5*dz*(s(m)^2);
end
diss(j) = ds;
ke(j) = kes;

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


%  vorticity budget over this volume 

avgvort(j) = s(nv1) - s(nv2);
topvflux(j) = -K*delsq(nv1)/dz^2;
botvflux(j) = -K*delsq(nv2)/dz^2;
netvflux(j) = botvflux(j) - topvflux(j);

%  plot profiles every nplot steps
nplot = 100;
if mod(j,nplot) == nplot - 1
   hold on
   znd = z/sqrt(K*time(j));
   Und = Uo;
   if stresBC == 1
      Und = tauBC*(tauBC/K);
end
   
 figure(1)  
 subplot(2,1,1)
 hold on
 plot(s, z)
    
 subplot(2,1,2)  
 hold on 
 if stresBC == 0
    plot(s/Und, znd)
 else  
    for m=1:ns-1
       shear(m) = -(s(m+1) - s(m))/dz;
       zshear(m) = 0.5*(z(m+1) + z(m));
    end 
 plot(shear/(tauBC/K), zshear/sqrt(K*time(j)))
end


%  some things done once to figure 1:
if j <= 2*nplot
% evaluate and plot the similarity solution, even if it
%  may not be appropriate!

deta = 0.01;
eta = 0:deta:5;
e2 = exp(-0.25*(eta.^2));
integ = cumsum(e2)*deta;
usim = (1 - integ/sqrt(pi));
plot(fliplr(usim), fliplr(-eta),'r.', 'LineWidth',1.6)
axis([-0.2 1. -6 0])
xxx = plot([-99 -90], [-93 -90], 'b', [-99 -97], [-93 -99], 'r', 'LineWidth',2.0)
legend(xxx, 'numerical soln', 'similarity soln',4)

%  add scales and labels
subplot(2,1,1)
xlabel('U, m/s')
ylabel('Z, m')
delt = nplot*dt/(L^2/K);
tit = strcat('Couette flow, numerically, dt = ',num2str(delt), ...
   ' (non-d)')
title(tit)

subplot(2,1,2)
if stresBC == 0
   xlabel('current, U/U_o')
else
   xlabel('shear, \partialU/\partialz/(U_{*}^{2}/K)')
end
ylabel('Z/sqrt(K*t)')
end   %  if on things done once to fig 1

drawnow   %  flush the print buffer
end   % if on whether to plot this time step

end   %  loop on j for the time step


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

totstr = cumsum(topstress-botstress)*dt;

figure(4)
clf reset

subplot(3,1,1)
plot(time, topstress, time, botstress,':')
%axis([ 0 4e6 0 2e-5])
ylabel('kin stress, m*m/s*s')
legend('top stress', 'bottom stress')
title('overall momentum budget')

subplot(3,1,2)
plot(time, totstr, time, trans, ':')
ylabel('stress*t m*m/s')
legend('integrated stress', 'transport')


subplot(3,1,3)
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

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')

%  vorticity budget
figure(7)
clf reset

dvortdt = diff(avgvort)/dt;
timedv = cumsum(diff(time));

plot(time, netvflux, '-', time, botvflux, '--', time, topvflux,':', ...
   timedv, dvortdt, '-.')
xlabel('time, s')
ylabel('vorticity fluxes')
legend('-div(flux)', 'bottom v flux', 'top v flux','dvort/dt')
title('Vorticity budget, Z = 20-30 m')