'... Decon.bas, Mark Brandon & Frank Pazzaglia, January 25, 1994, v. 1.1
'============================================================================
'    This program uses the erosion history provided by sediment yield data
' from offshore basins to calculate the tectonic history of a mountain belt.
' The basic assumption is that the average rate of mechanical erosion is
' proportional to the average elevation of the mountain belt. Based on this
' assumption, the mechanical erosion rate Ei at the ith time step can be
' related to  a set of incremental stimuli Bj and a response matrix Gij,
'
'                     Ei = Zi'/Kd = Gij Bj.
'
' The indices i and j represent elapse time (e.g. time_step * (i + 1/2)).
' The constant Kd indicates the ratio of mean erosion rate to mean basin
' elevation. The stimuli Bj = Vj - EChemdot - Ldotj where Vj and Ldotj are
' the source flux (tectonic signal) and rate of sea level change at time
' step j, and Echemdot is a constant rate of chemical erosion.
'    Zi' = Zi - Li is the relative elevation above Li contemporaneous sea
' level. The columns of Gij are response vectors describing the mechanical
' erosion rate at time i caused by an incremental stimulus Bj at time j.
'    The program uses a known set of mechanical erosion rates Ei to
' invert for parameters of the governing equations. There are two options:
' (1) estimate Bj and Vj assuming a constant Kd, and (2) estimate variable
' Kd's assuming that Vj = 0. The first attributes variations
' in mechanical erosion to variations in source flux, whereas the
' second assumes that variation in erosion are related to changes in
' erodability caused by climate or changing rock type.
'    Both inversions involve the solution of an evenly-determined system
' of equations. The first option uses a damped least-squares algorithm
' to simultaneously minimize the data residuals and the second
' derivatives (smoothness) of the bj time series.
'==========================================================================
DECLARE SUB DiskOpen (E$, DF$, FileNum%)
DECLARE SUB DiskOut (E$, F$, FileNum%)
DECLARE FUNCTION Func! (EdotKnown!, Wdot!, Zhi!, Dt!, KdGuess!)
DECLARE FUNCTION RchiSq! (Epsilon!)
DECLARE SUB Lubksb (A!(), N%, Indx%(), B!())
DECLARE SUB Ludcmp (A!(), N%, Indx%(), D!)
DECLARE FUNCTION MatMultS% (Alpha!(), Beta!(), Gamma!())
DECLARE FUNCTION MatTrans% (Alpha!(), Beta!())
DECLARE SUB RTbis (EdotKnown!, Wdot!, Zhi!, Dt!, Kd!)
DIM Head$(5), Temp!(200)
CONST C = .21   '... compensation ratio for emergent topography
'==========================================================================
'... Print opening screen
CLS
Head$(1) = "DECON PROGRAM v. 1.1  MTB & FJP 1/25/96"
Q$ = STRING$((79 - LEN(Head$(1))) \ 2, "=")
Head$(1) = Q$ + Head$(1) + Q$
PRINT Head$(1)
PRINT "This program uses the record of erosion provided by sediment yield data "
PRINT "to calculate the mean source flux for a mountain belt assuming a "
PRINT "elevation-dependent erosion law where the erosion rate is proportional "
PRINT "to elevation above contemporaneous sea level."
PRINT
PRINT "Format for input data :"
PRINT "LINE 1: Description of data set (less than 80 characters)"
PRINT "LINE 2: Time_step_(myr)"
PRINT "LINE 3: Estimated_rel_std_error_(%)_erosion_data"
PRINT "LINE 4: Chemical_erosion_rate(m/myr)     <assumed constant>"
PRINT "LINE 5: Present_mean_elevation_relative_to_present_short-term_SL"
PRINT "LINE 6: L0, B1, B2, B3  <3rd degree polynomial coef. for SL curve (m)> "
PRINT "LINES 7 to N+6: Mechanical_erosion_rate_(m/myr) <from youngest to oldest>"
PRINT
PRINT "The output will have the following format: "
PRINT "LINE 1-3: Header information"
PRINT "LINE 3-11: Description of solution"
PRINT "LINE 12: column labels for output data "
PRINT "LINE 13 to N+12: columns of output data as described in line 12"
PRINT
PRINT "Note that CONTROL-BREAK or CONTROL-C will abort the program at any time."
PRINT
PRINT "Press any key to continue:> ? ";
Q$ = INPUT$(1)
Start:
CLS
CLEAR
PRINT "--------------------------- OPEN INPUT FILE ---------------------------"
CALL DiskOpen("", FileName$, 1)
Head$(2) = "FILE: " + FileName$ + "  RUN TIME: " + DATE$ + "/" + TIME$
'... Read header information from input file
LINE INPUT #1, Head$(3)
INPUT #1, Dt
INPUT #1, RSEedot
RSEedot = RSEedot / 100
INPUT #1, EChemdot
INPUT #1, Z0
INPUT #1, L0, B1, B2, B3
Z0 = Z0 - L0
Num% = 0
DO WHILE NOT EOF(1)
	Num% = Num% + 1
	INPUT #1, Temp(Num%)
LOOP
CLOSE #1

