% geoadj_1d.m; adjustment under gravity, rotation and friction in 1-D % % A 1-D numerical ocean model of a single % shallow layer of rotating, incompressible fluid. % Useful for showing the numerical version of geostrophic % adjustment. The numerical methods used here are adequate % but not great: an A-grid (velocity and height at the % same grid points) and upwind differencing on advection % terms help preserve a shock, but at the expense of high % numerical viscosity. High frequency ripples on % ridge height should not be taken literally. Time-stepping is by the % leap-frog/trapezoidal time method. In all, this may be about the % simplest PE model that works reasonably well for non-linear % cases, i.e., you can recognize the major features of the % solution, including a fairly reasonable shock (or bore) in % the case of low latitude adjustment to a large initial ridge. % % by Jim Price, January, 2001. % Revised in June 2004 to refine the plotting. % In March, 2010 to add linear drag and gui interface. clear all close all format compact set(0,'DefaultLineLineWidth',1.5) set(0,'DefaultTextFontSize',14) set(0,'DefaultAxesLineWidth',1.4) set(0,'DefaultAxesFontSize',14) % text to display: str2(1) = {'geoadj\_1d.m' }; str2(2) = {'Geostrophic adjustment in one dimension '}; str2(3) = {'solved by a primitive equation numerical '}; str2(4) = {'model representing one (shallow) layer of'}; str2(5) = {'fluid having reduced gravity. Five parameters '}; str2(6) = {'may be modified below:'}; str2(7) = {' '}; str2(8) = {'latitude, 0 -- 90 degrees'}; str2(9) = {'layer thickness, 50 -- 5000 meters'}; str2(10)= {'half-width of the ridge, 20 -- 500 kilometers'}; str2(11)= {'height of the ridge, eta_0, 1 -- 500 meters'}; str2(12)= {'e-folding by Rayleigh drag, 1 -- 99999 days'}; str2(13)= {'ndays to integrate, 3 - 20 days. '}; hf3 = figure(9); clf set(hf3,'Position',[50 200 600 700]) set(gca,'Visible','off') text(-.05, 0.75, str2,'FontSize',14,'HorizontalAlignment','left'); latitude = 20.; hlat = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.37 0.40 0.06], ... 'FontSize', 14,'fontweight', 'demi', ... 'Callback', 'eval(get(hlat,''string''))', ... 'String', 'latitude = 20 '); thickness = 500.; hH = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.30 0.40 0.06], ... 'FontSize', 14, 'fontweight', 'demi',... 'Callback', 'eval(get(hH,''string''))', ... 'String', 'thickness = 500 '); width = 200.; hwid = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.23 0.40 0.06], ... 'FontSize', 14, 'fontweight', 'demi',... 'Callback', 'eval(get(hwid,''string''))', ... 'String', 'width = 200 '); eta0 = 20.; heta = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.16 0.40 0.06], ... 'FontSize', 14,'fontweight', 'demi', ... 'Callback', 'eval(get(heta,''string''))', ... 'String', 'etao = 20 '); efold = 999.; hefold = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.09 0.40 0.06], ... 'Callback', 'eval(get(hefold,''string''))', ... 'FontSize', 14, 'fontweight', 'demi',... 'String', 'efold = 999 '); ndays = 6.; hndays = uicontrol('Style', 'Edit', 'Units', 'Normalized', ... 'Position', [0.05 0.02 0.40 0.06], ... 'Callback', 'eval(get(hndays,''string''))', ... 'FontSize', 14,'fontweight', 'demi', ... 'String', 'ndays = 6 '); hgo = uicontrol('Style', 'Pushbutton', 'Units', 'Normalized', ... 'Position', [0.80 0.10 0.1 0.06], ... 'Callback', 'uiresume(9)', ... 'FontSize', 14,'fontweight', 'demi', ... 'String', 'start'); hquit = uicontrol('Style', 'Pushbutton', 'Units', 'Normalized', ... 'Position', [0.80 0.02 0.1 0.06], ... 'Callback', 'lquit = 1; uiresume(9)', ... 'FontSize', 14,'fontweight', 'demi', ... 'String', 'end'); lquit = 0; while (lquit == 0) disp (' ... waiting for control panel action.') uiwait(9) if lquit == 1; break; end; clear h u % define the basic statification H = thickness; g = 9.8; delr = 2/1030; % fractional density change c = sqrt(delr*g*H); % shallow water phase speed disp(' the non-disersive phase speed, m/s is ') disp(c) % ICtype = menu('What sort of IC?', ... % 'initial interface displacement', 'initial current') ICtype = 1; % initial interface displacement % lat = input(' Enter the latitude (degrees, 0 < lat < 90) '); f = 2*7.29e-5*sin(latitude*pi/180); IP = 2*pi/(f + 1.e-10); % an inertial period, seconds % set a flag to indicate whether rotation is 'on' rota = 1; if abs(latitude) <= 0.1 rota = 0; end % calculate the Radius of deformation Rd = 1.e10; if rota == 1 Rd = c/(f + 1.e-10); % Rossby radius of deformation disp(' the Rossby radius of deformation, km ') disp(Rd/1000) end % the velocity scale Ugeo = g*delr*eta0/(f*Rd); Usc = (eta0/H)*c; disp(' the velocity scale, m/s ') disp(Usc) disp(' the Mach number, Usc/c ') disp(Usc/c) % define the grid in x % dx = Rd/10; % grid interval, in units of Rd % L = 20.*Rd; % domain half-width, Rd (the larger, the better) % in general, probably better to do this..... % if rota == 0 dx = 5.e3; L = 1100.e3; Ld = 1000.e3; % end xu = -L:dx:L; % grid points for u,v,h nx = length(xu); n1e = round((L - Ld)/dx); n2e = nx - round((L - Ld)/dx); for k=1:nx Ldindex(k) = 1; if k <= n1e; Ldindex(k) = 0; end if k >= n2e; Ldindex(k) = 0; end end % define the Rayleigh drag from the given efolding time udecay = (1./(efold*8.64e4)); % Rayleigh drag Ekman = udecay/abs(f); disp(' The Ekman number ') disp(Ekman) % use CFL to set the time step CFL = 0.05; % set the CFL number, used to set dt dt = CFL*dx/c; % time step, seconds nstep = round(ndays*8.64e4/dt); % number of time steps to take % set a flag to make a movie (1) or not (0) mov = 1; if mov == 1 M = moviein(round(nstep/30)); end % define horizontal diffusion; helps reduce small scale numerical noise K = 1.e2; % horizontal diffusion; keep as small as possible disp(' the diffusion length scale (ndays), km ') Kd = sqrt(K*ndays*8.64e4)/1000.; disp(Kd) nonlinu = 1.; % set these both to 0 to get a linear model nonlinh = 1.; % having no advection of momentum and thickness [m nu] = size(xu); xh = ones(1,nu-1); xh = cumsum(xh)*dx; nh = nu - 1; % pre-allocate some arrays u = zeros(1,nu); unew = u; uold = u; uold2 = u; ut = u; v = u; vnew = u; vold = u; vold2 = u; vt = u; h = u; ht = u; hnew = u; holda = u; hold2 = u; na = 0; ha = u; ua = u; va = u; % store and plot the data at this interval nplot = 50; % pre-allocate some arrays for storage of diagncostics uc = zeros(1,nstep); hc = uc; vc = uc; he = uc; ue = uc; ve = vc; ts = uc; hf = he; uf = ue; vf = ue; ke = ue; pe = ue; fw = pe; rw = pe; ns9 = floor(nstep/nplot) + 1; % for data that are decimated in time ts9 = zeros(1, ns9); hs = zeros(nu, ns9); % vs = hs; ns = 0; % a counter used to keep track of how many stored % Specify the IC, first on h, the interface displacement above the % background. xl = width*1000.; % convert width to m for m = 1:nh xp = xh(m) - L; if abs(xp) <= xl; h(m) = eta0; % this makes a box-shaped ridge else h(m) = 0; end % if xp < 0; h(m) = delH*(0.5 + 0.5*tanh((xp + xl)/Hscale)); end; % if xp > 0; h(m) = delH*(0.5 + 0.5*tanh(-(xp - xl)/Hscale)); end; end % add an IC other than zero on the current ICtype = 0 if ICtype == 2 for j=1:nh v(j) = (eta0*f/H)*xu(j); if xu(j) <= -width*1000; v(j) = -eta0*f*width*1000/H; end; if xu(j) >= width*1000; v(j) = eta0*f*width*1000/H; end; vic(j) = v(j); end end % now smooth the ridge a little a4 = 0*h; a5 = a4; for m=2:nh a4(m) = 0.33*(h(m-1) + h(m+1) + h(m)); a5(m) = 0.33*(v(m-1) + v(m) + v(m+1)); end a4(1) = h(1); a4(nh+1) = h(nh+1); h = a4; v = a5; % initialize time levels holda = h; hold2 = h; h0 = h; vold = v; vold2 = v; v0 = v; % store the initial h and PV hi = h; PVi = f./(H + h); % % now we are ready to start time stepping % the first step has to be Euler forward du = 0.5*(diff([0 u]) + diff([u 0])); dv = 0.5*(diff([0 v]) + diff([v 0])); dh = 0.5*(diff([0 h]) + diff([h 0])); du(1) = 0.; du(nu) = 0.; dv(1) = 0.; dv(nu) = 0.; dh(1) = 0.; dh(nu) = 0.; ut = f*v - (delr*g*dh/dx + udecay*u + nonlinu*u.*du/dx); vt = -(f*u + udecay*v + nonlinu*u.*dv/dx); ht = -(H*du/dx + nonlinh*h.*du/dx + nonlinh*u.*dh/dx); uold = u; vold = v; holda = h; u = uold + dt*ut; v = vold + dt*vt; h = holda + dt*ht; hf1 = figure(1); clf reset % now for the real integration..... for m = 1:nstep % use a leap-frog method with an occasional % trapezoidal correction to supress time splitting trap = 20; mtype = mod(m,trap); % decide whether to do a trapezoidal step (20 - 100) if mtype >= 1 % come here for a simple leap frog du = 0.5*(diff([0 u]) + diff([u 0])); dv = 0.5*(diff([0 v]) + diff([v 0])); dh = 0.5*(diff([0 h]) + diff([h 0])); du(1) = 0.; du(nu) = 0.; dv(1) = 0.; dv(nu) = 0.; dh(1) = 0.; dh(nu) = 0.; % velocity and derivatives needed for an upwind difference ul = 0.5*(u + abs(u)); ur = 0.5*(u - abs(u)); dhl = diff([0 h]); dhl(1) = 0.; dhl(nu) = 0.; dhr = diff([h 0]); dhr(1) = 0.; dhr(nu) = 0.; dul = diff([0 u]); dul(1) = 0.; dul(nu) = 0.; dur = diff([u 0]); dur(1) = 0.; dur(nu) = 0.; dvl = diff([0 v]); dvl(1) = 0.; dvl(nu) = 0.; dvr = diff([v 0]); dvr(1) = 0.; dvr(nu) = 0.; du2 = [0 diff(u,2) 0]/dx^2; dv2 = [0 diff(v,2) 0]/dx^2; dh2 = [0 diff(h,2) 0]/dx^2; ut = f*v + K*du2 - (delr*g*dh/dx + udecay*u + ... nonlinu*(ur.*dur/dx + ul.*dul/dx)); vt = K*dv2 - (f*u + udecay*v + ... nonlinu*(ur.*dvr/dx + ul.*dvl/dx)); ht = K*dh2 - (H*du/dx + nonlinh*h.*du/dx + ... nonlinh*(ur.*dhr/dx + ul.*dhl/dx)); uold2 = uold; vold2 = vold; hold2 = holda; uold = u; vold = v; holda = h; u = uold2 + 2*dt*ut; v = vold2 + 2*dt*vt; h = hold2 + 2*dt*ht; else % do a leap-frog/trapezoidal time step du = 0.5*(diff([0 u]) + diff([u 0])); dv = 0.5*(diff([0 v]) + diff([v 0])); dh = 0.5*(diff([0 h]) + diff([h 0])); du(1) = 0.; du(nu) = 0.; dv(1) = 0.; dv(nu) = 0.; dh(1) = 0.; dh(nu) = 0.; % velocity and derivatives needed for an upwind difference ul = 0.5*(u + abs(u)); ur = 0.5*(u - abs(u)); dhl = diff([0 h]); dhl(1) = 0.; dhl(nu) = 0.; dhr = diff([h 0]); dhr(1) = 0.; dhr(nu) = 0.; dul = diff([0 u]); dul(1) = 0.; dul(nu) = 0.; dur = diff([u 0]); dur(1) = 0.; dur(nu) = 0.; dvl = diff([0 v]); dvl(1) = 0.; dvl(nu) = 0.; dvr = diff([v 0]); dvr(1) = 0.; dvr(nu) = 0.; du2 = [0 diff(u,2) 0]/dx^2; dv2 = [0 diff(v,2) 0]/dx^2; dh2 = [0 diff(h,2) 0]/dx^2; ut = f*v + K*du2 - (delr*g*dh/dx + udecay*u + ... nonlinu*(ur.*dur/dx + ul.*dul/dx)); vt = K*dv2 - (f*u + udecay*v + ... nonlinu*(ur.*dvr/dx + ul.*dvl/dx)); ht = K*dh2 - (H*du/dx + nonlinh*h.*du/dx + ... nonlinh*(ur.*dhr/dx + ul.*dhl/dx)); u1 = u + 0.5*dt*ut; v1 = v + 0.5*dt*vt; h1 = h + 0.5*dt*ht; du = 0.5*(diff([0 u1]) + diff([u1 0])); dv = 0.5*(diff([0 v1]) + diff([v1 0])); dh = 0.5*(diff([0 h1]) + diff([h1 0])); du(1) = 0.; du(nu) = 0.; dv(1) = 0.; dv(nu) = 0.; dh(1) = 0.; dh(nu) = 0.; % velocity and derivatives needed for an upwind difference ul = 0.5*(u1 + abs(u1)); ur = 0.5*(u1 - abs(u1)); dhl = diff([0 h1]); dhl(1) = 0.; dhl(nu) = 0.; dhr = diff([h1 0]); dhr(1) = 0.; dhr(nu) = 0.; dul = diff([0 u1]); dul(1) = 0.; dul(nu) = 0.; dur = diff([u1 0]); dur(1) = 0.; dur(nu) = 0.; dvl = diff([0 v1]); dvl(1) = 0.; dvl(nu) = 0.; dvr = diff([v1 0]); dvr(1) = 0.; dvr(nu) = 0.; du2 = [0 diff(u1,2) 0]/dx^2; dv2 = [0 diff(v1,2) 0]/dx^2; dh2 = [0 diff(h1,2) 0]/dx^2; ut = f*v + K*du2 - (delr*g*dh/dx + udecay*u + ... nonlinu*(ur.*dur/dx + ul.*dul/dx)); vt = K*dv2 - (f*u + udecay*v + ... nonlinu*(ur.*dvr/dx + ul.*dvl/dx)); ht = K*dh2 - (H*du/dx + nonlinh*h.*du/dx + ... nonlinh*(ur.*dhr/dx + ul.*dhl/dx)); u = uold + dt*ut; v = vold + dt*vt; h = holda + dt*ht; uold = u; vold = v; holda = h; ut = 0.*ut; vt = 0.*vt; ht = 0.*ht; end % if on mtype (whether leap or trap) % Apply a radiation BC to the grid ends. This is only moderately % succesful because the phase speed changes in time. bcfactor = 1.2; % an ad hoc 'correction' for the radiation BC rbc = bcfactor*c*dt/dx; u(1) = u(1) + rbc*(u(2) - u(1)); v(1) = v(1) + rbc*(v(2) - v(1)); h(1) = h(1) + rbc*(h(2) - h(1)); u(nu) = u(nu) + rbc*(u(nu-1) - u(nu)); v(nu) = v(nu) + rbc*(v(nu-1) - v(nu)); h(nu) = h(nu) + rbc*(h(nu-1) - h(nu)); % store some data for plotting later ts(m) = m*dt/8.64e4; % the time (days) np = length(u); nf = np - round(xl/dx); hf(m) = h(nf); % far field h,u,v uf(m) = u(nf); vf(m) = v(nf); ne = round(np/2 + round(xl/dx)); he(m) = h(ne); % h,u,v at a point near the edge of the front ue(m) = u(ne); ve(m) = v(ne); na = na + 1; % sum for averages ha = ha + h; ua = ua + u; va = va + v; % evaluate the energy balance at each time step over the % interval -Ld to Ld keall = 0.5*(H + h).*(u.^2 + v.^2); peall = 0.5*g*delr*(h.^2); ke(m) = dx*sum(Ldindex.*keall); pe(m) = dx*sum(Ldindex.*peall); fw3 = udecay*(u.^2 + v.^2).*(H + h); fs1 = -dt*dx*sum(Ldindex.*fw3); % the radiation (pressure work) evaluated at +- 1000 km. rad = dt*( u(n1e)*g*delr*h(n1e)*(H + h(n1e)) ... - u(n2e)*g*delr*h(n2e)*(H + h(n2e)) ); % sum the friction and radiation if m==1 fw(m) = fs1; rw(m) = rad; else fw(m) = fw(m-1) + fs1; rw(m) = rw(m-1) + rad; end % other diagnostic data are decimated in time if mod(m,nplot) == 1 % save the entire grid occasionally ns = ns+1; hs(:,ns) = h(:); % vs(:,ns) = sqrt(v(:).^2 + u(:).^2); ts9(ns) = ts(m); % plot the solution, occasionally figure(1) clf set(gcf,'DoubleBuffer','on') set(gcf, 'position', [250 250 800 400]) axis([-1000 1000 -0.5 3.5]) vv = axis; xyscale = (vv(2) - vv(1))/(vv(4) - vv(3)); for q=1:4:nu hold on % xa = xu(q)/1000.; if rota == 1 xa = xu(q)/1000.; else xa = xu(q)/1.e3; end Lplot = 1000.e3; % this is the scale of the plot ya = 2.5; xb = xa + (2*Lplot/(1000*4))*uold(q)/Usc; % xb = xa + 7.5*uold(q)/Usc; % this 5 depends upon the axis set below yb = ya + vold(q)/Usc; if rota == 0 % switch the components if there is no rotation xb = xa + 7.5*vold(q)/Usc; yb = ya + uold(q)/Usc; end plot([xa xb], [ya yb], 'LineWidth', 1.0) end plot(xu/1000., hi/eta0, 'g', 'linewidth', 1.0) % the IC if rota == 0 % title(['gravitational adjustment, latitude = ', num2str(latitude,2)]) else % title(['geostrophic adjustment; latitude = ', num2str(latitude,2), ... % ', Ekman = ', num2str(Ekman,2)]) end % add the Rd scale if abs(latitude) >= 0.1 plot([-500 (-500 + Rd/1000)], [1.5 1.5], 'k') text(-650, 1.5, 'Rd') end hold on if rota == 1 plot(xu/1000, holda/eta0); fill([xu(1)/1000 xu/1000 max(xu)/1000], [-9 holda/eta0 -9], [0.6 0.6 0.8]); else plot(xu/1.e3, holda/eta0, 'EraseMode', 'none'); fill([xu(1)/1.e3 xu/1.e3 max(xu)/1.e3], [-9 holda/eta0 -9], [0.6 0.6 0.8]); end ylabel('\eta/\eta_0 and U/(C \eta_0/H)') if rota == 1 ylabel('\eta/\eta_0 and V/(C \eta_0/H)') end xlabel('X distance, km', 'EraseMode', 'none') timed = (floor(100*m*dt/8.64e4))/100.; if rota == 1 text(250., 1.5, ['time = ', num2str(timed), ' days']) tIP = timed*8.64e4/IP; timeIP = (floor(100*tIP))/100.; text(250., 1.2, [' = ', num2str(timeIP), ' IP']) else text(250., 1.5, ['time = ', num2str(timed), ' days']) end % text(15., 1.3, ['latitude = ', num2str(lat)]) drawnow mov = 0; if mov == 1 M(:,ns) = getframe(gcf, [0 0 800 400]); end end % if on plot this time step ot not end % end of the time stepping loop % compute time averages ha = ha/na; ua = ua/na; va = va/na; % plot some things % plot some current components figure(2) clf reset tip = ts; subplot(2,1,1) plot(ts, uf/Usc, 'b', ts, vf/Usc, '--r') ylabel('U far/(C \eta_0/H)') legend('U','V') title('Currents in the far field (x=900 km) and on the ridge edge') subplot(2,1,2) plot(ts, ue/Usc, 'b', ts, ve/Usc, '--r') xlabel('time, days'); ylabel('Uedge/(C \eta_0/H)') legend('U','V') % 3-d plot of eta figure(3) clf reset set(gcf, 'paperposition', [1 1 6 6]) mesh(ts9, xu/1000., hs/eta0); xlabel('time, days') ylabel('X distance, km') zlabel('\eta/\eta0') axis([0 ndays -Ld/1000 Ld/1000 -0.5 1.5]) % title('geostrophic/gravitational adjustment in 1-d') view([-70 60]) % add the gravity wave speed tx = 0:0.1:5; cx = tx*8.64e4*c; zx = 0.55*ones(1,length(tx)); hold on plot3(tx, cx/1000., zx, 'r', tx, -cx/1000., zx, 'r') if rota == 1 % compute the final PV dv = 0.5*(diff([0 va]) + diff([va 0])); dv(1) = 0.; dv(nu) = 0.; rvort = dv/dx; PVf = (f + rvort)./(H + h); PVflin = (f + 0*rvort)./(H + h); % % figure(4) % clf reset % subplot(2,1,1) % plot(xu/1000., hi/eta0, 'g', xu/1000., h/eta0,'r'); % axis([-Ld/1000. Ld/1000 -1.2 1.2]); % ylabel('h/\eta_0') % legend('initial', 'final') % legend boxoff % title('Initial and final thickness and PV') % % subplot(2,1,2) % plot(xu/1000., PVi,'g', xu/1000., PVf, 'r', xu/1000., PVflin, 'b'); % axis([-L/1000 L/1000 0.95e-7 1.02e-7]); % xlabel('distance, km') % ylabel('Potential Vorticity') % legend('initial', 'final','f/ H + \eta') % legend boxoff % end % average current and rel vorticity figure(5) clf reset subplot(2,1,1) plot(xu/1000., ua/Usc, xu/1000., va/Usc, 'g') ylabel('U, V avg speeds/(C \eta_0/H)') legend('Uavg', 'Vavg') axis([-L/1000 L/1000 -0.5 0.5]) title('Time average currents, etc.') subplot(2,1,2) plot(xu/1000., rvort/f) xlabel('distance, km') ylabel('rel vorticity/f') axis([-L/1000 L/1000 -eta0/H eta0/H]) end % if on rotation or not % energy budget figure(8) clf reset set(gcf, 'paperposition', [0 0 5 5]) Esc = pe(1); plot(ts, ke/Esc, 'r', ts, pe/Esc, 'b', ts, (ke + pe)/Esc, 'g', ... ts, fw/Esc, 'r--', ts, rw/Esc, 'm') hl = legend('KE', 'PE','KE + PE', 'FW', 'PW', ... 'Location', 'SouthWest') legend boxoff set(hl, 'fontsize', 14) xlabel('time, days') ylabel('energy/PE_0') hold on plot([0 max(ts)], [0 0], 'k--') axis([0 max(ts) -1.2 1.2]) L = 1000000; figure(10) clf reset qscale = f/H; set(gcf,'paperposition', [0 0 5 8.5]) subplot(3,1,2) %plot(xu/1000., va/Usc) % , xu/1000., ua/Usc, 'r') ylabel('V/(C \eta_0/H)') %legend('Vavg', 'Uavg') legend boxoff axis([-L/2000 L/2000 -0.5 0.5]) grid on; box on hold on for k = 1:3:length(xu) plot([xu(k)/1000 xu(k)/1000], [0 va(k)/Usc], 'b') end subplot(3,1,3) plot(xu/1000., PVi/qscale,'g', xu/1000., PVf/qscale, 'b', ... xu/1000., PVflin/qscale, 'r'); axis([-L/2000 L/2000 0.95 1.02]); ylabel('q/(f/H)') hl = legend('q_0', 'q_{final}','q_{Sverdrup}') set(hl,'fontsize', 14) legend boxoff xlabel('distance, km') grid on % radius of deformation g = 9.8; dr = 2.; r0 = 1000.; H = 500.; lat = 20.; f = 2.*7.292e-5*sin(lat*pi/180); C = sqrt(g*dr*H/r0); Rd = C/f; Ly = 200.; y = -2.5*Ly:2:2.5*Ly; y = 1000*y; Ly = Ly*1000.; eta0 = 20. for j=1:length(y) if y(j) >= 0 eta(j) = (eta0/2)*(1 - exp((y(j) - Ly)/Rd)); if (y(j) - Ly) >= 0 eta(j) = eta0/2*(exp(-(y(j) - Ly)/Rd) - 1); end end if y(j) <= 0 eta(j) = (eta0/2)*(exp((y(j) + Ly)/Rd) - 1); if (y(j) + Ly) >= 0 eta(j) = eta0/2*(1 - exp(-(y(j) + Ly)/Rd)); end end end eta = eta + eta0/2 subplot(3,1,1) plot(xu/1000, hi/eta0, 'g', xu/1000, h/eta0, 'b', y/1000, eta/eta0, 'r') hl = legend('\eta_0', '\eta_{final}', '\eta_{q conserv}'); set(hl,'fontsize', 14) legend boxoff axis([-L/2000 L/2000 -0.2 1.2]) grid on ylabel('\eta/\eta_0') figure(9) % % if you want to see the movie again % replayn = 0; % replayn = input(' To replay the movie, enter 1, 2 or 3 (0 = no) '); % if replayn >= 1 % figure(1); % set(gcf, 'paperposition', [0 0 850 400]) % set(gca, 'position', [0 0 1 1]) % movie(M,replayn,6) % end end % end of the function % % % imov = input(' Shall we save an mpeg movie file? (1 = yes)'); % if imov == 1 % if exist('geoadjmovie.m') == 1; % delete 'geoadjmovie.mpg'; % end % mpgwrite(M, colormap, 'geoadjmovie.mpg') % end