C			ISLAND TRAPPED WAVES WITH STRATIFICATION AND
C			ARBITRARY BOTTOM TOPOGRAPHY  1/24/95
c				Modified 6/24/97 to include coupling
c				coefficients, etc.
c				modified 8/29/97 to have outputs
c				compatable with citwplt.m
c			Has output for Matlab (4.0 for Windows or better)
c				graphics
C			USE citwplt.m to grid and plot fields
C			FOR GRAPHICS OUTPUT (IF USED)
C				PGR.grd = PRESSURE
C				UGR.grd = U
C				VGR.grd = V
C				RGR.grd = RHO


		PROGRAM ISLTW
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
C			NN = NUMBER OF RADIAL GRID POINTS
C			MM = NUMBER OF VERICAL GRID POINTS
C			INPF = INPUT UNIT
C			IOPF = OUTPUT UNIT
C			IGPF = GRAPHICS OUTPUT UNIT
		PARAMETER (NM=NN*MM,NMX = 2*NN + 3)

	DOUBLE PRECISION AX,BX
	COMMON AX(NM,NMX),BX(NM)
	COMMON F,DX,DT
	COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
	DOUBLE PRECISION XX(NM)
	DIMENSION WW(3), RII(3)
	DIMENSION R(NN)
	DIMENSION BXB(NM)
	DIMENSION WH(10),NRH(10)
	CHARACTER*20 INP,IOP,IGQ
C
C		MAIN PROGRAM TO CALL ALL OF THE OTHER PIECES
C
C		DIMENSION(XL) = NM*(NN+2)
C		DIM(AX) = NM,2NN+3
C		


		WRITE(6,911)
		READ(5,*) INP
		WRITE(6,912)
		READ(5,*) IOP
		OPEN (7,FILE=INP)
		OPEN (8,FILE=IOP)


	READ(INPF,*) EPS,EST

C		EPS = NOMINAL FREQUENCY ACCURACY OF SOLUTION
C		EST = FRACTIONAL SEARCH INCREMENT FOR FREQUENCY

	READ(INPF,*) ICCM,NCALM,NITM,ISD

C		ICCM = NUMBER OF DISPERSION CURVES TO CALCULATE
C		NCALM = NUMBER OF POINTS ON EACH DISPERSION CURVE
C		NITM = MAXIMUM NUMBER OF ITERATIONS FOR A
C			FREQUENCY BEFORE GIVING UP         
C		ISD = 0 UNBIASED SEARCH
C		ISD = 1 SEARCH ONLY TOWARD LOWER FREQUENCIES
C		ISD = -1 SEARCH ONLY TOWARD HIGHER FREQUENCIES

	ICCC =1
	READ(INPF,*) IWS,IGG

C			IGG = 0, NO OUTPUT FOR GRAPHICS
C			IWS = 0  SHORT PRINTOUTS
C			IWS NOT = 0 FULL PRINTOUTS

3	READ(INPF,*) F,XMAX

C			F = CORIOLIS PARAMETER
C			XMAX = MAXIMUM RADIUS OF GRID

	XMAX = XMAX*1.0E+05
	F =F*1.0E-05

2	READ(INPF,*) NCAL,WH(1)

C			NCAL NOT 1 ONLY IF YOU ARE CONTINUING A
C				DISPERSION CURVE
C			SEE BIGLOAD DOCUMENTATION FOR USAGE

	READ(INPF,*) NRF,NRD

C			NRF = FIRST AZIMUTHAL WAVE NUMBER ON DISPERSION
C				CURVE
C			NRD = INCREMENT IN AZIMUTHAL WAVENUMBER FOR
C				DISPERSION CURVE

	WRITE(IOPF,914) NN,MM
	WH(1) = WH(1)*1.0E-06
	WRITE(IOPF,905) EPS,EST


	DT = 1./FLOAT(MM-1)
	DO 50 I = 1,NCALM
		NRH(I) = NRF + NRD*FLOAT(I-1)
50	Continue

	NR = NRH(NCAL)
	READ(INPF,*) (WW(J),J=1,3)

C			WW ARE THREE INITIAL GUESSES AT FIRST FREQUENCY
C				TO BE CALCULATED

	DO 6 J = 1,3
		WW(J) = WW(J)*1.0E-06
6	Continue

	CALL DEP(rz,xmax)
	WRITE(IOPF,907) F,rz,XMAX
	CALL NSQ
	Read(INPF,*) wms
c
c		WMS = wind decay length scale for distorted winds (km)
c		WMS < 0.5 km -> uniform winds
c
1	RIB = 0.
	I =1
15	W = WW(I)

	CALL MATS(NR,W,RI,rz)

	RII(I) = RI
	IF (RI.LT.RIB) GO TO 18
	RIB = RI
	WB = W
	IB = I
	DO 10 IJ = 1,NM
	 	BXB(IJ) = BX(IJ)
10	Continue
18	I = I+1
	IF (I.EQ.2) GO TO 15
	IF (I.EQ.4) GO TO 19
	IF (RII(1).LT.RII(2)) GO TO 15
	WW(3) = WW(1)*(1. -EST)
	GO TO 15
19	NIT = 3
	IGP = 0
20	CALL NGSW(WW,RII,WB,WN,ISUC,EPS,IN,IGP,EST,ISD)
	IF (WN.GT.F) GO TO 140
	IF (ISUC.EQ.2) GO TO 140
	IF (ISUC.EQ.1) GO TO 100

	CALL MATS(NR,WN,RI,rz)

	IF (IN.NE.0) GO TO 24
	IF (RI.GT.RIB) GO TO 21
	IF (WN.LT.WW(2)) GO TO 23
22	IN = 3
	GO TO 24
23	IN = 1
	GO TO 24
21	IF (WN.LT.WW(2)) GO TO 22
	GO TO 23
24	WW(IN) = WN
	RII(IN) = RI
	IF (RI.LT.RIB) GO TO 30
	WB = WN
	RIB = RI
	DO 25 I = 1,NM
		BXB(I) = BX(I)
