      SUBROUTINE MODES (NXTRA, IJFLAG)
      include 'param.inc'
      real*8 phasev(mdn), phasiv(mdn)
      INCLUDE 'COMMON.inc'
      DATA JIN / 10/   JOUT/ 11/
      PIE = 4 * ATAN (1.D0)
      WRITE (JOUT,*)
      DO 600 M = 1, MDN
      IF (IJ FLAG .EQ. 1) GO TO 53
 52   IF (NXTRA .GT. 0) GO TO 44
      READ (JIN,'(4f10.4, i10)') V,VI,STEP,STEPI,NXTRA
c      write (*,'(4f10.4, i10)') V,VI,STEP,STEPI,NXTRA
       IF (NXTRA .GE.0) GO TO 62
      V = V + VI * 1.D-10
      VI = STEP + STEPI * 1.D-10
      GO TO 85
 62   IF (V) 142,301,70
 142  IF(STEP) 44,44,143
C  SEARCH FOR MODE
 143  SIZE3 = -1.
      SIZE2=0
      IJ FLAG=1
      V=-V
      IF (NXTRA) 55,55,54
 55   NXTRA = 20
 54   XTRA = NXTRA
      HOP = (STEP - V) / XTRA
c      write (*,*) 'V, Step:', sngl(v), sngl(step) 
c      write (*,*) 'xtra, hop:', sngl(xtra), sngl(hop)
      HOPI=0.
      IF(STEPI.NE.0.) HOPI=(STEPI-VI)/XTRA
      DO 47 IJ = 1,NXTRA
      CALL SETUP
      DET = VEL
      DETI = VELI
      CALL DETNT(LX, VEL, VELI)
      DELTA = VEL
      DELTI = VELI
      SIZE = DELTA*DELTA + DELTI*DELTI
      if (K6 .LT. 2) WRITE (JOUT, 56) V, VI, SIZE, VEL, VELI
 56   FORMAT (2F12.3, 3D17.5)
      IF ((SIZE2.LT.SIZE3).AND.(SIZE.GT.SIZE2)) GO TO 45
   46 SIZE3=SIZE2
      SIZE2=SIZE
c      write (*,*) 'hop =', sngl(hop)
      V = V + HOP
      VI = VI + HOPI
      GO TO 47
 45   V = V - HOP
      TEMP = HOP / (SIZE - SIZE2)
      DELTI = TEMP * (DET * VELI - DETI * VEL)
      TEMP = .5D0 * (SIZE3 - SIZE) / (SIZE3 + SIZE - SIZE2 - SIZE2)
      DELTA = HOP * TEMP
      IF (HOPI .EQ. 0) GO TO 49
      VI=VI-HOPI
      DELTAI=HOPI*TEMP
      GO TO 49
 47   CONTINUE
      IJ FLAG=0
      NXTRA=0
      GO TO 52
   53 SIZE 2=-1.
      SIZE=0
      GO TO 46
 44   NXTRA = NXTRA - 1
      V = 3. * (phasev(M-1) - phasev(M-2)) + phasev(M-3)
      VI = 3.* (PHASIV(M-1) - phasiv(M-2)) + phasiv(M-3)
      STEP = (phasev(M-1) - phasev(M-2)) * .0001
 70   CALL SETUP
      CALL DETNT (LX, DET, DETI)
 80   FORMAT (/, F20.7, D20.11, 4D13.4)
 81   FORMAT (   F20.7, E20.11, 4D13.4)
      VEL = DET
      VELI = DETI
      DELTA = STEP
      DELTI = STEPI
      IF (DELTA .NE.0.) GO TO 49
      IF (DELTI .EQ.0.) DELTA = .01
 49   SIZE2 = 100.
      RX = DET**2 + DETI**2
      IF (K6 .EQ. 0) WRITE (JOUT, 80) V, VI, DET, DETI, SIZE, CNTR
      J = 0
 48   J = J + 1
      IF (J .GT. 15) GO TO 51
      V = V + DELTA
      VI = VI + DELTI
      IF (VI) 83,84,85
 83   DELTI = DELTI - VI
 84   VI = 1.D-18
 85   CALL SETUP
      DO 82 IA = 1, 2 * LX + 1
      DO 82 IB = 1,4
      BI(IA,IB) = Q(IA,IB)
 82   B(IA,IB) =  A(IA,IB)
      CALL DETNT (LX, DET, DETI)
      IF (K6 .GE. 2) GO TO 72
 71   WRITE (JOUT, 81) V, VI, DET, DETI, SIZE, CNTR
 72   IF (NXTRA .LT. 0) GO TO 51
      TEMNR =  DET * DELTA - DETI * DELTI
      TEMNI = DETI * DELTA + DET * DELTI
      TEMDR = VEL - DET
      TEMDI = VELI - DETI
      TEMDEN = TEMDR*TEMDR + TEMDI*TEMDI
      IF (TEMDEN .EQ. 0.) GO TO 51
      TEMRNU = TEMNR*TEMDR + TEMNI*TEMDI
      TEMINU = TEMNI*TEMDR - TEMNR*TEMDI
      DELTA = TEMRNU/TEMDEN
      DELTI = TEMINU/TEMDEN
