* PACKAGE LADMN  !" 線形力学固有値解析
*"                  ver 1.05    93/07/14 takepiro
*"                              95/06/21(takepiro)  出力固有関数計算数導入
**********************************************************************
      SUBROUTINE LANAL                      !"  線形力学解析
     O             ( GAU   , GAV   , GAW   , GAT   ,
     O               GATOR , GAPOR , 
     O               GRMAX , FREQ  , 
     I               MFIX  , 
     I               RAYL  , TAU   , PRND  , 
     I               ORSTR , MODE  , 
     C               ALON  , DLON  , ALAT  , DLAT  ,
     C               ARAD  , DRAD                    )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZLDIM)                      !" パラメター数, 係数行列
      INCLUDE   (ZHDIM)                      !" 文字列文字数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zldim.F"                    !" パラメター数, 係数行列
#include        "zhdim.F"                    !" 文字列文字数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT]
      REAL       GAU   ( IDIM, JDIM, 0:KDIM )  !" 西風   ｕ
      REAL       GAV   ( IDIM, JDIM, 0:KDIM )  !" 南風   ｖ
      REAL       GAW   ( IDIM, JDIM, 0:KDIM )  !" 鉛直風 ｗ
      REAL       GAT   ( IDIM, JDIM, 0:KDIM )  !" 温度   Ｔ
      REAL       GATOR ( IDIM, JDIM, 0:KDIM )  !" トロイダル Ψ
      REAL       GAPOR ( IDIM, JDIM, 0:KDIM )  !" ポロイダル Φ
*"       : 格子点データ(t) <DYNMCS>  格子点データ(t+Δt)
*
      REAL       GRMAX                       !" 成長率
      REAL       FREQ                        !" 振動数
*
*   [INPUT]
      REAL       ALON  ( IDIM )              !" 経度
      REAL       DLON  ( IDIM )              !" 経度荷重
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
      REAL       ARAD  ( 0:KDIM )            !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )            !" Δｒ(整数)
*"        : 座標格子
*
      REAL       RAYL                        !" レイリー数
      REAL       TAU                         !" √テイラー数
      REAL       PRND                        !" プランドル数
*
      INTEGER    MFIX                        !" 東西波数
*
      LOGICAL    ORSTR                       !" 固有関数出力するか否か
      INTEGER    MODE                        !" 出力固有関数のモード番号
*                                            !" ( 1 で最大成長率 )
*
*   [INTERNAL SAVE]
*
*   [INTERNAL ONCE]
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       EDEL  ( NMDIM  )            !" ζ，Ｄ→Ｕ，Ｖ
*
      REAL       UVFACT ( IDIM, JDIM )       !" u→U のファクター
*
      SAVE       NMO   , FLAPLA, EDEL  ,
     &           UVFACT
*
      LOGICAL    OSETC , OFIRST
      DATA       OSETC , OFIRST  / .FALSE., .TRUE.  /
      SAVE       OSETC , OFIRST
*
      LOGICAL    OATOR , OAPOR , OAT         !" 解析対象物理量
      SAVE       OATOR , OAPOR , OAT
*
*   [INTERNAL WORK]
      REAL       COFMTX( MATDMX , MATDMX )   !" 係数行列
      REAL       GROWTH( MATDMX )            !" 成長率(固有値実数部)
      REAL       FRQENC( MATDMX )            !" 振動数(固有値虚数部)
*
      REAL       EIGVCR( MATDMX, NVDIM )     !" 固有ベクトル実数部
      REAL       EIGVCI( MATDMX, NVDIM )     !" 固有ベクトル虚数部
*
      INTEGER    NMKMTX ( NMDIM, 0:KDIM, 3 ) !" 係数行列の添え字
      INTEGER    NMKMAX                      !" 係数行列の最大添え字
      INTEGER    NMTLST  ( NMDIM )           !" MFIX 成分リスト(トロイダル)
      INTEGER    NMPLST  ( NMDIM )           !" MFIX 成分リスト(ポロイダル)
      INTEGER    NMQLST  ( NMDIM )           !" MFIX 成分リスト(温度)
      INTEGER    NMTMAX                      !" NMTLST の最大添え字
      INTEGER    NMPMAX                      !" NMPLST の最大添え字
      INTEGER    NMQMAX                      !" NMQLST の最大添え字
      INTEGER    NUMEIG                      !" 固有値の数
*
      INTEGER    IERR                        !" サブルーチンエラーコード
*"       : 固有値計算用行列
*
      INTEGER    IUNIT                       !" 標準出力装置番号
*
c$$$      COMMON     /COMWRK/ EIGVCR , EIGVCI
*
*"  << SETCON : 定数設定 >>
*
      IF ( MODE .GT. NVDIM ) THEN
         CALL MSGDMP( 'E', 'LADMN', 'MODE NUMBER IS TOO LARGE' )
      ENDIF
*
*   [ONCE]
      IF ( .NOT. OSETC ) THEN
         OSETC = .TRUE.
         CALL    DSETC                 !" スペクトル定数設定
     O         ( NMO   , FLAPLA, EDEL  , UVFACT, 
     C           ALAT  , DLAT                      )
      ENDIF
*   [ONCE]
      IF ( OFIRST ) THEN
         OFIRST = .FALSE.
         CALL LSTPHS                         !" 解析対象物理量の設定
     O         ( OATOR  , OAPOR  , OAT  )
      ENDIF