PRINT
PRINT "--------------------------- OPEN OUTPUT FILE ---------------------------"
CALL DiskOut("", FileOut$, 2)
PRINT
PRINT
DO
	PRINT "Select type of solution:"
	PRINT "     1) Solve for variable source flux and constant Kd"
	PRINT "     2) Solve for variable Kd and zero source flux"
	PRINT "Enter option (1,2) >? ";
	Q$ = INPUT$(1)
	PRINT Q$
	PRINT
	Opt% = VAL(Q$)
LOOP WHILE Opt% <> 1 AND Opt% <> 2

REDIM Edot!(Num%), E!(Num%), EdotSE!(Num%), CovEdot!(Num%, Num%)
REDIM L!(Num%), Ldot!(Num%), Z!(Num%), ZSE!(Num%), Zdot!(Num%)
REDIM V!(Num%), Vdot!(Num%), CovVdot!(Num%, Num%), VdotSE!(Num%)
REDIM U!(Num%), Udot!(Num%), Kdv!(Num%), KdvSE!(Num%)
REDIM Gamma!(Num%), G!(Num%, Num%), Gtrans!(Num%, Num%)
REDIM We!(Num%, Num%), Wm!(Num%, Num%), EdotCalc!(Num%)
REDIM Dm!(Num% - 3, Num%), Dmtrans!(Num%, Num% - 3)
REDIM LHS!(Num%, Num%), M!(Num%, Num%), Mtrans!(Num%, Num%)
REDIM RHS!(Num%, Num%), Indx%(Num%), B!(Num%)

'... Invert order so that arrays go from oldest to youngest
FOR I% = 1 TO Num%
	Edot(I%) = Temp(Num% - I% + 1)
	Age = (Num% - I% + .5) * Dt    '... Age in Ma
	L(I%) = B1 * Age + B2 * Age ^ 2 + B3 * Age ^ 3
	Ldot(I%) = B1 + 2 * B2 * Age + 3 * B3 * Age ^ 2
	CovEdot(I%, I%) = (Edot(I%) * RSEedot) ^ 2
	EdotSE(I%) = Edot(I%) * RSEedot
NEXT I%

