Team 2 Solution

From Interactive System for Ice sheet Simulation
Revision as of 15:08, 5 August 2009 by Adamc (Talk | contribs)

Jump to: navigation, search
	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