Dynamics II: Sheet 7
Lecture: June 7 (Monday), 14:00 Prof. Dr. Gerrit Lohmann
Tutorial: June 7 (Monday), ca. 15:30 Justus Contzen, Lars Ackermann
Time required for Sheet 7: 9 h
June 7, 14:00: Lecture 7 (online G. Lohmann, 45 min)
7th lecture link to pdf
Content in the script: Shallow water equations and Waves
Topics: Rossby, Gravity, Kelvin waves, Plain waves, Scaling
Animations for Rossby waves (look some of them in advance)
Rossby waves naturally occur in rotating fluids. Within the Earth’s ocean and atmosphere, these planetary waves play a significant role in shaping weather:
Rossby wave animation,
Extremes and Rossby waves,
Jet streams,
experiment in a tank
Atmospheric Dynamics: Mountain Waves (12 min)
some additional background: Rossby wave in a tank, experiment in a tank
June 7, ca. 15:30: Tutorial (online 45 min)
Exercise 7 introduced, questions to the exercise (15 min)
Literature:
Equatorial waves: Theory of Matsuno (the code is below)
shallow2D_rossby.r for your R application
#This is just a definition of a function to plot vectorplots, you do not have to understand it...
par.uin<-function()
# determine scale of inches/userunits in x and y
# from http://tolstoy.newcastle.edu.au/R/help/01c/2714.html
# Brian Ripley Tue 20 Nov 2001 - 20:13:52 EST
{
u <- par("usr")
p <- par("pin")
c(p[1]/(u[2] - u[1]), p[2]/(u[4] - u[3]))
}
quiver<-function(lon,lat,u,v,scale=1,length=0.05,maxspeed=200, ...)
# first stab at matlab's quiver in R
# from http://tolstoy.newcastle.edu.au/R/help/01c/2711.html
# Robin Hankin Tue 20 Nov 2001 - 13:10:28 EST
{
ypos <- lat[col(u)]
xpos <- lon[row(u)]
speed <- sqrt(u*u+v*v)
u <- u*scale/maxspeed
v <- v*scale/maxspeed
matplot(xpos,ypos,type="p",cex=0,xlab="lon",ylab="lat", ...)
arrows(xpos,ypos,xpos+u,ypos+v,length=length*min(par.uin()))
}
#Program starts here
#Shallow water 2D,cyclic boundary conditions + Coriolis term
nn<- 50
ni<- 2*nn+1 #number of gridcells in one direction
nt<-10000 #number of timesteps
#The physical constants
g<-0.1 #low gravity, 0.1 m/s^2
dx<-1e5 #gridcell 10km
dy<-1e5
dt<-500 #timstep 1000 second
H<-1e3 #1km depth
Omega<-1e-4
#define three index vectors.. the middle one , one shifted one cell to the left, and one to the right
#(including the periodic boundary conditions)
ia.0<-1:ni
ia.m1<-c(ni,1:(ni-1))
ia.p1<-c(2:ni,1)
u<-matrix(0,ni,ni) #speed at each point
v<-matrix(0,ni,ni) #speed at each point
h<-matrix(0,ni,ni) #pertubation at each point
f<-matrix(0,ni,ni) #pertubation at each point
lat<-c(-nn:nn)*90/nn
weight<-sin(lat*pi/180)
lon<-c(-nn:nn)*180/nn
f<-rep(weight*2*Omega,each=ni)
dim(f)<-c(ni,ni)
#filled.contour(f)
u.new<-u
h.new<-h
v.new<-v
#Inital condition: One smooth blobs at each side of the "equator"(sin)
idit=nn/5*2
inix=ni-idit-1
iniy=ni-2*idit-1
endx=ni-1
endy=ni-1
endy2=2*idit+1
h[inix:endx,iniy:endy]<-sin(0:20/2*pi/10)%*%t(sin(0:40/2*pi/20))
h[inix:endx,1:endy2]<-sin(0:20/2*pi/10)%*%t(sin(0:40/2*pi/20))
#equator to study the Kelvin wave:
ii=idit+1
iy=nn-10
iy2=nn+10
h[1:ii,iy:iy2]<- -sin(0:20/2*pi/10)%*%t(sin(0:20/2*pi/10))
#1st step euler forward
u.new[ia.0,ia.0]<-u[ia.0,ia.0]-g*dt/2/dx*(h[ia.p1,ia.0]-h[ia.m1,ia.0])
v.new[ia.0,ia.0]<-v[ia.0,ia.0]-g*dt/2/dy*(h[ia.0,ia.p1]-h[ia.0,ia.m1])
h.new[ia.0,ia.0]<-h[ia.0,ia.0]-H*dt/2*((u[ia.p1,ia.0]-u[ia.m1,ia.0])/dx + (v[ia.0,ia.p1]-v[ia.0,ia.m1])/dy) #(du/dx + dv/dy)
#Divide the screen in two parts
#par(mfcol=c(2,1))
#Leapfrog from the third step on
for (n in 3:(nt-1))
{
u.old<-u
v.old<-v
h.old<-h
h<-h.new
u<-u.new
v<-v.new
u.new[ia.0,ia.0]<-u.old[ia.0,ia.0]-g*dt/dx*(h[ia.p1,ia.0]-h[ia.m1,ia.0])+dt*f*v
v.new[ia.0,ia.0]<-v.old[ia.0,ia.0]-g*dt/dy*(h[ia.0,ia.p1]-h[ia.0,ia.m1])-dt*f*u
h.new[ia.0,ia.0]<-h.old[ia.0,ia.0]-H*dt*((u[ia.p1,ia.0]-u[ia.m1,ia.0])/dx + (v[ia.0,ia.p1]-v[ia.0,ia.m1])/dy) #(du/dx + dv/dy)
#plot every 50th image
if ((n %% 50) == 0) {
#quiver(lon,lat,u,v,scale=200,maxspeed=1.5,length=3)
image(lon,lat,h,zlim=c(-1,1),col=rainbow(200)) #this time as color coated
#persp(h/3, theta = 0, phi = 40, scale = FALSE, ltheta = -120, shade = 0.6, border = NA, box = FALSE,zlim=c(-0.3,0.3))
}
}
