C******************************************************************************
C
C	***************************************
C	*******  > > > > L E M < < < <  *******
C	***************************************
C
C	Landscape Evolution Model
C
C	FORTRAN version 1.2
C	November, 1992
C
C	Programmed and translated from C by Greg Tucker
C
C******************************************************************************
	program main
	implicit none
C******************************************************************************
C	Grid size and other parameters. The arrays are dimensioned from
C	zero to XBND and zero to YBND, so that the boundaries (for a 
C	recantangular grid) are at x=0,x=XBND,y=0,and y=YBND. The area 
C	within the boundaries ranges from 1 to XMAX and 1 to YMAX.
C******************************************************************************
	integer XMAX, YMAX, XBND, YBND, GRIDSIZE 
	parameter (XMAX=60,YMAX=60,XBND=61,YBND=61,GRIDSIZE=3600)
C******************************************************************************
C
C	Variables in common blocks:
C
C		Dt		----	Size of a timestep (years)
C		kfl		----	Fluvial transport coefficient
C		mf,nf		----	Q and S exponents in sed transport func
C		threshold	----	Threshold value for channel initiation
C		mc,nc		----	Q and S exponents in channel init func
C		kc		----	Coefficient for channel init function
C		precip		----	Precipitation: read as m/yr, converted
C					   to m/yr/gridcell
C		uprate		----	Uplift rate: read as m/yr, converted
C					   to m/time-step
C		uptime		----	Duration of uplift, time steps
C		uptype		----	Code for uplift type
C		fault_loc	----	First column of nodes, measured from
C					   west boundary, that undergoes uplift
C		pivot_dist	----	Distance from fault to pivot line, km
C					   (used only for TILTBLOCK uplift)
C		dx		----	Grid cell spacing
C		dxx		----	Diagonal distance between cells
C		elev		----	Cell elevations
C		slope		----	Steepest downward slope @ each cell
C		Q		----	Discharge at each cell
C		load		----	Sediment flux into each cell
C		nbrx,nbry	----	Coordinates of each cell's steepest
C					   neighbor
C		flood		----	Code for flooding status of each cell
C		chan		----	TRUE=channel, FALSE=hillslope
C		cell		----	List of cells (coordinates of each)
C		uca,ucb,ucc	----	Tridiagonal matrix for ADI diff. scheme
C		vca,vcb,vcc	----	     "		"   "   "   "      "
C		ns_wrapped	----	Flag indicating whether north and
C					south boundaries connect (wrap around).
C		east_floats	----	Flag indicating indicating whether
C					east boundary floats (adjusts to 
C					neighboring cells) or is fixed.
C		outlet_north	----	Northern extent of outlet on west bnd.
C		outlet_south	----	Southern extent of outlet on west bnd.
C
C******************************************************************************
	integer ntsteps,dwrite
	integer Dt
	real*8 kfl, mf, nf
	real*8 diff,alpha
	real*8 threshold,mc,nc,kc
	real*8 precip
	real*8 uprate,pivot_dist
	integer uptime,uptype,fault_loc
	real*8 dx,dxx
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	real*8 slope
	dimension slope(0:XBND,0:YBND)
	real*8 Q
	dimension Q(0:XBND,0:YBND)
	real*8 load
	dimension load(0:XBND,0:YBND)
	integer nbrx,nbry
	dimension nbrx(0:XBND,0:YBND)
	dimension nbry(0:XBND,0:YBND)
	byte flood
	dimension flood(0:XBND,0:YBND)
	logical chan
	dimension chan(0:XBND,0:YBND)
	structure /cellcoord/
		integer x
		integer y
	end structure
	record /cellcoord/ cell(GRIDSIZE)
	real*8 uca,ucb,ucc,vca,vcb,vcc
	dimension uca(GRIDSIZE),ucb(GRIDSIZE),ucc(GRIDSIZE)
	dimension vca(GRIDSIZE),vcb(GRIDSIZE),vcc(GRIDSIZE)
	logical east_floats,ns_wrapped
	integer outlet_north,outlet_south
C******************************************************************************
C	Common Blocks
C******************************************************************************
	common /com1/ elev 
	common /com2/ slope
	common /com3/ Q
	common /com4/ load
	common /com5/ nbrx,nbry
	common /com6/ flood
	common /com7/ chan
	common /com8/ cell
	common /com11/ ntsteps,dwrite
	common /com12/ Dt
	common /com13/ kfl, mf, nf
	common /com14/ diff,alpha
	common /com15/ threshold,mc,nc,kc
	common /com16/ precip
	common /com17/ uprate,uptime,uptype,fault_loc,pivot_dist
	common /com18/ dx,dxx 
	common /com19/ uca,ucb,ucc,vca,vcb,vcc
	common /com20/ east_floats,ns_wrapped,outlet_north,outlet_south
C******************************************************************************
C	Variables local to MAIN
C******************************************************************************
	integer i, j, time, next_write
C******************************************************************************
C
C	MAIN	
C
C	The MAIN routine starts by calling the INITIALIZE subroutine, and 
C	calling WRITE_ELEV_DATA to record initial elevations. After that,
C	the routine counts time steps in a DO loop, calling subroutines to
C	compute drainage, sediment transport and uplift for each time step.
C
C		Calls: INITIALIZE, RESET_VARIABLES, FIND_NEXT_NODE,
C			DO_PONDS, STREAM_TRACE, CHANNEL_INIT, FLUVIAL,
C			HILLSLOPE, WRITE_ELEV_DATA, WRITE_NET_DATA, UPLIFT 
C
C******************************************************************************
	call initialize
	call write_elev_data( 0 )
	next_write = dwrite
	do time=1, ntsteps 
		print *, time
		call reset_variables
		do i=1,XMAX
			do j=1,YMAX
				call find_next_node(i,j)
			end do
		end do
		do i=1,XMAX
			do j=1,YMAX
				if( flood(i,j) .eq. 3 ) then
					call drain_basin( i, j )
				endif
			end do
		end do
		call stream_trace
		call channel_init( threshold, kc, mc, nc )
		if( diff .gt. 0.0 ) call hillslope
		call fluvial
		if( time .eq. next_write ) then
			call write_elev_data( time )
			call write_net_data( time )
			next_write = next_write + dwrite
		endif
		if( time .le. uptime ) call uplift()
	end do
	close(1)
	close(2)
	end
