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