Difference between revisions of "Team 2 Solution"

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
(deleted unused line)
(updated code and added python script)
 
Line 1: Line 1:
 +
* FORTRAN mountain glacier
 +
 
<source lang="fortran">
 
<source lang="fortran">
 
program Kees1
 
program Kees1
Line 11: Line 13:
 
real, parameter :: g = 9.2
 
real, parameter :: g = 9.2
 
real, parameter :: A = 1.0e-16
 
real, parameter :: A = 1.0e-16
 +
real, parameter :: Hunnecessary = 6378.1 * sqrt(9.8/9.2-1) !height where g =9.2m/s^2
 
real, parameter :: C = 2*A*(rho*g)**n/(n+2)*(1.e3)**n
 
real, parameter :: C = 2*A*(rho*g)**n/(n+2)*(1.e3)**n
 
integer :: nt = int(tend/dt)
 
integer :: nt = int(tend/dt)
Line 38: Line 41:
 
                 end where
 
                 end where
 
if (mod(t,100) .eq.0) then
 
if (mod(t,100) .eq.0) then
write (*,*) S(t,:)+1633.42 !height at which g=9.2
+
write (*,*) S(t,:)+Hunnecessary
 
end if
 
end if
 
end do
 
end do
 
end program
 
end program
 +
 +
</source>
 +
 +
*python plotting script
 +
 +
<source lang="python">
 +
#!/usr/bin/env python
 +
 +
# Import only what is needed
 +
from numpy import loadtxt,shape,linspace
 +
from pylab import plot,show,clf,show,ion,ylim
 +
 +
# Import data file, called 'data'
 +
d=loadtxt('data')
 +
 +
#Determine how much data came in
 +
dims = shape(d)
 +
 +
clf()    # Clears the screen
 +
ion()    # Interactive plot mode, critical for animation
 +
 +
# x data, note that this must correspond to program's domain
 +
x = linspace(0,30,dims[1]) 
 +
 +
# Initial plot, very Matlab(ish), note return of plot handle that allows plot to
 +
# be altered elsewhere in code.
 +
ph,=plot(x,d[0,:],'k')   
 +
ph.figure.show()          # matplot lib requires show to be called
 +
ylim(1629.7,1630.5)
 +
 +
# Loop to plot each time step
 +
for i in range(1,dims[0]):
 +
    ph.set_ydata(d[i,:])  # Only update y data (faster than replot)
 +
    ph.figure.show()
 +
   
 
</source>
 
</source>

Latest revision as of 15:21, 5 August 2009

  • FORTRAN mountain glacier
	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 :: Hunnecessary = 6378.1 * sqrt(9.8/9.2-1) !height where g =9.2m/s^2
	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)=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
	end program
  • python plotting script
#!/usr/bin/env python
 
# Import only what is needed
from numpy import loadtxt,shape,linspace
from pylab import plot,show,clf,show,ion,ylim
 
# Import data file, called 'data'
d=loadtxt('data')
 
#Determine how much data came in
dims = shape(d)
 
clf()    # Clears the screen
ion()    # Interactive plot mode, critical for animation
 
# x data, note that this must correspond to program's domain
x = linspace(0,30,dims[1])  
 
# Initial plot, very Matlab(ish), note return of plot handle that allows plot to
# be altered elsewhere in code.
ph,=plot(x,d[0,:],'k')     
ph.figure.show()           # matplot lib requires show to be called
ylim(1629.7,1630.5)
 
# Loop to plot each time step
for i in range(1,dims[0]):
    ph.set_ydata(d[i,:])   # Only update y data (faster than replot)
    ph.figure.show()