C******************************************************************************
C
C	Subroutine INITIALIZE
C
C	This subroutine reads model parameter values from a file called
C	'lem.in'. Also read are names for the initial topography file
C	and the output files. Once data have been read in, INITIALIZE
C	adjusts some of the parameter values and then sets up the arrays
C	used in the mass diffusion computation.
C
C		Calls: READ_TOPO, MAKE_SORT_LIST, OPEN_FILES
C		Called by: MAIN
C	
C	Local variables:
C		RHO		-	Density of water (kg/m^3)
C		SIGMA		-	Bulk density of sediment (kg/m^3)
C		SECPERHOUR	-	Number of seconds in an hour
C		HOURPERYEAR	-	Number of hours in a year
C
C******************************************************************************
	subroutine initialize()
	implicit none
	integer XMAX,YMAX,GRIDSIZE
	parameter(XMAX=60,YMAX=60,GRIDSIZE=3600)
	integer TILTBLOCK 
	parameter(TILTBLOCK=3)
	integer ntsteps,dwrite
	integer Dt
	real*8 kfl, mf, nf
	real*8 diff,alpha
	real*8 threshold,mc,nc,kc
	real*8 precip
	real*8 uprate,pivot_dist
	integer uptime,uptype,fault_loc
	real*8 dx,dxx
	real*8 uca,ucb,ucc,vca,vcb,vcc
	dimension uca(GRIDSIZE),ucb(GRIDSIZE),ucc(GRIDSIZE)
	dimension vca(GRIDSIZE),vcb(GRIDSIZE),vcc(GRIDSIZE)
	logical east_floats,ns_wrapped
	integer outlet_north,outlet_south
	common /com11/ ntsteps,dwrite
	common /com12/ Dt
	common /com13/ kfl, mf, nf
	common /com14/ diff,alpha
	common /com15/ threshold,mc,nc,kc
	common /com16/ precip
	common /com17/ uprate,uptime,uptype,fault_loc,pivot_dist
	common /com18/ dx,dxx 
	common /com19/ uca,ucb,ucc,vca,vcb,vcc
	common /com20/ east_floats,ns_wrapped,outlet_north,outlet_south
	character*16 topofile 
	character*16 outfile
	real*8 RHO,SIGMA
	real*8 SECPERHOUR
	real*8 HOURPERYEAR
	integer i,j,k
	open(1,file='lem.in')
	read(1,*) ntsteps
	read(1,*) Dt
	read(1,*) kfl 
	read(1,*) mf 
	read(1,*) nf 
	read(1,*) diff 
	read(1,*) threshold 
	read(1,*) kc 
	read(1,*) mc
	read(1,*) nc
	read(1,*) precip 
	read(1,*) uprate 
	read(1,*) uptime
	read(1,*) uptype
	read(1,*) fault_loc
	if( uptype .eq. TILTBLOCK ) then
		read(1,*) pivot_dist
	endif
	read(1,*) dwrite
	read(1,*) topofile
	read(1,*) outfile
	read(1,*) k
	if( k .eq. 1 ) then 
		east_floats = .TRUE.
	else 
		east_floats = .FALSE.
	endif
	read(1,*) k
	if( k .eq. 1 ) then
		ns_wrapped = .TRUE.
	else
		ns_wrapped = .FALSE.
	endif
	read(1,*) outlet_south
	read(1,*) outlet_north
	close(1)
C******************************************************************************
C	Include all the constant terms in the stream power coefficient,
C	and convert to SI units (the operation is broken up to avoid
C	round-off error).
C******************************************************************************
	dx = 500.0
	RHO = 1000.0
	SIGMA = 2100.0
	SECPERHOUR = 3600.0
	HOURPERYEAR = 8760.0
	kfl = ((kfl*Dt*RHO)/SIGMA) * (SECPERHOUR/dx) * (HOURPERYEAR/dx) 
C******************************************************************************
C	Calculate the diagonal distance between grid points
C******************************************************************************
	dxx = sqrt( 2*dx*dx )
C******************************************************************************
C	Convert precip from m/yr to m^3/s/gridcell (operation broken up
C	to avoid round-off error)
C******************************************************************************
	precip = precip * (dx/SECPERHOUR) * (dx/HOURPERYEAR)
C******************************************************************************
C	Convert uplift rate from m/yr to m/timestep
C******************************************************************************
	uprate = uprate * Dt
C******************************************************************************
C	Set up tridiagonal matrix (3 vectors) for ADI diffusion solution
C******************************************************************************
	if( diff .gt. 0 ) then
		alpha = (diff*Dt*0.5)/(dx*dx)
		do k=1,GRIDSIZE
			uca(k) = -1.0
			ucb(k) = 1.0/alpha + 2.0
			ucc(k) = -1.0
			vca(k) = -1.0
			vcb(k) = 1.0/alpha + 2.0
			vcc(k) = -1.0
		end do
		do j=1,YMAX
			k = (j-1)*XMAX + 1
			uca(k) = 0.0
			k = j*XMAX
			ucc(k) = 0.0
		end do
		do i=1,XMAX
			k = (i-1)*YMAX + 1
			vca(k) = 0.0
			k = i*YMAX
			vcc(k) = 0.0
		end do
	endif
C******************************************************************************
C	Call subroutines to initialize elevations, output, etc.
C******************************************************************************
	call read_topo( topofile )
	call make_sort_list
	call open_files( outfile )
	return
	end
C******************************************************************************
C
C	Subroutine READ_TOPO	
C
C	This subroutine reads the initial topography from the file named
C	in LEM.IN (see above).
C
C		Calls: (none)
C		Called by: INITIALIZE
C
C******************************************************************************
	subroutine read_topo( tfile )
	implicit none
	character*16 tfile
	integer XBND, YBND
	parameter (XBND=61,YBND=61)
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	common /com1/ elev
	integer i,j
	open(2,file=tfile)
	do j=0,YBND
		do i=0,XBND
			read(2,*,end=20) elev(i,j)
		end do
	end do
	close( 2 )
	return
20	print *, 'Cannot read topo file: End of file found'
	call exit(0)
	end
C******************************************************************************
C
C	Subroutine MAKE_SORT_LIST	
C
C	This subroutine enters the coordinates of each non-boundary cell 
C	into the CELL list, so that cells can conveniently be sorted.
C	(For FORTRAN versions that do not support structures, the list
C	can be broken into CELLX() and CELLY() arrays).
C
C		Calls: (none)
C		Called by: INITIALIZE
C
C******************************************************************************
	subroutine make_sort_list()
	parameter (XMAX=60,YMAX=60,GRIDSIZE=3600)
	structure /cellcoord/
		integer x
		integer y
	end structure
	record /cellcoord/ cell(GRIDSIZE)
	common /com8/ cell
	integer t
	t = 1
	do i=1,XMAX
		do j=1,YMAX
			cell(t).x = i
			cell(t).y = j
			t = t + 1
		end do
	end do		
	return
	end
C******************************************************************************
C
C	Subroutine OPEN_FILES	
C
C	This subroutine opens the output data files and writes preliminary
C	data. Each file is given a basic name that identifies the run 
C	(read from LEM.IN) along with an extension that identifies the
C	nature of the data (for example, RUN5.ELEV). This version of the
C	model maintains two output files that contain elevation and
C	drainage network data, respectively.
C
C		Calls: STRCAT
C		Called by: INITIALIZE
C
C******************************************************************************
	subroutine open_files( fname )
	implicit none
	character*16 fname
	integer XMAX, YMAX, XBND, YBND 
	parameter (XMAX=60,YMAX=60,XBND=61,YBND=61)
	character*24 fullname
	character*24 strcat
	integer ntsteps,dwrite
	real*8 dx,dxx
	common /com11/ ntsteps,dwrite
	common /com18/ dx,dxx 
	integer tts
	tts = ntsteps / dwrite
	fullname = strcat( fname, '.elev ' )
	open(1,status='new',file=fullname)
	write(1,*) XBND+1 
	write(1,*) YBND+1
	write(1,*) dx 
	fullname = strcat( fname, '.net ' )
	open(2,status='new',file=fullname)
	write(2,*) tts
	return
	end
