* PACKAGE LSLFRQ  !" 固有値解析 振動数選択
*
*"  [HIS] 93/02/19(takepiro) 
*
**********************************************************************
      SUBROUTINE LSLFRQ
     M           (  GROWTH , FRQENC ,
     M              EIGVCR , EIGVCI ,
     O              NUMEIG ,
     I              NMKMAX              )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZLDIM)                      !" パラメター数, 係数行列
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zldim.F"                    !" パラメター数, 係数行列
#endif
*
*   [MODIFY]
      REAL       GROWTH( MATDMX )            !" 成長率(固有値実数部)
      REAL       FRQENC( MATDMX )            !" 振動数(固有値虚数部)
      REAL       EIGVCR( MATDMX, MATDMX )    !" 固有ベクトル実数部
      REAL       EIGVCI( MATDMX, MATDMX )    !" 固有ベクトル虚数部
*
*   [OUTPUT]
      INTEGER    NUMEIG                      !" 固有値の数
*
*   [INPUT]
      INTEGER    NMKMAX                      !" 係数行列の最大添え字
*
*   [INTERNAL WORK]
      COMMON    /COMWRK/ WORK
      INTEGER    JORDER( MATDMX )            !" 入れ換え順序
      REAL       WORK  ( MATDMX, MATDMX )    !" 入れ換え用作業領域
      COMPLEX    ACS , ASN, EI
      DATA       EI/ ( 0.0, 1.0 )/
      INTEGER    I, J, JN, K
*
      JN = 0 
*
      DO 1000 J = 1, NMKMAX
         IF ( FRQENC( J ) . NE. 0 ) THEN
            DO 1100 I = 2, NMKMAX, 2
               ACS = EIGVCR( I-1,J ) + EI * EIGVCI( I-1,J )
               ASN = - ( EIGVCR( I,J )   + EI * EIGVCI( I,J ) )
               IF( ABS( ACS ) .NE. 0 )THEN
                  IF ( IMAG ( ASN/ACS ) .GT. 0 )THEN
                     JN = JN + 1
                     JORDER( JN ) = J 
                  ENDIF
                  GOTO 1000
               ENDIF
 1100       CONTINUE 
         ELSE
            JN = JN + 1
            JORDER( JN ) = J
         ENDIF
 1000 CONTINUE 
*
      NUMEIG = JN
*
*"  < 2. 固有値, 固有ベクトルのソート ; 固有値虚数部→振動数>
*
      CALL COPY( WORK, GROWTH, NMKMAX )
      DO 2100 J = 1, NUMEIG
         GROWTH( J ) = WORK( JORDER(J),1  )
 2100 CONTINUE
*
      CALL COPY ( WORK, FRQENC , NMKMAX )
      DO 2200 J = 1, NUMEIG
         FRQENC( J ) = - WORK( JORDER(J),1  )
 2200 CONTINUE
*
      CALL COPY( WORK, EIGVCR, MATDMX*NMKMAX )
      DO 2300 K = 1, NUMEIG
         DO 2300 J = 1, NMKMAX
            EIGVCR( J,K ) = WORK( J,JORDER(K) )
 2300 CONTINUE
*
      CALL COPY( WORK, EIGVCI, MATDMX*NMKMAX )
      DO 2400 K = 1, NUMEIG
         DO 2400 J = 1, NMKMAX
            EIGVCI( J,K ) = WORK( J,JORDER(K) )
 2400 CONTINUE
*
      RETURN
      END
