# Team 5 Solution

Jump to: navigation, search
```! A solutions to Kees'' problem - explicit non-linear diffusion equation
!	for time integration of mountain glacier flow

program main

implicit none

! x: vector along flowline
! b: bed elevation
! S: surface elevation
! H: ice thickness
! M: surface mass balance

real, allocatable :: x(:), b(:), S(:), H(:), M(:), dHdt(:) ! grid points on primary grid
real, allocatable :: Q_mid(:), D_mid(:), dSdx_mid(:)	! grid points for staggered grid
real, parameter :: dx = 1000, dt = 0.01  ! spatial (units: m) and time step (units: yr)
real, parameter :: grid_length = 50000 ! length of domain (m)
real, parameter :: S0 = 1000 ! elevation of first point on grid (m)
real, parameter :: M0 = 4, M1=-2e-4 ! parameters for linear mass balance
real, parameter :: dbdx = -0.1 ! bedslope
real, parameter :: g = 9.8  ! m/s2
real, parameter :: rho = 920 ! kg/m3
real, parameter :: A = 1e-16  ! kpa-3yr-1
real, parameter :: n = 3 ! flow law exponent
real, parameter :: t_int = 1  ! interval to output data
real, parameter :: tol = 1e-6  ! steady-state tolerance
real, parameter :: t_max = 500  ! maximum time
real :: t, t_old	! time variables
real :: C 			! diffusion constant

integer :: Nx, ii, errstat

Nx = floor(grid_length/dx)+1

allocate(x(Nx),stat=errstat)
call checkerr(errstat,"failed to allocate x")

allocate(b(Nx),stat=errstat)
call checkerr(errstat,"failed to allocate b")

allocate(S(Nx),stat=errstat)
call checkerr(errstat,"failed to allocate S")

allocate(H(Nx),stat=errstat)
call checkerr(errstat,"failed to allocate H")

allocate(M(Nx),stat=errstat)
call checkerr(errstat,"failed to allocate M")

allocate(dHdt(Nx),stat=errstat)
call checkerr(errstat,"failed to allocate dHdt")

allocate(Q_mid(Nx-1),stat=errstat)
call checkerr(errstat,"failed to allocate Q_mid")

allocate(D_mid(Nx-1),stat=errstat)
call checkerr(errstat,"failed to allocate D_mid")

allocate(dSdx_mid(Nx-1),stat=errstat)
call checkerr(errstat,"failed to allocate dSdx_mid")

open(unit=1,file='data')

! set up vectors
do ii=1,Nx
x(ii)=(ii-1)*dx
b(ii)=S0 + x(ii) * dbdx
S(ii)=b(ii)  ! start with zero ice thickness
H(ii)=S(ii)-b(ii)
dHdt(ii)=100 ! start with something silly for steady-state test
M(ii)=M0+M1*x(ii)
enddo
write(*,*) M

C = 2 * A/(n+2)*(rho*g)**n

! time loop
t=0
t_old=-t_int
time_loop: do while (maxval(abs(dHdt))>tol .and. t<t_max)
t=t+dt
!	write(*,*)t, maxval(abs(dHdt)), minval(M), maxval(M)
! surface slope
dSdx_mid=(S(2:Nx)-S(1:Nx-1))/dx
! diffusion
D_mid=C*((H(1:Nx-1)+H(2:Nx))/2)**(n+2)*(abs(dSdx_mid))**(n-1)
! flux
Q_mid=-D_mid*dSdx_mid
! change in ice thickness
dHdt(2:Nx-1)=-(Q_mid(2:Nx-1)-Q_mid(1:Nx-2))/dx+M(2:Nx-1)*dt
dHdt(1)=M(1)*dt
dHdt(Nx)=M(Nx)*dt

H=H+dHdt
S=b+H

! apply boundary conditions
H(1) = 0
do ii=1,Nx
if(H(ii).lt.0) then
H(ii)=0
endif
enddo

if(t>=t_old+t_int) then
write (1,*) H
t_old=t
endif

enddo time_loop
close(1)

contains

subroutine checkerr(errstat,msg)
implicit none
integer,      intent(in) :: errstat
character(*), intent(in) :: msg
if (errstat /= 0) then
write(*,*) "ERROR:", msg
stop
end if
end subroutine checkerr

end program```