'... Report information about input file
CLS
PRINT "Header recorder for data file:"
PRINT Head$(3)
PRINT
PRINT "Time Step (myr): "; Dt
PRINT "Relative standard error for erosion rate (%): "; RSEedot * 100
PRINT "Chemical erosion rate (assumed constant) (m/myr): "; EChemdot
PRINT "Compensation ratio (dZ/dT) for emergent topography: "; C
PRINT "Present mean elevation relative to modern short-term SL (m): "; Z0 + L0
PRINT "Present height of long-term SL (m): "; L0
PRINT "Present mean elevation relative to modern long-term SL (m): "; Z0
PRINT "Third-degree polynomial for sealevel (m) wrt present long-term SL:"
FF$ = "&##.###^^^^ & +#.###^^^^& +#.###^^^^"
PRINT USING FF$; "L ="; B1; "* t"; B2; "* t ^ 2"; B3; "* t ^ 3"
PRINT
PRINT "Press any key to continue:> ? ";
Q$ = INPUT$(1)
CLS
'... Solve erosion equations
SELECT CASE Opt%
	CASE 1    '... Solve for source flux and constant Kd
		PRINT
		PRINT "Estimate Kd using present relative elevation and average erosion rate for "
		PRINT "youngest time step."
		PRINT "First guess assumes no change in relative elevation  during the youngest "
		PRINT "time step."
		Kd1 = Edot(Num%) / Z0
		PRINT USING "&#.#####"; "First guess for Kd (1/myr): "; Kd1
		PRINT
		PRINT "Second guess assumes that Vdot = 0 for youngest time step. "
		Bdot = -Ldot(Num%) - EChemdot
		CALL RTbis(Edot(Num%), Bdot, Z0, Dt, Kd2)
		PRINT USING "&#.#####"; "Second guess for Kd (1/myr): "; Kd2
		PRINT "Implied elevation at youngest midpoint (m): "; Edot(Num%) / Kd2
		PRINT
		PRINT "Enter Kd (1/myr) <default:"; Kd1; ">: ";
		INPUT Q$
		IF Q$ = "" THEN Kd = Kd1 ELSE Kd = VAL(Q$)
		
		'... Construct response vector Gamma(i) of duration Num%
		G0 = (1 - EXP(-C * Kd * Dt)) / (C * Kd * Dt)
		Gamma(1) = 1 - G0
		FOR I% = 2 TO Num%
			Tau = Dt * (I% - 1)
			Gamma(I%) = G0 * (EXP(-C * Kd * (Tau - Dt)) - EXP(-C * Kd * Tau))
		NEXT I%

		'... Construct response matrix G(jk) where Edot(j) = G(jk) Bdot(k)
		' G(jk) contains the response vectors as column vectors with each vector
		' incremented downward one row for each column to the left.
		FOR K% = 1 TO Num%
			I% = 1
			FOR J% = 1 TO Num%
				IF J% < K% THEN
					G(J%, K%) = 0
				ELSE
					G(J%, K%) = Gamma(I%)
					I% = I% + 1
				END IF
			NEXT J%
		NEXT K%

		'... Constuct We() for weighting data
		FOR I% = 1 TO Num%
			FOR J% = 1 TO Num%
				IF I% = J% THEN
					We(I%, I%) = 1 / CovEdot(I%, I%)
				ELSE
					We(I%, J%) = 0
				END IF
			NEXT J%
		NEXT I%

		'... Find an optimal solution where RchiSq = 1 and Epsilon >0
		PRINT
		PRINT ">>> Search for optimal solution with smooth parameters <<<"
		Jmax = 40
		ACC = .0001
		Epsilonlw = 0
		Epsilonhi = .1
		Diff = RchiSq(Epsilonlw) - 1
		Diffmid = RchiSq(Epsilonhi) - 1
		IF Diff * Diffmid >= 0 THEN
			PRINT "  *** Root must be bracketed before bisection is started ***"
			STOP
		END IF
		IF Diff < 0 THEN
			Epsilon = Epsilonlw
			DEpsilon = Epsilonhi - Epsilonlw
		ELSE
			Epsilon = Epsilonhi
			DEpsilon = Epsilonlw - Epsilonhi
		END IF
		FOR J% = 1 TO Jmax
			DEpsilon = DEpsilon * .5
			EpsilonMid = Epsilon + DEpsilon
			Diffmid = RchiSq(EpsilonMid) - 1
			IF Diffmid <= 0 THEN Epsilon = EpsilonMid
			IF ABS(DEpsilon) < ACC OR Diffmid = 0 THEN EXIT FOR
			IF J% = Jmax THEN
				PRINT " *** Too many bisections in finding root to correction equation ***"
				STOP
			END IF
			FF$ = "& ##.##^^^^ ####.#######"
			PRINT USING FF$; "Epsilon and reduced chi-square: "; EpsilonMid; Diffmid + 1
		NEXT J%
		ReducedChiSq = Diffmid + 1
		Epsilon = EpsilonMid

		'... Calculate standard errors for calculated Vdot's assuming
		' that the main source of error is with the Edot's.
		' Method is from Meinke (1989, p. 58).

		'... Calculate inverse used for covariance calcuation
		FOR J% = 1 TO Num%
			FOR I% = 1 TO Num%
				B(I%) = RHS(I%, J%)
			NEXT I%
			CALL Lubksb(LHS(), Num%, Indx%(), B())
			FOR I% = 1 TO Num%
				M(I%, J%) = B(I%)
			NEXT I%
		NEXT J%
		ErrCode = MatTrans%(M(), Mtrans())
		ErrCode% = MatMultS%(M(), CovEdot(), LHS())
		ErrCode% = MatMultS%(LHS(), Mtrans(), CovVdot())

		'... Calculate source flux, surface uplift, and rock uplift
		FOR I% = 1 TO Num%
			Vdot(I%) = Vdot(I%) + Ldot(I%) / C + EChemdot
			VdotSE(I%) = SQR(CovVdot(I%, I%))
			Kdv(I%) = Kd
			KdvSE(I%) = 0
		NEXT I%

	CASE 2    '... Solve for variable Kd and zero source flux
		Zhi = Z0
		ZhiSE = 0
		FOR I% = Num% TO 1 STEP -1
			'... Calculate Kd for middle of timestep, starting with youngest
			' timestep and moving backward in time
			Bdot = -Ldot(I%) - EChemdot
			CALL RTbis(Edot(I%), Bdot, Zhi, Dt, Kdv(I%))
			'... Calculate elevation above contemporaneous sealevel at old limit
			' of timestep
			Zlw = Bdot / Kdv(I%) + (Zhi - Bdot / Kdv(I%)) * EXP(C * Kdv(I%) * Dt)
			EdotCalc(I%) = Bdot + (Kdv(I%) * Zhi - Bdot) * EXP(C * Kdv(I%) * Dt / 2)
			DEdotiDKi = (Zhi + C * Zhi * Kdv(I%) ^ 2 * Dt / 2 - C * Bdot * Kdv(I%) * Dt / 2)
			DEdotiDKi = DEdotiDKi * EXP(C * Kdv(I%) * Dt / 2)
			DZlwDKi = Bdot / Kdv(I%) ^ 2 + Zhi * C * Kdv(I%) * Dt - C * Bdot * Dt
			DZlwDKi = DZlwDKi * EXP(C * Kdv(I%) * Dt) - Bdot / Kdv(I%) ^ 2
			DZhiDKi = Bdot / Kdv(I%) ^ 2 - Zhi * C * Kdv(I%) * Dt + C * Bdot * Dt
			DZhiDKi = DZlwDKi * EXP(-C * Kdv(I%) * Dt) - Bdot / Kdv(I%) ^ 2
			DZlwDZhi = EXP(C * Kdv(I%) * Dt)
			DZlwDEdoti = DZlwDKi / DEdotiDKi
			KdvSE(I%) = SQR((ZhiSE2 / DZhiDKi ^ 2) + (EdotSE(I%) / DEdotiDKi) ^ 2)
			ZhiSE2 = (DZlwDKi * KdvSE(I%)) ^ 2 + (DZlwDZhi ^ 2 * ZhiSE2)
			Zhi = Zlw
		NEXT I%
		ReducedChiSq = 0
		Epsilon = 0
		FOR I% = 1 TO Num%
			Vdot(I%) = 0
			VdotSE(I%) = 0
		NEXT I%
