
c		newest version of program to compute free coastal
c		trapped waves. Modified for DOS use, 9/1/92
C		MODIFIED TO DO GRAPHICS OUTPUT FOR MATLAB USE 11/3/93
C			Updated for direct Matlab use (cctwplt.m) 9/1/97
C			Dtime updated 8/25/98
C			NEEDS MATLAB (4.0 or better) to do graphics 
C			Updated 12/1/2005 for more modern usage, deleted ICCM option
C			
C		OUTPUT (12) FILES FOR GRAPHICS:
C			PGR.grd = PRESSURE
C			VGR.grd = V
C			UGR.grd = U
C			RGR.grd = rho
C			DISP.grd = dispersion curve


	PROGRAM HYBD

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DIMENSION WW(3), RII(3)
	DIMENSION R(NN)
	DIMENSION BXB(NM)
	DIMENSION WH(30), CH(30),RLH(30),DTH(30)
	CHARACTER*15 FINP,FOT
C
C		MAIN PROGRAM TO CALL ALL OF THE OTHER PIECES
C
C		REVISED MARCH 1,1988
C		REVISED NOVEMBER 19, 1991 WITH MORE ACCURATE CONDITION
C			AT X = 0
C		Updated 6/30/97, 12/16/01
C	

C		Set up the input and output files
	WRITE(6,970)
	READ(5,*) FINP
	WRITE(6,971)
	READ(5,*) FOT
	OPEN (INPF,FILE=FINP)
	OPEN (IOPF,FILE=FOT)

c	To test this program, use input file "TESTBIGC.IN"



c	Input file (unit inpf) contents:
c
c	eps	est	dd1
c	ncalm	nitm	isd
c	iggg
c	f	xmax
c	beta	ilww
c	ncal	wh(1)
c	idiag
c	rlf	drl	
c	ww(1)	ww(2)	ww(3)
c	nrx						(read in from subroutine DEP)
c	x	dep
c	x	dep
c	:	:
c	:	:
c	nr 	dzr	alph			(read in from subroutine NSQ)
c	cmlt
c	N2
c	N2
c	:
c	:
c	nr						(read in from subroutine LGWH)
c	x	r
c	x	r
c	:	:
c	:	:
c
c
c
c	Where
c
c		EPS = nominal fractional accuracy of solution
c		EST = fractional search increment
c		DD1 = 1. for free surface
c		DD1 = 0. for rigid lid
c		NCALM = number of points on each curve
c		NITM =	maximum number of iterations to find a frequency
c		ISD		controls the search direction. 0 is unbiased
C			ISD = 0  NORMAL
C			ISD = 1  SEARCH DOWN in frequency
C			ISD = -1 SEARCH UP in frequency
c		IGGG = 1 for graphics output, = 0 for no graphics output
c		F =		Coriolis parameter (scaled by 1e-05 sec-1)
c		XMAX =	offshore distance to outer boundary (km)
c		Beta =	gradient of planetary vorticity (cm^-1 sec^-1)
c						perpendicular to the coast
c		ILWW = 1:	General frequency and wavenumber
c		ILWW = 0:	Long wave approximation made
c
c		[ If you are continuing a dispersion curve from an earlier
c			calculation, use these NCAL~=1, WH(1) = something
c		If not, use NCAL = 1, WH(1) anything
c		If continuing a curve, say NCAL = 2 and WH(1) = last
c			frequency computed scaled by 1e-06 sec-1 ]
c
c		IDIAG = 0  for no diagnostic calculations
c		RLF =	first scaled wavenumber to be used, 
c		DRL =	wavenumber increment
c					Both Scaled by 1e-07 cm-1
c		WW are first 3 guesses at scaled frequency
c					Scaled by 1e-06 sec-1
c		NRX is the number of (x,depth) pairs to read)
c		x,	depth 		pairs (NRX lines) in km, m
c		NR =	number of N2 values to read
c		dzr =	depth increment (m) for the N2 data
c		alph =	exponential tail length to extrapolate N2 data (km)
c		CMLT =	conversion factor to get into rad^2/sec^2
c		N2	N2 values (NR lines) in your units
c		NR = number of (x,r) pairs to read
c		x, r = offshore distance, bottom resistance coefficient 
c					(nr lines) (km, cm/sec)
c
c
c	Note: If you want to change the grid dimensions NN (offshore # of
c			 points), MM (vertical number of points), you have to change the
c			 PARAMETER statment in the main program and all of the subroutines.



	READ(INPF,*) EPS, EST, DD1
	READ(INPF,*) NCALM, NITM, ISD
	READ(INPF,*) IGGG


	READ(INPF,*) F,XMAX

	XMAX = XMAX*1.0E+05
	F =F*1.0E-05
	WRITE(IOPF,907) F, XMAX
	WRITE(IOPF,915) NN, MM
	WRITE(IOPF,905) EPS, EST
	READ(INPF,*) BETA, ILWW

	IF (ILWW.EQ.0) THEN
		NCALM = 1
		WRITE(IOPF,910)
	ELSE
		ILWW = 1
		WRITE(IOPF,911)
	END IF
	DDLW = FLOAT(ILWW)


	IF (DD1.GT.0.99) THEN
		WRITE(IOPF,914)
		DD1 = 1.
	ELSE
		WRITE(IOPF,913)
		DD1 = 0.
	END IF
	
		
	WRITE(IOPF,909) BETA

	READ(INPF,*) NCAL, WH(1)
	READ(INPF,*) IDIAG
	READ(INPF,*) RLF, DRL

	WH(1) = WH(1)*1.0E-06
	RLF = RLF*1.0E-07
	DRL = DRL*1.0E-07
	

	DX = XMAX/FLOAT(NN-1)
	DT = 1./FLOAT(MM-1)

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

	READ(INPF,*) (WW(J),J=1,3)
	DO J = 1,3
		WW(J) = WW(J)*1.0E-06
	END DO

	CALL DEP
	CALL NSQ
c		Both DEP (depth profile) and NSQ (stratification)
c			 read in information


10	RIB = 0.
	

C		Do the first three iterations
	DO I = 1,3
		W = WW(I)
		CALL MATS(RL,W,DD1,RI,BETA,DDLW)
		RII(I) = RI
	
		IF (RI.GT.RIB) THEN
			RIB = RI
			WB = W
			IB = I
			DO IJ = 1,NM
	 			BXB(IJ) = BX(IJ)
			END DO
		END IF


		IF ((I.EQ.2).AND.(RII(1).GT.RII(2))) THEN
				IF (ISD.NE.(-1)) THEN
					WW(3) = WW(1)*(1.-EST)
				END IF
		END IF
	
	END DO


	NIT = 3
	IGP = 0

