      PROGRAM EIGMOD
C
C---- FIND VERTICAL NORMAL MODES FOR A STEPWISE CONSTANT DENSITY PROFILE
C
      PARAMETER (II=6)
      REAL MATRIX(II,II),VECTOR(II,II),EIGVAL(II),WORK(90),DELP(II)
     .,SIGMA(II),RAD(II)
      DATA NDIM/6/,EPS1/-1./
C
C --- READ IN PROFILE (LAYER THICKNESS IN METERS, DENSITY IN SIGMA UNITS)
C
      READ 400,CORIO
 400  FORMAT(E12.4)
      PRINT 400,CORIO
1     READ 100,NLAY
      IF (NLAY.EQ.-1) STOP 
      PRINT 100,NLAY
      DO 201 N=1,NLAY
201   READ 200,DELP(N),SIGMA(N)
100   FORMAT (I2)
200   FORMAT (2F7.1)
      DO 20 I=1,NDIM
      DO 20 J=1,NDIM
20    MATRIX(I,J)=0.
C
      DO 2 I=1,NLAY
      DO 2 J=1,I
      MATRIX(I,J)=SQRT(DELP(I)*DELP(J))/(1.+AMAX1(SIGMA(I),SIGMA(J))
     .*1.E-3)
2     MATRIX(J,I)=MATRIX(I,J)
C --- CALL TO EIGENVALUE/EIGENVEVTOR ROUTINE
C     CALL EIGSFS(NDIM,NLAY,MATRIX,EIGVAL,VECTOR,EPS1,0.,6000.,NLAY,1
C    .,NO,IERR,WORK)
      CALL EVCSF(NLAY,MATRIX,NDIM,EIGVAL,VECTOR,NDIM)
      NO=NLAY
C --- NORMALIZE EIGENVECTORS TO UNIT KINETIC ENERGY
      DO 4 I=1,NO
      SUM=0.
      DO 5 N=1,NLAY
5     SUM=SUM+VECTOR(N,I)**2
      SUM=SIGN(SQRT(SUM),VECTOR(1,I))
      DO 4 N=1,NLAY
4     VECTOR(N,I)=VECTOR(N,I)/(SUM*SQRT(DELP(N)))
C
      DO 300 I=1,NO-1
 300     RAD(I)=1.E-3*SQRT(9.8*EIGVAL(I))/CORIO
      PRINT 103,(DELP(N),N=1,NLAY)
103   FORMAT (///'DENSITY PROFILE'/1X,10F12.1)
      PRINT 104,(SIGMA(N),N=1,NLAY)
104   FORMAT (1X,10F12.2)
      PRINT 101,(EIGVAL(I),I=1,NO)
101   FORMAT ('EIGENVALUES (EQUIVALENT DEPTHS)'/1X,1P10E12.4)
      PRINT 106,(RAD(I),I=1,NO-1)
 106  FORMAT ('ROSSBY RADII OF DEFORMATION (IN KM)'/1X,10E12.4)
      PRINT 105
105   FORMAT ('EIGENVECTORS (VELOCITY PROFILES)')
      DO 3 N=1,NLAY
3     PRINT 102,(VECTOR(N,I),I=1,NO)
102   FORMAT (1X,4P10F12.1)
      GO TO 1
      END
