# Team 3 Solution

• FORTRAN mountain glacier

Result for bed slope 0.025

Compile the program with gfortran -o simple_ice_model simple_ice_model_modules.f90 simple_ice_model.f90. The python script does the compiling automatically, so just running python run_and_plot.py will do it all.

Source: Kees_g3.tar.gz‎

### simple_ice_model.f90

```program simple_ice_model
use simple_ice_model_modules
implicit none

! local variables
real(8), parameter    :: pi = 3.1415926536

integer, parameter :: grid=51                    ! Number of nodes
real(8),    parameter :: dt = 0.1                ! length time step
integer, parameter :: t_final=1000                ! integration stop time
integer :: t_step, nt                                ! current time step, number of time steps
real(8), parameter :: xl=0                       ! start of domain
real(8), parameter :: xr=50e3                    ! end of domain
real(8) :: dx                                    ! node spacing
real(8) :: time                                  ! current time
real(8) :: const                                 ! a constant

!! unstaggered grid
real(8), dimension(:), allocatable :: xx         ! node positions
real(8), dimension(:), allocatable :: thick      ! ice thickness
real(8), dimension(:), allocatable :: dthick_dt  ! ice thickness time derivative
real(8), dimension(:), allocatable :: mass_b     ! mass balance
real(8), dimension(:), allocatable :: bed        ! bed elevation

!! staggered grid
real(8), dimension(:), allocatable :: diffus     ! diffusion
real(8), dimension(:), allocatable :: xx_s       ! xx on staggered

!! physical constants
real(8), parameter :: rho=920, g_grav=9.8, A_rate=1e-16, n_glen=3.
real(8), parameter :: dbed_dx=0.025, mass_b_0=4., mass_b_1=0.2e-3 !xxx

!! BC
real(8), parameter :: thickl=0, thickr=0 !xxx

!! counters and such
integer :: ii, jj
integer :: errstat                            ! for error checking

!! functions

! let us allocate some memory
allocate(xx(grid),stat=errstat)
call checkerr(errstat,"failed to allocate xx")

allocate(thick(grid),stat=errstat)
call checkerr(errstat,"failed to allocate thick")

allocate(dthick_dt(grid),stat=errstat)
call checkerr(errstat,"failed to allocate dthick_dt")

allocate(mass_b(grid),stat=errstat)
call checkerr(errstat,"failed to allocate mass_b")

allocate(bed(grid),stat=errstat)
call checkerr(errstat,"failed to allocate bed")

allocate(diffus(grid-1),stat=errstat)
call checkerr(errstat,"failed to allocate diffus")

allocate(xx_s(grid-1),stat=errstat)
call checkerr(errstat,"failed to allocate xx_s")

!! setup
dx = (xr-xl)/(grid-1)
do ii=1,grid
xx(ii) = xl + dx*(ii-1)
thick(ii) = 0   ! IC
mass_b(ii) = mass_b_0 - mass_b_1*xx(ii)
bed(ii) = -dbed_dx * xx(ii)
end do
const = 2*A_rate/(n_glen+2)*(rho*g_grav)**n_glen
nt = floor(t_final/dt) ! number of time steps

!! output, use a pipe to get into a file
print * , 'x'
print * , xx
print * , 'mass_b'
print * , mass_b
print * , 'bed'
print * , bed

print *, 'thick'
do ii=1,nt
! compute the diffusivity: diffus
diffus = calc_diffus(thick, bed, dx, const, n_glen)
! assemble dthick_dt
dthick_dt = pad(diff_1(diffus*diff_1((thick+bed), dx), dx), thickl, thickr) + mass_b
! time step
thick = time_step(thick, dthick_dt, dt)
! BC
thick(1) = thickl
do jj=1,grid
if (thick(jj)<0.) then
thick(jj) = 0.
end if
end do
print *, thick
end do

contains

subroutine checkerr(errstat,msg)
!! the error checker for allocation
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 simple_ice_model```

### simple_ice_model_modules.f90