20	CALL NGSW(WW,RII,WB,WN,ISUC,EPS,IN,IGP,EST,ISD)
	
C		Quit if frequency is superinertial
	IF (WN.GT.F) GO TO 140
C		Quit if search failure
	IF (ISUC.EQ.2) GO TO 140
C		Stop iterating if this has converged
	IF (ISUC.EQ.1) GO TO 100
	
	CALL MATS(RL,WN,DD1,RI,BETA,DDLW)
	
	IF (IN.EQ.0) THEN
		IF (RI.GT.RIB) THEN
			IF (WN.LT.WW(2)) THEN
				IN = 3
			ELSE
				IN = 1
			END IF
		ELSE
			IF (WN.LT.WW(2)) THEN
				IN = 1
			ELSE
				IN = 3
			END IF
		END IF	
	END IF
	

	WW(IN) = WN
	RII(IN) = RI
	IF (RI.GT.RIB) THEN
		WB = WN
		RIB = RI
c			Save the best estimate
		DO I = 1,NM
			BXB(I) = BX(I)
		END DO
	END IF

	NIT = NIT +1
	IF(NIT.LT.NITM) GO TO 20

C		Ran out of iterations
	WRITE(IOPF,902) NIT
	GO TO 140

100	RN = SQRT(RIB)
	RN = RN*BXB(MM)/ABS(BXB(MM))

	DO I = 1,NM
		BX(I) = BXB(I)/RN
	END DO

	CC = WB/RL
	WRITE(IOPF,903) WB,RL,CC,EPS,NIT

C		Assure sign convention
	IF (BX(MM).LT.0) THEN
		DO I = 1,NM
			BX(i) = -BX(i)
		END DO
	END IF

	CALL DTIME(WB,RL,DD1,DDLW,DTHH)

	IF (NCAL.EQ.1) THEN
		CALL LGWH(WB,RL,DD1,R,BETA,DDLW,IGGG,vke,spe)
c			LGWH reads in information about structure of bottom friction

		IF (IDIAG.NE.0) THEN
			CALL UCAL(WB,RL,BETA,DDLW,IGGG,uke)
			CALL RHOC(DD1,IGGG,rpe)
			RAT = (VKE + UKE)/(RPE/SPE)
			WRITE(IOPF,912) RAT
		END IF
	END IF	


	WRITE(IOPF,908)
C
C		DO P GRAPHICS OUTPUT
	IF (IGGG.NE.0) THEN
		WRITE(6,972)
		READ(5,*) IPX
		IF (IPX.NE.0) THEN
			OPEN (12, FILE = 'PGR.grd')

			DO N = 1,NN
				X = DX*FLOAT(N-1)		
				DO M = 1,MM			
					I = M + MM*(N-1)
					ZZZ = (-1.0 + DT*FLOAT(M-1))*RH(N)
					WRITE(12,*) X, ZZZ, BX(I)	
				END DO
			END DO
			CLOSE (12)
		END IF
	END IF


	DO N = 1,NN
		X = DX*FLOAT(N-1)
		DZ = RH(N)*DT
		WRITE(IOPF,904) X,RH(N),DZ
		ML = 1 +MM*(N-1)
		MH = MM*N
		WRITE(IOPF,901) (BX(M),M=ML,MH)
	END DO

	CALL DIAG(WB,RL,DDLW)

	WH(NCAL) = WB
	CH(NCAL) = CC
	DTH(NCAL) = DTHH

	NCAL = NCAL +1
C		Did all the required points on the dispersion curve
	IF (NCAL.GT.NCALM) GO TO 140
	
	RL = RLH(NCAL)
C		Estimate first guesses at next frequency on the
C			 dispersion curve
	IF (NCAL.EQ.2) 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.0 - EST)
	WW(3) = WW(2)*(1.0 + EST)

C		Check that the next guesses are OK
	IF (WW(1).LT.0.) GO TO 140
	IF (WW(3).GT.F) GO TO 140
	
	GO TO 10

140	IF (IGGG.NE.0) THEN
		WRITE(6,973)
		READ(5,*) IPX
		IF (IPX.NE.0) THEN 
			OPEN (12, FILE = 'DISP.grd')
			DO I = 1,NCALM
				WRITE(12,*) RLH(I), WH(I), CH(I), DTH(I)
			END DO
			CLOSE (12)
		END IF
	END IF

	WRITE(IOPF,916)
	DO I = 1,NCALM
		WRITE(IOPF,917) RLH(I), WH(I), CH(I), DTH(I)
	END DO


	CLOSE(INPF)
	CLOSE(IOPF)