END SELECT

'... Calculate Z(), ZSE(), Zdot(), Udot()
FOR I% = 1 TO Num%
	Age = (Num% - I% + .5) * Dt      '... Age in Ma
	Z(I%) = Edot(I%) / Kdv(I%)
	SELECT CASE Opt%
		CASE 1  '... variable Vdot, constant Kd
			ZSE(I%) = Z(I%) * EdotSE(I%) / Edot(I%)
		CASE 2  '... variable Kd, zero Vdot
			ZSE(I%) = Z(I%) * KdvSE(I%) / Kdv(I%)
	END SELECT
	Z(I%) = Z(I%) + L(I%)
	Zdot(I%) = C * (Vdot(I%) - Edot(I%) - EChemdot)
	Udot(I%) = (1 - C) * (Edot(I%) + EChemdot) + C * Vdot(I%)
	'PRINT USING "#####.#  "; Age; Vdot(I%); Edot(I%); EdotCalc(I%); : INPUT Q$
NEXT I%

'... Calculate total eroded section,
' and integrate Edot, Vdot, and Udot
FOR I% = 1 TO Num%
	E(I%) = 0
	'... From beginning to peneultimate time-step
	FOR J% = 1 TO I% - 1
		E(I%) = E(I%) + Edot(J%) * Dt
	NEXT J%
	'... Last half-step
	E(I%) = E(I%) + Edot(I%) * Dt / 2
NEXT I%
FOR I% = 1 TO Num%
	V(I%) = 0
	U(I%) = 0
	'... From beginning to peneultimate time-step
	' Skip Vdot(1) and Udot(1) because they are proxies for initial conditions
	FOR J% = 2 TO I% - 1
		V(I%) = V(I%) + Vdot(J%) * Dt
		U(I%) = U(I%) + Udot(J%) * Dt
	NEXT J%
	IF I% > 1 THEN
		'... Last half-step
		V(I%) = V(I%) + Vdot(I%) * Dt / 2
		U(I%) = U(I%) + Udot(I%) * Dt / 2
	END IF
NEXT I%

'... Write output file
FOR I% = 1 TO 3
	PRINT #2, Head$(I%)
NEXT I%
PRINT #2, "Relative standard error for erosion rate (%): "; RSEedot * 100
PRINT #2, "Chemical erosion rate (assumed constant) (m/myr): "; EChemdot
PRINT #2, "Compensation ratio (dZ/dT) for emergent topography: "; C
PRINT #2, "Present mean elevation relative to modern short-term SL (m): "; Z0 + L0
PRINT #2, "Present height of long-term SL (m): "; L0
PRINT #2, "Present mean elevation relative to modern long-term SL (m): "; Z0
PRINT #2, "Third-degree polynomial for sealevel (m) wrt present long-term SL:"
PRINT #2, "L = "; B1; "*t + "; B2; "*t^2 + "; B3; "*t^3"
PRINT #2, "Epsilon and reduced chi-square: "; Epsilon, ReducedChiSq
FF$ = "\         \"
PRINT #2, USING FF$; "Age_Ma"; "Time_myr"; "Edot_m/myr"; "Edot-SE"; "Edot+SE"; "EdotCalc"; "Ldot_m/myr";
PRINT #2, USING FF$; "Vdot_m/myr"; "Vdot-SE"; "Vdot+SE"; "Udot_m/myr"; "Zdot_m/myr";
PRINT #2, USING FF$; "Kdv_myr^-1"; "Kdv-SE"; "Kdv+SE"; "E_m"; "L_m"; "V_m"; "U_m"; "Z_m"; "Z-SE"; "Z+SE"
'Age,Time, Edot, Edot-SE, Edot+SE, Ldot
'Vdot, Vdot-SE, Vdot+SE, Udot, Zdot
'Kdv, Kdv-SE, Kdv+SE, E, L, V, U, Z, Z-SE, Z+SE
FOR I% = 1 TO Num%
	'... Middle of current timestep
	Age = (Num% - I% + .5) * Dt      '... Age in Ma
	time = (I% - .5) * Dt            '... Time since start of record in myr
	FF$ = "####.#     ####.#     ##.###^^^^ ##.###^^^^ ##.###^^^^ ##.###^^^^ "
	PRINT #2, USING FF$; Age; time; Edot(I%); Edot(I%) - EdotSE(I%); Edot(I%) + EdotSE(I%); EdotCalc(I%);
	IF I% = 1 THEN
		CLS
		FF$ = "##.###^^^^ \        \ \        \ \        \ "
		FF1$ = "    --    "
		PRINT #2, USING FF$; Ldot(I%); FF1$; FF1$; FF1$;
		FF$ = "\        \ \        \ ##.###^^^^ ##.###^^^^ ##.###^^^^ "
		PRINT #2, USING FF$; FF1$; FF1$; Kdv(I%); Kdv(I%) - KdvSE(I%); Kdv(I%) + KdvSE(I%);
		FF$ = "##.###^^^^ ##.###^^^^ ##.###^^^^ ##.###^^^^ ##.###^^^^ ##.###^^^^ "
		PRINT #2, USING FF$; E(I%); L(I%); V(I%); U(I%); Z(I%); Z(I%) - ZSE(I%); Z(I%) + ZSE(I%)
	ELSE
		FF$ = "##.###^^^^ ##.###^^^^ ##.###^^^^ ##.###^^^^ ##.###^^^^ ##.###^^^^ "
		PRINT #2, USING FF$; Ldot(I%); Vdot(I%); Vdot(I%) - VdotSE(I%); Vdot(I%) + VdotSE(I%);
		PRINT #2, USING FF$; Udot(I%); Zdot(I%); Kdv(I%); Kdv(I%) - KdvSE(I%); Kdv(I%) + KdvSE(I%);
		PRINT #2, USING FF$; E(I%); L(I%); V(I%); U(I%); Z(I%); Z(I%) - ZSE(I%); Z(I%) + ZSE(I%)
	END IF
