Global Index (all files) (short | long) | Local Index (files in subdir) (short | long)
Mout = sd(lat)
function sdriver(c,dt,mx,maxtime,lat,iplot,distribution) suggested values
| This function calls | |
|---|---|
function Mout = sd(lat)
c=40.;dt=1500;iplot=2;maxtime=76;mx=64;distribution=2
% c=40.;dt=1500;iplot=4;lat=0;maxtime=62;mx=64;distribution=2
% c=40.;dt=1500;iplot=4;lat=45;maxtime=62;mx=64;distribution=1
% lx=2.*pi*0.5e6
% sdriver(c,dt,mx,maxtime,lat,iplot,distribution)
%
% calculates the unforced, linear shallow water equation
% in two dimensions with rotation; n
% any text after the % is treated as a comment
%
%things you may want to change
%
% 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 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;
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
%new initial conditions
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);
%
% this one sets up the gausian phi
%
if distribution ==1;
[px,py] = meshgrid(-xmx2:1:xmx2,-ymy2:1:ymy2);
% [px,py] = meshdom(-xmx2-6:1:xmx2-6,-ymy2-6:1:ymy2-6);
phi = exp(-(.5*px).^2 -(.5*py).^2); px=[];py=[];
% phi = exp(-(1.*px).^2 -(1.*py).^2); px=[];py=[];
phi(:,[mx+1])=[];phi([my+1],:)=[];phio=phi;
% array=[ 'max' max(max(phi))]; mesh(phi);axisp=axis;xlabel(array);
%
elseif distribution == 2;
[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 the step phi intial condtion
%
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
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.;
nframes=floor(maxtime/iplot);
M= moviein(nframes);
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(2); 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])
M(:,jj/iplot) = getframe;
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.5 1])
drawnow
% 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
iiplot=iiplot+iplot;
else
end
% M(:,jj) = getframe(1);
% store the time averaged field for later viewing (start after t=1/f)
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