


c 		Barotropic coastal -trapped waves with a mean flow and
c		options for coastal and offshore boundary conditions.
c
c		January 5, 1995, April 16, 2004   2/7/2007
c
c		Changed to be consistent with DOS usage (1/4/95): "open"s files rather
c			than calling for assigning files before run time.
c		Program will ask for data file names once you start to run.
c		Program can also make flat files to use for matlab graphics.
c			Graphics options come up as you run the program.
c			Only one plot output is allowed per run.
c		Energy conserving norms are used for all computations:
c			this will mean that numerical values for ann and
c			bn differ from those in Brink and Chapman.
c
c		11/29/2005: Updated the Fortran usage and modified the outputs
c			slightly. (KHB)
c		2/7/2007:	Minor improvements made
c		3/9/2010:	Correct confusion about searching in u or p




	PROGRAM BTSW

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000), RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
	DIMENSION BH(4000),WW(3),RX(3),RLH(4000),WH(4000),CH(4000)
	DIMENSION WF(4000)
	DIMENSION TFH(4000)
	CHARACTER*12 FINP, FOT, FGP
C
C		MAIN PROGRAM CALLS ALL OF THE PIECES
C
C		NN LESS THAN OR EQUAL 4000
C



C		To test run this program, use input file "TESTBTCW.IN"

c			Find out input/output file names and open them
	WRITE(6,970)
	READ(5,*) FINP
	WRITE(6,971)
	READ(5,*) FOT
	OPEN (INPF,FILE=FINP)
	OPEN (IOPF,FILE=FOT)
	
	WRITE(6,950)
	READ(5,*) IGRP
	IF (IGRP.NE.0) THEN
		WRITE(6,972)
		READ(5,*) FGP
		OPEN (IGPF,FILE=FGP)
	END IF


	IMD = 1
	
c		Begin to model inputs	
	READ(INPF,*) NN
	WRITE(IOPF,929) NN

c		NN =      number of grid points to use in x
C		DD1 =  1. FREE SURFACE
C		DD1 =  0. RIGID LID
C		IDD3 = 1  U = 0 AT X = XMAX
C		IDD3 = 0  UX = 0 AT X = XMAX
C		IDD3 = 2  REAL B.C. AT X = XMAX FOR VX,HX = 0. THERE
C		IDD4 = 1  U = 0 AT X = 0
C		IDD4 = 0  UX = 0 AT X = 0
C		IDD4 = 2  DECAY B.C. AT X = 0 FOR VX,HX = 0 THERE
C		IPC =  0  REDUCED PRINT OUT
C		IUP =  0  SEARCH IN U
C		IUP =  1  SEARCH IN P
C		ILLW = 0  STRICTLY LONG WAVE
C		ILLW = 1  GENERAL FREQUENCY AND WAVENUMBER
C					(REVISED SEPTEMBER, 1988):
C					WFF (IN WFR) is now to be taken as the resistance
C					coefficient, R (cm/sec)
c		NITM =	  maximum number of iterations per frequency
c 		ISD =  0  unbiased search
c		ISD =  1  search only to lower frequencies
c		ISD = -1  search to higher frequencies only
c		EPS =     desired fractional accuracy for frequency
c		DEL =     fractioinal search increment (typically 0.05)

	READ(INPF,*) NITM,ISD,EPS,DEL
	READ(INPF,*) IUP,ILLW
	READ(INPF,*) IDD1,IDD3,IDD4
      DD1 = FLOAT(IDD1)
	NCAL = 1
	READ(INPF,*) NCALM,ILW
	READ(INPF,*) RLF,DRL
	READ(INPF,*) IPC

c		NCALM =	  Maximum number of frequencies per dispersion curve
c		ILW =  0  no long-wave parameters calculated
c		ILW ~= 0  long wave parameters calculated
c		RLF =     first wavenumber to calculate frequency
c					 for (*e+08cm^-1)
c		DRL =     wavenumber increment as in RLF
c		IPC =  0  short printout
c		IPC ~= 0  full printout



	DRL = DRL*1.0E-08
	RLF = RLF*1.0E-08
	READ(INPF,*) F,XMAX
c		F = coriolis parameter *e5 (sec^-1)
c		XMAX = distance (km) to offshore edge of grid


	XMAX = XMAX*1.0E+05
	F = F*1.0E-05
	DX = XMAX/FLOAT(NN-1)
	WRITE(IOPF,907) F,XMAX,DX


c		Sort through and print model assumptions and boundary
c			 conditions
	
	IF (IUP.EQ.1) THEN
c			search in p
		WRITE(IOPF,925)
	ELSE
c			search in U
		WRITE(IOPF,924)
		IUP = 0
	END IF

	IF (IDD4.EQ.2) THEN
c			decay boundary condition at x = 0
		WRITE(IOPF,928)
	ELSE IF (IDD4.EQ.0) THEN
c			Ux = 0 at x = 0
		WRITE(IOPF,921)
	ELSE
c			U = 0  at x = 0
		WRITE(IOPF,922)
		IDD4=1
	END IF 

	IF (IDD3.EQ.1) THEN
c			U = 0 at x = XMAX
		WRITE(IOPF,920)
	ELSE IF (IDD3.EQ.2) THEN
c			Decay (real) boundary condition at x = XMAX
		WRITE(IOPF,923)
	ELSE
c			Ux = 0 at x = XMAX
		WRITE(IOPF,919)
		IDD3 = 0
	END IF

	IF (IDD1.EQ.1) THEN