```module simple_ice_model_modules

contains

function on_stagger(xx)
!! interpolates a vector on the normal grid to the staggered grid
implicit none
real(8), intent(in), dimension(:):: xx
real(8), dimension(size(xx)-1):: on_stagger
integer :: ii
do ii=1,(size(xx)-1)
on_stagger(ii) = (xx(ii+1)+xx(ii))/2
end do
end function on_stagger

function calc_diffus(thick, bed, dx, const, n_glen)
! calculates the nonlinear diffusivity
implicit none
real(8), intent(in), dimension(:):: thick, bed
real(8), dimension(size(thick)-1):: calc_diffus
real(8), intent(in):: const, n_glen, dx
integer :: ii
calc_diffus = const * (on_stagger(thick))**(n_glen+2) * (diff_1((thick+bed), dx))**(n_glen-1)
end function calc_diffus

function diff_1(yy, dx)
! calculates the first staggered derivative of a normal vector
implicit none
real(8), intent(in), dimension(:):: yy
real(8), intent(in) :: dx
real(8), dimension(size(yy)-1):: diff_1
integer :: ii
do ii=1,(size(yy)-1)
diff_1(ii) = (yy(ii+1)-yy(ii))/dx
end do
end function diff_1

function time_step(yy, dyy_dt, dt)
! does a euler time step
implicit none
real(8), intent(in), dimension(:):: yy, dyy_dt
real(8), intent(in) :: dt
real(8), dimension(size(yy)):: time_step
time_step = yy + dt*dyy_dt
end function time_step

! pad a vector with left and right, used after double staggering
implicit none
real(8), intent(in), dimension(:):: on_stag
real(8), intent(in) :: left, right
integer :: ii
do ii=1,size(on_stag)
end do

end module simple_ice_model_modules```

### Python runner and plotter: run_and_plot.py

The runner and plotter. To make it work the fortran files need to be named as above

```import numpy as N
import pylab as P
import os
from py_fortran_tools import *

data_output_file = 'tmp.dat'
f_output_file = 'simple_ice_model'
f_input_files = ['simple_ice_model_modules.f90', 'simple_ice_model.f90']

compile_and_run(f_output_file, f_input_files, data_output_file)
# get the data
dict_ = parse_output(data_output_file)

# plot it
shape_ = dict_['thick'].shape
P.plot(dict_['x']/1e3, dict_['thick'][-1,:]+dict_['bed'], label='At \$t_{final}\$')
P.hold(True)
P.plot(dict_['x']/1e3, dict_['thick'][N.floor(shape_[0]/2),:]+dict_['bed']
, label='at \$t_{final}\$/2')
P.plot(dict_['x']/1e3, dict_['thick'][-100,:]+dict_['bed'],
label='100 time steps before \$t_{final}\$')
P.xlabel('x (km)')
P.ylabel('b+H (m)')
P.legend()
P.show()```

Using the module

```from __future__ import with_statement

import numpy as N
import pylab as P
import os

def parse_output(data_output_file):
"""
This function parses the output of a fortran program generated by
print *, 'x'
print *, x
etc.
statments.

input: data_output_file
output: dictionary with keys corresponding to the variable
names.
"""
# parse the file
with open(data_output_file, 'r') as fil:
dict_ = {}
# go through all the lines in the file
for line in fil:
ll = line.strip()
if ll[0].isalpha():
# if the first character in a line is a letter
# it's the beginning of a new variable
key = ll
# make a new dict_ entry and initialise as empty list
dict_[key] = []
else:
ll = ll.split() # split at white space
num = [float(el) for el in ll] # convert to float
dict_[key].append(num) # append to the list

# convert to numpy arrays
for key in dict_:
dict_[key] = N.array(dict_[key]).squeeze()
return dict_

def compile_and_run(f_output_file, f_input_files, data_output_file):
"""
Compiles a fortran program using
- f_output_file as the executable produced
- f_input_files as a list of the source code files (modules first)
- data_output_file as the data file into which the output to
standard output of the fortran program
is piped
"""

# compile the fortran program
command1 = 'gfortran -o %s %s' %(f_output_file, ' '.join(f_input_files))
exit_ = os.system(command1)
# some primitve error checking
if exit_!=0:
raise(TypeError('\nCompliling with command "%s" failed' %command1))
else:
print 'Compiling ... ok'
# run it
command2 = './%s > %s' %(f_output_file,  data_output_file)
exit_ = os.system(command2)
# some primitve error checking
if exit_!=0:
raise(TypeError('\nRunning of fortran program with command "%s" failed.' %command2))
else:
print 'Running ... finished'```