NEXT I%
CLOSE #2
PRINT
PRINT "Run the program again (y,n; default: n): > ? ";
Q$ = UCASE$(INPUT$(1))
IF Q$ = CHR$(13) THEN Q$ = "N"
PRINT Q$
IF Q$ = "Y" THEN GOTO Start
END

SUB DiskOpen (E$, DF$, FileNum%)
'===========================================================================
'... This routine helps open a disk file. The input arguement E$ allows
' the user to specify a listing of files with a specific file extension.
' Subroutine returns DF$ with full path and file name.
' 12/28/93: Parsing routine for path fixed and other minor changes.
' 11/21/94: Added Error routines to fix problems with empty directories
'===========================================================================
ON LOCAL ERROR GOTO ErrHandler1
IF E$ = "" THEN EXT$ = "*" ELSE EXT$ = E$
INDATA:
PRINT
PRINT "Current directory: "; CURDIR$
DO
	INPUT "Enter drive and path of directory"; DD$
	DD$ = UCASE$(RTRIM$(LTRIM$(DD$)))
	IF DD$ = "" OR RIGHT$(DD$, 1) = ":" THEN DD$ = CURDIR$(DD$)
	IF RIGHT$(DD$, 1) <> "\" THEN DD$ = DD$ + "\"
	Q$ = "SELECTED DIRECTORY: " + DD$ + "*." + EXT$
	L% = INT((80 - LEN(Q$)) / 2) - 1
	PRINT STRING$(L%, "-") + " " + Q$ + " " + STRING$(L%, "-")
	IF LEN(DIR$(DD$ + "*." + EXT$)) <> 0 THEN EXIT DO
	PRINT
	PRINT "    >>> Selected directory or selected files are not present: try again <<<   "
	PRINT
	PRINT
LOOP
Q$ = DIR$(DD$ + "*." + EXT$)
DO
	L% = LEN(Q$)
	Q$ = Q$ + SPACE$(20 - L%)
	PRINT Q$;
	IF POS(1) = 80 THEN PRINT
	Q$ = DIR$
LOOP WHILE Q$ <> ""
PRINT
INPUT ">>> File name"; DF$
PRINT
DF$ = UCASE$(DF$)
IF (LEN(DIR$(DD$ + DF$)) = 0) OR DF$ = "" THEN
	PRINT "    >>> File does not exist. Press any key to continue> ? ;"
	Q$ = INPUT$(1)
	CLS
	GOTO INDATA
END IF

'... Construct final output strings and change to new directory
IF RIGHT$(DD$, 1) = "\" THEN
	IF DD$ <> "\" AND RIGHT$(DD$, 2) <> ":\" THEN
		DD$ = LEFT$(DD$, LEN(DD$) - 1)
	END IF
END IF
IF MID$(DD$, 2, 1) = ":" THEN CHDRIVE DD$
CHDIR DD$
DD$ = CURDIR$
IF RIGHT$(CURDIR$, 1) <> "\" THEN DD$ = DD$ + "\"
DF$ = DD$ + DF$
CLOSE FileNum%
OPEN "INPUT", FileNum%, DF$
EXIT SUB

ErrHandler1:
Message$ = "    >>> "
SELECT CASE ERR
	CASE 52
		Message$ = Message$ + "Bad file name or number."
	CASE 53
		Message$ = Message$ + "File not found."
	CASE 57
		Message$ = Message$ + "Unformatted disk."
	CASE 64
		Message$ = Message$ + "Incorrect drive name."
	CASE 68
		Message$ = Message$ + "Device unavailable."
	CASE 71
		Message$ = Message$ + "Drive is not ready."
	CASE ELSE
		PRINT Message$ + "Unexpected FATAL error. Program halted. <<<"
		STOP
END SELECT
PRINT
PRINT Message$ + " Press any key to continue> ? ";
Q$ = INPUT$(1)
PRINT
RESUME INDATA
END SUB

