# Documentation of sd2

Global Index (all files) (short | long) | Local Index (files in subdir) (short | long)

## Function Synopsis

`Mout = sd2(lat)`

## Help text

```function sdriver(c,dt,mx,maxtime,lat,iplot,distribution)
```

## Cross-Reference Information

This function calls

## Listing of function sd2

```
function Mout = sd2(lat)

%  Mout = sd2(lat);
%
%  lat = latitude at which f is evaluated
%
%  This is a shallow water model of an f-plane in which the value
%  of f (among other things) can be specified.
%
% calculates the unforced, linear shallow water equation
% in two dimensions with rotation; n

% Input parameters:

c=40.;dt=1500;iplot=2;maxtime=76;mx=64;distribution=1;

% c       =  the phase speed;
% dt      =  time incr. for integration
% iplot   =  the number of time steps before plotting
% lat     =  latitude in degrees; lx=the dimensional lenght in x
% maxtime =  the maximum number of time steps
% mx, my  =  the maximum zonal and meridional wave
%            numbers allowed (they must be even and
%            should be in powers of two);
% inc     =  the number of centered forward time
%            steps before a forward time step
% if distribution is 1 or 2, then top hat; else step function

hold off; clg;

%to see short, long inertial gravity and inertial waves use the following
%  statement by removing the % character one the next line and
%  and inserting a % character in the front of the line beginning with
%    "c=40 ...." x lines below.

%c=160.;dt=1500;iplot=4;

%to see  inertial and long inertial gravity waves
%    only - short inertial are aliased out - use the followin line for
%  c, dt  etc., and comment out the line beginning "c=160...." above.

%c=40.;dt=3000;iplot=2;lat=45;maxtime=50;mx=32;

lx=2.*pi*0.20e6;
my=mx;
ly=lx;
beta=0.;

% below this line are found
% things you should never need to change.
%
hold off;clg;

%  Initialize integration variables

dtt=dt;cc=c*c;intime=1.5*24*3600;iiplot=iplot;
daysec=24*3600;inc=25;endtime=dt*maxtime/daysec
f=sin(lat*pi/180.)*1.4544e-4;

array=[c dt lx ];
format short e
disp('   c(m/s)      dt(secs)     domain size (m)  ');
disp(array);

array=[ lat beta mx iplot];
format
disp('latitude beta   m  plotting');
disp(array);

dxdy=sqrt((4*pi*lx/mx).^2+(4*pi*ly/my).^2);
format
disp('Ro/dx');disp(c/f/dxdy);

cfl=dt*c/dxdy;disp('cfl=');disp(cfl);pause(3);nstep=intime/dt;

%mx is the number of grid spaces in x (the number of modes is half that)
xmx=1./mx; xmax=1.-xmx/2.
xmx2=mx/2;xmx2m=xmx2-1;
x=[0:xmx:xmax]; x=x';

% now the y direction
ymy=1./my; ymax=1.-ymy/2.;
ymy2=my/2;ymy2m=ymy2-1;
y=[0:ymy:ymax];y=y';

dx=1./(mx-1);dy=1./(my-1);

%  Initial conditions:

if distribution ==1; % This sets up a sharp gaussian phi
[px,py] = meshgrid(-xmx2:1:xmx2,-ymy2:1:ymy2);
phi = exp(-(.5*px).^2 -(.5*py).^2); 		px=[];py=[];
phi(:,[mx+1])=[];phi([my+1],:)=[];phio=phi;
elseif distribution == 2; % This sets up a shallow gaussian phi
[px,py] = meshgrid(-xmx2:1:xmx2,-ymy2:1:ymy2);
phi = exp(-(0.25*px).^2 -(0.25*py).^2); 		px=[];py=[];
phi(:,[mx+1])=[];phi([my+1],:)=[];phio=phi;
else  %  This one sets up a step function
phi=[0:1:mx];phi=ones(size(phi));phi=[phi(1:xmx2m),phi(xmx2:mx)*0];
phi=(phi')*ones(size(phi));size(phi);              px=[];py=[];
phi(:,[mx+1])=[];phi([my+1],:)=[];phio=phi;
end

%  Draw initial conditions
u=0.*ones(mx,my);v=u;
mesh(phi);axisp=axis;phio=phi;
%    drawnow

%take out mean

mxy=mx*my;
mean=sum(sum(phi)/mxy);
phi=phi-mean;phisum=phi;
mean=sum(sum(u)/mxy);
u=u-mean;usum=u;
mean=sum(sum(v)/mxy);
v=v-mean;vsum=v;
uo=u; phio=phi;vo=v;
%
% set wave numbers
%
m=-i*[0:1:mx]/lx;
l=-i*[0:1:mx]/ly;
pemax=sum(sum(phi.^2))./cc;

% first need to truncate to get rid of half the modes

m=[m(1:xmx2m) m(xmx2:mx)*0];
m=m';
%plot(x,imag(m));pause(3);
l=[l(1:ymy2m) l(ymy2:my)*0];
l=l'; kke=0.; time=0.;ppe=1.;
%m=m*ones(length(m')); l=l*ones(length(l'));
m=m*ones(1,length(m')); l=l*ones(1,length(l'));

% inc is then number of time steps before forward time
% differencing (for stability reasons).

icrit=inc;

% start loop
%
%first time step special
%
dtt=dt/2.;

% Initialize movie out

nframes=floor(maxtime/iplot);
M= moviein(nframes);

%  Start integrations

for jj=1:maxtime
phif=fft2(phi); phiof=fft2(phio);
uf=fft2(u);uof=fft2(uo);vof=fft2(vo);vf=fft2(v);
% now time step
unf=2.*dtt*(m'.*phif+f*vf)+uof;
vnf=2.*dtt*(l.*phif-f*uf)+vof;
phinf=2*dtt*cc*(m'.*uf +l.*vf) + phiof;
%phinf=2*dtt*cc*(l.*vf) + phiof;
phin=real(ifft2(phinf));
un=real(ifft2(unf)); vn=real(ifft2(vnf));

% now switch registers
if jj ==  icrit,
uo=un;u=un;phio=phin;phi=phin;vo=vn;v=vn;
dtt=dt/2.;
icrit=inc+icrit;
else
uo=u;u=un;phio=phi;phi=phin;vo=v;v=vn;
dtt=dt;
end

% do we plot now ?
if jj == iiplot,
% plot phi in 3D
pphi=phi;
clg;
hold off;
subplot(221);
%     figure(1); subplot(1,1,1);

mesh(pphi);axisp=axis;
%array=[ 'max' max(max(phi))]; axis(axisp);mesh(phi);xlabel(array);
%  disp('max =');disp(max(max(pphi)));
hold off
% plot ke, pe
time=[time jj*dt];
ke=sum(sum(u.^2+v.^2))/pemax;kke=[kke ke];
pe=sum(sum(phi.^2))/(cc*pemax);ppe=[ppe pe];
array=[jj jj*dt/daysec ke pe pe+ke];
disp('      jj     time         ke          pe          tote');
disp(array);
axis([0 64 0 64 -0.5 1])
caxis([-0.25 0.25])

%     figure(1);

subplot(222)
axis([0 maxtime*dt/daysec  0 1])
hold on
plot(time/daysec,kke,'-',time/daysec,ppe,'--');
hold off
xlabel('time(days)');ylabel('Energy(m/s)**2');

% plot velocity and geopotential - overlay
%    wide screen - distorted
%     subplot(223);
%    screen not distorted
subplot(212);
vphi4(pphi,u,v);
caxis([-0.25 0.25])
drawnow
% plot geostrophic winds

%     if lat >  0
%     subplot(224);
%     vphi4(pphi,ug,vg);
%     else
%     end

M(:,jj/iplot) = getframe(gcf);
iiplot=iiplot+iplot;
else
end

if jj > jjcrit,
phisum=phisum+phi; vsum=vsum+v;usum=usum+u;
else
end

% go to top and integrate
end

if 0
pause
if f~=0,
phisum=phisum/(maxtime-1/(dt*f));
usum=usum/(maxtime-1/(dt*f));
vsum=vsum/(maxtime-1/(dt*f));
kesum=sum(sum(usum.^2+vsum.^2))
pesum=sum(sum(phisum.^2))
%
%  plot time mean fields; averaged from
%
figure(2)
subplot(111)
vphi4(phisum,usum,vsum)
else
end
end

if nargout == 1
Mout = M;
end
```