C******************************************************************************
C
C	Function STRCAT
C
C	This function concatenates the string ADDSTRING onto INSTRING and
C	returns the resulting string.
C
C		Calls: (none)
C		Called by: OPEN_FILES
C
C******************************************************************************
	character*24 function strcat( instring, addstring )
	implicit none
	character*16 instring
	character*8 addstring
	integer k1, k2
	do k1=1,16
		if( instring(k1:k1) .eq. ' ') goto 100
	end do
100	k1=k1-1
	do k2=1,8
		if( addstring(k2:k2) .eq. ' ') goto 110
	end do
110	k2=k2-1
	strcat=instring(1:k1)//addstring(1:k2)//char(0)
	return
	end
C******************************************************************************
C
C	Subroutine RESET_VARIABLES
C
C	This subroutine cleans up variables for a new time step.
C
C		Calls: (none)
C		Called by: MAIN
C
C******************************************************************************
	subroutine reset_variables()
	integer i, j
	integer XMAX, YMAX, XBND, YBND,TILTBLOCK 
	parameter (XMAX=60,YMAX=60,XBND=61,YBND=61)
	parameter(TILTBLOCK=3)
	real*8 elev 
	dimension elev(0:XBND,0:YBND)
	byte flood
	dimension flood(0:XBND,0:YBND)
	real*8 load
	dimension load(0:XBND,0:YBND)
	real*8 Q
	dimension Q(0:XBND,0:YBND)
	real*8 uprate,pivot_dist
	integer uptime,uptype,fault_loc
	logical east_floats,ns_wrapped
	integer outlet_north,outlet_south
	common /com1/ elev 
	common /com3/ Q
	common /com4/ load
	common /com6/ flood
	common /com17/ uprate,uptime,uptype,fault_loc,pivot_dist
	common /com20/ east_floats,ns_wrapped,outlet_north,outlet_south
	real*8 h,delh
	do i=0,XBND
		do j=0,YBND
			flood(i,j) = 0
			load(i,j) = 0.0
			Q(i,j) = 0.0
		end do
	end do
C******************************************************************************
C	Adjust north and south boundaries
C******************************************************************************
	if( ns_wrapped ) then
		do i=1,XMAX
			elev(i,0) = elev(i,YMAX)
			elev(i,YBND) = elev(i,1)
		end do
	else
		do i=1,XMAX
			elev(i,0) = elev(i,1)
			elev(i,YBND) = elev(i,YMAX)
		end do
	endif
C******************************************************************************
C	Adjust east boundary:
C	In the case of a tilt block, east boundary elevations are adjusted
C	so as to maintain the same slope at the boundary as exists between
C	the adjacent interior cells.
C******************************************************************************
	if( east_floats ) then
		if( uptype .eq. TILTBLOCK ) then
			h=0.0
			delh=0.0
			do j=1,YMAX
				h=h+elev(XMAX,j)
				delh=delh+0.5*(elev(XMAX-2,j)-elev(XMAX,j))
			end do
			h=(h-delh)/real(YMAX)
			if( h .lt. 0.0 ) h=0.0
			do j=0,YBND
				elev(XBND,j)=h
			end do
		else
			do j=0,YBND
				elev(XBND,j) = elev(XMAX,j)
			end do
		endif
	endif
C******************************************************************************
C	Adjust west boundary outside of outlet (if applicable)
C******************************************************************************
	if( outlet_south .gt. 0 ) then
		do j=0,outlet_south-1
			elev(0,j) = elev(1,j)
		end do
	endif
	if( outlet_north .lt. YBND ) then
		do j=outlet_north+1,YBND
			elev(0,j) = elev(1,j)
		end do
	endif
	return
	end
C******************************************************************************
C
C	Subroutine FIND_NEXT_NODE	
C
C	This subroutine identifies the direction of steepest downward slope
C	at a given cell, recording the coordinates of the adjacent (or
C	"neighbor") cell in the NBRX() and NBRY() arrays, and the slope in 
C	the SLOPE() array.
C
C		Calls: FIND_SLOPE
C		Called by: MAIN
C
C******************************************************************************
	subroutine find_next_node( x, y )
	implicit none
	integer x, y
	integer XMAX,YMAX,XBND,YBND 
	parameter (XMAX=60,YMAX=60,XBND=61,YBND=61)
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	byte flood
	dimension flood(0:XBND,0:YBND)
	integer nbrx,nbry
	dimension nbrx(0:XBND,0:YBND)
	dimension nbry(0:XBND,0:YBND)
	real*8 slope
	dimension slope(0:XBND,0:YBND)
	real*8 dx,dxx
	logical east_floats,ns_wrapped
	integer outlet_north,outlet_south
	common /com1/ elev
	common /com2/ slope
	common /com5/ nbrx,nbry
	common /com6/ flood
	common /com18/ dx,dxx
	common /com20/ east_floats,ns_wrapped,outlet_north,outlet_south
	real*8 slopen,slopes,slopee,slopew,slopene,slopenw,slopese,slopesw
	real*8 S
	real*8 find_slope
	integer ynorth,ysouth
C******************************************************************************
C	Calculate slope in each direction.
C******************************************************************************
	slopen = find_slope( elev(x,y), elev(x,y+1), dx )
	slopene = find_slope( elev(x,y), elev(x+1,y+1), dxx )
	slopee = find_slope( elev(x,y), elev(x+1,y), dx )
	slopese = find_slope( elev(x,y), elev(x+1,y-1), dxx )
	slopes = find_slope( elev(x,y), elev(x,y-1), dx )
	slopesw = find_slope( elev(x,y), elev(x-1,y-1), dxx )
	slopew = find_slope( elev(x,y), elev(x-1,y), dx )
	slopenw = find_slope( elev(x,y), elev(x-1,y+1), dxx )
C******************************************************************************
C	Set 'ynorth' and 'ysouth' cells according to NS boundary condition
C******************************************************************************
	ynorth = y+1
	ysouth = y-1
	if( ns_wrapped ) then
		if( y .eq. YMAX ) then
			ynorth = 1
		else if( y .eq. 1 ) then
			ysouth = YMAX 
		endif
	endif
C******************************************************************************
C	Check each direction to find the steepest, and assign neighbor cell.
C******************************************************************************
	S = slopen
	nbrx(x,y) = x
	nbry(x,y) = ynorth
	if( slopene .gt. S ) then
		S = slopene
		nbrx(x,y) = x+1
	endif
	if( slopee .gt. S ) then
		S = slopee
		nbrx(x,y) = x+1
		nbry(x,y) = y
	endif
	if( slopese .gt. S ) then
		S = slopese
		nbrx(x,y) = x+1
		nbry(x,y) = ysouth
	endif
	if( slopes .gt. S ) then
		S = slopes
		nbrx(x,y) = x
		nbry(x,y) = ysouth
	endif
	if( slopesw .gt. S ) then
		S = slopesw
		nbrx(x,y) = x-1
		nbry(x,y) = ysouth
	endif
	if( slopew .gt. S ) then
		S = slopew
		nbrx(x,y) = x-1
		nbry(x,y) = y
	endif
	if( slopenw .gt. S ) then
		S = slopenw
		nbrx(x,y) = x-1
		nbry(x,y) = ynorth
	endif
