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