Team 5 Solution

From Interactive System for Ice sheet Simulation
Revision as of 15:56, 8 August 2009 by Brian anderson (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Invalid language.

You need to specify a language like this: <source lang="html4strict">...</source>

Supported languages for syntax highlighting:

4cs, 6502acme, 6502kickass, 6502tasm, 68000devpac, abap, actionscript, actionscript3, ada, algol68, apache, applescript, apt_sources, arm, asm, asp, asymptote, autoconf, autohotkey, autoit, avisynth, awk, bascomavr, bash, basic4gl, bf, bibtex, blitzbasic, bnf, boo, c, c_loadrunner, c_mac, caddcl, cadlisp, cfdg, cfm, chaiscript, cil, clojure, cmake, cobol, coffeescript, cpp, cpp-qt, csharp, css, cuesheet, d, dcl, dcpu16, dcs, delphi, diff, div, dos, dot, e, ecmascript, eiffel, email, epc, erlang, euphoria, f1, falcon, fo, fortran, freebasic, freeswitch, fsharp, gambas, gdb, genero, genie, gettext, glsl, gml, gnuplot, go, groovy, gwbasic, haskell, haxe, hicest, hq9plus, html4strict, html5, icon, idl, ini, inno, intercal, io, j, java, java5, javascript, jquery, kixtart, klonec, klonecpp, latex, lb, ldif, lisp, llvm, locobasic, logtalk, lolcode, lotusformulas, lotusscript, lscript, lsl2, lua, m68k, magiksf, make, mapbasic, matlab, mirc, mmix, modula2, modula3, mpasm, mxml, mysql, nagios, netrexx, newlisp, nsis, oberon2, objc, objeck, ocaml, ocaml-brief, octave, oobas, oorexx, oracle11, oracle8, oxygene, oz, parasail, parigp, pascal, pcre, per, perl, perl6, pf, php, php-brief, pic16, pike, pixelbender, pli, plsql, postgresql, povray, powerbuilder, powershell, proftpd, progress, prolog, properties, providex, purebasic, pycon, pys60, python, q, qbasic, rails, rebol, reg, rexx, robots, rpmspec, rsplus, ruby, sas, scala, scheme, scilab, sdlbasic, smalltalk, smarty, spark, sparql, sql, stonescript, systemverilog, tcl, teraterm, text, thinbasic, tsql, typoscript, unicon, upc, urbi, uscript, vala, vb, vbnet, vedit, verilog, vhdl, vim, visualfoxpro, visualprolog, whitespace, whois, winbatch, xbasic, xml, xorg_conf, xpp, yaml, z80, zxbasic


! 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