25	Continue
30	NIT = NIT +1
	IF(NIT.LT.NITM) GO TO 20
	WRITE(IOPF,902) NIT
	GO TO 140
100	RN = SQRT(RIB)
	DO 110 I = 1,NM
		BX(I) = BXB(I)/RN
110	Continue
	
	Call Norm(rz)
	WRITE(IOPF,903) WB,NR,EPS,NIT
	CALL DCALC(NR,WB,rz)
	Call Coup(wb,rz,wms,nr)

	IF (NCAL.NE.1) GO TO 125

120		CALL VCAL(WB,NR,XX,RIV,IWS,IGG,rz)
		CALL UCAL(WB,NR,XX,RIU,IWS,IGG,rz)
		CALL RHOC(XX,WB,NR,RIR,IWS,IGG,rz)
		RKE = RIV + RIU
		RIT = RKE + RIR
		WRITE(IOPF,909) RIU,RIV,RIR
		RRR = RKE/RIR
		WRITE(IOPF,910) RRR

125		DO 230 N = 1,NN
		DO 230 M = 1,MM
			I = N + NN*(M-1)
			II = M + MM*(N-1)
			XX(II) = BX(I)
230		CONTINUE

		IF (IGG.EQ.0) GO TO 127
		write(6,935)
		read(5,*) ipx
		if (ipx.eq.0) go to 127
		open (12, file = 'pgr.grd')

		DO 126 N = 1,NN
			xxxx = rz + dx*float(n-1)
			do 126 j = 1,mm
				dept = (-1.0 +DT*float(j-1))*rh(n)
				m = mm*(n-1) + j		
				WRITE(12,*) xxxx,dept,xx(m)
126		continue

		close (12)
127		IF (IWS.EQ.0) GO TO 950
   	WRITE(IOPF,908)
	DO 130 N = 1,NN
		X = rz + 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) (XX(M),M=ML,MH)
130	CONTINUE

950	WH(NCAL) = WB
	NCAL = NCAL +1
	IF (NCAL.GT.NCALM) GO TO 140
		NR = NRH(NCAL)
	IF (NCAL.GE.3) GO TO 200
		WW(2) = WB*NR/NRH(1)
	GO TO 205
200	I1 = NCAL - 2
	I2 = NCAL - 1
	CG = (WH(I2) -WH(I1))/(NRH(I2) -NRH(I1))
	WW(2) = WH(I2) + CG*(NRH(NCAL) -NRH(I2))
205	WW(1) = WW(2)*(1.0 - EST)
	WW(3) = WW(2)*(1.0 + EST)
	IF (WW(1).LT.0.) GO TO 140
	IF (WW(3).GT.F) GO TO 140
	GO TO 1
140	ICCC = ICCC +1
	IF (ICCC.LE.ICCM) GO TO 3

	CLOSE(7)
	CLOSE(8)
	close(12)
	close(13)