SUB DiskOut (E$, F$, FileNum%)
'===========================================================================
'... Open disk file for output. Call includes specification (EXT$) for file
' extension. Subroutine returns F$ with full path and file name.
' 12/28/93: Parsing routine for path fixed and other minor changes.
' 11/21/94: Added Error routines to fix problems with empty directories
'===========================================================================
ON LOCAL ERROR GOTO ErrHandler2
IF E$ = "" THEN EXT$ = "*" ELSE EXT$ = E$
OUTDATA:
PRINT
PRINT "Current directory = "; CURDIR$
DO
	INPUT "Enter drive and path of directory"; DD$
	DD$ = UCASE$(RTRIM$(LTRIM$(DD$)))
	IF DD$ = "" OR RIGHT$(DD$, 1) = ":" THEN DD$ = CURDIR$(DD$)
	IF RIGHT$(DD$, 1) <> "\" THEN DD$ = DD$ + "\"
	Q$ = "SELECTED DIRECTORY: " + DD$ + "*." + EXT$
	L% = INT((80 - LEN(Q$)) / 2) - 1
	PRINT STRING$(L%, "-") + " " + Q$ + " " + STRING$(L%, "-")
	IF (LEN(DIR$(DD$ + "*.*")) <> 0) THEN EXIT DO
	PRINT
	PRINT "    >>> Selected directory or selected files are not present: try again <<<   "
	PRINT
LOOP
Q$ = DIR$(DD$ + "*." + EXT$)
DO
	L% = LEN(Q$)
	Q$ = Q$ + SPACE$(20 - L%)
	PRINT Q$;
	IF POS(1) = 80 THEN PRINT
	Q$ = DIR$
LOOP WHILE Q$ <> ""
PRINT
INPUT ">>> File name"; DF$
DF$ = UCASE$(DF$)
IF DF$ = "" THEN
	PRINT
	PRINT "   >>> Blank file name. Press any key to continue> ? "
	Q$ = INPUT$(1)
	CLS
	GOTO OUTDATA
END IF
IF (LEN(DIR$(DD$ + DF$)) <> 0) THEN
	PRINT
	DO
		PRINT "   >>> File already exists. Do you wish to overwrite (Y,N)> ? ";
		Q$ = UCASE$(INPUT$(1))
		PRINT Q$
	LOOP WHILE Q$ <> "Y" AND Q$ <> "N"
	IF Q$ = "N" THEN CLS : GOTO OUTDATA
END IF

'... Construct final output strings and change to new directory
IF RIGHT$(DD$, 1) = "\" THEN
	IF DD$ <> "\" AND RIGHT$(DD$, 2) <> ":\" THEN
		DD$ = LEFT$(DD$, LEN(DD$) - 1)
	END IF
END IF
IF MID$(DD$, 2, 1) = ":" THEN CHDRIVE DD$
CHDIR DD$
DD$ = CURDIR$
IF RIGHT$(CURDIR$, 1) <> "\" THEN DD$ = DD$ + "\"
F$ = DD$ + DF$
CLOSE FileNum%
OPEN "OUTPUT", FileNum%, DF$
EXIT SUB

ErrHandler2:
Message$ = "    >>> "
SELECT CASE ERR
	CASE 52
		Message$ = Message$ + "Bad file name or number."
	CASE 53
		Message$ = Message$ + "File not found."
	CASE 57
		Message$ = Message$ + "Unformatted disk."
	CASE 64
		Message$ = Message$ + "Incorrect drive name."
	CASE 68
		Message$ = Message$ + "Device unavailable."
	CASE 71
		Message$ = Message$ + "Drive is not ready."
	CASE ELSE
		PRINT Message$ + "Unexpected FATAL error. Program halted. <<<"
		STOP
END SELECT
PRINT
PRINT Message$ + " Press any key to continue> ? ";
Q$ = INPUT$(1)
PRINT
RESUME OUTDATA
END SUB

FUNCTION Func! (EdotKnown!, Bdot!, Zhi!, Dt!, KdGuess!)
'======================================================================
'... Used to house the function from which a root will be found for
' the RTbis subroutine. The objective is to find a root for a
' transcendental function, which is a function Edot=f(Kd) that cannot be
' inverted algebraically to the form Kd=g(Edot).
' The Func routine here uses the forward function and KdGuess to
' calculate Edotguess. Func is returned with the difference between
' Edotguess and Edotknown. When the correct root is found, Func = 0.
' C = compensation ratio for emergent topography, defined as constant
'======================================================================
'Edotguess = f(Kdguess)
EdotGuess = Bdot + (KdGuess * Zhi - Bdot) * EXP(C * KdGuess * Dt / 2)
Func = EdotKnown - EdotGuess
END FUNCTION

SUB Lubksb (A!(), N%, Indx%(), B!())
'=========================================================================
'... Solves the set of N linear equations AX = B. Here A is input, not as
' the matrix A(), but rather as its LU decomposition, as determined by
' the routine LUDCMP. Indx%() is input as the permutation vector returned
' by LUDCMP. B(1:N%) is input as the right-hand side vector B, and
' is returned with the solution vector X. A(), N%, and Indx%() are not
' modified by this routine and can be left in place for successive calls
' with different right-hand sides B(). This routine takes into account
' the possiblity that B() will begin with many zero elements, so that
' it is efficient for use in matrix inversion.
'=========================================================================
II% = 0
FOR I% = 1 TO N%
  LL% = Indx%(I%)
  Sum = B(LL%)
  B(LL%) = B(I%)
  IF II% <> 0 THEN
	FOR J% = II% TO I% - 1
	  Sum = Sum - A(I%, J%) * B(J%)
	NEXT J%
  ELSEIF Sum <> 0! THEN
	II% = I%
  END IF
  B(I%) = Sum
