c	Program for computing barotropic seamount-trapped waves with a mean
c		flow   2/1/93
c	Subroutine VRD corrected (l.c.) 8/21/93
c	Subroutine STABCK corrected 1/6/95
c	Subroutines VCAL, UCAL corrected 4/5/2002
c	Modified for graphics output 1/6/95
c		If you opt for graphics, you will prompted for whether you
c		want p, u, or v from the first mode computed.
c		This will give a flat file of (x,mode) that can easily be
C		used with MATLAB to plot it up.


		PROGRAM BTSTW
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		DIMENSION U(400), RX(3), WW(3), BH(400),WH(10),NH(10)
		CHARACTER*20 INP,IOP, IGRP

C
C			RIGID LID ONLY
C			VBAR = MEAN AZIMUTHAL FLOW
C			ASSUMES THAT
C				1) VBAR = 0 AT R = 0
C				2) VBAR = 0, HR= 0 AT OUTER LIMIT
C			SOLVES FOR PRESSURE USING RESONANCE ITERATION
C			THEN COMPUTES VELOCITY, DECAY TIME
C			ALL UNITS WITHIN PROGRAM ARE CGS
C			
C INPUT:
C		INPUT, OUTPUT, GRAPHICS FILE NAMES FROM TERMINAL
C	
C		FROM INPUT FILE:
C
C		NN, NITM, EPS, DEL
C		NFIRST, NMAX
C		IPC (PRINTOUT CONTROL)
C		F,RMAX	(F, W MULTIPLIED BY 1E5, RMAX IN KILOMETERS)
C		NDEPTHS
C			H AT R = 0
C			R, H PAIRS (NDEPTHS TIMES, AT LEAST ONCE)
C		NVELOCITIES (CAN BE ZERO)
C			R, V PAIRS (NVELOCITIES TIMES)
C		WW,WW,WW
C

C			FIND OUT THE INPUT, OUTPUT FILE NAMES
		WRITE(6,901)
		READ(5,*) INP
		WRITE(6,906)
		READ(5,*) IOP
		OPEN (7,FILE=INP)
		OPEN (8,FILE=IOP)
		WRITE(6,910)
		READ(5,*) IGRPX
		IF (IGRPX.EQ.0) GO TO 2
		WRITE(6,911)
		READ(5,*) IGRP
		OPEN (12,FILE=IGRP)


2		READ(7,*) NN,NITM,EPS,DEL

C			NN = NUMBER OF GRID POINTS (<400)
C			NITM = MAX NUMBER OF ITERATIONS
C			EPS = NOMINAL ACCURACY
C			DEL = SEARCH INCREMENT

C			USE ONLY UNBIASED SEARCH
		ISD = 0
		READ(7,*) NF,NMAX

C			NF = FIRST N (AZIMUTHAL NUMBER)
C			NMAX = LAST N

		READ(7,*) IPC

C			IPC = 0: NO PRINTOUT OF MODAL STRUCTURES
C			IPC NE 0: PRINTOUT

		READ(7,*) F,XMAX

C			F, OUTER RADIUS OF GRID (KM)

		F = F*1.0E-05
		XMAX = XMAX*1.0E+05
		DX = XMAX/FLOAT(NN)
		WRITE(8,902) F,XMAX,NN
		WRITE(8,905) DX

C			DEP READS THE DEPTH PROFILE
C			VRD READS THE MEAN VELOCITY PROFILE
C			STABCK CHECKS FOR BAROTROPIC STABILITY

		CALL DEP(NN,DX)
		CALL VRD(NN,DX)
		CALL STABCK(NN,DX,F)

		READ(7,*) WW(1),WW(2),WW(3)

C			FIRST GUESS FREQUENCIES

	DO 5 I = 1,3
		WW(I) = WW(I)*1.0E-05
5	Continue

		N = NF
		IX = 1

C			START BY LOOKING AT THE 3 FIRST GUESSES
10		RB = 0.
		WRITE(8,908)
	DO 20 J = 1,3
		W = WW(J)
		CALL MATS(W,N,F,DX,NN,RI)
		RX(J) = RI
 		WRITE(8,903) W,N,RI
		IF (RI.LT.RB) GO TO 20
C			SAVE THE BEST INFORMATION
		RB = RI
		WB = W
		DO 15 I = 1,NN
