Mountain glacier forward model

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
!-----------------------------
        module kees1_variables
!-----------------------------
 
	real(8), parameter :: dx = 1
	real(8), parameter :: xend = 30
 
        end module kees1_variables
 
!-----------------------------
        subroutine kees1(xx,fc)
!-----------------------------
 
        use kees1_variables
	implicit none
 
! routine arguments
real(8), intent(inout), dimension(int(xend/dx)+1) :: xx
real(8), intent(inout) :: fc
 
! local variables
	real(8), parameter :: dt = 1/12.0
	real(8), parameter :: tend = 5000
	real(8), parameter :: n = 3
	real(8), parameter :: rho = 920
	real(8), parameter :: g = 9.2
	real(8), parameter :: A = 1.0e-16
	real(8), parameter :: Hunnecessary = 6378.1 * sqrt(9.8/9.2-1) !height where g =9.2m/s^2
	real(8), parameter :: C = 2*A*(rho*g)**n/(n+2)*(1.e3)**n
	integer :: nt = int(tend/dt)
	integer :: nx = int(xend/dx)
	real(8), dimension(int(tend/dt)+1,int(xend/dx)+1)  :: S
	integer :: t,i
	real(8), dimension(int(xend/dx)+1)::xarr=(/((i-1)*dx,i=1,&
		int(xend/dx)+1)/), MB, H, B
	real(8), dimension(int(xend/dx)) :: Diff
 
	B = 1-0.0001*xarr		
	MB = 0.004 - xarr * 0.0002
	S(1,:) = B
 
! seed independent variable
        MB = MB + xx
 
	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,:)+Hunnecessary 
		end if
	end do
 
! compute dependent variable fc
        fc = 0.
        do i = 1, size(H)
           fc = fc + H(i)*dx
        end do
 
	end subroutine