c			Free surface
		WRITE(IOPF,918)
	ELSE
c			Rigid lid
		DD1 = 0.
		WRITE(IOPF,917)
	END IF

	IF (ILLW.EQ.1) THEN
c			general frequency and wavenumber
		WRITE(IOPF,927)
	ELSE
c			long wave approximation
		WRITE(IOPF,926)
		ILLW = 0
		NCALM = 1
	END IF
	REPS = FLOAT(ILLW)


c		read depth profile
	CALL DEP(NN,XMAX)

c		read mean velocity profile
	CALL VRD(NN,XMAX,F,NV)
	

C		Check to see that x = 0 boundary condition is consistent
	STEST = 1.0E-7
	AA = ABS(RHX(1))
	IF ((IDD4.EQ.2).AND.(AA.GT.STEST)) THEN
		WRITE(IOPF,930)
		GO TO 140
	END IF
	AA = ABS(RHX(NN))
	IF ((IDD3.EQ.2).AND.(AA.GT.STEST)) THEN
		WRITE(IOPF,931)
		GO TO 140
	END IF


	DO I = 1,NCALM
		RLH(I) = DRL*FLOAT(I-1) + RLF
	END DO
	RL = RLH(NCAL)


c		read scaled first 3 guesses at frequency
	READ(INPF,*) (WW(J),J=1,3)

	DO J = 1,3
		WW(J) = WW(J)*1.0E-05
	END DO

c		Start iterating frequency values
1	RB = 0.
	WRITE(IOPF,908)

c		Do the calculation for the first three estimates
	DO 18 J = 1,3
		W = WW(J)
		IF (IUP.EQ.0) THEN
			CALL MUTS(NN,XMAX,DD1,F,W,RL,IDD3,RI,IDD4,REPS)
		ELSE
			CALL MATS(NN,XMAX,DD1,F,W,RL,IDD3,RI,IDD4,REPS)
		END IF

		CP = W/RL
		CALL CRLC(W,RL,NN,ICR)
		WRITE(IOPF,905) W,RL,CP,RI,ICR
		RX(J) = RI
		IF (RI.GT.RB) THEN
			RB = RI
			ICRH = ICR
			WB = W
			DO I = 1,NN
				BH(I) = B(I)
			END DO
		END IF
18	CONTINUE
	
	NIT = 3
	IGP = 0

c		Now guess the next estimate for frequency
20	CALL NGW(WW,RX,WB,WN,ISUC,EPS,DEL,IN,ISD,IGP)

	IF (ISUC.EQ.1) GO TO 100

	IF (IUP.EQ.0) THEN
		CALL MUTS(NN,XMAX,DD1,F,WN,RL,IDD3,RI,IDD4,REPS)
	ELSE
		CALL MATS(NN,XMAX,DD1,F,WN,RL,IDD3,RI,IDD4,REPS)
	END IF

	IF (IN.EQ.0) THEN
		IF (RI.LT.RB) THEN
			IF (WN.GT.WW(2)) THEN
				IN = 3
			ELSE	
				IN = 1
			END IF
		ELSE
			IF (WN.LT.WW(2)) THEN
				IN = 3
			ELSE
				IN = 1
			END IF
		END IF
	END IF

	RX(IN) = RI
	WW(IN) = WN
	CP = WN/RL
	CALL CRLC(WN,RL,NN,ICR)
	WRITE(IOPF,905) WN,RL,CP,RI,ICR

	IF (RI.GT.RB) THEN
		WB = WN
		ICRH = ICR
		RB = RI
		DO I = 1,NN
			BH(I) = B(I)
		END DO
	END IF

	NIT = NIT + 1
c		check iteration count	
	IF (NIT.LT.NITM) GO TO 20
		WRITE(IOPF,902) NIT
	GO TO 250

100	DO I = 1,NN
		BH(I) = BH(I)/RB
	END DO
	
	CC = WB/RL
	IF (ICRH.NE.0) THEN
		WRITE(IOPF,911) ICRH
	END IF

	WRITE(IOPF,903) WB,RL,CC,EPS,DX,NIT
	CALL DTIM(WB,RL,NN,XMAX,F,IDD3,IDD4,ILLW,IDD1,RTF)
	TFH(NCAL) = RTF
	
120	IF ((IUP.NE.0).AND.(NCAL.EQ.1)) THEN
		CALL MUTS(NN,XMAX,DD1,F,WB,RL,IDD3,RI,IDD4,REPS)
		DO I = 1,NN
			B(I) = B(I)/RI
		END DO
		WRITE(IOPF,913)
		WRITE(IOPF,904) (B(I),I = 1,NN)
	ELSE IF ((RH(1).NE.0.).AND.(NCAL.EQ.1)) THEN
		CALL MATS(NN,XMAX,DD1,F,WB,RL,IDD3,RI,IDD4,REPS)
		DO I = 1,NN
			B(I) = B(I)/RI
		END DO
c		IF (IUP.EQ.1) GO TO 145
		WRITE(IOPF,912)
		WRITE(IOPF,904) (B(I),I=1,NN)
		IF ((ILW.EQ.1).AND.(NV.NE.0).AND.(IDD4.EQ.1)) THEN
			WRITE(IOPF,908)
			CALL LGWV(NN,XMAX,F,WF,WB,RL,DD1,IGRP)
		END IF
	END IF