*
*" << 係数行列の計算 >>
*
      CALL CLCSTR( 'LSTMTX' )
*
      CALL MKMLST
     O      ( NMKMTX , NMKMAX ,
     O        NMTLST , NMTMAX , 
     O        NMPLST , NMPMAX , 
     O        NMQLST , NMQMAX , 
     F        OATOR  , OAPOR  , OAT   , 
     I        NMO    , MFIX             )
*
      CALL LSTMTX
     O      ( COFMTX ,
     I        NMKMTX , NMKMAX , 
     I        NMTLST , NMTMAX , 
     I        NMPLST , NMPMAX , 
     I        NMQLST , NMQMAX , 
     I        RAYL   , TAU    , PRND   ,
     F        OATOR  , OAPOR  , OAT    , 
     C        NMO    , FLAPLA , 
     C        ARAD   , DRAD                       )
*
      CALL CLCEND( 'LSTMTX' )
*
*" << 固有値計算 >>
*
      CALL CLCSTR( 'LEIGEN' )
*
      CALL LEIGEN
     M        ( COFMTX , 
     O          GROWTH , FRQENC ,
     O          EIGVCR , EIGVCI ,
     O          IERR   , 
     I          NMKMAX , .TRUE. , NVDIM     )
*
      IF( IERR .NE. 0 )THEN
         CALL MSGDMP( 'W', 'LADMN', 'ERROR OCCURED AT ' )
         CALL GLPGET( 'MSGUNIT', IUNIT )
         WRITE ( IUNIT, * ) 'Rayleigh Number = ', RAYL
         WRITE ( IUNIT, * ) 'Taylor Number   = ', TAU**2
         WRITE ( IUNIT, * ) 'Prandtl Number  = ', PRND
         WRITE ( IUNIT, * ) 'E-W Wave Numver = ', MFIX
      ENDIF
*
      CALL CLCEND( 'LEIGEN' )
*
      IF ( MFIX .NE. 0 )THEN
         CALL LSLFRQ
     M         (  GROWTH , FRQENC ,
     M            EIGVCR , EIGVCI ,
     O            NUMEIG ,
     I            NMKMAX , NVDIM      )
      ELSE
         NUMEIG = NMKMAX
      ENDIF
*
*" << 固有値診断出力 >>
*
      CALL CLCSTR( 'LDIAG' )
*
      CALL LDIAG
     I      (  GROWTH , FRQENC , 
     I         EIGVCR , EIGVCI ,
     I         NMKMTX , NUMEIG , NVDIM  ,
     I         NMTLST , NMTMAX , 
     I         NMPLST , NMPMAX , 
     I         NMQLST , NMQMAX , 
     F         OATOR  , OAPOR  , OAT    ,
     C         NMO    , FLAPLA , UVFACT , 
     C         ARAD   , DRAD   ,
     W         GAU    , GAV    , GAW    , GAT   ,
     W         GATOR  , GAPOR  )
*
      CALL CLCEND( 'LDIAG' )
*
*" << 固有関数計算 >>
*
      CALL CLCSTR( 'LEIVEC' )
*
      GRMAX = GROWTH( MODE )
      FREQ  = FRQENC( MODE )
*
      IF ( ORSTR ) THEN
         CALL LEIVEC
     O       ( GAU    , GAV    , GAW    , GAT    ,
     O         GATOR  , GAPOR  , 
     I         EIGVCR(1,MODE) , EIGVCI(1,MODE) ,
     I         NMKMTX , 
     I         NMTLST , NMTMAX , 
     I         NMPLST , NMPMAX , 
     I         NMQLST , NMQMAX , 
     F         OATOR  , OAPOR  , OAT    ,
     C         NMO    , FLAPLA , UVFACT , 
     C         ARAD   , DRAD                       )
      ENDIF
*
      CALL CLCEND( 'LEIVEC' )
*
      RETURN
      END
**********************************************************************
      SUBROUTINE LSTPHS                         !" 解析対象物理量の設定
     O         ( OATOR  , OAPOR  , OAT  )
*
*"  [OUTPUT]
      LOGICAL    OATOR                       !" トロイダル解析スイッチ
      LOGICAL    OAPOR                       !" ポロイダル解析スイッチ
      LOGiCAL    OAT                         !" 温度解析スイッチ
*
*"  [INTERNAL ONCE]
      LOGICAL    OTOR                        !" トロイダル解析スイッチ
      LOGICAL    OPOR                        !" ポロイダル解析スイッチ
      LOGiCAL    OTEMP                       !" 温度解析スイッチ
      DATA       OTOR  ,  OPOR  ,  OTEMP 
     &        / .TRUE. , .TRUE. , .TRUE. /
*
      NAMELIST   /NMLPHS/  OTOR , OPOR , OTEMP
*
*   [INTERNAL WORK]
      INTEGER    IFPAR, JFPAR
*
*"     < 1. 解析対象物理量の設定 >
*
C      WRITE  ( JFPAR, NMLPHS )
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMLPHS , END=1190 )
 1190 WRITE  ( JFPAR, NMLPHS )
*
      OATOR = OTOR
      OAPOR = OPOR
      OAT   = OTEMP
*
      IF ( (.NOT.OATOR ) .OR. (.NOT. OAPOR ) .OR. (.NOT. OAT ) )THEN
         CALL MSGDMP( 'W', 'LSTPHS', 'ANALIZING RESTRICTED DYNAMICS' )
      ENDIF
*
      RETURN
      END