Team 2 Solution
From Interactive System for Ice sheet Simulation
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 :: 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,:)+1633.42 !height at which g=9.2 end if end do end program