% fredgui; geostrophic adjustment of a particle on a slope % % integrate the Fredholm problem for a single particle on a % rotating earth, numerically. This is useful for seeing how an % Ekman balance or a geostrophic balance appears as the time average % of a start up problem, or, as the instantaneous balance if there % is some friction and you can wait a bit. May also be useful as a % first look at some forms of numerical time integration of ODEs. % % By Jim Price, October, 1995; modified in June 2000 % Public domain for educational purposes. % % The equations of motion for a particle on a rotating earth are: % du/dt - f*v = Fx/rho - k*u, for the east component, and % dv/dt + f*u = Fy/rho - k*v, for the north component. % % d()/dt is an ordinary time derivative since this model is meant % to describe the motion of a single particle (and not a fluid!), % u and v are the east and north velocity components of the particle, % Fx and Fy are the east and north forces per unit volume applied % to the particle (in a fluid problem these could be boundary stress), % rho is the nominal density of sea water, % f is the usual Coriolis parameter, and % k is a linear decay coefficient that represents, say, bottom drag. clear all close all % display the text string str2 str2(1) = {'fred4:' }; str2(2) = {'' }; str2(3) = {'Solve for the motion of a particle sitting on a '}; str2(4) = {'sloping surface and subject to Coriolis acceleration'}; str2(5) = {'and bottom drag. This demonstrates a local version'}; str2(6) = {'of geostrophic adjustment, and the effect of '}; str2(7) = {'friction on a more or less geostrophic motion. '}; str2(8) = {'You may change the latitude (deg), and the e-folding'}; str2(9) = {'time (days) for linear (Rayleigh) drag, and the '}; str2(10) = {'the initial velocity (= 0 at rest; =1 in '}; str2(11) = {'geostrophic balance; >1 super geostrophic).'}; str2(12) = {''}; hf3 = figure(9); clf set(hf3,'Position',[50 200 500 500]) set(gca,'Visible','off') text(-.05, 0.70, str2,'FontSize',12,'HorizontalAlignment','left') latitude = 20. hlat = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.18 0.35 0.06], ... 'FontSize', 10, ... 'Callback', 'eval(get(hlat,''string''))', ... 'String', 'latitude = 20 '); efold = 5.; hefold = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.10 0.35 0.06], ... 'Callback', 'eval(get(hefold,''string''))', ... 'FontSize', 10, ... 'String', 'efold = 5 '); ufrac = 0.; hufrac = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.02 0.35 0.06], ... 'Callback', 'eval(get(hufrac,''string''))', ... 'FontSize', 10, ... 'String', 'ufrac = 0'); hgo = uicontrol('Style', 'Pushbutton', 'Units', 'Normalized', ... 'Position', [0.80 0.10 0.1 0.06], ... 'Callback', 'uiresume(9)', ... 'FontSize', 10, ... 'String', 'Start'); hquit = uicontrol('Style', 'Pushbutton', 'Units', 'Normalized', ... 'Position', [0.80 0.02 0.1 0.06], ... 'Callback', 'lquit = 1; uiresume(9)', ... 'FontSize', 10, ... 'String', 'Quit'); lquit = 0; while (lquit == 0) disp (' ... waiting for control panel action.') uiwait(9) if lquit == 1; break; end; day = 8.64e4; % lat = input('Enter the latitude (deg)'); f = 2.*7.292e-5*sin(latitude*pi/180.); if abs(f) <= 1.e-10; f = 1.e-8; end % don't let f = 0 % The force (Fx, Fy) is evaluated as if it were a buoyanct force. % As now set up, the force is northward only i.e., forcex = 0 and % forcey > 0. h = 100.; % a layer thickness, m. slope = 0.01; g = 10.; delrho = 0.5; rho = 1000.; % nominal sea water density. % There is a linear drag or frictional force, that will cause an % e-folding of an unforced current in time efold (days). % efold = input('Enter the e-folding time for drag (days, 1 - 9999999)' ); k = 1./(efold*day); % compute the linear drag consistent with this efold. spd = 100.; % time steps per day ndays = 10.; % the number of days to run (10 is plenty). dt = (2*pi/f)/spd; % the time step, sec. nstep = ndays*(2*pi/f)/dt; % the number of time steps needed to go ndays. u = 1:nstep; v=1:nstep; t=1:nstep; % set up some vectors to store the solution. % the (geostrophic) scales used to plot things usc = g*(delrho/rho)*slope/f; xsc = 0.001*usc/f; tsc = 2*pi/f; zsc = 1000.*xsc*slope; esc = h*usc^2; u(1) = ufrac*usc; v(1) = 0; t(1) = 0.; % define the initial condition. fricwork(1) = 0.; % the initial position of the particle x(1) = 10*xsc; y(1) = xsc; % Begin a loop that will step the equations ahead in time. for i=1:nstep-1 t(i+1) = t(i) + dt; % increment the time. Fx = 0.; Fy = g*delrho*slope; % A fourth order Runga-Kutta method; very accurate but relatively expensive. du1 = dt*( Fx/rho + f*v(i) - k*u(i) ); dv1 = dt*( Fy/rho - (f*u(i) + k*v(i)) ); du2 = dt*( Fx/rho + f*(v(i)+dv1/2.) - k*(u(i)+du1/2.) ); dv2 = dt*( Fy/rho - (f*(u(i)+du1/2.) + k*(v(i)+dv1/2.)) ); du3 = dt*( Fx/rho + f*(v(i)+dv2/2.) - k*(u(i)+du2/2.) ); dv3 = dt*( Fy/rho - (f*(u(i)+du2/2.) + k*(v(i)+dv2/2.)) ); du4 = dt*( Fx/rho + f*(v(i)+dv3) - k*(u(i)+du3) ); dv4 = dt*( Fy/rho - (f*(u(i)+du3) + k*(v(i)+dv3)) ); u(i+1) = u(i) + du1/6. + du2/3. + du3/3. + du4/6.; v(i+1) = v(i) + dv1/6. + dv2/3. + dv3/3. + dv4/6.; forx(1,i+1) = Fx/rho; forx(2,i+1) = f*v(i); forx(3,i+1) = -k*u(i); fory(1,i+1) = Fy/rho; fory(2,i+1) = -f*u(i); fory(3,i+1) = -k*v(i); x(i+1) = x(i) + dt*u(i+1)/1000.; y(i+1) = y(i) + dt*v(i+1)/1000.; z(i+1) = y(i)*1000*slope; ke(i+1) = 0.5*h*(u(i+1)^2 + v(i+1)^2); pe(i+1) = -g*delrho*h*(y(i+1)*slope); fricwork(i+1) = fricwork(i) - k*h*(u(i+1)^2 + v(i+1)^2)*dt; esum(i+1) = ke(i+1) + pe(i+1) - fricwork(i+1); end % the end of the time stepping loop. % set default graphics things set(0,'DefaultTextFontSize',12) set(0,'DefaultAxesFontSize',12) set(0,'DefaultAxesLineWidth',1.2) set(0,'DefaultLineLineWidth',1.6) % make a movie-like plot of position and depth x3 = 0:4:80; y3 = 0:1:4; nx3 = max(size(x3)); ny3 = max(size(y3)); for i=1:nx3 for j=1:ny3 z3(j,i) = -y3(j)*xsc*slope*1000./zsc; end end notthis = 1 if notthis == 0 figure(4) clf reset hsurf4 = surfl (x3, y3, z3); set (hsurf4, 'FaceColor', [.8 .9 1.]); depth = -y*slope*1000.; hold on hp9 = plot3(x/xsc, y/xsc, (depth+20)/zsc, 'r'); set(hp9, 'Linewidth', 1.6) view(-70, 60); xlabel(' x/(g\alpha/f^2)') ylabel(' y/(g\alpha/f^2)') zlabel(' z/(g\alpha^2/f^2)') axis([0 60 0 4 -4 0]) title('dense particle on a rotating slope w/ friction') end m = 0; for j=1:20:nstep m = m + 1; if m == 1 figure(1) clf reset hsurf1 = surfl(x3, y3, z3); set (hsurf1, 'FaceColor', [.8 .9 1.]); depth = -y*slope*1000.; hold on hp9 = plot3 (x(1)/xsc, y(1)/xsc, (depth(1)+10)/zsc, 'r*'); set(hp9, 'erasemode', 'none') % Do not give initial line plot. % plot3 (x/xsc, y/xsc, (depth+50)/zsc, 'm-'); view(-70, 60); xlabel(' x/(g\alpha/f^2)') ylabel(' y/(g\alpha/f^2)') zlabel(' z/(g\alpha^2/f^2)') axis([0 60 0 4 -4 0]) title('dense particle on a slope w/ rotation and friction') tip = num2str(t(j)/(2*pi/f)); tipstr = [tip, ' IP']; hth = text(50., 3.5, 0., tipstr); drawnow end if m > 1 hold on % This updated the point but caused refresh blinking. % set (hp9, 'XData', x(j)/xsc, 'YData', y(j)/xsc, ... % 'ZData', (depth(j)+50)/zsc ); plot3 (x(jl:j)/xsc, y(jl:j)/xsc, (depth(jl:j)+10)/zsc, 'r-', ... 'LineWidth', 1.6); if mod(m,4) == 1 tip = num2str(t(j)/(2*pi/f)) % set(htext, 'String', 'tip'; tipstr = [tip, ' IP']; set(hth, 'String', tipstr) end drawnow end jl = j; end % loop on j for plotting each time step figure(2) clf reset subplot(2,1,1) plot(t/tsc,u/usc,':',t/tsc,v/usc) grid xlabel('time/(2\pi/f)') ylabel('U/(g\alpha/f)') title('along and across slope current components') subplot(2,1,2) plot(x/xsc, y/xsc) grid axis('square') xlabel('along slope distance, x/(g\alpha/f^2)') ylabel('across slope, y/(g\alpha/f^2)') axis('equal') figure(3) clf reset plot(t/tsc, ke/esc, t/tsc, pe/esc, t/tsc, fricwork/esc, t/tsc, esum/esc) xlabel('time/(2\pi/f)') ylabel('energy/h(g\alpha/f)^2') legend('kinetic E', 'potential E', 'frictional W', 'E + W') title('Work and energy exchanges') format compact % Check the transport and a few other things. % Assumes that the external force is in the north direction only. eky = -Fx/(rho*f) % Ekman numbers ekx = Fy/(rho*f) ekman = k/f % Ekman number alpha = atan(ekman) % ufric = Fy/( rho*(f*cos(alpha) + k*sin(alpha)) ) urat = ufric/ekx vfric = ufric*sin(alpha) % the end used to stop the script with the quit button end % while close all