901	FORMAT(2X,10E12.5)
902	FORMAT(//'  USED UP ',I5,'  ITERATIONS'//)
903	FORMAT(//' CONVERGED W,N,EPS,NIT =',E15.5,I5,E15.5,I10)
904	FORMAT(/'  R,H,DZ =',3E15.5)
905	FORMAT(/'  EPS,EST =', 2E15.5/)
906	FORMAT(3F10.5)
907	FORMAT(//'  F, rz, RMAX =',3E15.5//)
908	FORMAT(//'  PRESSURE'/)
909	FORMAT(//' UKE,VKE,PE =',3E15.5)
910	FORMAT(/' KE/PE =',E15.5)
911	FORMAT(' INPUT FILE?')
912	FORMAT(' OUTPUT FILE?')
913	FORMAT(' GRAPHICS FILE?')
914	format(' number of radial, vertical grid points =',2i5)
935	format(' Do you want to plot pressure? 0 or 1')

	STOP
	END

	SUBROUTINE NGSW(WW,RII,WB,WN,ISUC,EPS,IN,IGP,EST,ISD)
	PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM =NN*MM, NMX = 2*NN +3)

	DOUBLE PRECISION AX,BX
	COMMON AX(NM,NMX),BX(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 10 I = 1,2
		I1 = 1 +I
		IF (WW(I1).GT.WW(I)) GO TO 10
		IC = 1
		RX = RII(I)
		WX = WW(I)
		WW(I) = WW(I1)
		WW(I1) = WX
		RII(I) = RII(I1)
		RII(I1) =RX
10	CONTINUE

	IF (IC.NE.0) GO TO 5
	IF (RII(3).GT.RII(2)) GO TO 150
	IF (RII(1).GT.RII(2)) GO TO 160
	ISD = 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)) GO TO 50
	IF (WN.GT.WW(3)) GO TO 45
	GO TO 20
50	WN = WB1
	GO TO 20
45	WN = WB2
20	EP = WN - WB
	EP = ABS(EP)
	EP = EP/WB
	IF (EP.LE.EPS) GO TO 100
	IGP = IGP +1
	ISUC = 0
	IN = 0
	GO TO 120
100	ISUC = 1
	IF (IGP.NE.0) GO TO 120
	IGP = IGP +1
	ISUC = 0
	IF (RII(3).GT.RII(1)) GO TO 110
	WN = WB2
	GO TO 20
110	WN = WB1
	GO TO 20
150	IF (RII(1).GT.RII(2)) GO TO 180
155	IF (ISD.EQ.1) GO TO 165
	WN = WW(3)*(1. + DEL)
	ISUC = 0
	IGP = 0
	IN = 1
	GO TO 120
160	IF (RII(3).GT.RII(2)) GO TO 180
165	IF (ISD.EQ.-1) GO TO 155
	WN = WW(1)*(1. - DEL)	
	ISUC = 0
	IGP = 0
	IN = 3
	GO TO 120
180	IF (RII(3).GT.RII(1)) GO TO 155
	GO TO 165
120	IF (WN.LE.0.) GO TO 125
	GO TO 130
125	ISUC = 2

130	RETURN
	END

	SUBROUTINE NSQ
	PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
	PARAMETER (NM = NN*MM , NMX = 2*NN +3)

	DOUBLE PRECISION AX,BX
	COMMON AX(NM,NMX),BX(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
	READ(INPF,*) NR,DZR,ALPH

C		NR = NUMBER OF NSQUARED VALUES TO READ
C		DZR = DEPTH INCREMENT
C		ALPH = EXPONENTIAL TAIL LENGTH IN KM.
C		SEE DOCUMENTATION FOR BIGLOAD

	ALPH = ALPH*1.0E+05
	DZR = DZR*100.
	MX = MM -1
	READ(INPF,*) CMLT

c		CMLT = constant to multiply Nsquared values by 
c			to get rad^2/sec^2

	DO 20 I = 1,NR
		READ(INPF,*) BX(I)
		BX(I) = BX(I)*CMLT
20	Continue
	TF = BX(NR)
	ZZZ = DZR*FLOAT(NR-1)
	DZI = RH(NN) - ZZZ
	NL = NR +1
	NH = 425
	DO 30 I = NL,NH
	Z = DZR*FLOAT(I-1)
30	BX(I) = TF*EXP((ZZZ -Z)/ALPH)

	DO 200 N = 1,NN
	DD = RH(N)
	DZ = DT*DD
	IF (DZ.GT.DZR) GO TO 110
	BV(N,MM) = BX(1)
	DO 50 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) = BX(IBXL) +(Z-ZS)*(BX(IBXH)-BX(IBXL))/DZR
50	Continue

	Z = DD - DZ
	IBXL = 1 + IFIX(Z/DZR)
	IBXH = IBXL + 1
	ZS = DZR*FLOAT(IBXL-1)
	AQ = BX(IBXL) +(Z-ZS)*(BX(IBXH)-BX(IBXL))/DZR
	GO TO 145
110	ZD = DZ/2.
	XBB = 0.
	NAVG = 0

	DO 120 I = 1,NR
		ZC = DZR*FLOAT(I-1)
		IF (ZC.GT.ZD) GO TO 120
		XBB = XBB + BX(I)
		NAVG = NAVG +1
120	CONTINUE

	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
125	XBB = 0.
	NAVG = 0

	DO 130 I = 1,NH
		ZC = DZR*FLOAT(I-1)
		IF ( ZC.LT.ZS) GO TO 130
		IF (ZC.GT.ZD) GO TO 130
		XBB = XBB + BX(I)
		NAVG = NAVG + 1
130	CONTINUE

	IF (NAVG.NE.0) GO TO 135
	IBXL = 1 + IFIX(Z/DZR)
	IBXH = IBXL + 1
	ZSS = DZR*FLOAT(IBXL-1)
	IF (M.EQ.0) GO TO 133
	BV(N,M) = BX(IBXL) +(Z-ZSS)*(BX(IBXH)-BX(IBXL))/DZR
	GO TO 140
133	AQ = BX(IBXL) +(Z-ZSS)*(BX(IBXH)-BX(IBXL))/DZR
	GO TO 140
135	IF (M.EQ.0) GO TO 138
	BV(N,M) = XBB/FLOAT(NAVG)
	GO TO 140
138	AQ = XBB/FLOAT(NAVG)
140	CONTINUE

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

	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 RMAX '/)
902	FORMAT(2X,10E12.5)
250	RETURN
	END

		SUBROUTINE MATS(NR,W,RI,rz)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX = 2*NN + 3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
C
C			SUBROUTINE TO SET UP AND SOLVE MATRIX EQUATION
C
		DO 10 I = 1,NM
			B(I) = 0.
			DO 10 J = 1,NMX
				AX(I,J) = 0.
10		Continue

		B(39) = 1.0E-10
		DRT4 = 4.0*DR*DT
		DR2 = 2.0*DR
		DRR = DR*DR
		DT2 = 2.0*DT
		DTT = DT*DT
		DTDR = DT/DR
		FN2 = FLOAT(NR*NR)
		FW = F*F - W*W
		FNW = FLOAT(NR)*F/W
		RMAX = DR*FLOAT(NN)
		BB3 = (1.0 + FNW)/RMAX
		X = 1.0/DR + BB3/2.0
		CC1 = 2.0/(DR*X)
		CC2 = -(1.0/DR-0.5*BB3)/X

   		J = NN + 2
		IP1 = J + 1
		IM1 = J - 1
		IPNM = 2*NN + 1
		IPN = IPNM + 1
		IPN1 = IPN + 1
		IPP1 = J + 2
		IMM1 = J - 2


		DO 100 N = 1,NN
		R = rz +DR*FLOAT(N-1)
		AA1 = R
		AA2 = -FN2/R

		DO 100 M = 1,MM
		I = N + NN*(M-1)
		CALL AS(A1,A2,A3,N,M)
		AA3 = 2.0*A1*R
		AA4 = R*A2 + A1 -R*BVZ(N,M)*A3*FW/(BV(N,M)**2)
		AA5 = R*A1*A1 + R*FW*A3*A3/BV(N,M)
		
		AX(I,1) = AA3/DRT4
		AX(I,2) = -AA4/DT2 + AA5/DTT
		AX(I,3) = -AA3/DRT4
		AX(I,IM1) = AA1/DRR -1.0/DR2
		AX(I,J) = -2.0*AA1/DRR + AA2 -2.0*AA5/DTT
		AX(I,IP1) = AA1/DRR + 1.0/DR2
		AX(I,IPNM) = -AA3/DRT4
		AX(I,IPN) = AA4/DT2 + AA5/DTT
		AX(I,IPN1) = AA3/DRT4

		IF (M.EQ.1) GO TO 20
		IF (M.EQ.MM) GO TO 30
		IF (N.EQ.1)  GO TO 40
		IF (N.EQ.NN) GO TO 50
		GO TO 100

C			BOTTOM B.C.
20		CALL BCAL(BB1,BB2,N,FNW,FW,rz)
		IF (N.EQ.NN) GO TO 22
		if (n.eq.1) go to 26
		AX(I,IPN) = AX(I,IPN) + AX(I,2)
		AX(I,IP1) = AX(I,IP1) -DTDR*BB1*AX(I,2)
		AX(I,IM1) = AX(I,IM1) + DTDR*BB1*AX(I,2)
		AX(I,J) = AX(I,J) -DT2*BB2*AX(I,2)
		AX(I,2) = 0.
		NX = N + 1
		CALL BCAL(BB1,BB2,NX,FNW,FW,rz)
		AX(I,IPN1) = AX(I,IPN1) + AX(I,3)
		AX(I,IPP1) = AX(I,IPP1) -DTDR*BB1*AX(I,3)
		AX(I,J) = AX(I,J) +DTDR*BB1*AX(I,3)
		AX(I,IP1) = AX(I,IP1) -DT2*BB2*AX(I,3)
		AX(I,3) = 0.

		NX = N -1
		CALL BCAL(BB1,BB2,NX,FNW,FW,rz)
		AX(I,IPNM) = AX(I,IPNM) + AX(I,1)
		AX(I,J) = AX(I,J) -DTDR*BB1*AX(I,1)
		AX(I,IMM1) = AX(I,IMM1) + DTDR*BB1*AX(I,1)
		AX(I,IM1) = AX(I,IM1) -DT2*BB2*AX(I,1)
		AX(I,1) = 0.
		IF (N.NE.2) GO TO 100
		ax(i,j) = ax(I,j) + ax(i,imm1)
		axx = 2.*dr*fnw/rz
		ax(i,im1) = ax(I,im1) + axx*ax(i,imm1)
		AX(I,IMM1) = 0.
		GO TO 100
C			BOTTOM OFFSHORE
22		AX(I,J) = AX(I,J) + CC1*AX(I,IP1)
		AX(I,IM1) = AX(I,IM1) + CC2*AX(I,IP1)
		AX(I,IPN) = AX(I,IPN) + AX(I,2)
		AX(I,1) = 0.
		AX(I,2) = 0.
		AX(I,3) = 0.
		AX(I,IP1)= 0.
		AX(I,IPN1) = 0.
		GO TO 100
C			BOTTOM ,INNER
26		a1t = -rhx(1)/rh(1)
		c1 = 2.*dr*(a1t + fnw/r)
		c2 = 4.*a1/dtdr
		c3 = 8.*a1/dtdr
		ccc1 = c1 + c2
		ccc2 = c1 - c2
		ax(i,ipn1) = ax(i,ipn1) - ax(i,1)
		ax(i,3) = ax(i,3) + ax(i,1)
		ax(i,ipn) = ax(i,ipn) -ccc1*ax(i,1)
		ax(i,2) = ax(i,2) + ccc2*ax(i,1)
		ax(i,j) = ax(i,j) + c3*ax(i,1)
		c1 = 2.*dr*fnw/r
		ax(i,ip1) = ax(i,ip1) + ax(i,im1)
		ax(i,j) = ax(i,j) +c1*ax(i,im1)
		ax(i,ipn) = ax(i,ipn) + ax(i,2)
		ax(i,ipn1) = ax(i,ipn1) + ax(i,3)
		ax(i,2) = 0.
		ax(i,3) = 0.

		AX(I,1) = 0.
		AX(I,IM1) = 0.
		AX(I,IPNM) = 0.
		AX(I,IMM1) = 0.
		GO TO 100
			
C			SURFACE B.C.
30		AX(I,1) = AX(I,1) + AX(I,IPNM)
		AX(I,2) = AX(I,2) + AX(I,IPN)
		AX(I,3) = AX(I,3) + AX(I,IPN1)
		AX(I,IPNM) = 0.
		AX(I,IPN) = 0.
		AX(I,IPN1) = 0.
		IF (N.EQ.1) GO TO 32
		IF (N.EQ.NN) GO TO 36
		GO TO 100
C			TOP,INNER
32		ax(i,ip1) = ax(i,ip1) + ax(i,im1)
		c1 = 2.*dr*fnw/r
		ax(i,j) = ax(i,j) + c1*ax(i,im1)
		AX(I,1) = 0.
		AX(I,IM1) = 0.
		GO TO 100
C			TOP,OUTER
36		AX(I,J) = AX(I,J) + CC1*AX(I,IP1)
		AX(I,IM1) = AX(I,IM1) + CC2*AX(I,IP1)
		AX(I,IP1) = 0.
		AX(I,3) = 0.
		GO TO 100

C			INNER B.C.
40		a1t = -rhx(1)/rh(1)
		c1 = 2.*dr*(a1t + fnw/r)
		c2 = 4.*a1/dtdr
		c3 = 8.*a1/dtdr
		ccc1 = c1 + c2
		ccc2 = c1 - c2
		ax(i,ipn1) = ax(i,ipn1) -ax(i,1)
		ax(i,3) = ax(i,3) + ax(i,1)
		ax(i,ipn) = ax(i,ipn) -ccc1*ax(i,1)
		ax(i,2) = ax(i,2) + ccc2*ax(i,1)
		ax(i,j) = ax(i,j) +c3*ax(i,1)
		c1 = a1/dtdr
		c2 = 2.*dr*fnw/r
		ax(i,ip1) = ax(i,ip1) + ax(i,im1)
		ax(i,ipn) = ax(i,ipn) + c1*ax(i,im1)
		ax(i,2) = ax(i,2) -c1*ax(i,im1)
		ax(i,j) = ax(i,j) + c2*ax(i,im1)

		AX(I,IM1) = 0.
		AX(I,IPNM) =0.
		AX(I,1) = 0.
		GO TO 100

C			OUTER B.C.
50		AX(I,J) = AX(I,J) + CC1*AX(I,IP1)
		AX(I,IM1) = AX(I,IM1) + CC2*AX(I,IP1)
		AX(I,3) = 0.
		AX(I,IP1) = 0.
		AX(I,IPN1) = 0.
100		CONTINUE

		CALL BANDG


		CALL CALI(RI,rz)
		WRITE(IOPF,901) W,NR,RI

901		FORMAT(/' W,N,RI =',E15.5,I5,E15.5)
		RETURN
		END

		SUBROUTINE AS(A1,A2,A3,N,M)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX = 2*NN + 3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
C			SUBROUTINE TO COMPUTE COEFFICIENTS FOR MATS
		TH = -1.0 + DT*FLOAT(M-1)
		A3 = 1.0/RH(N)
		A1 = -RHX(N)*TH/RH(N)
		A2 = -RHXX(N)*TH/RH(N) + 2.0*TH*(RHX(N)/RH(N))**2

		RETURN
		END

		SUBROUTINE BCAL(BB1,BB2,N,FNW,FW,rz)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX = 2*NN +3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
C			ANOTHER SUBROUTINE TO COMPUTE COEFFICIENTS FOR 
C				MATS
		CALL AS(A1,A2,A3,N,1)
		R = rz +DR*FLOAT(N-1)
		X = BV(N,1)*RHX(N)/FW
		BB1 = -X/(A3 +X*A1)
		BB2 = FNW*BB1/R

		RETURN
		END

		SUBROUTINE CALI(RI,rz)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,igpf=12)
		PARAMETER (NM = NN*MM, NMX = 2*NN + 3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
		
C			SUBROUTINE TO COMPUTE DOUBLE INTEGRAL OF
C				PSQUARED AS A MEASURE OF WHETHER
C				YOU ARE NEAR RESONANCE
		RI = 0.

		DO 50 M = 1,MM
			DTT = DT*DR
			IF (M.EQ.1) GO TO 5
			IF (M.NE.MM) GO TO 10
5			DTT = DR*DT/2.0
10			DO 50 N = 1,NN
				R = rz +DR*FLOAT(N-1)
				DZ = DTT*RH(N)
				if (n.ne.1) go to 40
				dz = dz/2.
40				I = N + NN*(M-1)
				RI = RI + B(I)*B(I)*DZ*R
50		CONTINUE

		RETURN
		END

		SUBROUTINE VCAL(W,NR,XX,RI,IWS,IGG,rz)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX =2*NN + 3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
		DOUBLE PRECISION XX(NM)
C			SUBROUTINE TO COMPUTE AZIMUTHAL VELOCITY V
		
		DR2 = 2.0*DR
		DT2 = 2.0*DT
		WRITE(IOPF,901)
		DAA = DT*DR
		RI = 0.
		FW = F*F - W*W
		FNW = FLOAT(NR)*F/W
		WN = W*FLOAT(NR)
		BB3 = (1.0 + FNW)/(DR*FLOAT(NN))
		X = 1.0/DR + 0.5*BB3
		CC1 = 2.0/(X*DR)
		CC2 = -(1.0/DR -BB3*0.5)/X

		IF (IGG.EQ.0) GO TO 5
		WRITE(6,910)
		READ(5,*) IPX
		IF (IPX.EQ.0) GO TO 5
		OPEN (12, FILE = 'VGR.grd')		

5		DO 60 N = 1,NN
		R = rz +DR*FLOAT(N-1)
		WNR = WN/R
		
		DO 50 M = 1,MM
		CALL AS(A1,A2,A3,N,M)
		I = N + NN*(M-1)
		IP1 = I + 1
		IM1 = I - 1
		IPN = I + NN
		IMN = I - NN
		
		RP = B(I)
		
		DA = DAA*RH(N)
		IF (M.EQ.1) GO TO 10
		IF (M.EQ.MM) GO TO 15

		RPT = (B(IPN) - B(IMN))/DT2
		
		IF (N.EQ.1) GO TO 20
		IF (N.EQ.NN) GO TO 25
		RPR = (B(IP1) - B(IM1))/DR2
		GO TO 45

C			BOTTOM
10		CALL BCAL(BB1,BB2,N,FNW,FW,rz)
		IF (N.EQ.1) GO TO 12
		IF (N.EQ.NN) GO TO 13
		RPR = (B(IP1) - B(IM1))/DR2
		GO TO 14
C			BOTTOM,INNER
12		RPR = -fnw*rp/r
		GO TO 14
C			BOTTOM,OUTER
13		X = CC1*B(I) + CC2*B(IM1)
		RPR = (X-B(IM1))/DR2
14		RPT = BB1*RPR + BB2*RP
		DA = DA/2.0
		GO TO 45

C			TOP
15		RPT = 0.
		IF (N.EQ.1) GO TO 17
		IF (N.EQ.NN) GO TO 18
		RPR = (B(IP1) - B(IM1))/DR2
		GO TO 19
C			TOP,INNER
17		RPR = -fnw*rp/r
		GO TO 19
C			TOP,OUTER
18		X = CC1*B(I) + CC2*B(IM1)
		RPR = (X - B(IM1))/DR2
19		DA = DA/2.0
		GO TO 45

C			INNER
20		RPR = -a1*rpt - fnw*rp/r
		GO TO 45

C			OUTER
25		X = CC1*B(I) + CC2*B(IM1)
		RPR = (X -B(IM1))/DR2
45		RPR = RPR + A1*RPT
		XX(M) = (F*RPR + WNR*RP)/FW
		RI = RI + XX(M)*XX(M)*0.5*DA*R
50		CONTINUE

		IF (IGG.EQ.0) GO TO 55
		IF (IPX.EQ.0) GO TO 55
		DO 52 M = 1,MM
			ZZ = (-1.0 + DT*FLOAT(M-1))*RH(N)
			WRITE(12,*) R, ZZ, XX(M)
52		Continue

55		IF (IWS.EQ.0) GO TO 60
		DZ = RH(N)*DT
		WRITE(IOPF,902) R,DZ,RH(N)
		WRITE(IOPF,903) (XX(J),J=1,MM)
60		CONTINUE

		CLOSE (12)

901		FORMAT(////' V STRUCTURE')
902		FORMAT(/' R,DZ,H =',3E12.5/)
903		FORMAT(2X,10E12.5)
910		FORMAT(' DO YOU WANT V GRAPHICS? 0 OR 1')
		RETURN
		END

		SUBROUTINE UCAL(W,NR,XX,RI,IWS,IGG,rz)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM =NN*MM,NMX =2*NN +3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
		DOUBLE PRECISION XX(NM)

C			SUBROUTINE TO COMPUTE RADIAL VELOCITY U
		
		DR2 = 2.0*DR
		DT2 = 2.0*DT
		WRITE(IOPF,901)
		DAA = DT*DR
		RI = 0.
		FW = F*F - W*W
		FNW = FLOAT(NR)*F/W
		FN = F*FLOAT(NR)
		WN = W*FLOAT(NR)
		BB3 = (1.0 + FNW)/(DR*FLOAT(NN))
		X = 1.0/DR + 0.5*BB3
		CC1 = 2.0/(X*DR)
		CC2 = -(1.0/DR -BB3*0.5)/X

		IF (IGG.EQ.0) GO TO 5
		WRITE(6,910)
		READ(5,*) IPX
		IF (IPX.EQ.0) GO TO 5
		OPEN (12, FILE = 'UGR.grd')

5		DO 60 N = 1,NN
		R = rz +DR*FLOAT(N-1)
		WNR = WN/R
		FNR = FN/R
		
		DO 50 M = 1,MM
		CALL AS(A1,A2,A3,N,M)
		I = N + NN*(M-1)
		IP1 = I + 1
		IM1 = I - 1
		IPN = I + NN
		IMN = I - NN
		
		RP = B(I)
		
		DA = DAA*RH(N)
		IF (M.EQ.1) GO TO 10
		IF (M.EQ.MM) GO TO 15

		RPT = (B(IPN) - B(IMN))/DT2
		
		IF (N.EQ.1) GO TO 20
		IF (N.EQ.NN) GO TO 25
		RPR = (B(IP1) - B(IM1))/DR2
		GO TO 45

C			BOTTOM
10		CALL BCAL(BB1,BB2,N,FNW,FW,rz)
		IF (N.EQ.1) GO TO 12
		IF (N.EQ.NN) GO TO 13
		RPR = (B(IP1) - B(IM1))/DR2
		GO TO 14
C			BOTTOM,INNER
12		RPR = -fnw*rp/r
		GO TO 14
C			BOTTOM,OUTER
13		X = CC1*B(I) + CC2*B(IM1)
		RPR = (X-B(IM1))/DR2
14		RPT = BB1*RPR + BB2*RP
		DA = DA/2.0
		GO TO 45

C			TOP
15		RPT = 0.
		IF (N.EQ.1) GO TO 17
		IF (N.EQ.NN) GO TO 18
		RPR = (B(IP1) - B(IM1))/DR2
		GO TO 19
C			TOP,INNER
17		RPR = -fnw*rp/r
		GO TO 19
C			TOP,OUTER
18		X = CC1*B(I) + CC2*B(IM1)
		RPR = (X - B(IM1))/DR2
19		DA = DA/2.0
		GO TO 45

C			INNER
20		RPR = -a1*rpt - fnw*rp/r
		GO TO 45

C			OUTER
25		X = CC1*B(I) + CC2*B(IM1)
		RPR = (X -B(IM1))/DR2
45		RPR = RPR + A1*RPT
		XX(M) = (-W*RPR - FNR*RP)/FW
		RI = RI + XX(M)*XX(M)*0.5*DA*R
50		CONTINUE

		IF (IGG.EQ.0) GO TO 55
		IF (IPX.EQ.0) GO TO 55
		DO 52 I = 1,MM
			ZZZ = (-1.0 + DT*FLOAT(I-1))*RH(N)
			WRITE(12,*) R, ZZZ, XX(I)
52		Continue

55		IF (IWS.EQ.0) GO TO 60
  		DZ = RH(N)*DT
		WRITE(IOPF,902) R,DZ,RH(N)
		WRITE(IOPF,903) (XX(J),J=1,MM)
60		CONTINUE

		CLOSE (12)

901		FORMAT(////' U STRUCTURE')
902		FORMAT(/' R,DZ,H =',3E12.5/)
903		FORMAT(2X,10E12.5)
910		FORMAT(' DO YOU WANT U GRAPHICS? 0 OR 1')
		RETURN
		END
		
		SUBROUTINE RHOC(XX,W,NR,RI,IWS,IGG,rz)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX =2*NN + 3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)
		DOUBLE PRECISION XX(NM)

C			SUBROUTINE TO COMPUTE DENSITY PERTURBATIONS

		WRITE(IOPF,901)

		DR2 = 2.0*DR
		DAA = DT*DR
		RI = 0.
		FW = F*F - W*W
		FNW = FLOAT(NR)*F/W
		G2 = (980.0)**2
		G = 1.0/980.

		IF (IGG.EQ.0) GO TO 5
		WRITE(6,910) 
		READ(5,*) IPX
		IF (IPX.EQ.0) GO TO 5
		OPEN (12, FILE = 'RGR.grd')

5		DO 60 N = 1,NN
		R = rz + DR*FLOAT(N-1)
		DZ =RH(N)*DT
		DZ2 = 2.0*DZ

		DO 50 M = 1,MM
		I = N + NN*(M-1)
		IP1 = I + 1
		IM1 = I- 1
		IPN = I + NN
		IMN = I - NN

		RP = B(I)

		IF (M.EQ.1) GO TO 10
		IF (M.EQ.MM) GO TO 20
		RPZ = (B(IPN) - B(IMN))/DZ2
		DA = DAA*RH(N)*R
		GO TO 40

C			BOTTOM
10		CALL BCAL(BB1,BB2,N,FNW,FW,rz)
		DA = DAA*RH(N)*R/2.0
		IF (N.EQ.1) GO TO 12
		IF (N.EQ.NN) GO TO 14
		RPR = (B(IP1) - B(IM1))/DR2
		GO TO 16
C			INNER BOTTOM
12		RPR = -fnw*rp/r
		GO TO 16
C			OUTER BOTTOM
14		RPZ = 0.
		GO TO 40
16		RPT = BB1*RPR + BB2*RP
		RPZ = RPT/RH(N)
		GO TO 40

C			TOP
20		RPZ = 0.
		DA = 0.
40		XX(M) = -RPZ*G
		RI = RI + XX(M)*XX(M)*0.5*DA*G2/BV(N,M)
50		CONTINUE

		IF (IGG.EQ.0) GO TO 55
		IF (IPX.EQ.0) GO TO 55
		DO 52 I = 1,MM
			ZZZ = (-1.0 + DT*FLOAT(I-1))*RH(N)
			WRITE(12,*) R, ZZZ, XX(I)
52		Continue

55		IF (IWS.EQ.0) GO TO 60
		WRITE(IOPF,902) R,DZ,RH(N)
		WRITE(IOPF,903) (XX(J),J=1,MM)
60		CONTINUE

		CLOSE (12)

901		FORMAT(////' RHO STRUCTURE')
902		FORMAT(/' R,DZ,H =',3E15.5/)
903		FORMAT(2X,10E12.5)
910		FORMAT(' DO YOU WANT RHO GRAPHICS? 0 OR 1')
		RETURN
		END

		SUBROUTINE DEP(rz,xmax)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM=NN*MM,NMX = 2*NN + 3)

		DOUBLE PRECISION AX,BX
		COMMON AX(NM,NMX),BX(NM)
		COMMON F,DX,DT
		COMMON RH(NN),RHR(NN),RHRR(NN),BV(NN,MM),BVZ(NN,MM)

c			island-trapped waves for arbitrary depth profile

C			SUBROUTINE TO COMPUTE DEPTH PROFILE, GIVEN 
C				INPUTS

		READ(INPF,*) NRD

C			NRD = NUMBER OF R,H PAIRS TO READ
C				R IN KM, H IN METERS
C			The first radius read is the coast
c			Depth at coast must be nonzero

		DO 10 I = 1,NRD
			READ(INPF,*) RHR(I),RHRR(I)
C				RHR = RADIUS
C				RHRR = DEPTH
C				SPACING IS ARBITRARY
			RHR(I) = RHR(I)*1.0E+05
			RHRR(I) = RHRR(I)*100.
10		Continue

		rz = rhr(1)
		rh(1) = rhrr(1)
		dx = (xmax-rz)/float(nn-1)
		
		DO 40 I = 2,NN
			X = rz +DX*FLOAT(I-1)

15			IC = 2
20			IF (X.GT.RHR(IC)) GO TO 25
			IL = IC - 1
			XX = (X -RHR(IL))/(RHR(IC) - RHR(IL))
			RH(I) = RHRR(IL) + XX*(RHRR(IC) -RHRR(IL))
			GO TO 40

25			IC = IC + 1
			IF (IC.GT.NRD) GO TO 30
			GO TO 20

30			RH(I) = RHRR(NRD)

40		CONTINUE

		RHR(NN) = 0.
		RHRR(NN) = 0.
		DXX = DX*DX
		DX2 = DX*2.0
		RHR(1) = (RH(2)-RH(1))/DX
		RHRR(1) = 0.

		NX = NN - 1
		DO 60 I = 2,NX
			IL = I - 1
			IH = I + 1
			RHR(I) = (RH(IH) -RH(IL))/DX2
			RHRR(I) = (RH(IH) - 2.0*RH(I) + RH(IL))/DXX
60		CONTINUE

C		Now,	RH = DEPTH
C			RHR = D(DEPTH)/DR
C			RHRR = D2(DEPTH)/DR2

100		RETURN
		END


		SUBROUTINE BANDG
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX = 2*NN + 3)

		DOUBLE PRECISION A,BB
		COMMON A(NM,NMX),BB(NM)
		COMMON F,DR,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.LE.NM) GO TO 5
		MH = NM
5		DO 50 IR = IP,MH
			ICD = MBW + N + 1 -IR
			R = A(IR,ICD)/A(N,MDIAG)
			BB(IR) = BB(IR) -R*BB(N)
			DO 50 IC = IP,MH
				ICD = MBW + IC + 1 -IR
				ICB = MBW + IC + 1 -N
				A(IR,ICD) = A(IR,ICD) -R*A(N,ICB)
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.GE.1) GO TO 60
		IL = 1
60		DO 100 IR = IL,IH
			ICD = MBW + N + 1 -IR
			BB(IR) = BB(IR) -A(IR,ICD)*BB(N)
100	CONTINUE

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

	RETURN
	END

		SUBROUTINE DCALC(NR,W,rz)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX = 2*NN + 3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)