C******************************************************************************
C	Record steepest slope in SLOPE(). Negative slope means a closed basin.	
C******************************************************************************
	if( S .le. 0.0 ) then
		flood(x,y) = 3
		slope(x,y) = 0.0
	else
		slope(x,y) = S
	endif
	return
	end
C******************************************************************************
C
C	Subroutine DRAIN_BASIN
C
C	This subroutine calculates the drainage of closed basins (lakes).
C	To do this, it uses an iterative algorithm that maintains a list of
C	all cells that are in the flooded region. During each iteration, the
C	lowest cell along the perimeter of the flooded region is identified.
C	If this cell can drain downwards into a cell that is not part of the
C	flooded region (but which can be part of a previously calculated
C	lake), then that cell is flagged as the outlet for the lake. If not,
C	then the cell is added to the list and the process repeats.
C	Boundary cells are always considered to be outlets.
C
C		Calls: IS_BOUNDARY, FIND_DRY_NEIGHBOR 
C		Called by: MAIN
C
C******************************************************************************
	subroutine drain_basin( x, y )
	implicit none
	integer x, y
	integer XMAX,YMAX,XBND,YBND
	parameter(XMAX=60,YMAX=60,XBND=61,YBND=61)
	byte UNFLOODED,NEWPOND,OLDPOND,BASIN
	real*8 VERYHIGH
	parameter(UNFLOODED=0,NEWPOND=1,OLDPOND=2,BASIN=3)
	parameter(VERYHIGH=500000.0)
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	byte flood
	dimension flood(0:XBND,0:YBND)
	integer nbrx,nbry
	dimension nbrx(0:XBND,0:YBND)
	dimension nbry(0:XBND,0:YBND)
	logical east_floats,ns_wrapped
	integer outlet_north,outlet_south
	common /com1/ elev
	common /com5/ nbrx,nbry 
	common /com6/ flood
	common /com20/ east_floats,ns_wrapped,outlet_north,outlet_south
	integer n,a,b,bb,nnodes
	logical done
	integer LISTMAX
	parameter(LISTMAX=3600)
	structure /cellcoord/
		integer x
		integer y
	end structure
	record /cellcoord/ fcell
	dimension fcell(LISTMAX)
	integer lowx,lowy
	real*8 lowelev
	logical is_boundary,find_dry_neighbor	
C******************************************************************************
C	At first, the flooded region consists of only one cell (x,y)
C******************************************************************************
	nnodes = 1
	flood(x,y) = NEWPOND
	fcell(1).x = x
	fcell(1).y = y
	done = .FALSE.
C******************************************************************************
C	Iterate until an outlet for the flooded region is found
C******************************************************************************
	do while( .not. done )
		lowelev = VERYHIGH
C******************************************************************************
C	Find the lowest cell on the perimeter of the flooded area
C	(A WHILE loop is needed because NNODES can change within the loop)
C******************************************************************************
		n=1
		do while( n .le. nnodes )
			do a=fcell(n).x-1,fcell(n).x+1
				do bb=fcell(n).y-1,fcell(n).y+1
					b = bb
					if( ns_wrapped ) then
						if( bb .eq. 0 ) then
							b = YMAX
						else if( bb .eq. YBND ) then
							b = 1
						endif
					endif
					if(flood(a,b).eq.UNFLOODED
     &  				  .or.flood(a,b).eq.BASIN) then
						if(elev(a,b).lt.lowelev) then
							lowx = a
							lowy = b
							lowelev = elev(a,b)
						endif
					else if( flood(a,b) .eq. OLDPOND ) then
						flood(a,b) = NEWPOND
						nnodes = nnodes+1
						fcell(nnodes).x = a
						fcell(nnodes).y = b
					endif
				end do
			end do
			n = n+1
		end do
C******************************************************************************
C	Now we've found the lowest cell. If it's a boundary cell, then treat 
C	it as an outlet for the flooded area and we're done iterating. If not,
C	scan all the cells surrounding the "low" cell. If the "low" cell can
C	drain to another cell that isn't part of this flooded region, then
C	the "low" cell is the outlet, and we're done iterating.
C******************************************************************************
		if( is_boundary(lowx,lowy) ) then
			done = .TRUE.
		else if( find_dry_neighbor( lowx, lowy ) ) then
			done = .TRUE.
C******************************************************************************
C	If these conditions aren't met, the "low" cell becomes part of the
C	flooded region and we continue iterating to find an outlet.
C******************************************************************************
		else
			flood(lowx,lowy) = NEWPOND
			nnodes = nnodes+1
			fcell(nnodes).x = lowx
			fcell(nnodes).y = lowy
			done = .FALSE.
		endif
	end do
C******************************************************************************
C	Now we've found an outlet at LOWX,LOWY. This cell now becomes the
C	neighbor cell for all the cells in this flooded region. Each flooded
C	cell is recorded as an OLDPOND, to distinguish it from other flooded
C	regions that might be calculated later during this time step.
C******************************************************************************
	do n=1,nnodes
		nbrx(fcell(n).x,fcell(n).y) = lowx
		nbry(fcell(n).x,fcell(n).y) = lowy
		flood(fcell(n).x,fcell(n).y) = OLDPOND 
	end do
	return
	end
C******************************************************************************
C
C	Function FIND_DRY_NEIGHBOR 
C
C	This function is nearly identical to the subroutine FIND_NEXT_NODE,
C	with this difference: in order to be considered as a neighbor cell,
C	a cell must be "dry." This means that it must not be part of the
C	flooded region currently being processed by DRAIN_BASIN, and must
C	not drain into the current cell (see ISDRY). The function evaluates
C	to TRUE when a downslope neighbor that meets these criteria is found.
C
C		Calls: FIND_SLOPE, ISDRY
C		Called by: DRAIN_BASIN
C
C******************************************************************************
	logical function find_dry_neighbor( x, y )
	implicit none
	integer x, y
	integer XMAX,YMAX,XBND, YBND 
	parameter (XMAX=60,YMAX=60,XBND=61,YBND=61)
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	integer nbrx,nbry
	dimension nbrx(0:XBND,0:YBND)
	dimension nbry(0:XBND,0:YBND)
	real*8 slope
	dimension slope(0:XBND,0:YBND)
	real*8 dx,dxx
	logical east_floats,ns_wrapped
	integer outlet_north,outlet_south
	common /com1/ elev
	common /com2/ slope
	common /com5/ nbrx,nbry
	common /com18/ dx,dxx
	common /com20/ east_floats,ns_wrapped,outlet_north,outlet_south
	real*8 slopen,slopes,slopee,slopew,slopene,slopenw,slopese,slopesw
	real*8 S
	real*8 find_slope
	logical isdry 
	integer ynorth,ysouth
C******************************************************************************
C	Calculate slope in each direction.
C******************************************************************************
	slopen = find_slope( elev(x,y), elev(x,y+1), dx )
	slopene = find_slope( elev(x,y), elev(x+1,y+1), dxx )
	slopee = find_slope( elev(x,y), elev(x+1,y), dx )
	slopese = find_slope( elev(x,y), elev(x+1,y-1), dxx )
	slopes = find_slope( elev(x,y), elev(x,y-1), dx )
	slopesw = find_slope( elev(x,y), elev(x-1,y-1), dxx )
	slopew = find_slope( elev(x,y), elev(x-1,y), dx )
	slopenw = find_slope( elev(x,y), elev(x-1,y+1), dxx )
