Team 2 Solution
From Interactive System for Ice sheet Simulation
- 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()