C			CALCULATE DAMPING TIME

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

	DAA = DT*DR
	DR2 = 2.0*DR
	DT2 = 2.0*DT
	FW = F*F - W*W
	FPW = F*F + W*W
	FNW = (F*FLOAT(NR))/W
	BB3 = (1.0 + FNW)/(DR*FLOAT(NN))
	X =  1.0/DR + BB3*0.5
	CC1 = 2.0/(X*DR)
	CC2 = -(1.0/DR -BB3*0.5)/X
	fn = float(nr)
	FN2 = FLOAT(NR*NR)
	W2 = W*W
	FX = FPW/FW
	FX2 = 2.0*FLOAT(NR)*W*F/FW
	MX = MM-1

	DO 100 N = 1,NN
		R = rz + DR*FLOAT(N-1)
		IB = N
		dxaa = 1.0
		
		IPN = N + NN
		IP1 = N + 1
		IM1 = N - 1

		CALL BCAL(BB1,BB2,N,FNW,FW,rz)
		CALL AS(A1,A2,A3,N,1)

		IF (N.EQ.1) GO TO 10
		IF (N.EQ.NN) GO TO 20
		
		RPBR = (B(IP1) -B(IM1))/DR2
		GO TO 30
C			INNER BOTTOM
10		RPBR = -fnw*b(ib)/r
		dxaa = 0.5
		GO TO 30