C******************************************************************************
C	Set 'ynorth' and 'ysouth' cells according to NS boundary condition
C******************************************************************************
	ynorth = y+1
	ysouth = y-1
	if( ns_wrapped ) then
		if( y .eq. YMAX ) then
			ynorth = 1
		else if( y .eq. 1 ) then
			ysouth = YMAX 
		endif
	endif
C******************************************************************************
C	Check each direction to find the steepest, and assign neighbor cell.
C******************************************************************************
	S = -10.0 
	if( slopen .gt. S .and. isdry(x,ynorth,x,y) ) then
		S = slopen
		nbrx(x,y) = x
		nbry(x,y) = ynorth
	endif
	if( slopene .gt. S .and. isdry(x+1,ynorth,x,y) ) then
		S = slopene
		nbrx(x,y) = x+1
		nbry(x,y) = ynorth 
	endif
	if( slopee .gt. S .and. isdry(x+1,y,x,y) ) then
		S = slopee
		nbrx(x,y) = x+1
		nbry(x,y) = y
	endif
	if( slopese .gt. S .and. isdry(x+1,ysouth,x,y) ) then
		S = slopese
		nbrx(x,y) = x+1
		nbry(x,y) = ysouth 
	endif
	if( slopes .gt. S .and. isdry(x,ysouth,x,y) ) then
		S = slopes
		nbrx(x,y) = x
		nbry(x,y) = ysouth 
	endif
	if( slopesw .gt. S .and. isdry(x-1,ysouth,x,y) ) then
		S = slopesw
		nbrx(x,y) = x-1
		nbry(x,y) = ysouth 
	endif
	if( slopew .gt. S .and. isdry(x-1,y,x,y) ) then
		S = slopew
		nbrx(x,y) = x-1
		nbry(x,y) = y
	endif
	if( slopenw .gt. S .and. isdry(x-1,ynorth,x,y) ) then
		S = slopenw
		nbrx(x,y) = x-1
		nbry(x,y) = ynorth 
	endif
C******************************************************************************
C	If a "dry" cell with a positive downward slope exists, the fn is TRUE	
C******************************************************************************
	if( S .le. 0.0 ) then
		slope(x,y) = 0.0
		find_dry_neighbor = .FALSE.
	else
		slope(x,y) = S
		find_dry_neighbor = .TRUE.
	endif
	return
	end
C******************************************************************************
C
C	Function ISDRY
C
C	This function evaluates whether a given cell is "dry." Dry means
C	that it is not part of a flooded region currently processed by
C	the DRAIN_BASIN subroutine (in which case it is coded as a NEWPOND),
C	and does not drain back to the cell LOWX, LOWY.
C
C		Calls: (none)
C		Called by: FIND_DRY_NEIGHBOR
C
C******************************************************************************
	logical function isdry( x, y, lowx, lowy )
	implicit none
	integer x,y,lowx,lowy
	byte UNFLOODED,NEWPOND,BASIN
	parameter(UNFLOODED=0,NEWPOND=1,BASIN=3)
	integer XBND,YBND
	parameter(XBND=61,YBND=61)
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	byte flood
	dimension flood(0:XBND,0:YBND)
	integer nbrx,nbry
	dimension nbrx(0:XBND,0:YBND),nbry(0:XBND,0:YBND)
	common /com1/ elev 
	common /com5/ nbrx,nbry
	common /com6/ flood
	integer outx,outy
	if( flood(x,y) .eq. UNFLOODED .or. flood(x,y) .eq. BASIN ) then
		isdry = .TRUE.
	else if( flood(x,y) .eq. NEWPOND ) then
		isdry = .FALSE.
C******************************************************************************
C	If it passes these tests, it's an OLDPOND. Occasionally, the outlet
C	may be to the LOWCELL itself, or else the OLDPOND cell may have an
C	outlet whose elevation is equal to the LOWCELL. In these situations,
C	the test fails, since an endless loop would result.
C******************************************************************************
	else
		outx = nbrx(x,y)
		outy = nbry(x,y)
		if( elev(outx,outy) .eq. elev(lowx,lowy) ) then
			isdry = .FALSE.
		else
			isdry = .TRUE.
		endif
	endif
	return
	end
C******************************************************************************
C
C	Subroutine STREAM_TRACE	
C
C	This function calculates the discharge at each cell by tracing the
C	route of each incoming "parcel" of precipitation, until that "parcel"
C	flows into a boundary cell.
C
C		Calls: (none)
C		Called by: MAIN
C
C******************************************************************************
	subroutine stream_trace()
	implicit none
	integer XMAX, YMAX, XBND, YBND 
	parameter (XMAX=60,YMAX=60,XBND=61,YBND=61)
	integer i, j, m, n, newm
	logical is_boundary
	real*8 Q
	dimension Q(0:XBND,0:YBND)
	integer nbrx,nbry
	dimension nbrx(0:XBND,0:YBND)
	dimension nbry(0:XBND,0:YBND)
	real*8 precip
	common /com3/ Q
	common /com5/ nbrx,nbry
	common /com16/ precip
	do i=1,XMAX
		do j=1,YMAX
			m = i
			n = j
			Q(i,j) = Q(i,j) + precip
			do while( .not. is_boundary(nbrx(m,n),nbry(m,n) ) ) 
				newm = nbrx(m,n)
				n = nbry(m,n)
				m = newm
				Q(m,n) = Q(m,n) + precip
			end do
		end do
	end do
	return
	end
C******************************************************************************
C
C	Subroutine CHANNEL_INIT	
C
C	This subroutine tests for the initiation of channels, based on a
C	threshold function of slope and discharge. Channels are defined to
C	exist wherever this function exceeds a given threshold value.
C	Elsewhere, cells are considered to be hillslopes and subject only to
C	hillslope processes. Channels are therefore not permanent, but can
C	appear and disappear as tectonics and stream piracy alter the form
C	of the drainage network.
C
C		Calls: (none)
C		Called by: MAIN
C
C******************************************************************************
	subroutine channel_init( threshold, kc, mc, nc )
	implicit none
	real*8 threshold,kc,mc,nc
	integer XMAX, YMAX, XBND, YBND 
	parameter (XMAX=60,YMAX=60,XBND=61,YBND=61)
	real*8 slope
	dimension slope(0:XBND,0:YBND)
	real*8 Q
	dimension Q(0:XBND,0:YBND)
	byte flood	
	dimension flood(0:XBND,0:YBND)
	logical chan
	dimension chan(0:XBND,0:YBND)
	common /com2/ slope
	common /com3/ Q
	common /com6/ flood 
	common /com7/ chan
	real*8 cinit
	integer i,j
	do i=1,XMAX
		do j=1,YMAX
			if( flood(i,j) .eq. 0 ) then
				cinit = kc*(Q(i,j)**mc)*(slope(i,j)**nc)
				if( cinit .gt. threshold ) then
					chan(i,j) = .true. 
				else
					chan(i,j) = .false. 
				endif
			else
				chan(i,j) = .false. 
			endif
		end do
	end do
	return
	end
