Difference between revisions of "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 = 0.1                  ! length time step (years)
  integer, parameter :: nt = 10000                 ! 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.1                ! 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.0*A)/real(n+2)*(rho*g)**n
 
!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
		d(ii) = Const * ((H(ii)+H(ii+1))/2.0)**(n+2) * ((elev(ii+1)-elev(ii))/(dx))**(n-1)
		mflux(ii)= -d(ii) * (elev(ii+1)-elev(ii) )/dx
	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)-mflux(jj-1))/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 ?)
	if (mod(t,10)==0) then
		write (*,*) elev
	endif
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