Team 6 Solution

!> 1D Convection Diffusion equations solver in Fortran !! !! Solves the equation: !! !!\f[ !!\frac{du}{dt}=\frac{d}{dx}D(x)\frac{du}{dx} + C(x)\frac{du}{dx}+F(x)u-S(x) !!\f] !! for \f$u\f$, given functions for \f$D\f$, \f$C\f$, \f$F\f$, and \f$S\f$, defined in this program !! !! Explicit methods are used !! !! \author Matt & Erin & Sophie (jvj) !! \date 8-5-09

program OurCode

 implicit none

 ! local variables

 integer  :: nx                                 ! Number of nodes
 real,    parameter :: dt = 1                  ! length time step (years)
 integer, parameter :: nt = 1000                 ! number of time steps
 integer :: t                                  ! current time step
 real :: xl                                    ! start of domain
 real :: xr                                    ! end of domain (m)
 real :: Const					! (2*A)*(rho*grav)^n/(n+2)
 real, parameter :: dx = 1000                  ! node spacing (m)
 real, parameter :: dbdx = -0.0                ! bedslope (m/m)
 real, parameter :: g = 9.8                    ! gravity (m/s2)
 real, parameter :: rho = 917                ! density of ice (g/cm3)
 real, parameter :: A = 1e-16                   ! Glen rate factor (kPa-3 a-1)
 real, parameter :: n = 3                      ! Glen Flow Exponent (unitless)
 real, parameter :: M0 = 4.0                   ! m/yr
 real, parameter :: M1 = 2.0/10000.0           ! m/yr/m
 
 real, dimension(:), allocatable :: elev       ! surface elevation (m)
 real, dimension(:), allocatable :: bedelev    ! bed elevation (m), y origin is at the bed elev in the left of the domain. up is up!
 real, dimension(:), allocatable :: H          ! thicknesss (m)
 real, dimension(:), allocatable :: Mb         ! Mass Balance (m/yr)
 real, dimension(:), allocatable :: dhdt_store ! space to store du/dx
 real, dimension(:), allocatable :: xref       ! refence distance

 real, dimension(:), allocatable :: d          ! diffusivity coeff
 real, dimension(:), allocatable :: mflux      ! mass flux between grid points

 integer :: ii                                 ! a counter
 integer :: jj  
 integer :: errstat                            ! for error checking
 	

 ! Set up grid 
 ! Space
 xl = 0.0
 xr = 60000.0                 		        !(m) our guess at how large of a domain we need (started with 60km)
 nx = int(   ((xr - xl) / dx) +1   )
 
 ! let's allocate some memory
 allocate(elev(nx),stat=errstat)
 call checkerr(errstat,"failed to allocate elev")

 allocate(xref(nx),stat=errstat)
 call checkerr(errstat,"failed to allocate xref")

 allocate(bedelev(nx),stat=errstat)
 call checkerr(errstat,"failed to allocate bedelev")

 allocate(H(nx),stat=errstat)
 call checkerr(errstat,"failed to allocate H")

 allocate(Mb(nx),stat=errstat)
 call checkerr(errstat,"failed to allocate Mb")

 allocate(d(nx-1),stat=errstat)
 call checkerr(errstat,"failed to allocate d")

 allocate(mflux(nx-1),stat=errstat)
 call checkerr(errstat,"failed to allocate mflux")

 allocate(dhdt_store(nx),stat=errstat)
 call checkerr(errstat,"failed to allocate dhdt_store")

!We could've done this as a function, Matt and Erin revolted vs. Sophie. Said No.

 do ii=1,nx  

bedelev(ii) = dx*dbdx*real(ii-1)  !reminder, at x = 0, bedelev = 0 m Mb(ii) = M0 - M1*(dx*real(ii-1)) !Mass balance equation (m/yr) xref(ii)= real(ii-1)*dx

 enddo

!Constant C

Const = (2.*A)*(rho*g)**n/real(n+2) 

!Initial conditions

 H  =  0.0					!set thickness everywhere in x as 0.0 m
 elev = bedelev + H				!No glacier yet

time_loop: do t=1,nt

spatial_midpoint_loop: do ii=1,nx-1 !Calculate flux midway between elevation points mflux(ii) = Const * ((H(ii)+H(ii+1))/2.0)**(n+2) * ((elev(ii+1)-elev(ii))/(dx))**n enddo spatial_midpoint_loop

H(1) = 0.0  ! left boundary condition (this could be moved above time_loop)

spatial_gridpoint_loop: do jj=2,nx-1 dhdt_store(jj) = - (mflux(jj+1)-mflux(jj))/dx + Mb(jj) H(jj) = H(jj) + dhdt_store(jj) * dt  !new thickness !search for terminus location... if (H(jj)<0) then H(jj) = 0 !write (*,*) jj, (jj-1)*dx

          endif

enddo spatial_gridpoint_loop

!to update 'elev'

	elev = bedelev + H

!output geometry (also need to output geom at t=0 ?) write (*,*) H end do time_loop

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 OurCode