901	FORMAT(2X,10E12.5)
902	FORMAT(//'  USED UP ',I5,'  ITERATIONS'//)
903	FORMAT(///'   CONVERGED:  W,L,C,EPS,NIT =',4E15.5,I10//)
904	FORMAT(/'  X,H,DZ =',3E15.5)
905	FORMAT(/'  EPS,EST =', 2E15.5/)
906	FORMAT(3F10.5)
907	FORMAT(/'  F,XMAX =',2E15.5//)
908	FORMAT(//'  PRESSURE'/)
909	FORMAT(/'  BETA= ',E15.5)
910	FORMAT('  LONG WAVE LIMIT')
911	FORMAT('  GENERAL FREQUENCY AND WAVENUMBER')
912	FORMAT(/' KE/PE =',e15.5)
913	FORMAT('  Rigid lid')
914	FORMAT('  Free surface')
915	FORMAT('  NN, MM =',2i5)
916	FORMAT(//'  RL           W           C           rTf'/)
917	FORMAT(4E12.5)
970	FORMAT(' Input file?')
971	FORMAT(' Output file?')
972	FORMAT(' DO YOU WANT PRESSURE GRAPHICS? 0 (for no) OR 1 (for
     & yes)')
973	FORMAT(' Do you want a dispersion curve for graphics? 0 (for no)
     & or 1 (for yes)')

	STOP
	END




	SUBROUTINE NGSW(WW,RII,WB,WN,ISUC,EPS,IN,IGP,EST,ISD)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DIMENSION  WW(3),RII(3)
C
C		SUBROUTINE TO COMPUTE THE NEXT GUESS AT W
C
C		ISD = 0  NORMAL
C		ISD = 1  SEARCH DOWN
C		ISD = -1 SEARCH UP

	DEL = EST

5	IC = 0

	DO I = 1,2
		I1 = 1 +I
		IF (WW(I1).LT.WW(I)) THEN
			IC = 1
			RX = RII(I)
			WX = WW(I)
			WW(I) = WW(I1)
			WW(I1) = WX
			RII(I) = RII(I1)
			RII(I1) =RX
		END IF
	END DO

	IF (IC.NE.0) GO TO 5


	IF (RII(3).GT.RII(2)) THEN
		ISUC = 0
		IGP = 0
		IF (ISD.EQ.1) THEN
			WN = WW(1)*(1-DEL)
			IN = 3
		ELSE IF (ISD.EQ.-1) THEN
			WN = WW(3)*(1 + DEL)
			IN = 1
		ELSE IF (RII(2).GT.RII(1)) THEN
			WN = WW(3)*(1 + DEL)
			IN = 1
		ELSE IF	(RII(1).GT.RII(3)) THEN
			WN = WW(1)*(1 - DEL)
			IN = 3
		ELSE 
			WN = WW(3)*(1 + DEL)
			IN = 1
		END IF
	END IF

	IF (RII(1).GT.RII(2)) THEN
		ISUC = 0
		IGP = 0
		IF (ISD.EQ.-1) THEN
			WN = WW(3)*(1+DEL)
			IN = 3
		ELSE IF (ISD.EQ.-1) THEN
			WN = WW(3)*(1 + DEL)
			IN = 1
		ELSE IF (RII(2).GT.RII(3)) THEN
			WN = WW(3)*(1 - DEL)
			IN = 3
		ELSE IF	(RII(1).GT.RII(3)) THEN
			WN = WW(1)*(1 - DEL)
			IN = 3
		ELSE 
			WN = WW(3)*(1 + DEL)
			IN = 1
		END IF
	END IF
	

	IF ((RII(2).GE.RII(1)).AND.(RII(2).GE.RII(3))) THEN

C		If you get here, the peak is localized
		ISD = 0
		IN = 0
		WB1 = (WW(1) +WW(2))/2.
		WB2 = (WW(2) + WW(3))/2.
		RW1 = (RII(2) - RII(1))/(WW(2) -WW(1))
		RW2 = ( RII(3) - RII(2))/(WW(3) - WW(2))
		A = (RW2 -RW1)/(WB2 -WB1)
		B = RW1 - A*WB1
		WN = -B/A

		IF (WN.LE.WW(1)) THEN
			WN = WB1
		ELSE IF (WN.GT.WW(3)) THEN
			WN = WB2
		END IF

		EP = ABS(WN - WB)/WB
		IF ((EP.LE.EPS).AND.(IGP.GT.0)) THEN
			ISUC = 1
		ELSE IF ((EP.LE.EPS).AND.(IGP.LT.1)) THEN
			IGP = IGP + 1
			ISUC = 0
			IF (RII(3).GT.RII(1)) THEN
				WN = WB1
			ELSE
				WN = WB2
			END IF
		END IF
	END IF


	IF (WN.LE.0.) THEN
		ISUC = 2
	END IF

	RETURN
	END




	SUBROUTINE NSQ

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
C
C		SUBROUTINE TO READ AND INTERPOLATE BUOYANCY FREQUENCY
C			SQUARED PROFILE BV, AND Z DERIVATIVE BVZ
C

c		Read in the Stratification information
	READ(INPF,*) NR,DZR,ALPH

c		NR = number of Nsq values to read
c		DZR = depth increment in m
c		ALPH = exponential decay scale for Nsq below depth of last
c			value (km)
c		CMLT = value to mutiply read Nsq by to get rad^2/sec^2
	ALPH = ALPH*1.0E+05
	DZR = DZR*100.
	MX = MM -1
	READ(INPF,*) CMLT

	DO I = 1,NR
		READ(INPF,*) XL(I)
		XL(I) = XL(I)*CMLT
	END DO

c		XL are the values of Nsq (user units)
c		Can have no more than NM lines of input

C		Extrapolate the N^2 information with the exponential tail
	TF = XL(NR)
	ZZZ = DZR*FLOAT(NR-1)
	DZI = RH(NN) - ZZZ
	NL = NR +1
	NH = 600
	DO I = NL,NH
		Z = DZR*FLOAT(I-1)
		XL(I) = TF*EXP((ZZZ -Z)/ALPH)
	END DO


	DO 200 N = 1,NN
		DD = RH(N)
		DZ = DT*DD
		IF (DZ.LE.DZR) THEN
C				If the N^2 information is coaser than the grid, fill in the 
C					model N^2 by interpolation

			BV(N,MM) = XL(1)
			DO M = 1,MX
				Z = DD - DZ*FLOAT(M-1)
				IBXL = 1 + IFIX(Z/DZR)
				IBXH = IBXL +1
				ZS = DZR*FLOAT(IBXL -1)
				BV(N,M) = XL(IBXL) +(Z-ZS)*(XL(IBXH)-XL(IBXL))/DZR
			END DO
			Z = DD - DZ
			IBXL = 1 + IFIX(Z/DZR)
			IBXH = IBXL + 1
			ZS = DZR*FLOAT(IBXL-1)
			AQ = XL(IBXL) +(Z-ZS)*(XL(IBXH)-XL(IBXL))/DZR
			
		ELSE
C				If the N^2 information is finer than the model grid, fill in by
C					averaging the input information		
			ZD = DZ/2.
			XBB = 0.
			NAVG = 0
			DO I = 1,NR
				ZC = DZR*FLOAT(I-1)
				IF (ZC.LE.ZD) THEN
					XBB = XBB + XL(I)
					NAVG = NAVG +1
				END IF
			END DO
			BV(N,MM) = XBB/FLOAT(NAVG)

			DO 140 MQ = 1,MM
				M = MQ -1
				Z = DD - DZ*FLOAT(M-1)
				ZS = Z - DZ/2.
				ZD = ZS + DZ
				XBB = 0.
				NAVG = 0
				DO I = 1,NH
					ZC = DZR*FLOAT(I-1)
					IF ((ZC.GE.ZS).AND.(ZC.LE.ZD)) THEN
						XBB = XBB + XL(I)
						NAVG = NAVG + 1
					END IF
				END DO
		
				IF (NAVG.EQ.0) THEN
					IBXL = 1 + IFIX(Z/DZR)
					IBXH = IBXL + 1
					ZSS = DZR*FLOAT(IBXL-1)
			
					IF (M.NE.0) THEN
						BV(N,M) = XL(IBXL) +(Z-ZSS)*(XL(IBXH)-XL(IBXL))/DZR
					ELSE
						AQ = XL(IBXL) +(Z-ZSS)*(XL(IBXH)-XL(IBXL))/DZR
					END IF
			
				ELSE
					IF (M.EQ.0) THEN
						AQ = XBB/FLOAT(NAVG)
					ELSE
						BV(N,M) = XBB/FLOAT(NAVG)
					END IF
				END IF

140			CONTINUE
		END IF

		DO M = 2,MX
			IP = M +1
			IM = M -1
			BVZ(N,M) = (BV(N,IP) - BV(N,IM))/(2.*DZ)
		END DO

		BVZ(N,1) = (BV(N,2) -AQ)/(2.0*DZ)
		BVZ(N,MM) = (BV(N,MM) - BV(N,MX))/DZ
200	CONTINUE

	WRITE(IOPF,901)
	WRITE(IOPF,902) (BV(NN,J),J=1,MM)

901	FORMAT(//'   NSQUARED AT XMAX '/)
902	FORMAT(2X,10E12.5)

	RETURN
	END





	SUBROUTINE MATS(RL,W,DD1,RI,BETA,DDLW)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)

C
C		SUBROUTINE TO SET UP AND SOLVE THE PRESSURE EQUATION
C


c		Calculate a bunch of constants
	GG = 1./980.
	DXDT = DX/(2.*DT)
	DXX = DX*DX
	DDD = DXX/(DT*DT)
	RLL = RL*RL*DDLW
	DDX = DX*DXDT
	DQ = 2.0*DX
	FLW = F*RL/W
	BLW=BETA*RL/W
	CCC1 = W*(F*F - DDLW*W*W)/DXX
	CCC2 = 0.5*F*(2.0*BETA*W + RL*(F*F - DDLW*W*W))/DX
	CCC3 = BETA*RL*(F*F + DDLW*W*W)
	AA5 = 2.0*F*BETA/(F*F -W*W*DDLW)
	RLQ = DDLW*RLL - BLW*(F*F + DDLW*W*W)/(F*F -DDLW*W*W)


C		Initialize matrix to zero
	DO I = 1,NM
		BX(I) = 0.
	END DO

	DO J = 1,NMX
		DO I = 1,NM
			AX(I,J) = 0.
		END DO
	END DO
	BX(39) = 1.



	J = NN + 2
	IP1 = J + 1
	IM1 = J - 1
	IPM1 = 2*NN +1
	IPM = 2*NN + 2
	IPMM = 2*NN + 3

	DO 100 M = 1,MM
		DO 90 N = 1,NN
			I = N + NN*(M-1)
			CALL AS(A1,A2,A3,W,N,M,C1,C2,C3,DDLW)
			AX(I,1) = A2*DXDT
			AX(I,IM1) = 1.-AA5*DX/2.0
			AX(I,IPM1) = -A2*DXDT
			AX(I,2) = A1*DDD - DDX*(A3+A2*AA5)
			AX(I,J) = -2. -2.*DDD*A1 - RLQ*DXX
			AX(I,IPM) = A1*DDD + DDX*(A3+A2*AA5)
			AX(I,3) = -A2*DXDT
			AX(I,IP1) = 1.+AA5*DX/2.0
			AX(I,IPMM) = A2*DXDT
		
		
			IF ( N.EQ.1) THEN
			
C 					X = 0 BOUNDARY
			
				AX(I,IPM1) = 0.
				AX(I,3) = 0.
				AX(I,IPMM) = 0.
				AX(I,IP1) = 2.
C	
				B11 = DQ*(-FLW + RHX(N)/RH(N))
				B12 = -A2*2.0*DQ/DT
				AX(I,IPM) = AX(I,IPM) + AX(I,IM1)*A2*2.0*DXDT +AX(I,1)*(B11+B12)
				AX(I,J) = AX(I,J) + AX(I,IM1)*(DQ*FLW) + AX(I,1)*(A2*16.*DXDT)
				AX(I,2) = AX(I,2) -AX(I,IM1)*A2*2.*DXDT + AX(I,1)*(-B11 + B12)
				AX(I,IM1) = 0.
				AX(I,1) = 0.
C
				IF (M.EQ.MM) THEN
					D5 = DD1*RH(N)*BV(N,M)*GG
					AX(I,2) = AX(I,2) + AX(I,IPM)
					AX(I,J) = AX(I,J) - 2.0*DT*D5*AX(I,IPM)
					AX(I,IPM) = 0.
				ELSE IF (M.EQ.1) THEN
					AX(I,IPM) = AX(I,IPM) + AX(I,2)
					AX(I,2) = 0.
				END IF
		
			ELSE IF (N.EQ.NN) THEN 
				
C					X = XMAX BOUNDARY
			
				DZ = DT*RH(N)
				AX(I,3) = 0.
				AX(I,IPMM) = 0.
				AX(I,1) = 0.
				AX(I,IPM1) = 0.
				AX(I,J) = AX(I,J) +AX(I,IP1)*(2.0*CCC1-CCC3)/(CCC1+CCC2)
				AX(I,IM1) = AX(I,IM1) +AX(I,IP1)*(CCC2-CCC1)/(CCC1+CCC2)
				AX(I,IP1) = 0.
		
				IF (M.EQ.MM) THEN
					AX(I,2) = AX(I,2) +AX(I,IPM)
					AX(I,J) = AX(I,J) -2.0*DZ*DD1*GG*BV(N,M)*AX(I,IPM)
					AX(I,IPM) = 0.
				ELSE IF (M.EQ.1) THEN
					AX(I,IPM) = AX(I,IPM) + AX(I,2)
					AX(I,2) = 0.
				END IF
		
	
			ELSE IF (M.EQ.1) THEN

C					Z = -H(X) BOUNDARY
			
				D1 = 2.*C2*DT/(C1 +C2*C3)
				AX(I,IPM) = AX(I,IPM) + AX(I,2)
				AX(I,IP1) = AX(I,IP1) + D1*AX(I,2)/DQ
				AX(I,IM1) = AX(I,IM1) -D1*AX(I,2)/DQ
				AX(I,J) = AX(I,J) +D1*FLW*AX(I,2)
				NNN = N +1
				IJJP = J + 2
				CALL AS(A1,A2,A3,W,NNN,M,C1,C2,C3,DDLW)
				D1 = 2.0*C2*DT/(C1 + C2*C3)
				AX(I,IPMM) = AX(I,IPMM) + AX(I,3)
				AX(I,IJJP) = AX(I,IJJP) + D1*AX(I,3)/DQ
				AX(I,J) = AX(I,J) -D1*AX(I,3)/DQ
				AX(I,IP1) = AX(I,IP1) + D1*FLW*AX(I,3)
		
				IF (N.EQ.2) THEN
					AX(I,IPM1) = AX(I,IPM1) + AX(I,1)
				ELSE
			
					IJJM = J - 2
					NNN = N - 1
					CALL AS(A1,A2,A3,W,NNN,M,C1,C2,C3,DDLW)
					D1 = 2.0*C2*DT/(C1 + C2*C3)
					AX(I,IPM1) = AX(I,IPM1) + AX(I,1)
					AX(I,J) = AX(I,J) + D1*AX(I,1)/DQ
					AX(I,IJJM) = AX(I,IJJM) -D1*AX(I,1)/DQ
					AX(I,IM1) = AX(I,IM1) + D1*FLW*AX(I,1)
				END IF

				AX(I,1) = 0.
				AX(I,2) = 0.
				AX(I,3) = 0.


			ELSE IF (M.EQ.MM) THEN
			
C					Z = 0 BOUNDARY
	
				AX(I,1) = 0.
				AX(I,3) = 0.
				D5 = RH(N)*DD1*BV(N,M)*GG
				AX(I,2) = AX(I,2) + AX(I,IPM)
				AX(I,J) = AX(I,J) - D5*2.0*DT*AX(I,IPM)
				AX(I,IM1) = AX(I,IM1) -D5*2.0*DT*AX(I,IPM1)
				AX(I,IP1) = AX(I,IP1) -D5*DT*2.0*AX(I,IPMM)
				AX(I,IPM1) = 0.0
				AX(I,IPMM) = 0.0
				AX(I,IPM) = 0.0
			END IF

90		CONTINUE

100	CONTINUE

C		Now solve the matrix equation and evaluate it's resonance	
	CALL BANDG
	CALL CALI(RI)

	CC = W/RL
	WRITE(IOPF,901) W,RL,CC,RI

901	FORMAT(/'  W,L,C,RI =',4E15.5)

	RETURN
	END



	SUBROUTINE AS(A1,A2,A3,W,N,M,C1,C2,C3,DDLW)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
C
C		SUBROUTINE TO COMPUTE COEFFICIENTS FOR MATRIX
C

	F2 = F*F
	W2 = W*W*DDLW
	TZ = 1./RH(N)
	TH = -1. +DT*FLOAT(M-1)
	TX = -RHX(N)*TH*TZ
	AA = (RHX(N)/RH(N))**2
	TXX = (2.*AA -RHXX(N)/RH(N))*TH
	BW = BV(N,M)
	BWW = BW*BW
	A1 = TX*TX + TZ*TZ*(F2-W2)/BW
	A2 = TX
	A3 = TXX -TZ*BVZ(N,M)*(F2-W2)/BWW
	C1 = TZ
	C2 = RHX(N)*BV(N,1)/(F2-W2)
	C3 = C1*RHX(N)

	RETURN
	END




	SUBROUTINE DEP

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
C
C		SUBROUTINE TO READ AND INTERPOLATE DEPTH PROFILE
C			RH = DEPTH
C			RHX = X DERIVATIVE OF DEPTH
C			RHXX = SECOND X DERIVATIVE OF DEPTH
C

	READ(INPF,*) NRX

c 		NRX = number of depth values to read
c		XL = x in km
c		BX = depth in m
c		Can not read in more than NM lines

	DO I =1,NRX
		READ(INPF,*) XL(I),BX(I)
		BX(I) = BX(I)*100.
		XL(I) = XL(I)*1.0E+05
	END DO
	RHMA = BX(NRX)
	
		
	RH(1) = BX(1)

	DO 20 N = 2,NN
		X = DX*FLOAT(N-1)
		IF (X.GT.XL(NRX)) THEN
			RH(N) = RHMA
		ELSE
			IC = 0
			DO J = 2,NRX
				IF (IC.EQ.0) THEN
					I = J
					IF (X.LE.XL(I)) THEN
						IC = I
					END IF
				END IF
			END DO
			IM = I-1
			RHX(N) = (BX(I) -BX(IM))/(XL(I)-XL(IM))
			XX = X - XL(IM)
			RH(N) = BX(IM) + RHX(N)*XX
		END IF
20	CONTINUE


	RHX(1) = (RH(2) - RH(1))/DX
	RHXX(1) = 0.
	D2 = 2.*DX
	DXX = DX*DX
	NX = NN -1
	DO N = 2,NX
		IP = N +1
		IM = N -1
		RHX(N) =(RH(IP) - RH(IM))/D2
		RHXX(N) = (RH(IP) -2.*RH(N) + RH(IM))/DXX
	END DO
	RHX(NN) = 0.
	RHXX(NN) = 0.

	RETURN
	END




	SUBROUTINE CALI(RIX)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DOUBLE PRECISION RI,RZ
C
C		SUBROUTINE TO CALCULATE THE INTEGRAL OF RESPONSE SQUARED
C
C		FIRST, REORDER THE MATRIX

	DO M = 1,MM
		DO N = 1,NN
			IN = N + NN*(M -1)
			II = M + MM*(N-1)
			XL(II) = BX(IN)
		END DO
	END DO


	DO I = 1,NM
		BX(I) = XL(I)
	END DO
	RI = 0.
	RR = 0.
	RZ = 0.

	DO N = 1,NN
		IF ((N.EQ.1).OR.(N.EQ.NN)) THEN
			DXX = DX/2.
		ELSE
			DXX = DX
		END IF

		DO  M = 1,MM
			DZ = DT*RH(N)
			IF ((M.EQ.1).OR.(M.EQ.MM)) THEN
				DZZ = DZ/2.
			ELSE
				DZZ = DZ
			END IF
			I = (N-1)*MM +M
			RZ = RZ + DZZ*DXX
			RI = RI + DZZ*DXX*BX(I)**2
		END DO
	END DO

	RI = RI/RZ
	RIX = RI

	RETURN
	END



	SUBROUTINE DIAG(W,RL,DDLW)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
C
C		SUBROUTINE TO COMPUTE MOMENTUM DIAGNOSTICS
C

	N = 5
	I = N*MM
	IP = I + MM
	IM = I -MM
	PX = (BX(IP) -BX(IM))/(2.*DX)
	F2 = F*F - W*W*DDLW
	VT = W*(F*PX + DDLW*RL*W*BX(I))/F2
	FU = -F*(RL*F*BX(I) + W*PX)/F2
	PY = RL*BX(I)
	X = DX*FLOAT(N-1)
	WRITE(IOPF,901) X
	WRITE(IOPF,902) VT,FU,PY

901	FORMAT(//'  Y MOMENTUM TERMS AT Z = 0, X =',E15.5)
902	FORMAT(/'  VT,FU,PY =',3E15.5/)

	RETURN
	END



	SUBROUTINE LGWH(W,RL,DD1,R,BETA,DDLW,IGGG,vke,spe)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DOUBLE PRECISION RD,RP,RI,XINT,XP
	DIMENSION R(NN)
C
C		SUBROUTINE TO COMPUTE LONG WAVE PARAMETERS BN,ANN
C

	XINT = 0.
	RHO = 1.03
	F2 = F*F
	MX = MM - 1

C		NORMALIZE
	NX = NN -1
	RI = (BX(1)**2 + BX(MM)**2)/2.
	RD = (RHX(1)*BX(1)**2)/2.
	DO I = 2,MX
		RI = RI + BX(I)**2
	END DO

	DO I = 2,NX
		IJ = MM*(I-1) +1
		RD = RD + RHX(I)*BX(IJ)**2
	END DO

	RD = RD*DX
	RI = RI*RH(1)*DT
	RD = (RD + RI)/F
	RD = DSQRT(RD*RD)
	RD = DSQRT(RD)

	DO I = 1,NM
		BX(I) = BX(I)/RD
	END DO



	CALL VCAL(W,RL,BETA,DDLW,IGGG,vke)

C		COMPUTE WIND COUPLING

	BN = BX(MM)/F

C		NOW GET BOTTOM FRICTION
C		READ R(X) (bottom resistance coefficient)
	READ(INPF,*) NR

	DO I = 1,NR
		READ(INPF,*) AX(I,1), AX(I,2)

c			NR = number of r values to read (< NM)
c			AX(i,1) =  x in km
c			AX(i,2) =  r in cm/sec

		AX(I,1) = AX(I,1)*1.0E+05
	END DO

	R(1) = AX(1,2)
	RMA = AX(NR,2)

C		Interpolate bottom friction
	DO 50 N = 2,NN
		X = DX*FLOAT(N-1)
		IF (X.GT.AX(NR,1)) THEN
			R(N) = RMA
		ELSE
			IC = 0
			DO J = 2,NR
				IF (IC.EQ.0) THEN
					I = J
					IF (X.LE.AX(I,1)) THEN
						IC = I
					END IF
				END IF
			END DO
			IM = I -1
			RX = (AX(I,2) -AX(IM,2))/(AX(I,1)-AX(IM,1))
			XX = X - AX(IM,1)
			R(N) = AX(IM,2) + RX*XX
		END IF
50	CONTINUE

C		Now compute ann
C			X=0 CONTRIBUTION
	RI = 0.5*R(1)*(RL*F*BX(1)/W)**2

C			X = L CONTRIBUTION
	I = 1 + MM*(NN-1)
	BXP = F*XL(I)
	RI = RI + 0.5*R(NN)*(BXP)**2

C			INTERMEDIATE X
	DO N = 2,NX
		I = 1 + MM*(N-1)
		IP = I + MM
		IM = I - MM
		RX = 0.5*(BX(IP)-BX(IM))/DX
		RI = RI + R(N)*RX*F*XL(I)
	END DO

	ANN = RI*DX/F2
C
	WRITE(IOPF,902) BN,ANN
	WRITE(IOPF,901)
	WRITE(IOPF,905) (R(I),I=1,NN)

C		Compute free surface contribution to PE
	CC = 0.5/(980.*RHO)
	DO N = 1,NN
		IF ((N.EQ.1).OR.(N.EQ.NN)) THEN
			XX = DX/2
		ELSE
			XX = DX
		END IF
		I = N*MM
		XINT = XINT + CC*XX*BX(I)*BX(I)
	END DO
	WRITE(IOPF,904) XINT
	SPE = XINT

901	FORMAT(/'   R (cm/sec) ='/)
902	FORMAT('   BN ([cm sec]^-1/2), ANN (cm^-1) =',2E15.5///)
903	FORMAT(2F10.5)
904	FORMAT(//'  FREE SURFACE CONTRIBUTION TO PE =',E15.5/)
905	FORMAT(2X,10E12.5)

	RETURN
	END	



	SUBROUTINE VCAL(W,RL,BETA,DDLW,IGGG,VKE)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DOUBLE PRECISION XINT
C
C		SUBROUTINE TO COMPUTE THE V FIELD OF THE WAVE
C

	RHO = 0.515
	XINT = 0.
	WRITE(IOPF,905)
	FW = F*F - DDLW*W*W
	WL = DDLW*W*RL
	RLW = RL/W
	FLW = F*RL/W
	MX = MM -1
	DXX = 2.*DX
	W2 = W*W*DDLW
	
C		Do you want graphics printout?
	IF (IGGG.NE.0) THEN
		WRITE(6,970)
		READ(5,*) IPX
		IF (IPX.NE.0) THEN
			OPEN (IGPF, FILE = 'VGR.grd')
		END IF
	END IF

	DO 50 N = 1,NN
		XX = DX
		X = DX*FLOAT(N-1)
		D = RH(N)
		DD = RHX(N)/D
		DZ = DT*D
		RHZ = RHX(N)**2
	
		IF (N.EQ.1) THEN
C				Coastal boundary
			DO M = 1,MM
				I = (N-1)*MM +M
				XL(I) = -RLW*BX(I)
			END DO
			XX = DX/2.
		
		ELSE IF (N.EQ.NN) THEN
C				Offshore boundary	
			CC1 = W*FW/(DX*DX)
			CC2 = 0.5*(2.0*F*BETA*W +RL*F*FW)/DX
			CC3 = BETA*RL*(F*F +DDLW*W*W)
			DQ = 1./(F*DX)
			DO M = 1,MM
				I = (N-1)*MM +M
				IM = I -MM
				XP = BX(I)*(2.0*CC1 -CC3)/(CC1+CC2)
				XP = XP + BX(IM)*(CC2-CC1)/(CC1+CC2)
				XPX = (XP - BX(IM))/(2.0*DX)
				XL(I) = (F*XPX + WL*BX(I))/FW
			END DO
			XX = DX/2.
		
		ELSE
C				Interior points
			I = (N-1)*MM +1
			IP = I + MM
			IM = I - MM
			XX = DX
			CALL AS(A1,A2,A3,W,N,1,C1,C2,C3,DDLW)
			CQ = C2*C3/(C1+ C2*C3)
			XPX = (1.0 -CQ)*(BX(IP) -BX(IM))/DXX
			XPX = XPX -CQ*FLW*BX(I)
			XL(I) = (WL*BX(I) + F*XPX)/FW
			I = N*MM
			IP = I + MM
			IM = I -MM
			XPX = (BX(IP) - BX(IM))/DXX
			XL(I) =(WL*BX(I) + F*XPX)/FW
			DO M = 2,MX
				I = (N-1)*MM +M
				IP = I + MM
				IM = I -MM
				I2 = I +1
				I1 = I -1
				XPX = (BX(IP) - BX(IM))/DXX
				T = -1. +DT*FLOAT(M-1)
				XPT = (BX(I2) -BX(I1))/(2.*DT)
				XPX = XPX - T*DD*XPT
				XL(I) = (WL*BX(I) + F*XPX)/FW
			END DO
		END IF	

C			Compute VKE integral
		DO M = 1,MM
			I = (N-1)*MM + M
			IF ((M.EQ.1).OR.(M.EQ.MM)) THEN
				ZZ = DZ/2.
			ELSE
				ZZ = DZ
			END IF
			XINT = XINT + RHO*XL(I)*XL(I)*XX*ZZ
		END DO

		WRITE(IOPF,904) X,D,DZ
		IL = (N-1)*MM +1
		IH = N*MM
		WRITE(IOPF,901) (XL(IJ),IJ=IL,IH)

C			Graphics output if desired
		IF ((IGGG.NE.0).AND.(IPX.NE.0)) THEN
			X = DX*FLOAT(N-1)
			DO M = 1,MM
				ZZZ = (-1.0 + DT*FLOAT(M-1))*RH(N)
				I = (N-1)*MM + M
				WRITE(12,*) X, ZZZ, XL(I)		
			END DO
		END IF
	
50	CONTINUE

	WRITE(IOPF,902) XINT
	VKE = XINT

	IF ((IGGG.NE.0).AND.(IPX.NE.0)) THEN
		CLOSE (IGPF)
	END IF

901	FORMAT(2X,10E12.5)
902	FORMAT(//'  V CONTRIBUTION TO KE =',E15.5/)
905	FORMAT(////'  V'//)
904	FORMAT(/'    X, H, DZ =',3E15.5)
970	FORMAT(' DO YOU WANT V GRAPHICS? 0 (no) OR 1 (yes)')

	RETURN
	END



	SUBROUTINE RHOC(DD1,IGGG,RPE)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DOUBLE PRECISION XINT
C
C		SUBROUTINE TO COMPUTE THE DENSITY FIELD OF THE WAVE
C

	G2 = 980.*980./2.06
	XINT = 0.
	WRITE(IOPF,903)
	G = 980.

	IF (IGGG.NE.0) THEN
		WRITE(6,970)
		READ(5,*) IPX
		IF (IPX.NE.0) THEN
			OPEN (IGPF, FILE = 'RGR.grd')		
		END IF
	END IF

	DO 50 N = 1,NN
	
		IF ((N.EQ.1).OR.(N.EQ.NN)) THEN
			DXX = DX/2.
		ELSE
			DXX = DX
		END IF
		X = DX*FLOAT(N-1)
		DZ = DT*RH(N)
		GDZ = G*DZ

		DO M = 1,MM
			DZZ = DZ
			I = (N-1)*MM +M
			IF (M.EQ.1) THEN
				IP = I + 1
				DZZ = DZ/2.
				XL(M) = -(BX(IP) - BX(I))/GDZ
			ELSE IF (M.EQ.MM) THEN
				DZZ = DZ/2.
				XL(M) = BV(N,M)*BX(I)*DD1/(G*G)
			ELSE
				IP = I +1
				IM = I -1
				XL(M) = -(BX(IP) -BX(IM))/(2.*GDZ)
			END IF
			XX = G2*XL(M)*XL(M)/BV(N,M)
			XINT = XINT + XX*DZZ*DXX
		END DO

		WRITE(IOPF,901) X,RH(N),DZ
		WRITE(IOPF,902) (XL(M),M=1,MM)

		IF ((IGGG.NE.0).AND.(IPX.NE.0)) THEN
			X = DX*FLOAT(N-1)
			DO M = 1,MM
				ZZZ = (-1.0 + DT*FLOAT(M-1))*RH(N)
				WRITE(IGPF,*) X, ZZZ, XL(M)	
			END DO
		END IF

50	CONTINUE	

	WRITE(IOPF,904) XINT
	RPE = XINT

	IF ((IGGG.NE.0).AND.(IPX.NE.0)) THEN
		CLOSE(12)
	END IF

901	FORMAT(/'  X,D,DZ =',3E15.5)
902	FORMAT(2X,10E12.5)
903	FORMAT(///'  RHO'/)
904	FORMAT(//'  RHO CONTRIBUTION TO PE =',E15.5/)
970	FORMAT(' DO YOU WANT RHO GRAPHICS? 0 (for no) OR 1 (for yes)')

  	RETURN
	END



	SUBROUTINE UCAL(W,RL,BETA,DDLW,IGGG,UKE)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM, NMX =2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DOUBLE PRECISION XINT
C
C		SUBROUTINE TO COMPUTE THE U FIELD OF THE WAVE
C

	RHO = 0.515
	XINT = 0.
	FW = F*F - DDLW*W*W
	FL = F*RL

	WRITE(IOPF,903)

	IF (IGGG.NE.0) THEN
		WRITE(6,970)
		READ(5,*) IPX
		IF (IPX.NE.0) THEN
			OPEN (IGPF, FILE = 'UGR.grd')		
		END IF
	END IF

	DO 100 N = 1,NN
		X = DX*FLOAT(N-1)
		DZ = DT*RH(N)
	
		IF (N.EQ.1) THEN
C				Onshore
			DO I = 1,MM
				XL(I) = 0.
			END DO
			DXX = DX/2.
		ELSE IF (N.EQ.NN) THEN
C				Offshore	
			CC1 = W*FW/(DX*DX)
			CC2 = 0.5*(2.0*F*BETA*W +RL*F*FW)/DX
			CC3 = BETA*RL*(F*F +DDLW*W*W)
			DXX = DX/2.

			DO M = 1,MM
				I = (NN-1)*MM +M
				IM = I - MM
				XPIM = BX(I)*(2.0*CC1-CC3)/(CC1+CC2)
				XPIM = XPIM + BX(IM)*(CC2-CC1)/(CC1+CC2)
				XL(M) =-(FL*BX(I) + W*0.5*(XPIM - BX(IM))/DX)/FW
			END DO
		ELSE 
C				Interior points
			DD = RHX(N)/RH(N)
			DXX = DX

			DO M = 1,MM
				I = (N-1)*MM +M
				IP = I + MM
				IM = I -MM
				T = -1. + DT*FLOAT(M-1)
				IF (M.EQ.1) THEN
					I1 = I +1
					XPT = (BX(I1) - BX(I))/DT
				ELSE IF (M.EQ.MM) THEN
					XPT = 0.
				ELSE
					I2 = I+1
					I1 = I-1
					XPT = (BX(I2) - BX(I1))/(2.*DT)
				END IF
				XPX = (BX(IP) - BX(IM))/(2.*DX)
				XPX = XPX -T*DD*XPT
				XL(M) = -(W*XPX + FL*BX(I))/FW
			END DO
		END IF
		
C			Contribute to UKE integral
		DO M = 1,MM
			IF ((M.EQ.1).OR.(M.EQ.MM)) THEN
				DZZ = DZ/2.
			ELSE
				DZZ = DZ
			END IF
			XINT = XINT + RHO*XL(M)*XL(M)*DXX*DZZ
		END DO

		WRITE(IOPF,901) X,RH(N),DZ
		WRITE(IOPF,902) (XL(M),M=1,MM)
	
	
		IF ((IGGG.NE.0).AND.(IPX.NE.0)) THEN
			X = DX*FLOAT(N-1)
			DO M = 1,MM
				ZZZ = (-1.0 + DT*FLOAT(M-1))*RH(N)
				WRITE(12,*) X, ZZZ, XL(M)
			END DO
		END IF


100	CONTINUE

	WRITE(IOPF,904) XINT
	UKE = XINT

	IF ((IGGG.NE.0).AND.(IPX.NE.0)) THEN
		CLOSE(IGPF)
	END IF

901	FORMAT(/'  X,D,DZ =',3E15.5)
902	FORMAT(2X,10E12.5)
903	FORMAT(///'   U'/)
904	FORMAT(//'  U CONTRIBUTION TO KE =',E15.5/)
970	FORMAT(' DO YOU WANT U GRAPHICS? 0 (for no) OR 1 (for yes)')

	RETURN
	END


	SUBROUTINE BANDG

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM,NMX =2*NN + 3)
	DOUBLE PRECISION A,BB,XL
	COMMON A(NM,NMX),BB(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DOUBLE PRECISION R
C
C			BAND GAUSSIAN ELIMINATION
C

	MBW = NN + 1
	NH = NM -1
	MDIAG =MBW + 1

	DO 50 N = 1,NH
		IP = N +1
		MH = N + MBW
		IF (MH.GT.NM) THEN
			MH = NM
		END IF
		DO IR = IP,MH
			ICD = MBW + N + 1 -IR
			R = A(IR,ICD)/A(N,MDIAG)
			BB(IR) = BB(IR) -R*BB(N)
			DO IC = IP,MH
				ICD = MBW + IC + 1 -IR
				ICB = MBW + IC + 1 - N
				A(IR,ICD) = A(IR,ICD) -R*A(N,ICB)
			END DO
		END DO
50	CONTINUE
		
	DO 100 I = 1,NH
		N = NM - I + 1
		BB(N) =BB(N)/A(N,MDIAG)
		IL = N - MBW
		IH = N - 1
		IF (IL.LT.1) THEN
			IL = 1
		END IF
		DO IR = IL,IH
			ICD = MBW + N + 1 - IR
			BB(IR) = BB(IR) -A(IR,ICD)*BB(N)
		END DO
100	CONTINUE

	BB(1) = BB(1)/A(1,MDIAG)

	RETURN
	END



	SUBROUTINE DTIME(W,RL,DD1,DDLW,RTF)

	PARAMETER (NN=55,MM=50,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM=NN*MM,NMX=2*NN +3)
	DOUBLE PRECISION AX,BX,XL
	COMMON AX(NM,NMX),BX(NM),XL(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DOUBLE PRECISION RI1,RI2,RI3,RI4,RI5


c		Calculates perturbation decay time assuming r = constant

	RI1 = 0.
	RI2 = 0.
	RI3 = 0.
	RI4 = 0.
	RI5 = 0.
	RI6 = 0.

	FLW = F*RL/W
	WW2 = -2.0*W*W*DDLW
	WW2G = (WW2/980.)*DD1
	FW = F*F - DDLW*W*W
	FPW = F*F + DDLW*W*W
	RLL = RL*RL
	DT2 = 2.0*DT
	DX2 = DX*2.0

	DO 100 N = 1,NN
		DXX = DX
		RH2 = RH(N)*RH(N)
		DO M   = 1,MM
			DTT = DT*RH(N)
			I = M  + MM*(N-1)
			IP1 = I + 1
			IM1 = I - 1
			IPP = I + MM
			IMM = I - MM

			IF (M.EQ.1) THEN 
				DTT = DT*RH(N)/2.0
				IF (N.EQ.1) THEN
					RPT = 0.
					RPBX = -FLW*BX(I)
					DXX = DX/2.
					RPX = RPBX
				ELSE IF (N.EQ.NN) THEN
					CC1 = W*FW/(DX*DX)
					CC2 = 0.5*RL*F*FW/DX
					DQ = 1.0/(F*DX)
					XP = BX(I)*(2.0*CC1)/(CC1 + CC2)
					XP = XP + BX(IMM)*(CC2 -CC1)/(CC1 + CC2)
					RPX = (XP -BX(IMM))/DX2
					RPBX = RPX
					RPT = 0.
				ELSE
					RPBX = (BX(IPP) -BX(IMM))/DX2
					XX = -RHX(N)*(FLW*BX(I) + RPBX)
					YY = (FW/BV(N,M) + RHX(N)*RHX(N))/RH(N)
					RPT = XX/YY
					RPX = RPBX + RHX(N)*RPT/RH(N)
				END IF

				RI1 = RI1 + DXX*RHX(N)*BX(I)*BX(I)
				RI5 = RI5 + (RPBX*RPX + DDLW*RLL*BX(I)*BX(I))*DXX
				RI6 = RI6 + BX(I)*(RPX + RPBX)*DXX
			

			ELSE IF (M.EQ.MM) THEN
				DTT = DT*RH(N)/2.0
				RPT = -DD1*RH(N)*BV(N,M)*BX(I)/980.
				IF (N.EQ.1) THEN
					DXX = DX/2.
				END IF
				RI3 = RI3 + BX(I)*BX(I)*DXX
			ELSE
				RPT = (BX(IP1)-BX(IM1))/DT2	
				IF (N.EQ.1) THEN
					DXX = DX/2.
					RI2 = RI2 +DTT*BX(I)*BX(I)
				END IF			
			END IF

			RI4 = RI4 + (RPT*RPT/(BV(N,M)*RH2))*DXX*DTT
		END DO
100	CONTINUE

	RI1 = RI1*FLW
	RI2 = RI2*FLW
	RI3 = RI3*WW2G
	RI4 = RI4*WW2
	RI5 = RI5*FPW/FW
	RI6 = RI6*DDLW*2.0*W*RL*F/FW

	RTF = (RI1 + RI2 + RI3+ RI4)/(RI5 +RI6)
	WRITE(IOPF,901) RTF

901	FORMAT(//' Assuming r is constant,  r*Tf (cm)=',E15.5//)

	RETURN
	END