% dogyre_new.m % a procedure to run a two dimensional gyre initially % with zero concentration, and whose % "southern" boundary (row 1) is set at the beginning % of the run to be 10 units. The domain is 50 nodes by 50 nodes, % with equal spacing of 50 km (50,000 m). The velocity field is % sinsusoidal (see text of problem set) such that it is rotating % counter-clockwise. The maximum velocity is 4 cm/s (.04 m/s). % % The code as written does four simulations, with four different turbulent % diffusivities (0, 500, 1000 and 2000 m^2/s), and plots the evolving % property fields at one year intervals, running for a total of 10 years. % The time step is adjusted for each simulation conservatively for % courant/diffusion stability % There is some stuff in here that is not absolutely necessary, but % makes the code a little more general for possible later use. % % This routine takes less than 7 minutes to run on a 100 MHz 486 % % Modif'd 2004:12:03 D. Glover to get contour to work with zero matrix % Modif'd 2008:09:17 N. Levine to fix boundary conditions % Modif'd 2008:11:17 DMG to clean up excessive tabbing in code, etc. % Modif'd 2008:11:18 DMG a few minor labeling changes to clarify % Modif'd 2008:11:19 DMG to add documentation about the Stommel gyre warning off % to take care of v4 vs. v5 warnings tyr=3600*24*365.25; nx=50; ny=50; % various constants dx=50000; dy=dx; % node spacing X=[0:nx-1]*dx/nx/1000;Y=[0:ny-1]*dy/ny/1000; % distance vectors for plotting (in km) %Coordinates y=meshgrid([0:ny]/(ny)); x=meshgrid([0:nx]/(nx)); y=y'; % Stommel gyre to calculate A so that the average advective velocity is correct % % the magnitude of A controls the average velocity, epsilon controls % western intensification (a value of 0.1 yields Stommel's gyre) A=.02; epsilon=1; lamda1=(-1+sqrt(1+4*(pi^2)*(epsilon^2)))/(2*epsilon); lamda2=(-1-sqrt(1+4*(pi^2)*(epsilon^2)))/(2*epsilon); c1=(1-exp(lamda2))/(exp(lamda2)-exp(lamda1)); c2=-(1+c1); for j=1:nx+1 for i=1:ny+1 psi(i,j)=A.*sin(pi*y(i,j)).*(c1.*exp(lamda1.*x(i,j)) + c2.*exp(lamda2.*x(i,j)) +1); % the stream function (not used here) end end u=zeros(ny+1, nx+1); v=u; for j=1:nx+1 for i=1:ny+1 yy = y(i,j)+1/(2*ny); xx = x(i,j)+1/(2*nx); u(i,j)=pi.*cos(pi.*yy).*A.*(c1.*exp(lamda1.*x(i,j)) + c2.*exp(lamda2.*x(i,j)) +1); v(i,j)=-sin(pi.*y(i,j)).*A.*(lamda1.*c1.*exp(lamda1.*xx)+ lamda2.*c2.*exp(lamda2.*xx)); end end km=[0 500 1000 2000]; % vector of diffusivities to try V=[0:10]; % vector of contour line values vmax=max(max(v)); umax=max(max(u)); % maximum possible velocities colour=['r' 'g' 'b' 'c']; % colour code the curves figure(1); clg; figure(2); clg; % start with a clean slate for np=1:length(km) % now do the simulation for various diffusivities dt=0.5*(dx/(umax+vmax+4*km(np)/dx)); % stable time step nt=10*tyr/dt; % number of steps in 10 years kx=km(np)*ones(ny+1,nx+1); ky=kx; % kx, ky over all domain c=zeros(ny,nx); cold=c; % initialize concentrations to zero % now calculate weight matrices % now calculate weight matrices wxm=zeros(ny,nx); wym=wxm; w0=wxm; w0=1-kx(1:ny,1:nx)*dt/dx/dx -kx(1:ny,2:nx+1)*dt/dx/dx -ky(1:ny,1:nx)*dt/dy/dy- ky(2:ny+1,1:nx)*dt/dy/dy; wxm=kx(1:ny,1:nx)*dt/dx/dx; wxp=kx(1:ny,2:nx+1)*dt/dx/dx; wym=ky(1:ny,1:nx)*dt/dy/dy; wyp=ky(2:ny+1,1:nx)*dt/dy/dy; for i=1:ny for j=1:nx if u(i,j)>=0 & u(i,j+1)>=0 w0(i,j)=w0(i,j)-u(i,j+1)*dt/dx; wxm(i,j)=wxm(i,j)+u(i,j)*dt/dx; elseif u(i,j)<=0 & u(i,j+1)<=0 w0(i,j)=w0(i,j)+u(i,j)*dt/dx; wxp(i,j)=wxp(i,j)-u(i,j+1)*dt/dx; elseif u(i,j)>0 & u(i,j+1)<0 wxm(i,j)=wxm(i,j)+u(i,j)*dt/dx; wxp(i,j)=wxp(i,j)-u(i,j+1)*dt/dx; elseif u(i,j)<0 & u(i,j+1)>0 w0(i,j)=w0(i,j)-u(i,j+1)*dt/dx+u(i,j)*dt/dx; end if v(i,j)>=0 & v(i+1,j)>=0 w0(i,j)=w0(i,j)-v(i+1,j)*dt/dy; wym(i,j)=wym(i,j)+v(i,j)*dt/dy; elseif v(i,j)<=0 & v(i+1,j)<=0 w0(i,j)=w0(i,j)+v(i,j)*dt/dy; wyp(i,j)=wyp(i,j)-v(i+1,j)*dt/dy; elseif v(i,j)>0 & v(i+1,j)<0 wym(i,j)=wym(i,j)+v(i,j)*dt/dy; wyp(i,j)=wyp(i,j)-v(i+1,j)*dt/dy; elseif v(i,j)<0 & v(i+1,j)>0 w0(i,j)=w0(i,j) - v(i+1,j)*dt/dy + v(i,j)*dt/dy; end end end %Boundary conditions wxm(:,1)=0; wxp(:,nx)=0; wym(1,:)=0; wyp(ny,:)=0; w0(:,1)=w0(:,1)+ kx(1:ny,1)*dt/dx/dx; %Western w0(:,nx)=w0(:,nx)+ kx(1:ny,nx)*dt/dx/dx; %Eastern w0(1,:)=w0(1,:)+ ky(1,1:nx)*dt/dy/dy; %South w0(ny,:)=w0(ny,:)+ ky(ny,1:nx)*dt/dy/dy; %North oldyear=0; % our trusty year counter figure(1) % the main screen for t=1:nt % your main time step loop here cold(1,:)=10*ones(1,nx); % set boundary conditions (others are zero) newyear=ceil(t*dt/tyr); % time counter if newyear ~= oldyear % happy new year? subplot(2,2,np) % go to the correct plot if t > 1 % [Ccon,Hcon]=contour(X,Y,c,V); clabel(Ccon,Hcon); title(sprintf('%d %d',newyear-1,km(np)));pause(1); clabel(contour(X,Y,c,V)); title(sprintf('%d %d',newyear-1,km(np)));pause(1); end%if oldyear=newyear; % only do once a year end %The main event c(2:ny-1,2:nx-1) = w0(2:ny-1,2:nx-1).*cold(2:ny-1,2:nx-1) +... wxm(2:ny-1,2:nx-1).*cold(2:ny-1,1:nx-2) +... wxp(2:ny-1,2:nx-1).*cold(2:ny-1,3:nx) +... wym(2:ny-1,2:nx-1).*cold(1:ny-2,2:nx-1) +... wyp(2:ny-1,2:nx-1).*cold(3:ny,2:nx-1); %Western Boundary c(2:ny-1,1) = w0(2:ny-1,1).*cold(2:ny-1,1) +... wxp(2:ny-1,1).*cold(2:ny-1,2) +... wym(2:ny-1,1).*cold(1:ny-2,1) +... wyp(2:ny-1,1).*cold(3:ny, 1); %Eastern Boundary c(2:ny-1,nx) = w0(2:ny-1,nx).*cold(2:ny-1,nx) +... wxm(2:ny-1,nx).*cold(2:ny-1,nx-1) +... wym(2:ny-1,nx).*cold(1:ny-2,nx) +... wyp(2:ny-1,nx).*cold(3:ny, nx); %Southern Boundary c(1,2:nx-1) = w0(1,2:nx-1).*cold(1,2:nx-1) +... wxm(1,2:nx-1).*cold(1,1:nx-2) +... wxp(1,2:nx-1).*cold(1,3:nx) +... wyp(1,2:nx-1).*cold(2,2:nx-1); %Northern Boundary c(ny,2:nx-1) = w0(ny,2:nx-1).*cold(ny,2:nx-1) +... wxm(ny,2:nx-1).*cold(ny,1:nx-2) +... wxp(ny,2:nx-1).*cold(ny,3:nx) +... wym(ny,2:nx-1).*cold(ny-1,2:nx-1); % The 4 corners %southwestern boundary c(1,1) = w0(1,1).*cold(1,1) +... wxp(1,1).*cold(1,2) +... wyp(1,1).*cold(2,1); %southeastern boundary c(1,nx) = w0(1,nx).*cold(1,nx) +... wxm(1,nx).*cold(1,nx-1) +... wyp(1,nx).*cold(2,nx); %northwestern boundary c(ny,1) = w0(ny,1).*cold(ny,1) +... wxp(ny,1).*cold(ny,2) +... wym(ny,1).*cold(ny-1,1); %southeastern boundary c(ny,nx) = w0(ny,nx).*cold(ny,nx) +... wxm(ny,nx).*cold(ny,nx-1) +... wym(ny,nx).*cold(ny-1,nx); cold=c; % update concentration matrix end % now must plot one final time clabel(contour(X,Y,c,V)); title(sprintf('K = %d',km(np)));pause(1); csum=sum(sum(c)); ccent=c(ny/2,nx/2); figure(2) subplot(2,2,1); plot(km(np),csum,[colour(np) '*']); xlabel('Diff'); ylabel('Sum'); % axis([0 2000 0 6000]); hold on axis([0 2000 0 2e4]); hold on subplot(2,2,2); plot(km(np),ccent,[colour(np) '*']); xlabel('Diff');ylabel('Central Value'); % axis([0 2000 0 2.5]); hold on axis([0 2000 0 6]); hold on legend('K_{H}=0','K_{H}=500','K_{H}=1000','K_{H}=2000','location','northwest'); legend('boxoff'); subplot(2,2,3); plot(X,c(ny/2,:),colour(np)); xlabel('Grid pt');ylabel('C along x-transect'); axis([0 50 0 6]); hold on % an "east-west" line across the center subplot(2,2,4); plot(c(:,nx/2),Y,colour(np)); xlabel('C along y-transect');ylabel('Grid pt'); axis([0 10 0 50]); hold on % a "north-south" line up the center end