%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SLAM: Channel Version.
%
% Single Level Atmospheric Model
% based on the Shallow Water Equations.
%
% u is the zonal velocity component
% v is the meridional velocity component
% z is the geopotential height.
%
% Shallow Water Equations:
% du/dt - fv + g dz/dx = 0
% dv/dt + fu + g dz/dy = 0
% dz/dt + z(du/dx+dv/dy) = 0
%
% A-grid: (u,v,z) given at all points.
% (This causes noise, which must be controlled)
%
% Timestep: Leapfrog (first step forward)
%
% Initial Conditions:
% z0 = zRef + zBar + zWave
% zRef = 5500 metres
% zBar for mean zonal flow
% zWave for perturbation
%
% Initial winds (u0,v0) are geostrophic.
%
% Boundary Conditions: Periodic in x
% Walls at N and S
% (v = 0, dzdx = 0)
%
% Damping:
% Second-order diffusion, equivalent to
% additional terms in the equations, e.g.
% du/dt - fv + g dz/dx = nu del^2 u
%
% Fourier Filtering:
% Even with second-order diffusion, there is
% a severe problem with two-grid-interval noise.
% This is rectified by Fourier-clipping the
% smallest scales. The extent of the clipping
% is determined by the parameter Fband.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear %%% Clear the memory
clf %%% Clear the plot window
% Define basic parameters.
a = (4e7)/(2*pi); % Radius of the Earth
Omega = 2*pi / ( 60*60*24); % Rotation rate of the Earth
g = 9.8066; % Acceleration of gravity
nu = 2.5*10^6; % Damping coefficient
% nu = 1.25*10^6; % Damping coefficient
nu = 0.00; % Damping coefficient
% nu = 5.0*10^7; %%%%%%%%% BIG !!!
Filter = 1; % (1) Fourier Filtering (0) No Filtering.
Fband = 1; % Parameter for Filtering band-width.
% Parameters specifying the domain (4:1 aspect ratio).
Lx = 32*10^6; % Channel length in x-direction
Ly = 8*10^6; % Channel width in y-direction
Nx = 81; % Number of points in x-direction
Ny = 21; % Number of points in y-direction
Nx = 41; Ny = 11;
x = linspace(0,Lx,Nx); % Independent variable x
y = linspace(0,Ly,Ny); % Independent variable y
[X, Y] = meshgrid(x,y); % 2D arrays, mainly for plotting
X = X'; Y = Y'; % MatLab is strange!
dx = Lx/(Nx-1); % Grid step in x-direction
dy = Ly/(Ny-1); % Grid step in y-direction
% Set the forcast length and time step
Tdays = 2; % Forecast length in days
Thours = 24*Tdays; % Forecast length in hours
Tsecs = Thours*3600; % Forecast length in seconds
dt = 600; % Timestep in seconds
dt = 1800; % Timestep in seconds
% dt = 2400; % Timestep in seconds
Nsteps = Tsecs / dt; % Total number of timesteps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Coriolis parameter (beta-plane version).
f0 = 2*Omega*sin(pi/4); % Mean value of f
beta = 2*Omega*cos(pi/4)/a; % Constant beta-parameter
y0 = Ly/2; % Central value of y
f = f0 + beta*(Y-y0); % Coriolis parameter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Define initial conditions
mx = 1; % Wave number in x-direction
my = 1; % Wave number in y-direction
zRef = 5500; % Constant reference geopotential height
% Initial geopotential height: cos-squared in y-direction.
y0 = Ly/2;
zAmp = 500;
zWave = zAmp*sin(2*pi*mx*X/Lx).*cos(my*pi*(Y-y0)/Ly);
% cos-squared in y-direction.
% zWave = zAmp*sin(2*pi*mx*X/Lx).*(cos(my*pi*(Y-y0)/Ly)).^2;
% Add a mean flow (if uBar non-zero).
uBar = 100; uBar = 0;
% Mean height from geopotential balance with uBar.
zBar = -(uBar/g)*(f0*(Y-y0)+0.5*beta*(Y-y0).^2);
% Construct the initial geopotential height
z0 = zRef + zWave + zBar;
% Store the geopotential height.
z = z0;
% Plot initial geopotential height z
SURF = 1;
if (SURF==1 )
surf(X,Y,z); shading('interp');
axis([0,Lx,0,Ly]); % ,5000,6000]);
view(0,90);
xlabel('X'); ylabel('Y'); zlabel('Z')
title('Initial Geopotential Height')
drawnow
else
lines = 4800:40:6200;
contourf(X,Y,z,lines)
drawnow
end
pause
% Calculate the x-gradient of z
for jj=1:Ny
for ii=2:Nx-1
dzdx(ii,jj) = (z(ii+1,jj)-z(ii-1,jj))/(2*dx);
end
dzdx( 1,jj) = (z(2,jj)-z(Nx-1,jj))/(2*dx);
dzdx(Nx,jj) = dzdx(1,jj);
end
% Set N and S boundary x-gradients to zero explicitly.
dzdx(:,1) = 0; dzdx(:,Ny) = 0;
% Calculate the y-gradient of z
for ii=1:Nx
for jj=2:Ny-1
dzdy(ii,jj) = (z(ii,jj+1)-z(ii,jj-1))/(2*dy);
end
dzdy(ii,1) = (z(ii,2) -z(ii,1) )/(dy);
dzdy(ii,Ny) = (z(ii,Ny)-z(ii,Ny-1))/(dy);
end
% Define geostrophic initial winds.
u0 = -g*dzdy./f;
v0 = +g*dzdx./f;
% Store the initial winds
u = u0;
v = v0;
% Modify ICs to accentuate Gravity Waves.
GW = 0;
if(GW==0) end % No action.
if(GW==1) u = 0*u; v = 0*v; end % Zero winds.
if(GW==2) u = -u; v = -v; end % Reverse winds.
if(GW==3) z = zRef+0*zWave; end % Winds only.
% Plot the wind field as arrows.
quiver(u',v')
pause
% Store the fields for start-up of leapfrog scheme.
uOld = u;
vOld = v;
zOld = z;
% Compute the Rossby wave-speed
kx = (2*pi)/(Lx/mx); ky = (2*pi)/(Ly/my);
cR = uBar - beta/(kx^2+ky^2);
Period = Lx/(mx*cR)/(24*3600);
% Define the sizes of the second derivative fields.
d2zdx2 = zeros(Nx,Ny);
d2udx2 = zeros(Nx,Ny);
d2vdx2 = zeros(Nx,Ny);
d2zdy2 = zeros(Nx,Ny);
d2udy2 = zeros(Nx,Ny);
d2vdy2 = zeros(Nx,Ny);
% Define the sizes of the Laplacians
lapz = zeros(Nx,Ny);
lapu = zeros(Nx,Ny);
lapv = zeros(Nx,Ny);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main time-stepping loop.
for nn=1:Nsteps
% Calculate the x-gradients
for jj=1:Ny
for ii=2:Nx-1
dzdx(ii,jj) = (z(ii+1,jj)-z(ii-1,jj))/(2*dx);
dudx(ii,jj) = (u(ii+1,jj)-u(ii-1,jj))/(2*dx);
dvdx(ii,jj) = (v(ii+1,jj)-v(ii-1,jj))/(2*dx);
end
dzdx( 1,jj) = (z(2,jj)-z(Nx-1,jj))/(2*dx);
dudx( 1,jj) = (u(2,jj)-u(Nx-1,jj))/(2*dx);
dvdx( 1,jj) = (v(2,jj)-v(Nx-1,jj))/(2*dx);
dzdx(Nx,jj) = dzdx(1,jj);
dudx(Nx,jj) = dudx(1,jj);
dvdx(Nx,jj) = dvdx(1,jj);
end
% Set N and S boundary values to zero explicitly.
dzdx(:,1) = 0; dzdx(:,Ny) = 0;
dvdx(:,1) = 0; dvdx(:,Ny) = 0;
% Calculate the y-gradients
for ii=1:Nx
for jj=2:Ny-1
dzdy(ii,jj) = (z(ii,jj+1)-z(ii,jj-1))/(2*dy);
dudy(ii,jj) = (u(ii,jj+1)-u(ii,jj-1))/(2*dy);
dvdy(ii,jj) = (v(ii,jj+1)-v(ii,jj-1))/(2*dy);
end
dzdy(ii, 1) = ( z(ii, 2) - z0(ii, 1))/(dy);
dudy(ii, 1) = ( u(ii, 2) - u(ii, 1))/(dy);
dvdy(ii, 1) = ( v(ii, 2) - 0 )/(dy);
dzdy(ii,Ny) = (z0(ii,Ny) - z(ii,Ny-1))/(dy);
dudy(ii,Ny) = ( u(ii,Ny) - u(ii,Ny-1))/(dy);
dvdy(ii,Ny) = ( 0 - v(ii,Ny-1))/(dy);
end
% Compute and plot the divergence.
div = dudx + dvdy;
% fprintf('%g\n',max(max(abs(div))));
% if(nn==1)
% contourf(X,Y,div); title('Divergence'); pause
% end
%% Calculate the Right-hand sides of the equations.
Ru = -(u.*dudx + v.*dudy - f.*v + g*dzdx);
Rv = -(u.*dvdx + v.*dvdy + f.*u + g*dzdy);
Rz = -(u.*dzdx + v.*dzdy + z.*div );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Add the damping term to the Right-hand sides
Ru = Ru + nu*lapu;
Rv = Rv + nu*lapv;
Rz = Rz + nu*lapz;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP FORWARD IN TIME
% Set the time step
tstep = 2*dt;
if(nn==1)
tstep = dt; % Forward step to begin.
end
% Step the equations forward:
% Compute the values at (nn+1)*dt.
uNew = uOld + tstep*Ru;
vNew = vOld + tstep*Rv;
zNew = zOld + tstep*Rz;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute damping using the fields at (n-1)dt.
%
% Calculate the Laplacians of z, u and v
for ii=2:Nx-1
for jj=2:Ny-1
d2zdx2(ii,jj) = (z(ii+1,jj)+z(ii-1,jj)-2*z(ii,jj))/(dx^2);
d2udx2(ii,jj) = (u(ii+1,jj)+u(ii-1,jj)-2*u(ii,jj))/(dx^2);
d2vdx2(ii,jj) = (v(ii+1,jj)+v(ii-1,jj)-2*v(ii,jj))/(dx^2);
d2zdy2(ii,jj) = (z(ii,jj+1)+z(ii,jj-1)-2*z(ii,jj))/(dy^2);
d2udy2(ii,jj) = (u(ii,jj+1)+u(ii,jj-1)-2*u(ii,jj))/(dy^2);
d2vdy2(ii,jj) = (v(ii,jj+1)+v(ii,jj-1)-2*v(ii,jj))/(dy^2);
end
end
for jj=2:Ny-1
d2zdx2(1,jj) = (z(2,jj)+z(Nx-1,jj)-2*z(1,jj))/(dx^2);
d2udx2(1,jj) = (u(2,jj)+u(Nx-1,jj)-2*u(1,jj))/(dx^2);
d2vdx2(1,jj) = (v(2,jj)+v(Nx-1,jj)-2*v(1,jj))/(dx^2);
d2zdy2(1,jj) = (z(1,jj+1)+z(1,jj-1)-2*z(1,jj))/(dy^2);
d2udy2(1,jj) = (u(1,jj+1)+u(1,jj-1)-2*u(1,jj))/(dy^2);
d2vdy2(1,jj) = (v(1,jj+1)+v(1,jj-1)-2*v(1,jj))/(dy^2);
d2zdx2(Nx,jj) = d2zdx2(1,jj); d2zdy2(Nx,jj) = d2zdy2(1,jj);
d2udx2(Nx,jj) = d2udx2(1,jj); d2udy2(Nx,jj) = d2udy2(1,jj);
d2vdx2(Nx,jj) = d2vdx2(1,jj); d2vdy2(Nx,jj) = d2vdy2(1,jj);
end
% Now get the Laplacians (for use next time round)
lapz = d2zdx2 + d2zdy2;
lapu = d2udx2 + d2udy2;
lapv = d2vdx2 + d2vdy2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Move the nn-fields to (nn-1)-slots
uOld = u;
vOld = v;
zOld = z;
% Move the (nn+1)-fields to (nn)-slots
u = uNew;
v = vNew;
z = zNew;
% Store some central values
uMid(nn) = uNew((Nx+1)/2,(Ny+1)/2);
vMid(nn) = vNew((Nx+1)/2,(Ny+1)/2);
zMid(nn) = zNew((Nx+1)/2,(Ny+1)/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Force geostrophic u at N and S boundaries.
% u(:, 1) = -( (Fi(:, 2)-Fi(:, 1))/dy )./f(:, 1);
% u(:,Ny) = -( (Fi(:,Ny)-Fi(:,Ny-1))/dy )./f(:,Ny);
% Force zero divergence at the N and S boundaries.
% First, apply consistency condition, that the integral
% of v along the second row up or down vanishes.
vSodd = mean(v(1:2:Nx,2));
vSevn = mean(v(2:2:Nx,2));
v(1:2:Nx,2) = v(1:2:Nx,2) - vSodd;
v(2:2:Nx,2) = v(2:2:Nx,2) - vSevn;
vNodd = mean(v(1:2:Nx,Ny-1));
vNevn = mean(v(2:2:Nx,Ny-1));
v(1:2:Nx,Ny-1) = v(1:2:Nx,Ny-1) - vNodd;
v(2:2:Nx,Ny-1) = v(2:2:Nx,Ny-1) - vNevn;
% Now integrate along boundaries.
% u(2, 1) = u(Nx-1, 1) - (2*dx/dy)*v(1, 2);
% u(2,Ny) = u(Nx-1,Ny) + (2*dx/dy)*v(1,Ny-1);
for ii=2:Nx-1
u(ii+1, 1) = u(ii-1, 1) - (2*dx/dy)*v(ii, 2);
u(ii+1,Ny) = u(ii-1,Ny) + (2*dx/dy)*v(ii,Ny-1);
end
% u(1, 1) = u(Nx, 1);
% u(1,Ny) = u(Nx,Ny);
% Adjust u to maintain average value on boundaries.
uS = mean(u0(:, 1));
u(:, 1) = u(:, 1)+uS;
uN = mean(u0(:,Ny));
u(:,Ny) = u(:,Ny)+uN;
% Re-set S and N boundary values for z and v.
z(:,1) = z0(:,1); z(:,Ny) = z0(:,Ny);
v(:,1) = 0 ; v(:,Ny) = 0 ;
% Re-set S and N boundary values for u.
u(:,1) = u0(:,1) ; u(:,Ny) = u0(:,Ny) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the z field each timestep.
SURF = 1 ;
if (SURF==1 )
surf(X,Y,z); shading('interp');
axis([0,Lx,0,Ly]); % ,5000,6000]);
view(0,90);
xlabel('X'); ylabel('Y'); zlabel('Z');
title('Geopotential Height')
colorbar
drawnow
else
contourf(X,Y,z,lines)
drawnow
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fourier Filtering in space.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(Filter==1)
% compute the s-to-N trend and remove it.
for nx=1:Nx
zTrend(nx,:) = z(nx,1) + (y/Ly).*(z(nx,Ny)-z(nx,1));
uTrend(nx,:) = u(nx,1) + (y/Ly).*(u(nx,Ny)-u(nx,1));
vTrend(nx,:) = v(nx,1) + (y/Ly).*(v(nx,Ny)-v(nx,1));
end
zDash = z - zTrend;
uDash = u - uTrend;
vDash = v - vTrend;
% Define filtering masks in x and y directions.
NxHalf = (Nx+1)/2;
NyHalf = (Ny+1)/2;
xMask = ones(Nx,Ny);
yMask = ones(Nx,Ny);
xMask(NxHalf-Fband:NxHalf+Fband,:) = 0;
yMask(:,NyHalf-Fband:NyHalf+Fband) = 0;
Mask = xMask.*yMask;
zHat = fft2(zDash); zHatFilt = zHat.*Mask;
uHat = fft2(uDash); uHatFilt = uHat.*Mask;
vHat = fft2(vDash); vHatFilt = vHat.*Mask;
zDashFilt = real(ifft2(zHatFilt));
uDashFilt = real(ifft2(uHatFilt));
vDashFilt = real(ifft2(vHatFilt));
for nx=1:Nx
z(nx,:) = zDashFilt(nx,:)+zTrend(nx,:);
u(nx,:) = uDashFilt(nx,:)+uTrend(nx,:);
v(nx,:) = vDashFilt(nx,:)+vTrend(nx,:);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % End of time-stepping loop.
pause
time = 1:Nsteps;
plot(time,zMid)
title('Geopotential Height at mid-point')
% subplot(1,3,1); plot(time,zMid); legend('z')
% subplot(1,3,2); plot(time,uMid); legend('u')
% subplot(1,3,3); plot(time,vMid); legend('v')
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%