# Team 2 Solution

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
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)=S(t-1,2)
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