15		BH(I) = B(I)
20	CONTINUE

		NIT = 3
		IGP = 0
C			NOW START TO SEARCH, GIVEN INFORMATION FROM THE
C				FIRST 3 GUESSES

30		CALL NGW(WW,RX,WB,WN,ISUC,EPS,DEL,IN,ISD,IGP)
		IF (ISUC.EQ.1) GO TO 100

		CALL MATS(WN,N,F,DX,NN,RI)
		IF (IN.NE.0) GO TO 50
		IF (RI.GT.RB) GO TO 40
 		IF (WN.LT.WW(2)) GO TO 25
22		IN = 3
		GO TO 50
25		IN = 1
		GO TO 50
40		IF (WN.GT.WW(2)) GO TO 25
		GO TO 22

50		RX(IN) = RI
		WW(IN) = WN
		WRITE(8,903) WN,N,RI
		IF (RI.LT.RB) GO TO 70
C			SAVE THE BEST INFORMATION TO DATE
		WB = WN
		RB = RI
	DO 60 I = 1,NN
		BH(I) = B(I)
60	Continue


70		NIT = NIT + 1
C			CHECK TO SEE IF YOU RAN OUT OF ITERATIONS
		IF (NIT.LT.NITM) GO TO 30
		WRITE(8,904) NIT
		GO TO 200

C			SOLUTION HAS BEEN FOUND
C			NORMALIZE THE BEST SOLUTION

100		RBB = SQRT(RB)
	DO 110 I = 1,NN
		B(I) = BH(I)/RBB
110	Continue

		WRITE(8,907) WB,N,NIT
C			CLCK CHECKS FOR CRITICAL LAYERS
		CALL CLCK(W,N,NN,DX)
C			SAVE THE RESULTS

		IF (IGRPX.EQ.0) GO TO 120
		WRITE(6,912)
		READ(5,*) IXXX
		IF (IXXX.EQ.0) GO TO 120
		IGRPX = 0
	DO 115 I = 1,NN
		XXX = DX*FLOAT(I)
		WRITE(12,*) XXX, B(I)
115	Continue

120		WH(IX) = WB
		NH(IX) = N

C			UCAL COMPUTES RADIAL VELOCITY MODAL STRUCTURE
C			DTIME COMPUTES FRICTIONAL DECAY TIME
C			VCAL COMPUTES AZIMUTHAL VELOCITY MODAL STRUCTURE

		CALL UCAL(WB,N,F,DX,U,NN,IPC)
		IF (IGRPX.EQ.0) GO TO 135
		WRITE(6,913)
		READ(5,*) IXXX
		IF (IXXX.EQ.0) GO TO 135
		IGRPX = 0
	DO 130 I = 1,NN
		XXX = DX*FLOAT(I)
		WRITE(12,*) XXX,U(I)
130	Continue

135		CALL DTIME(WB,N,F,NN,DX,U)
		IF (IPC.EQ.0) GO TO 140
		CALL VCAL(WB,N,F,DX,U,NN)
		IF (IGRPX.EQ.0) GO TO 140
		WRITE(6,914) 
		READ(5,*) IXXX
		IF (IXXX.EQ.0) GO TO 140
		IGRPX = 0
	DO 138 I = 1,NN
		XXX = DX*FLOAT(I)
		WRITE(12,*) XXX,U(I)
138	Continue

C			GO TO THE NEXT POINT ON THE 'DISPERSION CURVE'
140 		N = N + 1
		IF (N.GT.NMAX) GO TO 200
C			GENERATE THE FIRST GUESS FOR THE NEXT POINT
		IXX = IX
		IX = IX + 1
		IF (IX.GT.2) GO TO 150
		WW(2) = W*FLOAT(N)/FLOAT(NH(IXX))
		GO TO 160
150		IXXX = IXX - 1
		WW(2) = WH(IXX) + (WH(IXX) -WH(IXXX))
160		WW(1) = WW(2)*(1.0 -DEL)
		WW(3) = WW(2)*(1.0 + DEL)
		GO TO 10

200		CLOSE(7)
		CLOSE(8)
		CLOSE(12)