NEXT I%
FOR I% = N% TO 1 STEP -1
  Sum = B(I%)
  FOR J% = I% + 1 TO N%
	Sum = Sum - A(I%, J%) * B(J%)
  NEXT J%
  B(I%) = Sum / A(I%, I%)
NEXT I%
END SUB

SUB Ludcmp (A!(), N%, Indx%(), D!)
'=========================================================================
'... Given a matrix A(1:N%, 1:N%), this routine replaces it by the LU
' decomposition of a row-wise permutation of itself. Indx%(1:N%) is
' output as a vector containing a record of the row permutation effected
' by the partial pivoting. D is output as +/-1, depending on whether the
' number of row interchanges was even or odd, respectively. This routine
' is used in combination with LUBKSB to solve linear equations or to
' invert a matrix.
'=========================================================================
TINY = 1E-20
DIM VV(N%)
D = 1!
FOR I% = 1 TO N%
  AAMAX = 0!
  FOR J% = 1 TO N%
	IF ABS(A(I%, J%)) > AAMAX THEN AAMAX = ABS(A(I%, J%))
  NEXT J%
  IF AAMAX = 0! THEN PRINT ">>> Singular matrix <<<": STOP
  VV(I%) = 1! / AAMAX
NEXT I%
FOR J% = 1 TO N%
  FOR I% = 1 TO J% - 1
	Sum = A(I%, J%)
	FOR K% = 1 TO I% - 1
	  Sum = Sum - A(I%, K%) * A(K%, J%)
	NEXT K%
	A(I%, J%) = Sum
  NEXT I%
  AAMAX = 0!
  FOR I% = J% TO N%
	Sum = A(I%, J%)
	FOR K% = 1 TO J% - 1
	  Sum = Sum - A(I%, K%) * A(K%, J%)
	NEXT K%
	A(I%, J%) = Sum
	Dum = VV(I%) * ABS(Sum)
	IF Dum >= AAMAX THEN
	  IMAX% = I%
	  AAMAX = Dum
	END IF
  NEXT I%
  IF J% <> IMAX% THEN
	FOR K% = 1 TO N%
	  Dum = A(IMAX%, K%)
	  A(IMAX%, K%) = A(J%, K%)
	  A(J%, K%) = Dum
	NEXT K%
	D = -D
	VV(IMAX%) = VV(J%)
  END IF
  Indx%(J%) = IMAX%
  IF A(J%, J%) = 0! THEN A(J%, J%) = TINY
  IF J% <> N% THEN
	Dum = 1! / A(J%, J%)
	FOR I% = J% + 1 TO N%
	  A(I%, J%) = A(I%, J%) * Dum
	NEXT I%
  END IF
NEXT J%
ERASE VV
END SUB

FUNCTION MatMultS% (Alpha!(), Beta!(), Gamma!())
'============================================================================
'... MatMultS% multiplies two single precision matrices and places the
' product in a result matrix. All arrays must be two-dimensional.
' A one-dimensional array is converted to two-dimensions by making the
' second dimension equal to 1.
' Parameters: matrices Alpha,Beta,Gamma

' Returns: Gamma() = Alpha() * Beta()
'============================================================================
ON LOCAL ERROR GOTO smulterr: MatMultS% = 0
' Check inside dimensions
IF (LBOUND(Alpha, 2) <> LBOUND(Beta, 1)) THEN ERROR 197
IF (UBOUND(Alpha, 2) <> UBOUND(Beta, 1)) THEN ERROR 197

' Check dimensions of result matrix
IF (LBOUND(Alpha, 1) <> LBOUND(Gamma, 1)) THEN ERROR 195
IF (UBOUND(Alpha, 1) <> UBOUND(Gamma, 1)) THEN ERROR 195
IF (LBOUND(Beta, 2) <> LBOUND(Gamma, 2)) THEN ERROR 195
IF (UBOUND(Beta, 2) <> UBOUND(Gamma, 2)) THEN ERROR 195

'... Loop through, Gamma(row,col)=inner product of Alpha(row,*)
'    and Beta(*,col)
FOR row% = LBOUND(Gamma, 1) TO UBOUND(Gamma, 1)
   FOR col% = LBOUND(Gamma, 2) TO UBOUND(Gamma, 2)
	  Gamma(row%, col%) = 0!
	  FOR inside% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
		 Gamma(row%, col%) = Gamma(row%, col%) + Alpha(row%, inside%) * Beta(inside%, col%)
	  NEXT inside%
   NEXT col%
NEXT row%

smultexit:
EXIT FUNCTION

smulterr:
   ErrCode% = ERR
   PRINT "   *** Error in matrix calculation ErrCode% = "; ErrCode%; " ***"
   MatMultS% = ERR
   RESUME smultexit
END FUNCTION

FUNCTION MatTrans% (Alpha!(), Beta!())
'=========================================================================
'... uses Alpha!() to create a tranposed matrix Beta!()
'    MatTranS% returns ErrCode%
'=========================================================================
ON LOCAL ERROR GOTO TranErr: MatTrans% = 0
'check if Alpha, Beta have same dimensions if not, exit and send back error
IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 2)) THEN ERROR 196
IF (UBOUND(Alpha, 1) <> UBOUND(Beta, 2)) THEN ERROR 196
IF (LBOUND(Alpha, 2) <> LBOUND(Beta, 1)) THEN ERROR 196
IF (UBOUND(Alpha, 2) <> UBOUND(Beta, 1)) THEN ERROR 196

