Team 5 Solution

From Interactive System for Ice sheet Simulation
Jump to: navigation, search

Mountain flow.png

! 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