C******************************************************************************
C
C	Subroutine FLUVIAL	
C
C	This subroutine computes sediment transport by streams. The carrying
C	capacity of a stream segment is defined to be proportional to local
C	stream power. In general, some of this capacity will be taken up by
C	sediment passed from upstream. For this reason, cells must be processed
C	in upstream-to-downstream order (see SORT_BY_ORDER).
C
C	The subroutine uses the BISECT function to find an elevation change
C	DH such that the resultant stream power is exactly sufficient to
C	transport the resulting sediment load (this method is more stable
C	than computing stream power based on the initial, rather than the
C	final slope). The LOW and HIGH values are the bounds for the bisection
C	algorithm. They represent the cases of maximum erosion and maximum
C	deposition, respectively. Once DH has been found, the new sediment
C	load at the cell is added to the sediment load of the downstream cell.
C
C	Flooded cells are treated as follows: the sediment load passed to 
C	each flooded cell is deposited in that cell, up to just above the 
C	water level (which is defined by the elevation of the outlet cell;
C	evaporation, etc, are not considered). Any excess sediment is then
C	passed directly to the outlet cell. This treatment is considered
C	reasonable as long as flooded areas and sediment loads are small.
C
C	For hillslope cells, any sediment loads passed from discontinuous
C	channel segments upstream are simply dumped. 
C
C		Calls: SORT_BY_ORDER, OMEGA, BISECT
C		Called by: MAIN
C
C******************************************************************************
	subroutine fluvial()
	implicit none
	integer XBND,YBND,GRIDSIZE
	parameter (XBND=61,YBND=61,GRIDSIZE=3600)
	structure /cellcoord/
		integer x
		integer y
	end structure
	record /cellcoord/ cell(GRIDSIZE)
	real*8 Q
	dimension Q(0:XBND,0:YBND)
	logical chan
	dimension chan(0:XBND,0:YBND)
	real*8 omega
	real*8 bisect 
	integer a,x,y,nx,ny
	real*8 dh,dist,low,high
	real*8 slope
	dimension slope(0:XBND,0:YBND)
	real*8 elev 
	dimension elev(0:XBND,0:YBND)
	integer nbrx,nbry
	dimension nbrx(0:XBND,0:YBND)
	dimension nbry(0:XBND,0:YBND)
	real*8 load
	dimension load(0:XBND,0:YBND)
	byte flood
	dimension flood(0:XBND,0:YBND)
	real*8 dx,dxx
	common /com1/ elev 
	common /com2/ slope
	common /com3/ Q
	common /com4/ load
	common /com5/ nbrx,nbry 
	common /com6/ flood 
	common /com7/ chan
	common /com8/ cell
	common /com18/ dx,dxx 
	real*8 depth,SMALL
	SMALL = 0.001
	call sort_by_order
	do a=1,GRIDSIZE
		x = cell(a).x
		y = cell(a).y
		nx = nbrx(x,y)
		ny = nbry(x,y)
C******************************************************************************
C	Channel cells.
C******************************************************************************
		if( chan(x,y) ) then
			if( abs(x-nx) .eq. 0 .or. abs(y-ny) .eq. 0 ) then
				dist = dx
			else
				dist = dxx
			endif
			low = -1 * omega(slope(x,y),Q(x,y),dist,0.0D0)
			high = load(x,y)
			dh = bisect(low,high,slope(x,y),Q(x,y),dist,load(x,y))
			elev(x,y) = elev(x,y) + dh
			load(nx,ny) = load(nx,ny) + load(x,y) - dh
C******************************************************************************
C	Flooded cells.
C******************************************************************************
		else if( flood(x,y) .eq. 2 ) then
			depth = elev(nx,ny) + SMALL - elev(x,y)
			if( load(x,y) .le. depth ) then
				elev(x,y) = elev(x,y)+load(x,y)
			else
				elev(x,y) = elev(x,y)+depth
				load(nx,ny)=load(nx,ny)+load(x,y)-depth
			endif
C******************************************************************************
C	Hillslope cells.
C******************************************************************************
		else
			elev(x,y) = elev(x,y) + load(x,y)
		endif
	end do
	return
	end
C******************************************************************************
C
C	Subroutine SORT_BY_ORDER	
C
C	This subroutine sorts the list of cell coordinates CELL() in 
C	ascending order of discharge. Since this version of the model
C	assumes a uniform precipitation field across the landscape,
C	discharge is directly proportional to drainage area. Sorting 
C	the cells by drainage area enables the FLUVIAL subroutine to 
C	process cells in upstream-to-downstream	order. The insertion
C	sort is adapted from D.E. Knuth (1973), volume 3.
C
C		Calls: (none)
C		Called by: FLUVIAL
C
C******************************************************************************
	subroutine sort_by_order()
	parameter (XBND=61,YBND=61,GRIDSIZE=3600)
	integer a,b,ax,ay,bx,by
	structure /cellcoord/
		integer x
		integer y
	end structure
	real*8 k
	record /cellcoord/ tmpcell
	record /cellcoord/ cell(GRIDSIZE)
	real*8 Q
	dimension Q(0:XBND,0:YBND)
	common /com3/ Q
	common /com8/ cell
	do a=2,GRIDSIZE
		b = a - 1
		ax = cell(a).x
		ay = cell(a).y
		bx = cell(b).x
		by = cell(b).y
		k = Q(ax,ay)
		tmpcell.x = ax
		tmpcell.y = ay
		do while( k .lt. Q(bx,by) .and. b .gt. 0 )
			cell(b+1).x = cell(b).x
			cell(b+1).y = cell(b).y
			b = b - 1
			if( b .gt. 0 ) then
				bx = cell(b).x
				by = cell(b).y
			endif
		end do
		cell(b+1).x = tmpcell.x
		cell(b+1).y = tmpcell.y
	end do
	return
	end
C******************************************************************************
C
C	Function BISECT	
C
C	This function uses the method of bisection to find the root of an
C	equation that balances carrying capacity (stream power) at a cell
C	with the total sediment load carried in the cell's channel segment.
C	The operative variable is a predicted change in cell elevation.
C	In other words, BISECT finds some elevation change Dz such that
C	the resulting stream power is just sufficient to carry sediment from
C	upstream, plus the change in load dictated by Dz (positive for
C	erosion, and negative for deposition). See TRANSPORT.
C
C		Calls: TRANSPORT
C		Called by: FLUVIAL
C
C******************************************************************************
	real*8 function bisect( a,b,ss,qq,dd,ld )
	implicit none
	real*8 a,b,ss,qq,dd,ld
	real*8 TOL
	parameter( TOL = 0.0001 )
	integer j
	real*8 delx,u,v,w,c
	real*8 transport
	u = transport( ss,qq,dd,ld,a )
	v = transport( ss,qq,dd,ld,b )
	if( u*v .ge. 0.0 ) then
		print *, 'Root must be bracketed for bisection in BISECT'
		call exit(0)
	endif
	delx = b - a
	do j=1,100
		c = 0.5*(a+b)
		w = transport( ss,qq,dd,ld,c )
		if( w .eq. 0.0 .or. delx .le. TOL ) goto 89
		delx = delx * 0.5
		if( w*u .lt. 0.0 ) then
			b = c
			v = w
		else
			a = c
			u = w
		endif
	end do
	print *, 'BISECT: Solution not found after 100 bisections'
	call exit( 0 )