'... duplicate matrix Alpha()
FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
		FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
				Beta(col%, row%) = Alpha(row%, col%)
		NEXT col%
NEXT row%

TranExit:
EXIT FUNCTION

TranErr:
   ErrCode% = ERR
   PRINT "   *** Error in matrix calculation ErrCode% = "; ErrCode%; " ***"
   MatTrans% = ErrCode%
   RESUME TranExit
END FUNCTION

FUNCTION RchiSq! (Epsilon!)
'============================================================================
SHARED Dm!(), Dmtrans!(), Edot!(), EdotCalc!(), Vdot!(), Wm!(), We!()
SHARED G!(), Gtrans!(), LHS!(), RHS!(), Indx%(), Num%
'============================================================================
'... Construct Dm() which forces a smooth solution. The smoothness
' constrain is not enforced for the oldest parameter which is
' a proxy for the initial elevation.
FOR I% = 3 TO Num% - 1
	FOR J% = 1 TO Num%
		SELECT CASE J%
			CASE IS = (I% - 1)
				Dm(I% - 2, J%) = Epsilon
			CASE IS = I%
				Dm(I% - 2, J%) = -2 * Epsilon
			CASE IS = (I% + 1)
				Dm(I% - 2, J%) = Epsilon
			CASE ELSE
				Dm(I% - 2, J%) = 0
		END SELECT
	NEXT J%
NEXT I%

'... Weighted Damped Least Squares (Meinke, 1989, p. 54)
' Equation:  Vdot() = [Gtrans() We() G() + Wm()]^-1 Gtrans() We() Edot()
' Note that Wm() here is equal to epsilon^2 Wm() in Meinke

'... Calculate RHS matrices: Gtrans() We() Edot()
ErrCode = MatTrans%(G(), Gtrans())
ErrCode% = MatMultS%(Gtrans(), We(), RHS())
FOR I% = 1 TO Num%
	Vdot(I%) = 0
	FOR J% = 1 TO Num%
		Vdot(I%) = Vdot(I%) + RHS(I%, J%) * Edot(J%)
	NEXT J%
NEXT I%

'... Calculate LHS matrices: Gtrans() We() G() + Wm()
ErrCode = MatTrans%(Dm(), Dmtrans())
ErrCode% = MatMultS%(Dmtrans(), Dm(), Wm())
ErrCode% = MatMultS%(RHS(), G(), LHS())
FOR I% = 1 TO Num%
	FOR J% = 1 TO Num%
		LHS(I%, J%) = LHS(I%, J%) + Wm(I%, J%)
	NEXT J%
NEXT I%

'... Use LU decomposition to find solution. On input, Vdot() is equal to
' the product of the RHS matrices. On output, Vdot() is set to the solution.
CALL Ludcmp(LHS(), Num%, Indx%(), D)
CALL Lubksb(LHS(), Num%, Indx%(), Vdot())

Sum = 0
FOR I% = 1 TO Num%
	EdotCalc(I%) = 0
	FOR J% = 1 TO Num%
		EdotCalc(I%) = EdotCalc(I%) + G(I%, J%) * Vdot(J%)
	NEXT J%
   Sum = Sum + We(I%, I%) * (Edot(I%) - EdotCalc(I%)) ^ 2
NEXT I%
RchiSq = Sum / Num%
END FUNCTION

SUB RTbis (EdotKnown!, Bdot!, Zhi!, Dt!, Kd!)
'==========================================================================
'... Using bisection, find for Kd for a function Edot=f(Kd) where Edot
' is known. Kdlw and Kdhi bracket the possible range for Kd. Edotknow is
' the known value for f(Kd). The subroutine iterates towards Kd until the
' estimate Kd has an accuracy of +/- KdACC.
' From Press et al., 1986, p. 247. C is defined in the main program.
'==========================================================================
CONST Jmax = 40
CONST KdACC = .0001
'==========================================================================
Kdlw = 0
Kdhi = 10
Diff = Func(EdotKnown, Bdot, Zhi, Dt, Kdlw)
Diffmid = Func(EdotKnown, Bdot, Zhi, Dt, Kdhi)
IF Diff * Diffmid >= 0 THEN
	PRINT "  *** Root must be bracketed before bisection is started ***"
	STOP
END IF
IF Diff < 0 THEN
	Kd = Kdlw
	DKd = Kdhi - Kdlw
ELSE
	Kd = Kdhi
	DKd = Kdlw - Kdhi
END IF
FOR J% = 1 TO Jmax
	DKd = DKd * .5
	Kdmid = Kd + DKd
	Diffmid = Func(EdotKnown, Bdot, Zhi, Dt, Kdmid)
	IF Diffmid <= 0 THEN Kd = Kdmid
	IF ABS(DKd) < KdACC OR Diffmid = 0 THEN EXIT SUB
NEXT J%
PRINT " *** Too many bisections in finding root to correction equation ***"
STOP
END SUB