901		FORMAT(' INPUT FILE?')
902		FORMAT(' F, RMAX, NN = ',2E15.5, I5)
903		FORMAT(' W, N, RI =',E12.5,I5,3X,E12.5)
904		FORMAT(' USED UP',I3,' ITERATIONS')
905		FORMAT(' DX = ',E12.5)
906		FORMAT(' OUTPUT FILE?')
907		FORMAT(/' CONVERGED: W, N, NITERATIONS =',E12.5,2I5)
908		FORMAT(///)
910 		FORMAT(' DO YOU WANT OUTPUT FOR GRAPHICS? 0 OR 1')
911		FORMAT(' GRAPHICS FILE NAME?')
912		FORMAT(' DO YOU WANT PRESSURE GRAPHICS OUTPUT? 0 OR 1')
913		FORMAT(' DO YOU WANT u GRAPHICS OUTPUT? 0 OR 1')
914		FORMAT(' DO YOU WANT v GRAPHICS OUTPUT? 0 OR 1')

250		STOP
      		END


		SUBROUTINE MATS(W,N,F,DX,NN,RI)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		DOUBLE PRECISION A1,A2,A3

C 			SET UP PRESSURE MATRIX EQUATION AND SOLVE IT

	DX2 = DX/2.
	DXX = DX*DX

	DO 50 I = 1,NN
		II = I + NN
		III = I + 2*NN

		CALL COCAL(W,N,F,DX,I,A1,A2,A3)
		A(I) = -2.0*A1 + DXX*A3
		A(II) = A1 - A2*DX2
		A(III) = A1 + A2*DX2

		IF (I.EQ.1) GO TO 5
		IF (I.EQ.NN) GO TO 10
		GO TO 50

C			R = R BOUNDARY: P = 0
5		A(II) = 0.
		GO TO 50

C			R = RMAX BOUNDARY: PR + N*P/R = 0
10		X = DX*FLOAT(I)
		XX = 2.0*DX*FLOAT(N)/X
		A(I) = A(I) -A(III)*XX
		A(II) = A(II) + A(III)
		A(III) = 0.

50	CONTINUE

	DO 60 I = 1,NN
		B(I) = 0.
60	Continue
		B(10) = 1.

		CALL TRI(NN)
		CALL CALI(NN,DX,RI)

	RETURN
	END


		SUBROUTINE CALI(NN,DX,RI)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		DOUBLE PRECISION RII

C			 COMPUTE MEASURE OF RESONANCE = RI

	RII = 0.

	DO 50 I = 1,NN
		R = DX*FLOAT(I)
		RII = RII + R*B(I)*B(I)*DX*RH(I)
50	Continue
	RI = RII
		
		RETURN 
		END


		SUBROUTINE COCAL(W,N,F,DX,I,A1,A2,A3)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		DOUBLE PRECISION A1,A2,A3

C			COMPUTE COEFFICIENTS OF MATRIX EQUATION

		X = DX*FLOAT(I)
		XX = X*X
		FN = FLOAT(N)
		WP = W + FN*V(I)/X
		FP = F + V(I)/X
		FPR = (VX(I) - V(I)/X)/X
		WPR = FN*FPR

		FF = (FP + V(I)/X)*(FP + VX(I)) -WP*WP
		FR = 2.0*FPR*(FP + VX(I))
		FR = FR + (FP+V(I)/X)*(VX(I)/X -V(I)/XX + VXX(I))
		FR = FR - 2.0*WP*WPR

		A1 = WP*RH(I)
		A2 = WP*RH(I)/X + WP*RHX(I) -WP*RH(I)*FR/FF
		A2 = A2 + RH(I)*WPR + RH(I)*FN*(V(I)/X - VX(I))/X
		A3 = FN*RHX(I)*(FP + V(I)/X)/X
		A3 = A3 -FN*RH(I)*(FP +V(I)/X)*FR/(X*FF)
		A3 = A3 + FN*RH(I)*2.0*FPR/X
		A3 = A3 - FN*FN*RH(I)*WP/XX

		RETURN
		END


		SUBROUTINE UCAL(W,N,F,DX,U,NN,IPC)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		DIMENSION U(400)

C			COMPUTE AND PRINT U (RADIAL VELOCITY)

		FN = FLOAT(N)
		DX2 = DX*2.0

		DO 50 I = 1,NN
		X = DX*FLOAT(I)
		WP = W + FN*V(I)/X
		FP = F + V(I)/X
		FF = (FP + V(I)/X)*(FP + VX(I)) -WP*WP

		IP = I +1
		IM = I - 1
		IF (I.EQ.1) GO TO 5
		IF (I.EQ.NN) GO TO 10
		RPR = (B(IP) - B(IM))/DX2
		GO TO 20

C			R = 0 BOUNDARY
5		RPR = B(IP)/DX2
		GO TO 20

C			R = RMAX BOUNDARY
10		RPR = -FN*B(I)/X

20		U(I) = (-WP*RPR -FN*(FP +V(I)/X)*B(I)/X)/ff

50		CONTINUE

		IF (IPC.EQ.0) GO TO 100
		WRITE(8,902)
		WRITE(8,901) (U(I),I=1,NN)

C			PRINTOUT STARTS AT R = DX, INCREMENTING BY DX

901		FORMAT(1X,5E12.5)
902		FORMAT(//' U'/)

100		RETURN
		END


		SUBROUTINE VCAL(W,N,F,DX,U,NN)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		DIMENSION U(400)

C 			COMPUTE AND PRINT V (AZIMUTHAL VELOCITY)

	FN = FLOAT(N)
	DX2 = DX*2.0

	DO 50 I = 1,NN
		X = DX*FLOAT(I)
		WP = W + FN*V(I)/X
		FP = F + V(I)/X
		FF = (FP + V(I)/X)*(FP + VX(I)) -WP*WP

		IP = I +1
		IM = I - 1
		IF (I.EQ.1) GO TO 5
		IF (I.EQ.NN) GO TO 10
		RPR = (B(IP) - B(IM))/DX2
		GO TO 20

C			R = 0 BOUNDARY
5		RPR = B(IP)/DX2
		GO TO 20

C			R = RMAX BOUNDARY
10		RPR = -FN*B(I)/X

20		U(I) = ((FP+VX(I))*RPR + WP*FN*B(I)/X)/ff

50	CONTINUE
		
	WRITE(8,902)
	WRITE(8,901) (U(I),I=1,NN)

901	FORMAT(1X,5E12.5)
902	FORMAT(//' V'/)
	RETURN
	END


		SUBROUTINE VRD(NN,DX)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)

C			READ VBAR AND COMPUTE VX, VXX

	READ(7,*) NRD

C			NUMBER OF R,VBAR PAIRS

	IF (NRD.NE.0) GO TO 20

C			VBAR = 0 EVERYWHERE

	WRITE(8,901)
	DO 10 I = 1,NN
		V(I) = 0.
		VX(I) = 0.
		VXX(I) = 0.
10	Continue
		GO TO 100

20	DO 25 I = 1,NRD

C			READ R IN KM, VBAR IN CM/SEC
C			R MUST BE IN INCREASING ORDER
C			LAST VBAR MUST BE 0

		READ(7,*) A(I),B(I)
		A(I) = A(I)*1.0E+05
25	Continue

C 			NOW INTERPOLATE

	DO 60 N = 1,NN
		X = DX*FLOAT(N)
		IF (X.GE.A(NRD)) GO TO 50
		IF (X.LT.A(1)) GO TO 40

		J = 0
		DO 30 I = 1,NRD
			IF (J.NE.0) GO TO 30
			IF (a(i).lt.x) GO TO 30
			J = I
30		CONTINUE

		jj = j
		j = j-1
		XX = X -A(J)
		XY = A(JJ) - A(J)
		V(N) = B(J) + XX*(B(JJ)-B(J))/XY
		GO TO 60

40		V(N) = X*B(1)/A(1)
		GO TO 60

50		V(N) = 0.
60	CONTINUE

C			NOW COMPUTE DERIVATIVES

		DX2 = DX*2.0
		DXX = DX*DX
		VX(1) = V(2)/DX2
		VXX(1) = 0.
		VXX(NN) = 0.
		NQ = NN - 1
		VX(NN) = (V(NN) - V(NQ))/DX

	DO 70 I = 2,NQ
		IP = I + 1
		IM = I - 1
		VX(I) = (V(IP) - V(IM))/DX2
		VXX(I) = (V(IP) - 2.0*V(I) + V(IM))/DXX
70	Continue
	WRITE(8,902)
	WRITE(8,903) (V(I),I=1,NN)

901	FORMAT(/' VBAR = 0 EVERYWHERE')
902	FORMAT(/' V (CM/SEC)')
903	FORMAT(1X,5E12.5)

100	RETURN
	END


                SUBROUTINE DEP(NN,DX)
      		COMMON RH(400),RHR(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		DIMENSION RHRR(400)

C                       SUBROUTINE TO COMPUTE DEPTH PROFILE, GIVEN 
C                               INPUTS

        READ(7,*) NRD
        READ(7,*) RHZ
	
C                       NRD = NUMBER OF R,H PAIRS TO READ
C                               R IN KM, H IN METERS
C			FOR R > (LAST READ), USES LAST DEPTH READ
C                       RHZ = DEPTH AT R=0

        RHZ = RHZ*100.
                
        DO 10 I = 1,NRD
                READ(7,*) RHR(I),RHRR(I)

C                       RHR = RADIUS IN KM
C                       RHRR = DEPTH IN M
C                       SPACING IS ARBITRARY, BUT R MUST BE INCREASING

                RHR(I) = RHR(I)*1.0E+05
                RHRR(I) = RHRR(I)*100.
10	Continue

C			NOW INTERPOLATE
                
        DO 40 I = 1,NN
                X = DX*FLOAT(I)
                IF (X.GT.RHR(1)) GO TO 15
                XX = X/RHR(1)
                RH(I) = RHZ + XX*(RHRR(1) - RHZ)
                GO TO 40

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

C			NOW COMPUTE DERIVATIVE: ASSUMES A FLAT SUMMIT:
C				Hr = 0 AT R = 0

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

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


                WRITE(8,901)
		WRITE(8,903) RHZ
                WRITE(8,902) (RH(I),I=1,NN)

C                       NOW RH = DEPTH IN CM
C                       RHR = D(DEPTH)/DR

901		FORMAT(/' DEPTH')
902		FORMAT(1X,5E12.5)
903		FORMAT(' CENTRAL DEPTH =',E15.5)

100             RETURN
                END



        SUBROUTINE TRI(N)
        COMMON RH(400),RHX(400),A(1200),B(400),V(400),VX(400),VXX(400)
        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
       
	DO 5 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
5	Continue
C			NOW BACKSOLVE
        DO 10 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
10	Continue
        B(1) = B(1)/A(1)

        RETURN
        END



        SUBROUTINE NGW(WW,RX,WB,WN,ISUC,EPS,DEL,IN,ISD,IGP)
		COMMON RH(400),RHX(400), AQQ(1200), BQQ(400), V(400), VX(400)
		COMMON VXX(400)
        DIMENSION WW(3),RX(3)

C               SUBROUTINE TO GUESS THE NEXT W

C		ORDER W BY FREQUENCY

5       IC = 0
        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

C		YOU HAVE THE RESONANCE SURROUNDED
        IL = 1
        IH = 3
        IF (RX(3).GT.RX(1)) GO TO 9
        IL = 3
        IH = 1
C		FIRST TRY NEWTON'S RULE
9       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.LT.WW(1)) GO TO 15
        IF (WN.LT.WW(3)) GO TO 20
C		IF YOU ARE HERE, NEWTON'S RULE DID NOT WORK: TRY BISECTION
15      WN = (WW(2) +WW(IH))/2.
		EP = ABS((WN-WB)/WB)
		IF (EP.LT.EPS) GO TO 100
20      IN = 0
        GO TO 130

100      ISUC = 1
        IF (IGP.NE.0) GO TO 130
        IGP = IGP + 1
        ISUC = 0
        WN = (WW(2) + WW(IL))/2.
        GO TO 20

C		YOU HAVE NOT SURROUNDED THE RESONANCE YET
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

C		ISUC = 0 MEANS TRY AGAIN
C		ISUC = 1 MEANS YOU HAVE SUCCEEDED TO WITHIN EPS

130     RETURN
        END

		SUBROUTINE CLCK(W,N,NN,DX)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		
C			SEARCH FOR CRITICAL LAYERS

	FN = FLOAT(N)
		
	IC = 0
	II = 0
	DO 10 I = 1,NN
		X = DX*FLOAT(I)
		WP = W + FN*V(I)/X
		IF (WP.GT.0.) GO TO 10
		IC = 1
		II = I
10	CONTINUE

	IF (IC.EQ.0) GO TO 20
	WRITE(8,901) II

901		FORMAT(//' SINGULARITY BEGINNING AT I = ',I3//)
20		RETURN
		END	

		SUBROUTINE DTIME(W,N,F,NN,DX,U)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)
		DIMENSION U(400)
		DOUBLE PRECISION R1,R2,R3

C			SUBROUTINE TO COMPUTE DECAY TIME OF WAVES

	FN = FLOAT(N)

C			COMPUTE STREAM FUNCTION

	DO 10 I = 1,NN
		X = DX*FLOAT(I)
		A(I) = -X*RH(I)*U(I)/FN
10	Continue

C			COMPUTE RADIAL DERIVATIVE

		DX2 = DX*2.0
		II = 1 + NN
		A(II) = A(2)/DX2
		II = 2*NN
		IM = NN - 1
		A(II) = (A(NN) - A(IM))/DX

	DO 20 I = 2,IM
		IPP = I + 1
		IMM = I - 1
		II = I + NN
		A(II) = (A(IPP)-A(IMM))/DX2
20	Continue
C			COMPUTE INTEGRALS

	R1 = 0.
	R2 = 0.
	R3 = 0.

	DO 30 I = 1,NN
		II = I + NN
		X = DX*FLOAT(I)
		WP = W + V(I)*FN/X
		RQ = (F + V(I)/X + VX(I))*(-RHX(I)/RH(I))
		RQ = RQ + VX(I)/X + VXX(I) - V(I)/(X*X)	
		RQ = RQ/RH(I)
		WPR = FN*(VX(I)/X -V(I)/(X*X))

		R1 = R1 + RQ*A(I)*A(I)/(WP*WP)
		XX = X*A(II)/(WP*RH(I)*RH(I))
		R2 = R2 + XX*(A(II) -WPR*A(I)/WP)
		R3 = R3 + A(I)*A(I)/(WP*X*RH(I)*RH(I))
30	CONTINUE

	R1 = R1*FN*DX
	R2 = R2*DX
	R3 = R3*DX*FN*FN

	DTF = -R1/(R2+R3)

	WRITE(8,901) DTF
	WRITE(8,903) R1,R2,R3

C			DEL = BOTTOM RESISTANCE COEFFICIENT (CM/SEC)
C			TF = FRICTIONAL DECAY TIME (SEC)

901		FORMAT(/' DEL*TF (CM) =',E15.5/)
903		FORMAT('R1, R2, R3 =',3E15.5)
		RETURN
		END




		SUBROUTINE STABCK(NN,DX,F)
		COMMON RH(400),RHX(400), A(1200), B(400), V(400), VX(400)
		COMMON VXX(400)

C			CHECK FOR BAROTROPIC STABILITY OF THE MEAN FLOW

C			COMPUTE VORTICITY GRADIENT

	DO 50 I = 1,NN
		X = DX*FLOAT(I)
		RQ = F + VX(I) + V(I)/X
		RQ = -RQ*RHX(I)/(RH(I)*RH(I))
		RQ = RQ + (VXX(I) + VX(I)/X - V(I)/(X*X))/RH(I)
		A(I) = RQ
50	Continue

C			CHECK FOR SIGN CHANGES
C			START BY LOOKING FOR NON-ZERO VALUES

		I = 1
		EPSS = 1.0E-20
55		AA = ABS(A(I))
		IF (AA.GT.EPSS) GO TO 60
		IF (I.GE.NN) GO TO 100
		I = I + 1
		GO TO 55

60		XX = A(I)
C			XX IS THE FIRST NONZERO VALUE OF THE GRADIENT
C			NOW LOOK FOR A SIGN CHANGE

	ITT = 0
	DO 70 I = 1,NN
		IF (ABS(A(I)).LT.EPSS) GO TO 70
		YY = A(I)/XX
		IF (YY.GE.0.0) GO TO 70
		ITT = 1
70	CONTINUE

		IF (ITT.NE.0) GO TO 80
		WRITE(8,901)
		GO TO 100

80		WRITE(8,902)
		WRITE(8,903) (A(I),I=1,NN)

C			POTENTIAL VORTICITY GRADIENT IN (SEC*CM*CM)-1

901		FORMAT(/' FLOW IS STABLE')
902		FORMAT(/' FLOW IS POTENTIALLY UNSTABLE, PV GRADIENT ='/)
903		FORMAT(1X,5E12.5)

100		RETURN
		END