C * * * THE NEXT CONSTANT DEPENDS ON WORD LENGTH AND SIZE OF PHASE VELOCITY * 
      IF (ABS(DELTA) .LT. 1.D-14) GO TO 51
      SIZE = DELTA*DELTA + DELTI*DELTI
      IF ((SIZE.GT.SIZE2).AND.(J.GT.3)) GO TO 51
 92   SIZE2 = SIZE * 2.
      VEL = DET
      VELI = DETI
      GO TO 48
C  Find Coefficients of Depth Functions:
   51 AF(1,M) = A(1,3)
      AG(1,M) = Q(1,3)
      BF(1,M) =-A(1,4)
      BG(1,M) =-Q(1,4)
      phasev(M) = V
      phasiv(M) = VI
      IF (K6 .GE. 2) GO TO 73
 74   WRITE (JOUT, 81) V, VI, DET, DETI, SIZE, CNTR
  73  IF (LX - 2) 95,96,97
 96   I = 0
      GO TO 98
 97   DO 110 L = 2, LX-1
      I = L + L - 2
      TEMNR = A(I,2)*AF( L-1 ,M) - Q(I,2)*AG( L-1 ,M)
      TEMNI = Q(I,2)*AF( L-1 ,M) + A(I,2)*AG( L-1 ,M)
      TEMDR = A(I,3)*A( I+1 ,4) - Q(I,3)*Q( I+1 ,4) -
     :   A(I,4)*A( I+1 ,3) + Q(I,4)*Q( I+1 ,3)
      TEMDI = Q(I,3)*A( I+1 ,4) + A(I,3)*Q( I+1 ,4)  -
     :    Q(I,4)*A( I+1 ,3) - A(I,4)*Q( I+1 ,3)
      TEMDEN = TEMDR*TEMDR + TEMDI*TEMDI
      TEMRNU = TEMNR*TEMDR + TEMNI*TEMDI
      TEMINU = TEMNI*TEMDR - TEMNR*TEMDI
      TEMP = TEMRNU / TEMDEN
      TEMPI = TEMINU / TEMDEN
      BF(L,M) = -(TEMP*A( I+1 ,4) - TEMPI*Q( I+1 ,4))
      BG(L,M) = -(TEMPI*A( I+1 ,4) + TEMP*Q( I+1 ,4))
      AG(L,M) = TEMPI*A( I+1 ,3) + TEMP*Q( I+1 ,3)
  110 AF(L,M) = TEMP*A( I+1 ,3) - TEMPI*Q( I+1 ,3)
 98   TEMNR = -(A(I+2,2) * AF(LX-1,M) -Q(I+2,2) * AG(LX-1,M))
      TEMNI = -(Q( I+2 ,2)*AF(LX-1,M) + A( I+2 ,2)*AG(LX-1,M))
      TEMDEN = A( I+2 ,3)*A( I+2 ,3) + Q( I+2 ,3)*Q( I+2 ,3)
      TEMRNU = TEMNR*A( I+2 ,3) + TEMNI*Q( I+2 ,3)
      TEMINU = TEMNI*A( I+2 ,3) - TEMNR*Q( I+2 ,3)
      BF(LX,M) = TEMRNU / TEMDEN
      BG(LX,M) = TEMINU / TEMDEN
 95   AF(LX,M) = 0.
      AG(LX,M) = 0.
