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;

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(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);
%     [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