Global Index (all files) (short | long) | Local Index (files in subdir) (short | long)
Mout = sd2(lat)
function sdriver(c,dt,mx,maxtime,lat,iplot,distribution)
| This function calls | |
|---|---|
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=100;mx=32;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;
disp('rossby radius=');disp(c/f);jjcrit=2*pi/(f*dt);
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(211);
% 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 mx 0 mx -0.5 1])
caxis([-0.25 0.25])
view([-37.5 50]);
% 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
axis square
% plot geostrophic winds
% if lat > 0
% subplot(224);
% [px,py]=gradient(pphi,dx,dy);ug=-py/f; vg=px./f; py=[];px=[];
% 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