Team 2 Solution

• FORTRAN mountain glacier

	program Kees1
 
	implicit none
	real, parameter :: dt = 1/12.0
	real, parameter :: dx = 1
	real, parameter :: tend = 5000
	real, parameter :: xend = 30
	real, parameter :: n = 3
	real, parameter :: rho = 920
	real, parameter :: g = 9.2
	real, parameter :: A = 1.0e-16
	real, parameter :: Hunnecessary = 6378.1 * sqrt(9.8/9.2-1) !height where g =9.2m/s^2
	real, parameter :: C = 2*A*(rho*g)**n/(n+2)*(1.e3)**n
	integer :: nt = int(tend/dt)
	integer :: nx = int(xend/dx)
	real, dimension(int(tend/dt)+1,int(xend/dx)+1)  :: S
	integer :: t,i
	real, dimension(int(xend/dx)+1)::xarr=(/((i-1)*dx,i=1,&
		int(xend/dx)+1)/), MB, H, B
	real, dimension(int(xend/dx)) :: Diff
 
	B = 1-0.0001*xarr		
	MB = 0.004 - xarr * 0.0002
 
	S(1,:) = B
 
	do t=2, nt+1
                S(t,1)=0
                H = S(t-1,:) - B
                Diff=C*((H(1:nx)+H(2:nx+1))/2)**(n+2)&
                	*abs((S(t-1,2:nx+1)-S(t-1,1:nx))/dx)**(n-1)
                S(t,2:nx)=S(t-1,2:nx)+(dt/dx**2)&
                	*(Diff(2:nx)*(S(t-1,3:nx+1)-S(t-1,2:nx))&
                	- Diff(1:nx-1)*(S(t-1,2:nx)-S(t-1,1:nx-1)))&
                	+MB(2:nx)*dt
                where(S(t,:)<B)
                	S(t,:)=B
                end where
		if (mod(t,100) .eq.0) then
			write (*,*) S(t,:)+Hunnecessary 
		end if
	end do
	end program

• python plotting script

#!/usr/bin/env python
 
# Import only what is needed
from numpy import loadtxt,shape,linspace
from pylab import plot,show,clf,show,ion,ylim
 
# Import data file, called 'data'
d=loadtxt('data')
 
#Determine how much data came in
dims = shape(d)
 
clf()    # Clears the screen
ion()    # Interactive plot mode, critical for animation
 
# x data, note that this must correspond to program's domain
x = linspace(0,30,dims[1])  
 
# Initial plot, very Matlab(ish), note return of plot handle that allows plot to
# be altered elsewhere in code.
ph,=plot(x,d[0,:],'k')     
ph.figure.show()           # matplot lib requires show to be called
ylim(1629.7,1630.5)
 
# Loop to plot each time step
for i in range(1,dims[0]):
    ph.set_ydata(d[i,:])   # Only update y data (faster than replot)
    ph.figure.show()