c		Print flat file for graphics if desired
	IF (IGRP.NE.0) THEN
		WRITE(6,955) NCAL
		READ(5,*) IXXX
		IF (IXXX.NE.0) THEN
			IGRP = 0
			DO I = 1,NN
				XXX = DX*FLOAT(I-1)
				WRITE(IGPF,*) XXX,B(I)
			END DO
		END IF
	END IF

	WH(NCAL) = WB
	CH(NCAL) = CC
	NCAL = NCAL + 1
	IF ( NCAL.GT.NCALM) GO TO 250

c		Make a first guess at the next frequency on the curve
	RL = RLH(NCAL)
	IF (NCAL.LT.3) THEN
		WW(2) = CC*RL
	ELSE
		I1 = NCAL - 2
		I2 = NCAL - 1
		CG = (WH(I2) - WH(I1))/(RLH(I2) - RLH(I1))
		WW(2)=  WH(I2) + CG*(RLH(NCAL)-RLH(I2))
	END IF

	WW(1) = WW(2)*(1. - DEL)
	WW(3) = WW(2)*(1. + DEL)
	GO TO 1


c		Reached maximum number of wavenumbers for this curve
250	WRITE(IOPF,909)
	DO I = 1,NCALM
		WRITE(IOPF,910) WH(I),RLH(I),CH(I),TFH(I)
	END DO

	IF (IGRP.NE.0) THEN
c			Print out dispersion curve for graphics, if desired		
		WRITE(6,951)
		READ(5,*) IXX
		IF (IXX.NE.0) THEN
			WRITE(IGPF,*) 0., 0., 0.
			DO I = 1,NCALM
				WRITE(IGPF,*) RLH(I),WH(I),TFH(I)
			END DO
		END IF
	END IF