89	bisect = c
	return
	end
C******************************************************************************
C
C	Function TRANSPORT 
C
C	This function finds the difference between stream power at a cell,
C	which is the same as the stream's carrying capacity, and the total
C	sediment carried after material has been eroded or deposited. Both
C	are functions of a predicted change in elevation DELZ. The correct
C	DELZ results when carrying capacity (stream power) exactly equals
C	the sediment load from upstream, plus new sediment entrained (or
C	minus sediment deposited) at the current cell.
C
C		Calls: OMEGA
C		Called by: BISECT
C
C******************************************************************************
	real*8 function transport( slp, dsch, dst, sedload, delz )
	implicit none	
	real*8 slp,dsch,dst,sedload,delz
	real*8 rhs,lhs
	real*8 omega
	rhs = omega( slp,dsch,dst,delz )
	lhs = sedload - delz
	transport = rhs - lhs
	return
	end
C******************************************************************************
C
C	Function OMEGA 
C
C	This function calculates transport capacity at a cell, given initial 
C	slope S, predicted change in elevation L, discharge DISCH, and 
C	length of channel segment DLX.
C
C		Calls: (none)
C		Called by: FLUVIAL, TRANSPORT
C
C******************************************************************************
	real*8 function omega( S, disch, dlx, L )
	implicit none
	real*8 S,disch,dlx,L
	real*8 kfl, mf, nf
	common /com13/ kfl, mf, nf
	omega = kfl*disch**mf*( S + L/dlx )**nf
	return
	end
C******************************************************************************
C
C	Subroutine HILLSLOPE
C
C	This subroutine calculates the diffusion of mass in two dimensions
C	using the Alternating-Direction Implicit (ADI) method of
C	Peaceman and Rachford. The boundaries are configured to represent
C	a landscape surface that slopes downward to the west, striking
C	north-south. The west boundary (x=0) is fixed at its
C	initial value. The east boundary follows the adjacent column of
C	cells. For the north and south boundaries, a 3-cell gradient is main-
C	tained such that dh/dy = 0. This introduces a small amount of 
C	inaccuracy near	these boundaries.
C
C	This subroutine can be easily modified to fix the east boundary 
C	as well, in order to simulate a rising orogen that strikes north-
C	south and is bounded on the east and west by regions of fixed 
C	elevation.
C
C		Calls: TRIDAG
C		Called by: MAIN
C
C******************************************************************************
	subroutine hillslope()
	implicit none
	integer XMAX,YMAX,XBND,YBND,GRIDSIZE,TILTBLOCK
	parameter(XMAX=60,YMAX=60,XBND=61,YBND=61,GRIDSIZE=3600)
	parameter(TILTBLOCK=3)
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	real*8 diff,alpha
	real*8 uprate,pivot_dist
	integer uptime,uptype,fault_loc
	real*8 uca,ucb,ucc,vca,vcb,vcc
	dimension uca(GRIDSIZE),ucb(GRIDSIZE),ucc(GRIDSIZE)
	dimension vca(GRIDSIZE),vcb(GRIDSIZE),vcc(GRIDSIZE)
	logical east_floats,ns_wrapped
	integer outlet_north,outlet_south
	common /com1/ elev
	common /com14/ diff,alpha
	common /com17/ uprate,uptime,uptype,fault_loc,pivot_dist
	common /com19/ uca,ucb,ucc,vca,vcb,vcc
	common /com20/ east_floats,ns_wrapped,outlet_north,outlet_south
	integer i,j,k 
	real*8 uu,rr,u,rhs
	dimension uu(GRIDSIZE),rr(GRIDSIZE)
	dimension u(0:XBND,0:YBND)
	dimension rhs(XMAX,YMAX)
	real*8 h,delh
C******************************************************************************
C	Compute right hand side (RHS) for odd traverse
C******************************************************************************
	do j=1,YMAX
		do i=1,XMAX
			rhs(i,j) = elev(i,j-1)+(1.0/alpha-2.0)*elev(i,j)
     &					+ elev(i,j+1)
		end do
		rhs(1,j) = rhs(1,j) + elev(0,j)
		rhs(XMAX,j) = rhs(XMAX,j) + elev(XBND,j)
	end do
C******************************************************************************
C	Translate RHS() to vector form and solve for odd traverse
C******************************************************************************
	k = 1
	do j=1,YMAX
		do i=1,XMAX
			rr(k) = rhs(i,j)
			k = k+1
		end do
	end do
	call tridag(uca,ucb,ucc,rr,uu)
C******************************************************************************
C	Translate odd solution vector to matrix form
C******************************************************************************
	k = 1
	do j=1,YMAX
		do i=1,XMAX
			u(i,j) = uu(k)
			k = k+1
		end do
	end do
C******************************************************************************
C	Adjust east boundary
C******************************************************************************
	if( east_floats ) then
		if( uptype .eq. TILTBLOCK ) then
			h=0.0
			delh=0.0
			do j=1,YMAX
				h=h+u(XMAX,j)
				delh=delh+0.5*(u(XMAX-2,j)-u(XMAX,j))
			end do
			h=(h-delh)/real(YMAX)
			if( h .lt. 0.0 ) h=0.0
			do j=0,YBND
				u(XBND,j)=h
			end do
		else
			do j=0,YBND
				u(XBND,j) = u(XMAX,j)
			end do
		endif
	endif
C******************************************************************************
C	Adjust north and south boundaries: if not "wrapped", they are
C	adjusted for no gradient across 3 cells.
C******************************************************************************
	if( ns_wrapped ) then
		do i=1,XMAX
			u(i,0) = u(i,YMAX)
			u(i,YBND) = u(i,1)
		end do
	else
		do i=1,XMAX
			u(i,1) = u(i,2)
			u(i,0) = u(i,1)
			u(i,YMAX) = u(i,YMAX-1)
			u(i,YBND) = u(i,YMAX)
		end do
	endif
C******************************************************************************
C	Adjust fixed west boundary and floating east boundary 
C******************************************************************************
	do j=0,YBND
		u(0,j) = elev(0,j)
		u(XBND,j) = u(XMAX,j)
	end do
C******************************************************************************
C	Adjust west boundary outside of outlet (if applicable)
C******************************************************************************
	if( outlet_south .gt. 0 ) then
		do j=0,outlet_south-1
			u(0,j) = u(1,j)
		end do
	endif
	if( outlet_north .lt. YBND ) then
		do j=outlet_north+1,YBND
			u(0,j) = u(1,j)
		end do
	endif
C******************************************************************************
C	Compute RHS for even traverse
C******************************************************************************
	do i=1,XMAX
		do j=1,YMAX
			rhs(i,j) = u(i-1,j)+(1.0/alpha-2.0)*u(i,j) +u(i+1,j)
		end do
		rhs(i,1) = rhs(i,1) + elev(i,0)
		rhs(i,YMAX) = rhs(i,YMAX) + elev(i,YBND)
	end do
C******************************************************************************
C	Translate RHS to vector and solve for even traverse
C******************************************************************************
	k = 1
	do i=1,XMAX
		do j=1,YMAX
			rr(k) = rhs(i,j)
			k = k+1
		end do
	end do
	call tridag( vca,vcb,vcc,rr,uu )