C  FIND NORMALIZING FACTOR
      D(M) = 2.12429296D0 * rho(1)**3 / G(1)
      DI(M) = 0.
      DO 111 L = 2, LX
      TEMRSP = AF( L-1 ,M)*B( 2*L-2 ,2) - AG( L-1 ,M)*BI( 2*L-2 ,2)  +
     :   BF( L-1 ,M)*B( 2*L-2 ,1) - BG( L-1 ,M)*BI( 2*L-2 ,1)
      TEMISP = AG( L-1 ,M)*B( 2*L-2 ,2) + AF( L-1 ,M)*BI( 2*L-2 ,2)  +
     :   BG( L-1 ,M)*B( 2*L-2 ,1) + BF( L-1 ,M)*BI( 2*L-2 ,1)
      AX1 = TEMRSP*TEMRSP - TEMISP*TEMISP
      AX1I = TEMRSP * TEMISP
      AX1I = AX1I + AX1I
      TEMDR = (G(L-1)**2 + GI(L-1)**2)
      TEMDI = G(L)**2 + GI(L)**2
      TEMP = (RHO(L-1) / RHO(L)) / TEMDI
      TEM1  =(ZB (L-1) * G(L-1) + ZBI(L-1) * GI(L-1)) / TEMDR
     : -(ZT (L) * G(L) + ZTI(L) * GI(L)) * TEMP
      TEM1I =(ZBI(L-1) * G(L-1) - ZB (L-1) * GI(L-1)) / TEMDR
     : -(ZTI(L) * G(L) - ZT (L) * GI(L)) * TEMP
      TEMRSP = AF( L-1 ,M)*B( 2*L-1 ,2) - AG( L-1 ,M)*BI( 2*L-1 ,2)  +
     :   BF( L-1 ,M)*B( 2*L-1 ,1) - BG( L-1 ,M)*BI( 2*L-1 ,1)
      TEMISP = AG( L-1 ,M)*B( 2*L-1 ,2) + AF( L-1 ,M)*BI( 2*L-1 ,2)  +
     :   BG( L-1 ,M)*B( 2*L-1 ,1) + BF( L-1 ,M)*BI( 2*L-1 ,1)
      AX2 = TEMRSP * TEMRSP - TEMISP * TEMISP
      AX2I = TEMRSP * TEMISP
      AX2I = AX2I + AX2I
      TEMDR = rho(L-1) / (G CU(L-1)**2 + G CUI(L-1)**2)
      TEMDI = RHO(L  ) / (G CU(L  )**2 + G CUI(L )**2)
      TEM2 = G CU(L-1) * TEMDR - G CU(L) * TEMDI
      TEM2I = G CUI(L) * TEMDI - G CUI(L-1) * TEMDR
      TEMR1 = AX1*TEM1 - AX1I*TEM1I
      TEMI1 = AX1I*TEM1 + AX1*TEM1I
      TEMR2 = AX2 * TEM2 - AX2I * TEM2I
      TEMI2 = AX2I * TEM2 + AX2 * TEM2I
      D(M) =  D(M)  + TEMR1 /rho(L-1) + TEMR2
      DI(M) = DI(M) + TEMI1 /rho(L-1) + TEMI2
c      WD(M) = CMPLX(D(M), DI(M))
  111  CONTINUE
      EIGEN(M) = XI
      EIGENI(M) = XII
      IF (M .GT. 1) THEN
       FREQ = OMEGA /2 /PIE
       PHINV = V * phasev(M-1) /((V - phasev(M-1))* FREQ)
      END IF
      WRITE (JOUT,109) M, V, VI, D(M), DI(M), PHINV, - XII * 8686.D0 
  109 FORMAT (' Mode:',I4,'  Vel:',F12.4, '  Imag Vel:', 1p E12.5,
     : '  Norm:', 2E12.4, '  Int. Range:', 0p F8.0, '  L/K:', 1p e9.2)
  600  CONTINUE
  301 MX = M - 1
      RETURN
      END