C			OUTER BOTTOM
20		X = CC1*B(IB) + CC2*B(IM1)
		RPBR = (X-B(IM1))/DR2
30		RPT = BB1*RPBR + BB2*B(IB)
		RPZ = RPT/RH(N)
		RPRB = RPBR + A1*RPT

		RI1 = RI1 + RHX(N)*B(IB)*B(IB)*dxaa
		RI3 = RI3 + (B(IB)*B(IB)/R)*dxaa
		RI4 = RI4 + R*RPRB*RPBR*dxaa
		RI5 = RI5 + B(IB)*(RPRB +RPBR)*DXAA
		
		DA = DAA*RH(N)*R/2.0

		RI2 =  RI2 + DXAA*RPZ*RPZ*DA/BV(N,1)

		DA = 2.0*DA
		DZ2 = DT2*RH(N)

		DO 100 M = 2,MX
			I = N + NN*(M-1)
			IPN = I + NN
			IMN = I - NN

			RPZ = (B(IPN) -B(IMN))/DZ2
			RI2 = RI2 + dxaa*RPZ*RPZ*DA/BV(N,M)
100		CONTINUE

		dz = dt*rh(1)
		do 140 m = 1,mm
			i = 1 + nn*(m-1)
			da = dz
			if (m.eq.1) go to 110
			if (m.ne.mm) go to 120