C******************************************************************************
C	Translate solution vector to matrix
C******************************************************************************
	k = 1
	do i=1,XMAX
		do j=1,YMAX
			u(i,j) = uu(k)
			k = k+1
		end do
	end do
C******************************************************************************
C	Apply new surface elevations to ELEV array, except for fixed west bdy
C	(and fixed east boundary, if applicable).	
C******************************************************************************
	if( east_floats ) then
		k = XBND
	else
		k = XMAX
	endif
	do i=1,k
		do j=0,YBND
			elev(i,j) = u(i,j)
		end do
	end do
C******************************************************************************
C	Adjust north and south boundaries: if not "wrapped", they are
C	adjusted for no gradient across 3 cells.
C******************************************************************************
	if( ns_wrapped ) then
		do i=1,XMAX
			elev(i,0) = elev(i,YMAX)
			elev(i,YBND) = elev(i,1)
		end do
	else
		do i=1,XMAX
			elev(i,1) = elev(i,2)
			elev(i,0) = elev(i,1)
			elev(i,YMAX) = elev(i,YMAX-1)
			elev(i,YBND) = elev(i,YMAX)
		end do
	endif
C******************************************************************************
C	Adjust west boundary outside of outlet (if applicable)
C******************************************************************************
	if( outlet_south .gt. 0 ) then
		do j=0,outlet_south-1
			elev(0,j) = elev(1,j)
		end do
	endif
	if( outlet_north .lt. YBND ) then
		do j=outlet_north+1,YBND
			elev(0,j) = elev(1,j)
		end do
	endif
	return
	end
C******************************************************************************
C
C	Subroutine TRIDAG	
C
C	This is a modified version of the tri-diagonal matrix solving
C	subroutine given in NUMERICAL RECIPES (Press et al, 1986).
C
C		Calls: (none)
C		Called by: HILLSLOPE
C
C******************************************************************************
	subroutine tridag(a,b,c,r,u)
	implicit none
	integer GRIDSIZE
	parameter(GRIDSIZE=3600)
	real*8 bet,gam,a,b,c,r,u
	dimension a(GRIDSIZE),b(GRIDSIZE),c(GRIDSIZE)
	dimension r(GRIDSIZE),u(GRIDSIZE),gam(GRIDSIZE)
	integer j
	bet = b(1)
	u(1) = r(1)/bet
	do j=2,GRIDSIZE
		gam(j) = c(j-1)/bet
		bet = b(j) - a(j)*gam(j)
		u(j) = (r(j)-a(j)*u(j-1))/bet
	end do
	do j=GRIDSIZE-1,1,-1
		u(j) = u(j)-gam(j+1)*u(j+1)
	end do
	return
	end
C******************************************************************************
C
C	Subroutine UPLIFT
C
C	This subroutine simulates spatially uniform uplift by incrementing
C	the elevation of all cells except for those along the western
C	boundary, which remains fixed.
C
C		Calls: (none)
C		Called by: MAIN
C
C******************************************************************************
	subroutine uplift()
	implicit none
	integer XBND,YBND
	parameter (XBND=61,YBND=61)
	integer UNIFORM,TILTBLOCK 
	parameter (UNIFORM=1,TILTBLOCK=3)
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	real*8 uprate,pivot_dist
	integer uptime,uptype,fault_loc
	real*8 dx,dxx
	common /com1/ elev
	common /com17/ uprate,uptime,uptype,fault_loc,pivot_dist
	common /com18/ dx,dxx 
	integer i,j
	real*8 d
	if( uptype .eq. UNIFORM ) then
		do i=fault_loc,XBND
			do j=0,YBND
				elev(i,j) = elev(i,j) + uprate
			end do
		end do
	else if( uptype .eq. TILTBLOCK ) then
		do i=fault_loc,XBND
			do j=0,YBND
				d = pivot_dist*1000.0 - (i-fault_loc)*dx
				elev(i,j)=elev(i,j)+(d*uprate)/(pivot_dist*1000)
			end do
		end do
	endif
	return
	end
C******************************************************************************
C
C	Subroutine WRITE_ELEV_DATA	
C
C	This subroutine records the elevation of each cell (including
C	boundary cells) in a file called '<name>.elev'.
C
C		Calls: (none)
C		Called by: MAIN
C
C******************************************************************************
	subroutine write_elev_data( t )
	integer t
	integer i,j
	parameter (XBND=61,YBND=61)
	real*8 elev
	dimension elev(0:XBND,0:YBND)
	common /com1/ elev
90	format(I6)	
91	format(F8.3)	
	write( 1, 90 ) t
	do j=0,YBND
		write( 1, 91 ) (elev(i,j), i=0,XBND) 
	end do
	return
	end
C******************************************************************************
C
C	Subroutine WRITE_NET_DATA
C
C	This subroutine writes drainage network data to a file <name>.net.
C	Network data consists of the channel/hillslope status (1 or 0) and
C	the flow direction (x and y coordinates) for each cell.
C
C		Calls: (none)
C		Called by: MAIN
C
C******************************************************************************
	subroutine write_net_data( t )
	implicit none
	integer t
	integer i,j
	integer XBND,YBND
	parameter (XBND=61,YBND=61)
	logical chan
	real*8 elev
	integer nbrx,nbry
	byte flood
	dimension nbrx(0:XBND,0:YBND)
	dimension nbry(0:XBND,0:YBND)
	dimension chan(0:XBND,0:YBND)
	dimension elev(0:XBND,0:YBND)
	dimension flood(0:XBND,0:YBND)
	common /com1/ elev
	common /com5/ nbrx,nbry 
	common /com6/ flood 
	common /com7/ chan 
	byte ch
99	format(I1,1X,I3,1X,I3)	
	write( 2, * ) t
	do j=0,YBND
		do i=0,XBND
			if( chan(i,J) ) then
				ch = 1
			else if( flood(i,j) .gt. 0 ) then
				ch = 2
			else
				ch = 0
			endif
			write( 2, 99 ) ch, nbrx(i,j), nbry(i,j)
		end do
	end do
	return
	end
C******************************************************************************
C
C	Function FIND_SLOPE	
C
C	This function calculates the slope between two points, given the
C	elevation of each point and the distance between them. Slope is
C	considered to be positive downward.
C
C		Calls: (none)
C		Called by: FIND_NEXT_NODE, FIND_DRY_NEIGHBOR
C
C******************************************************************************
	real*8 function find_slope( h1, h2, delx )
	real*8 h1,h2,delx
	find_slope = (h1-h2) / delx
	return
	end
C******************************************************************************
C
C	Function IS_BOUNDARY 
C
C	This function determines whether a given cell is a boundary cell,
C	returning a value of TRUE if it is and FALSE if it is not.
C
C******************************************************************************
	logical function is_boundary( cx, cy )
	integer cx,cy
	integer XBND, YBND 
	parameter (XBND=61,YBND=61)
	if(cx .eq. XBND .or. cx .eq. 0 ) then
		is_boundary = .TRUE.
	else if( cy .eq. 0 .or. cy .eq. YBND ) then
		is_boundary = .TRUE.
	else
		is_boundary = .FALSE.
	endif
	return
	end