902	FORMAT(//'  USED UP',I3,' ITERATIONS'/)
903	FORMAT(//'  CONVERGED: W,L,C,EPS,DX,NIT =',5E15.5,I5)
904	FORMAT(10E13.5)
905	FORMAT(/' W,RL,CP,RI,ICR =',4E15.5,I10)
907	FORMAT(/' F,XMAX,DX =',3E15.5/)
908	FORMAT(///)
909	FORMAT(///'     W              L                C             RTF'/)
910	FORMAT(4E15.5)
911	FORMAT(/'  SOLUTION HAS ',I3,' CRITICAL LAYERS'/)
912	FORMAT(/'   ZETA'/)
913	FORMAT(/'   U'/)
917	FORMAT(' RIGID LID'/)
918	FORMAT(' FREE SURFACE'/)
919	FORMAT(' UX = 0 AT X = XMAX'/)
920	FORMAT(' U = 0 AT X = XMAX'/)
921	FORMAT(' UX = 0 AT X = 0'/)
922	FORMAT(' U = 0 AT X = 0'/)
923	FORMAT(' REAL B.C.AT X = XMAX'/)
924	FORMAT(/' SEARCH IN U'/)
925	FORMAT(/' SEARCH IN P'/)
926	FORMAT(' LONG WAVE EXACTLY')
927	FORMAT(' GENERAL FREQUENCY AND WAVENUMBER')
928	FORMAT(' DECAY CONDITION AT X = 0'/)
929	FORMAT(' NN=',I5)
930	FORMAT(//' ERROR! Must use IDD4=2 with RHX(1) = 0 !!!!'/)
931	FORMAT(//' ERROR! Must use IDD3=2 with RHX(NN) = 0 !!!!'/)
950 	FORMAT(' Do you want graphics output? 0 (for no) or 1 (for yes)')
951	FORMAT(' Do you want a dispersion curve? 0 (for no) or 1 (for
     & yes)')
955	FORMAT(' Do you want an unnormalized modal structure for ncal ='   
     & , I3,'?:  0 (for no) or 1 (for yes)')
970	FORMAT(' Input file name?')
971	FORMAT(' Output file name?')
972	FORMAT(' Graphics file name?')

140	CLOSE(INPF)
	CLOSE(IOPF)	
	IF (IGRP.NE.0) THEN
		CLOSE(IGPF)
	END IF

	STOP
	END




	SUBROUTINE NGW(WW,RX,WB,WN,ISUC,EPS,DEL,IN,ISD,IGP)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	DIMENSION WW(3),RX(3)
C
C		SUBROUTINE TO GUESS THE NEXT W
C
5	IC = 0
c			Put the frequency estimates in order
	DO 10 I = 1,2
		I1 = I +1
		AA = ABS(WW(I1))
		BB = ABS(WW(I))
		IF (AA.GT.BB) GO TO 10
			IC = 1
			RI = RX(I)
			WX = WW(I)
			RX(I) = RX(I1)
			WW(I) = WW(I1)
			WW(I1) = WX
			RX(I1) = RI
10	CONTINUE
	IF (IC.NE.0) GO TO 5


	ISUC = 0
	IF (RX(3).GT.RX(2)) GO TO 150
	IF (RX(1).GT.RX(2)) GO TO 160
	IF (RX(3).GT.RX(1)) THEN
		IL = 1
		IH = 3
	ELSE
		IL = 3
		IH = 1
	END IF

c		Try Newton's rule
	WB1 = (WW(1) +WW(2))/2.
	WB2 = (WW(2) + WW(3))/2.
	RW1 = (RX(2) - RX(1))/(WW(2) - WW(1))
	RW2 = (RX(3) - RX(2))/(WW(3) -WW(2))
	A = (RW2 -RW1)/(WB2 - WB1)
	B = RW1  - A*WB1
	WN = -B/A
	EP= (WN -WB)/WB
	EP = ABS(EP)
	IF (EP.LT.EPS) GO TO 100

	IGP = IGP + 1
	IF ((WN.LE.WW(1)).OR.(WN.GE.WW(3))) THEN
c			Newton's rule gave an out of bounds guess 		
		WN = (WW(2) + WW(IH))/2
	END IF
20	IN = 0
	GO TO 130
	

100	ISUC = 1
c		Convergence?
	IF (IGP.NE.0) GO TO 130

	IGP = IGP + 1
	IN = 0
	ISUC = 0
	WN = (WW(2) + WW(IL))/2.
	GO TO 130

c		New estimate has not yet been surrounded
150	IF (RX(1).GT.RX(2)) GO TO 180
155	IF (ISD.EQ.1) GO TO 165
		WN = WW(3)*(1. + DEL)
		WW(1) = WN
		RX(1) = 0.
		IN = 1
	GO TO 130

160	IF (RX(3).GT.RX(2)) GO TO 180
165	IF (ISD.EQ.-1) GO TO 155
		WN = WW(1)*(1. - DEL)
		WW(3) = WN
		RX(3) = 0.
		IN = 3
	GO TO 130

180	IF (RX(3).GT.RX(1)) GO TO 155
	GO TO 165

130	RETURN
	END



	SUBROUTINE CALC(NN,RX)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
	DOUBLE PRECISION RI,DSQRT
C
C		SUBROUTINE TO CALCULATE THE INTEGRAL OF RESPONSE SQUARED
C
	RI = (B(1)**2 + B(NN)**2)/2.
	NX = NN - 1
	DO I = 2,NX
		RI = RI + B(I)**2
	END DO
	RI = DSQRT(RI)
	RX = RI

	RETURN
	END



	SUBROUTINE DEP(NN,XMAX)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)	
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
	
C
C		SUBROUTINE TO READ AND INTERPOLATE THE DEPTH PROFILE
C			RH = DEPTH
C			RHX = X DERIVATIVE OF DEPTH
C
	
c		Read in depth profile
	READ(INPF,*) NRX
	DO I = 1,NRX
		READ(INPF,*) A(I),B(I)
		A(I) = A(I)*1.0E+05
		B(I) = B(I)*100.
	END DO
	
c		Interpolate depths from input information
	DX = XMAX/FLOAT(NN-1)
	RH(1) = B(1)
	DO 20 N = 2,NN
		X = DX*FLOAT(N-1)
		IF (X.GT.A(NRX)) THEN 
			RH(N) = B(NRX)
		ELSE
			IC = 0
			DO J = 2,NRX
				IF (IC.EQ.0) THEN
					I = J
					IF (X.LE.A(i)) THEN
						IC = I
					END IF
				END IF 
		    END DO
			IM = I -1
			AA = (B(I) -B(IM))/(A(I) - A(IM))
			XX = X - A(IM)
			RH(N) = B(IM) +AA*XX
		END IF
20	CONTINUE


c		Now compute the derivative
	RHX(1) = (RH(2) - RH(1))/DX
	NM = NN -1
	RHX(NN) = (RH(NN) - RH(NM))/DX
	D2 = 2.*DX
	DO N = 2,NM
		IP = N +1
		IM = N -1
		RHX(N) = (RH(IP) - RH(IM))/D2
	END DO

c		write out depth profile
	WRITE(IOPF,903)
	WRITE(IOPF,902) (RH(N),N=1,NN)

902	FORMAT(10E13.5)
903	FORMAT(//'  DEPTH IN CM'/)

	RETURN
	END





	SUBROUTINE LGWV(NN,XMAX,F,WFF,WB,RL,DD1,igrp)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
	DIMENSION WFF(4000)
C
C		SUBROUTINE TO CALCULATE THE BN AND ANN IN THE LONG WAVE
C			LIMIT
C
	DX = XMAX/FLOAT(NN -1)
	DXX = 2.*DX
	GF = 980./F
	CALL WFR(NN,XMAX,WFF)

	IF (DD1.NE.0.) GO TO 200
c		The following is for rigid lid only

c		Compute stream function from the pressure
	A(1) = 0.
	DO 100 N = 2,NN
		A(N) = (RHX(N)*B(N) + RHX(1)*B(1))/2.
		NX = N - 1
					
		IF (NX.GT.1) THEN	
			DO I = 2,NX
				A(N) = A(N) + RHX(I)*B(I)
			END DO
		END IF
		
		A(N) = A(N)*DX
		A(N) = A(N) + RH(1)*B(1) -RH(N)*B(N)
		A(N) = -GF*A(N)
100	CONTINUE

c		Now normalize stream function mode
	NX = NN - 1
	WF = RHX(1)*(A(1)/RH(1))**2
	WF = WF + RHX(NN)*(A(NN)/RH(NN))**2
	WF = WF/2.
	DO N = 2,NX
		WF = WF + RHX(N)*(A(N)/RH(N))**2
	END DO
	
	WF = WF*DX*F
	WF = SQRT(ABS(WF))
	DO N = 1,NN
		A(N) = A(N)/WF
	END DO

c		Now compute wind coupling BN and damping WF
	PX = (A(2) -A(1))/(DX*RH(1))
	WF = 0.5*PX*PX*WFF(1)
	NH = NN-1
	DO N = 2,NH
		BN = BN + RHX(N)*A(N)/(RH(N)**2)
		I1 = N - 1
		I2 = N + 1
		PX = (A(I2) -A(I1))/DXX
		PX = PX/RH(N)
		WF = WF+ PX*PX*WFF(N)
	END DO
	BN = BN*DX
	WF = WF*DX

	WRITE(IOPF,905)
	WRITE(IOPF,901) WF
	WRITE(IOPF,903) BN
	WRITE(IOPF,904)
	WRITE(IOPF,902) (A(I),I=1,NN)

	IF (IGRP.NE.0) THEN
c			option of printing flat file of streamfunction modal
c				 structure	
		WRITE(6,950) 
		READ(5,*) IXXX
		IF (IXXX.NE.0) THEN
			IGRP= 0
			DO I = 1,NN
				XXX = DX*FLOAT(I-1)
				WRITE(IGPF,*) XXX,A(i)
			END DO
		END IF
	END IF



200	C = WB/RL
c		Now work with pressure mode
c		
c		First, normalize the mode
	NX = NN - 1
	WF = RHX(1)*B(1)*B(1)*0.5
	DO I = 2,NX
		WF = WF +RHX(I)*B(I)*B(I)
	END DO
	WF = WF + 0.5*RHX(NN)*B(NN)*B(NN)
	WF = WF*DX
	WF = WF + RH(1)*B(1)*B(1)

	WF =SQRT(WF/F)

	DO I = 1,NN
		B(I) = B(I)/WF
	END DO

c		Now compute wind coupling BN and damping WF
	BX = -F*B(1)/C
	WF = 0.5*WFF(1)*BX*BX
	DO I = 2,NX
		IP = I + 1
		IM = I - 1
		BX = (B(IP) - B(IM))/DXX
		WF = WF + WFF(I)*BX*BX
	END DO
	WF = WF*DX/F
	ALX =DD1/(980.0*RH(NN))
	ALX = SQRT(ALX)
	WF = WF + 0.5*ALX*B(NN)*B(NN)*WFF(NN)
	BN = -B(1)/F
	WF = WF/F

	WRITE(IOPF,906)
	WRITE(IOPF,907) WF
	WRITE(IOPF,908) BN
	WRITE(IOPF,909)
	WRITE(IOPF,902) (B(I),I = 1,NN)


	IF (IGRP.NE.0) THEN
c		Option of printing a flat file for graphics
		WRITE(6,951)
		READ(5,*) IXXX
		IF (IXXX.NE.0) THEN
			IGRP = 0
			DO I = 1,NN
				XXX = DX*FLOAT(I-1)
				WRITE(IGPF,*) XXX,B(I)
			END DO
		END IF
	END IF
	
901	FORMAT(' ANN =',E15.5,' 1/CM')
902	FORMAT(10E13.5)
903	FORMAT(' BN =',E15.5,' CM-1/2 S1/2')
904	FORMAT(//'  PHI (CM1/2 S1/2) ='/)
905	FORMAT(//' USING STREAM FUNCTION')
906	FORMAT(//' USING PRESSURE')
907	FORMAT(' ANN = ',E15.5,' CM-1')
908	FORMAT(' BN =',E15.5,' CM-1/2 S1/2')
909	FORMAT(//' NORMALIZED PRESSURE (CM-1/2 S-1/2) ='/)
950	FORMAT(' Want a normalized streamfunction modal plot? 0 (for no)
     &or 1 (for yes)')
951	FORMAT(' Want a normalized pressure modal plot? 0 (for no) or 1 
     &(for yes)')

500	RETURN
	END





	SUBROUTINE MUTS(NN,XMAX,DD1,F,W,RL,IDD3,RI,IDD4,REPS)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
	DOUBLE PRECISION A1,A2,A3,AL,DSQRT
C
C		SUBROUTINE TO SET UP AND SOLVE THE U EQUATION
C
	NU = NN
	IF (IDD3.EQ.1) THEN
		NU = NN-1
	END IF

	IF (IDD4.EQ.1) THEN
		NU = NU - 1
	END IF 

	DX = XMAX/FLOAT(NN-1)
	GG = DD1/980.
	RLL = RL*RL
	DXX = DX*DX
	DDD = DX/2.

	DO 100 N = 1,NU
		N1 = N
		IF (IDD4.EQ.1) THEN
			N1 = N+1
		END IF

		B(N) = 0.
		WP = W + RL*V(N1)
		W2 = WP*WP
		RR = RH(N1)
		RX = RHX(N1)
		WX = RL*VX(N1)
		FP = F + VX(N1)
		A1 = RR*(GG*W2 - RR*RLL)
		A2 = RR*(RX*RLL -2.0*WP*WX*GG)
		A3 = -GG*GG*W2*(F*FP -REPS*W2)
		A3 = A3 + REPS*RR*RR*RLL*RLL
		A3 = A3 + GG*RR*RLL*(-2.0*W2*REPS + F*FP + 2.0*FP*VX(N1))
		A3 = A3 - RR*RLL*RL*FP*RX/WP
		A3 = A3 - RR*RL*VXX(N1)*(GG*W2 - RR*RLL)/WP
		N1 = N + NU
		N2 = N + 2*NU
		A(N) = -2.*A1 + A3*DXX
		A(N1) = A1 -A2*DDD
		A(N2) = A1 + A2*DDD
	

		IF (N.EQ.1) THEN
c				x = 0 boundary conditions
			IF (IDD4.EQ.0) THEN
				A(N2) = A(N1) + A(N2)
			ELSE IF (IDD4.EQ.2) THEN
				AL = DSQRT(-A3/A1) 
				A(N2) = A(N2) + A(N1)
				A(N) = A(N) -2.0*DX*AL*A(N1)
			END IF
			A(N1) = 0.

		ELSE IF (N.EQ.NU) THEN
c				x = XMAX  boundary conditions
			IF (IDD3.EQ.0) THEN
				A(N1) = A(N1) + A(N2)
			ELSE IF (IDD3.EQ.2) THEN
				AL = DSQRT(-A3/A1)
				A(N1) = A(N1) +A(N2)
				A(N) = A(N) -2.0*DX*AL*A(N2)
			END IF
			A(N2) = 0.

		END IF
100	CONTINUE

	B(5) = 1.
	CALL TRI(NU)
	
	IF (IDD4.EQ.1) THEN	
		DO I = 1,NU
			N = NU -I +1
			N1 = N + 1
			B(N1) = B(N)
		END DO
		B(1) = 0.
	END IF

	IF (IDD3.EQ.1) THEN
		B(NN) = 0.
	END IF

	CALL CALC(NN,RI)

	RETURN
	END




	SUBROUTINE VRD(NN,XMAX,F,NRD)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
C
C		SUBROUTINE TO READ AND INTERPOLATE THE MEAN V PROFILE
C			V = ALONGSHORE VELOCITY
C			VX = X DERIVATIVE OF V
C			VXX = SECOND X DERIVATIVE OF V
C

	DX = XMAX/FLOAT(NN-1)

c		Read mean velocity information
	READ(INPF,*) NRD
	IF (NRD.EQ.0) GO TO  200

	DO I = 1,NRD
		READ(INPF,*) A(I),B(I)
C			A = DISTANCE FROM SHORE IN KM
C			B = V IN CM/SEC
		A(I) = A(I)*1.0E+05
	END DO

c		Interpolate velocity field from inputs
	DO 100 N = 1,NN
		X = DX*FLOAT(N-1)
		IF (X.GT.A(NRD)) THEN
			V(N) = B(NRD)
		ELSE IF (X.LT.A(1)) THEN
			V(N) = B(1)
		ELSE	
			J = 0
			DO I = 2,NRD
				IF (J.EQ.0) THEN
					IF (X.LE.A(I)) THEN
						J = I
					END IF
				END IF
			END DO
			JJ = J - 1
			XX = A(J) - X
			XY = A(J) - A(JJ)
			V(N) = B(J) -XX*(B(J) -B(JJ))/XY
		END IF

100	CONTINUE


c		Now compute first and second derivatives of mean velocity
	VX(1) = (V(2) - V(1))/DX
	NQ = NN - 1
	VX(NN) = (V(NN) - V(NQ))/DX
	DXX = DX*DX
	DZ = 2.*DX
	VXX(1) = 0.
	VXX(NN) = 0.
	DO N = 2,NQ
		I1 = N - 1
		I2 = N + 1
		VX(N) = (V(I2) - V(I1))/DZ
		VXX(N) = (V(I2) - 2.*V(N) + V(I1))/DXX
	END DO

	WRITE(IOPF,904)
	WRITE(IOPF,903) (V(N),N=1,NN)

c		Check for potential instability
	IST = 0
	DO N = 1,NN
		Q = -RHX(N)*(F+VX(N)) + RH(N)*VXX(N)
		A(N) = Q
		IF (IST.EQ.0) THEN
			IF (N.EQ.1) THEN
				PQ = Q
			ELSE IF (PQ.NE.0) THEN	
				RQ = Q/PQ
				IF (RQ.LT.0) THEN
					IST = 1
				END IF
				IF (Q.NE.0) THEN
					PQ = Q
				END IF
			ELSE 
				PQ = Q
			END IF
		END IF
	END DO

	IF (IST.NE.0) THEN
C			It is potentially unstable if you are here
		WRITE(IOPF,905)
		WRITE(IOPF,906)
		WRITE(IOPF,905)
		WRITE(IOPF,907)
		WRITE(IOPF,903) (A(N),N=1,NN)
	END IF
	GO TO 250

200	DO N = 1,NN
		V(N) = 0.
		VX(N) = 0.
		VXX(N) = 0.
	END DO

903	FORMAT(10E13.5)
904	FORMAT(//'  V IN CM/SEC'/)
905	FORMAT(/' *************************************************'/)
906	FORMAT('  POSSIBILITY OF UNSTABLE MODES')
907	FORMAT(//'  PBARX'/)

250	RETURN
	END



	SUBROUTINE MATS(NN,XMAX,DD1,F,W,RL,IDD3,RI,IDD4,REPS)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
	DOUBLE PRECISION A1,A2,A3,B1,B2,B3,BB1,BB2,BX,DSQRT,AL
C
C		SUBROUTINE TO SET UP AND SOLVE THE PRESSURE EQUATION
C
	DX = XMAX/FLOAT(NN - 1)
	GG = DD1/980.
	RLL = RL*RL*REPS
	DXX = DX*DX
	DDX = 2.*DX
	DDD = DX/2.

	DO N = 1,NN
		B(N) = 0.
		WP = W + RL*V(N)
		FP = F + VX(N)
		FW = F*FP - WP*WP*REPS
		RR = RH(N)
		RX = RHX(N)
		WX = RL*VX(N)
		FLW = F*RL/WP
		A1 = FW*RR
		A2 = FW*RX - RR*(F*VXX(N) -2.0*WP*WX*REPS)
		A3 = -GG*FW*FW +RX*FLW*FW -RR*FLW*(F*VXX(N) -2.*WP*WX*REPS)
		A3 = A3 -RR*RLL*FW
		N1 = N + NN
		N2 = N + 2*NN
		A(N) = -2.*A1 + A3*DXX
		A(N1) = A1 - A2*DDD
		A(N2) = A1 + A2*DDD
		
		IF (N.EQ.1) THEN
			IF (IDD4.EQ.0) THEN
				B1 = RL*RR*FP
				B2 = RLL*WP*RR +GG*WP*FW
				A(N) = A(N) + 2.0*DX*B2*A(N1)/B1
				A(N2) = A(N2) + A(N1)
			ELSE IF (IDD4.EQ.1) THEN
				A(N2) = A(N2) + A(N1)
				A(N) = A(N) + DDX*FLW*A(N1)
			ELSE
				AL = DSQRT(-A3/A1)
				A(N2) = A(N2) + A(N1)
				A(N) = A(N) - 2.0*DX*AL*A(N1)
			END IF
			A(N1) = 0

		ELSE IF (N.EQ.NN) THEN
			IF (IDD3.EQ.0) THEN 
				B1 = -RR*WP*FW
				B2 = -FW*(RR*RL*F + RX*WP + RR*WX)
				B2 = B2 + RR*WP*(F*VXX(N) -2.0*WP*WX*REPS)
				B3 = RL*F*(-RX*FW + RR*(F*VXX(N) -2.*WP*WX*REPS))
				BX = 1. + DDD*B2/B1
				BB1 = -(-2. +B3*DXX/B1)/BX
				BB2 = -(1. -B2*DDD/B1)/BX
				A(N) = A(N) + A(N2)*BB1
				A(N1) = A(N1) + A(N2)*BB2
			ELSE IF (IDD3.EQ.1) THEN
				A(N1) = A(N1) + A(N2)
				A(N) = A(N) -DDX*FLW*A(N2)
			ELSE
				AL = DSQRT(-A3/A1)
				A(N1) = A(N1) + A(N2)
				A(N) = A(N) -2.0*DX*AL*A(N2)
			END IF
			A(N2) = 0.

		END IF

	END DO

C		An arbitrary forcing
	B(5) = 0.00001

C		Check for dangerously small diagonal elements
	EQQ = 1.0E-36
	IKK = 0
	DO I = 1,NN
		AQ = ABS(A(I))
		IF (AQ.LT.EQQ) THEN
			IKK = I
		END IF
	END DO
	IF (IKK.NE.0) THEN  
		WRITE(IOPF,901)
		RI = 0
	ELSE
C			Solve it if it is OK
		CALL TRI(NN)
		CALL CALC(NN,RI)
	END IF	

901	FORMAT(' NUMERICAL PROBLEM: SMALL DIAGONAL ELEMENT')

	RETURN
	END





	SUBROUTINE CRLC(W,RL,NN,ICR)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
C
C		SUBROUTINE TO CHECK FOR CRITICAL LAYERS
C
	ICR = 0
	DO N = 1,NN
		WP = W + RL*V(N)
		IF (N.EQ.1) THEN
			PQ = WP
		ELSE
			PQ = WP/PQ
			IF (PQ.LT.0.) THEN
				ICR = ICR + 1
				PQ = WP
			END IF
		END IF
	END DO

	RETURN
	END
	




	SUBROUTINE TRI(N)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
	DOUBLE PRECISION AA,BB
C
C
C		SUBROUTINE TO SOLVE A TRIDIAGONAL MATRIX BY GAUSSIAN
C			ELIMINATION
C
C		STORE MAIN DIAGONAL FIRST
C		0,LOWER DIAGONAL
C		UPPER DIAGONAL,0
	N3 = N +1
	N2 = 2*N
	NN = N - 1

C		Eliminate the lower diagonal
	DO I = 1,NN
		II = N3 + I
		AA = A(II)/A(I)
		I2 = I + 1
		I3 =N2 + I
		BB = A(I2) -AA*A(I3)
		A(I2) = BB
		BB = B(I2) -AA*B(I)
		B(I2) = BB
	END DO

C		Backsolve
	DO IJ = 1,NN
		I = N3 - IJ
		B(I) = B(I)/A(I)
		II = I + N2 -1
		IP = I -1
		BB = B(IP) -A(II)*B(I)
		B(IP) = BB
	END DO

	B(1) = B(1)/A(1)

	RETURN
	END





	SUBROUTINE WFR(NN,XMAX,WF)

	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000),VXX(4000)
	DIMENSION WF(4000)
C
C		SUBROUTINE TO READ AND INTERPOLATE FRICTIONAL WEIGHT
C			FUNCTIONS, WF
C
	WRITE(IOPF,904)

C		Read in frictional weight functions
	READ(INPF,*) NW 
	IF (NW.EQ.0) GO TO 100

	DO I = 1,NW
		II = 100 + I
		READ(INPF,*) A(I),A(II)
		A(I) = A(I)*1.0E+05
	END DO

	DX = XMAX/FLOAT(NN-1)
	WF(1) = A(101)
	NMM = 100 + NW
	WFM = A(NMM)
	
	DO 50 N = 2,NN
		X = DX*FLOAT(N-1)
		IF (X.GT.A(NW)) THEN
			WF(N) = WFM
		ELSE IF (X.LT.A(1)) THEN
			WF(N) = WF(1)
		ELSE
			IC = 0
			DO J = 2,NW
				IF (IC.EQ.0) THEN
					I = J
					IF (X.LE.A(I)) THEN
						IC = I
					END IF
				END IF
			END DO
			II = 100 + I
			IIM = II - 1
			IM = I -1
			WX = (A(II) - A(IIM))/(A(I) -A(IM))
			XX =X -A(IM)
			WF(N) = A(IIM) + WX*XX
		END IF
50	CONTINUE

	WRITE(IOPF,901)
	WRITE(IOPF,902) (WF(I),I=1,NN)
	GO TO 300

	
100	DO I = 1,NN
		WF(I) = 1.0
	END DO
	WRITE(IOPF,903)
	

901	FORMAT(/' FRICTION WEIGHT FUNCTION')
902	FORMAT(10E12.5)
903	FORMAT(/' FRICTION WEIGHT FUNCTION = 1.0 EVERYWHERE')
904	format(' Assumed that the friction weight function is in cm/sec')

300	RETURN
	END




	SUBROUTINE DTIM(W,RL,NN,XMAX,F,IDD3,IDD4,ILLW,IDD1,RTF)
	
	PARAMETER (INPF = 7, IOPF = 8, IGPF = 12)
	COMMON A(12000),RH(4000),RHX(4000),B(4000),V(4000),VX(4000)
	COMMON VXX(4000)
	DOUBLE PRECISION AA,BB,RI(6)
C
C			COMPUTE FRICTIONAL DECAY TIME FOR ANY WAVE
C				FOR CONSTANT R
C			SEPTEMBER 30, 1988
C
C			Only does the calculation for a closed coast and
C				 open offshore BC
C			Assumes that bottom fricxtion is constant across
C				 the domain
		IF (IDD3.EQ.1) GO TO 200
		IF (IDD4.NE.1) GO TO 200

        	DD1 = FLOAT(IDD1)
		REPS = FLOAT(ILLW)

C			Get the pressure field for the mode
		CALL MATS(NN,XMAX,DD1,F,W,RL,IDD3,RI,IDD4,REPS)
		
		FL = F*RL
		DX = XMAX/FLOAT(NN-1)
		RLL = REPS*RL*RL
		DX2 = 2.0*DX
		GG = DD1/980.
		WP = W + RL*V(NN)
		FW = F*F - REPS*WP*WP
		RLB2 = RLL + GG*FW/RH(NN)
		RLB = SQRT(RLB2)

C		Compute the various integrals
	DO I = 1,6
		RI(I) = 0.
	END DO

	DO I = 1,NN
		DXX = DX
		WP = W + RL*V(I)
		FP = F + VX(I)
		FF = F*FP
		W2 = WP*WP
		FW = FF -REPS*W2
		FPW = FF + REPS*W2
		RP2 = B(I)*B(I)
		FW2 = FW*FW

		IF (I.EQ.1) THEN
			RPX2 = RP2*(FL/WP)**2
			RPP = -FL*RP2/WP
			DXX = DX/2.0
			AA = -RH(I)*FL*RP2*(FF -REPS*3.0*W2)
			AA = AA/(W2*FW2)
		ELSE IF (I.EQ.NN) THEN
			DXX = DX/2.0
			RPX2= RLB2*RP2
			RPP = -RLB*RP2
		ELSE
			IP = I + 1
			IM = I - 1
			RPX = (B(IP) - B(IM))/DX2
			RPP = RPX*B(I)
			RPX2 = RPX*RPX
		END IF

		X = RH(I)*WP*(RPX2 + RLL*RP2)/FW2
		RI(1) = RI(1) -2.0*DXX*X

		X = -RHX(I)*(FF-3.0*REPS*W2)
		X = (X+RH(I)*F*VXX(I))/(W2*FW2)
		RI(2) = RI(2) + FL*DXX*RP2*X

		X = -4.0*RH(I)*(F*VXX(I)  -2.0*WP*REPS*RL*VX(I))
		X = X/(FW*FW2)
		RI(3) = RI(3) + FL*DXX*X*RP2

		X = (FPW -2.0*F*VX(I))/(WP*FW2)
		RI(4) = RI(4) + RLL*X*RP2*DXX

		X = (-FW*VX(I) + 4.0*REPS*F*W2)/(W2*FW2)
		RI(5) = RI(5) + RL*RPP*X*DXX

		X = FPW/(WP*FW2)
		RI(6) = RI(6) + RPX2*X*DXX
	END DO

	IF (ILLW.NE.0) THEN
		X = -2.0*RH(NN)*WP/FW2
		RI(1) = RI(1) + X*0.5*RLB*RP2
		RI(1) = RI(1) + X*RLL*RP2*0.5/RLB
	END IF
	RI(1) = REPS*RI(1)

		RI(3) = RI(3)*REPS

	IF (ILLW.NE.0) THEN
		X = FPW/(WP*FW2)
		RI(4) = RI(4) + RLL*X*0.5*RP2/RLB
	END IF
	RI(4) = RI(4)*REPS

	X = RL*4.0*REPS*F/FW2
	RI(5) = RI(5) -X*RP2*0.5

	X = FPW/(WP*FW2)
	RI(6) = RI(6) + X*RP2*0.5*RLB
		
	DO I = 1,3
		AA = AA + RI(I)
	END DO
	
	BB = 0.
	DO I = 4,6
		BB = BB + RI(I)
	END DO
	RTF = -AA/BB
	WRITE(IOPF,901) RTF
	GO TO 300

200	RTF = 0.

901	FORMAT(//' R*TF (CM)=',E15.5//)

300	RETURN
	END