110			da = dz/2.
120			ri6 = ri6 + DA*B(I)*B(I)	
140		continue

		anm = (ri4 + fn2*ri3)*dr
		bnm = dr*ri5
		anm = w*anm/(fn*f)
		bnm = -2*w*bnm/(f*f)


		RI2 = RI2*2.0*W2
		RI1 = RI1 *DR*(-FNW)
		RI3 = RI3*DR*FN2*FX
		RI4 = RI4*DR*FX
		RI5 = RI5*DR*FX2
		ri6 = ri6*(-fnw)

		ETF = -(RI1 + RI2 +ri6)/(RI3 +RI4 +RI5)
		WRITE(IOPF,902) ETF
		write(iopf,903)
		write(Iopf,904)
		write(iopf,905) anm
		write(iopf,906) bnm


902		FORMAT(/' EPSILON*TF (CM) =',E15.5)
903		format('(Valid for any subinertial frequency)')
904		format(//' Friction coefficents for f^2 >> w^2 :')
905		format(' ann/eps (cm^-1) =',e15.5)
906		format(' bnn/eps (sec/cm) =',e15.5//)
		RETURN
		END


		Subroutine norm(rz)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX = 2*NN + 3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)

c			Normalize modes by low frequency (f^2 >> w^2) norm

		mt = 1 + nn*(mm-1)
		if (b(mt).ge.0) go to 8
		
		do 5 i = 1,nm
			b(i) = -b(i)
5		Continue

8		r1 = rhx(1)*0.5*b(1)*b(1)

		do 10 n = 2,nn
			r1 = r1 + rhx(n)*b(n)*b(n)
10		Continue
		r1 = r1*dr

		r2 = 0.5*(b(1)*b(1) + b(mt)*b(mt))
		mmax = mm-1

		do 20 m = 2,mmax
			i = 1 + nn*(m-1)
			r2 = r2 + b(i)*b(i)
20		Continue
		dz = dt*rh(1)
	
		r = r1 + dz*r2
		r = sqrt(r)

		do 30 i = 1,nm
			b(i) = b(i)/r
30		Continue

		return
		end		


		Subroutine coup(w,rz,wms,nr)
		PARAMETER (NN=67,MM=30,INPF=7,IOPF=8,IGPF=12)
		PARAMETER (NM = NN*MM, NMX = 2*NN + 3)

		DOUBLE PRECISION AX,B
		COMMON AX(NM,NMX),B(NM)
		COMMON F,DR,DT
		COMMON RH(NN),RHX(NN),RHXX(NN),BV(NN,MM),BVZ(NN,MM)

c			Compute wind coupling coefficents
c			Valid for f^2 >> w^2
c		
c			Taur = cos(theta)*(1 - exp(-zeta))
c			Tautheta = -sin(theta)*(1 + exp(-zeta))
c	
c			where     zeta = (r - rz)/wms
c
c			If wms < 0.5 km, winds are strictly uniform
c
c			gn,-1, fn,-1 reflect the fact that w < 0 for m < 0
c
c			This only makes sense for nr = +/- 1

		if (nr.ne.1) go to 120
		write(iopf,901)

		i = 1 + nn*(mm-1)
		pnz = b(i)

c			Uniform winds

		fn1 = -0.5*rz*pnz*w
		fnm1 = -fn1
		gn1 = -w*rz*pnz*0.5/f
		gnm1 = gn1
		write(iopf,906)
		write(iopf,902) fn1
		write(iopf,903) fnm1
		write(iopf,904) gn1
		write(iopf,905) gnm1

		If (wms.lt.0.5) go to 100

50		wms = wms*1.0e+05

c			Nonuniform winds

		ri = 0.5*pnz*(2.0 - rz/wms)

		do 60 n = 2,nn
			i = n + nn*(mm-1)
			z = dr*float(n)/wms
			r = rz + dr*float(n-1)
			ee = exp(-z)
			ri = ri + b(i)*ee*(2.0 - r/wms)
60		Continue

		ri = ri *dr

		fn1 = 2*rz*pnz + ri
		fn1 = -0.5*w*fn1
		fnm1 = -fn1
		gn1 = w*0.5*ri/f
		gnm1 = gn1
		wms = wms*1.0e-05
		write(iopf,907) wms

100		write(iopf,902) fn1
		write(iopf,903) fnm1
		write(iopf,904) gn1
		write(iopf,905) gnm1

901	format(//' Wind coupling coefficients with T and Tt removed')
902	format(/' fn,1  (cm^0.5 sec^-1) =',e15.5)
903	format(' fn,-1 (cm^0.5 sec^-1) =',e15.5)
904	format(' gn,1 (cm^0.5) = ',e15.5)
905	format(' gn,-1 (cm^0.5) =',e15.5//)
906	format(' Computed for uniform winds')
907	format(//' for non-uniform winds with decay scale (km)',e15.5)
120		return
		end
