* PACKAGE DSTMTX  !" 線形力学固有値解析 : 行列設定
*
*"[HIS]  93/01/19(takepiro)
*"       93/02/19(takepiro) 行 - 列 の設定変更
*"       93/04/30(takepiro) トロイダル, ポロイダルの波数設定分離
*"       93/06/27(takepiro) 繰り返し行列計算高速化
*"       93/07/10(takepiro) 波数リストチェック
*"       93/07/13(takepiro) ポロイダル, 温度の波数設定分離
*"       93/07/29(takepiro) 
*
**********************************************************************
      SUBROUTINE 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                       )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZLDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zldim.F"                    !" 係数行列の大きさ
#endif
*
*"  [OUTPUT]
      REAL       COFMTX( MATDMX , MATDMX )   !" 係数行列
*
*"  [INPUT]
      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 の最大添え字
*"       : 固有値計算用行列
*
      REAL       RAYL                        !" レイリー数
      REAL       TAU                         !" √テイラー数
      REAL       PRND                        !" プランドル数
*
      LOGICAL    OATOR , OAPOR , OAT         !" 解析対象物理量
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
*
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
*"  [INTERNAL SAVE]
      REAL       RATAMX( MATDMX , MATDMX ) !" レイリー/テイラー数依存係数行列
      REAL       PRNDMX( MATDMX , MATDMX ) !" プランドル数依存係数行列
      REAL       TAUS                      !" √テイラー数
*
      INTEGER    NMTS  ( NMDIM )           !" MFIX 成分リスト(トロイダル)
      INTEGER    NMPS  ( NMDIM )           !" MFIX 成分リスト(ポロイダル)
      INTEGER    NMQS  ( NMDIM )           !" MFIX 成分リスト(温度)
      INTEGER    NMTMX                     !" NMTLST の最大添え字
      INTEGER    NMPMX                     !" NMPLST の最大添え字
      INTEGER    NMQMX                     !" NMQLST の最大添え字
      DATA       NMTMX, NMPMX , NMQMX/ 0, 0, 0/
*
      LOGICAL    OMTRIX                    !" 行列計算スイッチ
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
*
      SAVE       RATAMX , PRNDMX , 
     &           NMTS   , NMPS   , NMQS , 
     &           NMTMX  , NMPMX  , NMQMX , 
     &           OFIRST , OMTRIX
*
*"  [INTERNAL WORK]
      REAL       WDTOR ( NMDIM , 0:KDIM     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , 0:KDIM     )  !" ポロイダル Φ
      REAL       WDDPOR( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
      REAL       WDT   ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
*
      REAL       WTTOR ( NMDIM, 0:KDIM      )  !" トロイダル Ψ
      REAL       WTPOR ( NMDIM, 0:KDIM      )  !" ポロイダル Φ
      REAL       WTT   ( NMDIM, 0:KDIM      )  !" 温度  Ｔ
*
      INTEGER    NML , NMLJ , K , KJ
*
      IF ( OFIRST )THEN
         OFIRST = .FALSE. 
         WRITE ( 6,* ) ' INTERACTION MATRIX FAST : DATE=93/07/29'
      ENDIF
*
*======================================================================
*"  << 0. 物理量間の相互作用行列設定 >>
*
*"    < 0.01 波数リストチェック : 前と異なっていたら行列を計算 >
*
      IF(       (NMTMX.EQ.NMTMAX) 
     &    .AND. (NMPMX.EQ.NMPMAX) 
     &    .AND. (NMQMX.EQ.NMQMAX)    )THEN
*
         DO 10 NML = 1, NMTMAX
            IF( NMTS(NML) .NE. NMTLST(NML) ) THEN
               OMTRIX = .TRUE.
               GOTO 99
            ENDIF
 10      CONTINUE 
         DO 20 NML = 1, NMPMAX
            IF( NMPS(NML) .NE. NMPLST(NML) ) THEN
               OMTRIX = .TRUE.
               GOTO 99
            ENDIF
 20      CONTINUE 
         DO 30 NML = 1, NMQMAX
            IF( NMQS(NML) .NE. NMQLST(NML) ) THEN
               OMTRIX = .TRUE.
               GOTO 99
            ENDIF
 30      CONTINUE 
         OMTRIX = .FALSE.
      ELSE
         OMTRIX = .TRUE.
      ENDIF
 99   CONTINUE 
*
      IF( OMTRIX )THEN  !" <------------------------------- OMTRIX 開始
*
*----------------------------------------------------------------------
*"    < 0.02 波数情報セーブ >
*
         CALL MSGDMP( 'M','LSTMTX', 'EVALUATE INTERACTION MATRIX' )
*
         DO 51 NML = 1, NMDIM
            NMTS( NML ) = NMTLST( NML )
            NMPS( NML ) = NMPLST( NML )
            NMQS( NML ) = NMQLST( NML )
 51      CONTINUE 
*
         NMTMX = NMTMAX
         NMPMX = NMPMAX
         NMQMX = NMQMAX
*
         IF ( TAU .GT. 1.0 )THEN
            TAUS = TAU
         ELSE
            TAUS = 1.0
         ENDIF
*
*----------------------------------------------------------------------
*"    < 0.1 トロイダル >
*
         IF ( OATOR ) THEN
            DO 100 K = 1, KDIM-1
            DO 100 NML = 1, NMTMAX
*
*"          < 0.1.1 トロイダル拡散 >
*
               CALL RESET ( WDTOR , NMDIM*KMAX )
               CALL RESET ( WDPOR , NMDIM*KMAX )
               CALL RESET ( WDDPOR, NMDIM*KMAX )
               CALL RESET ( WDT   , NMDIM*KMAX )
*
               WDTOR( NMTLST(NML), K ) = 1.0
*
               CALL LANDYN      !" 力学 線型項(回転なし)
     O               ( WTTOR , WTPOR , WTT   ,
     M                 WDTOR , WDPOR , WDT   , WDDPOR ,  
     C                 0.0   , 0.0   , 1.0   , 
     C                 NMO   , FLAPLA, 
     C                 ARAD  , DRAD                      )
*
               DO 101 KJ = 1, KDIM-1
               DO 101 NMLJ = 1, NMTMAX
                  PRNDMX( NMKMTX( NMLJ,KJ,1 ), NMKMTX( NML,K,1 ) )
     &               = WTTOR( NMTLST(NMLJ), KJ )
 101           CONTINUE 
*
*"          < 0.1.2 トロイダル回転依存 >
*
               CALL RESET ( WDTOR , NMDIM*KMAX )
               CALL RESET ( WDPOR , NMDIM*KMAX )
               CALL RESET ( WDDPOR, NMDIM*KMAX )
               CALL RESET ( WDT   , NMDIM*KMAX )
*
               WDTOR( NMTLST(NML), K ) = 1.0
*
               CALL LANDYN                 !" 力学線型項(回転あり)
     O                ( WTTOR , WTPOR , WTT   ,
     M                  WDTOR , WDPOR , WDT   , WDDPOR ,  
     C                  0.0   , TAUS  , 1.0   , 
     C                  NMO   , FLAPLA, 
     C                  ARAD  , DRAD                      )
*
               DO 110 KJ = 1, KDIM-1
               DO 110 NMLJ = 1, NMTMAX
                  RATAMX( NMKMTX( NMLJ,KJ,1 ), NMKMTX( NML,K,1 ) )
     &               = WTTOR( NMTLST(NMLJ), KJ )
     &                - PRNDMX( NMKMTX( NMLJ,KJ,1 ), NMKMTX( NML,K,1 ) )
 110           CONTINUE 

               IF ( OAPOR ) THEN
                  DO 120 KJ = 2, KDIM-2
                  DO 120 NMLJ = 1, NMPMAX
                     RATAMX( NMKMTX( NMLJ,KJ,2 ), NMKMTX( NML,K,1 ) )
     &                 = WTPOR( NMPLST(NMLJ), KJ )
 120              CONTINUE 
               ENDIF
*
 100        CONTINUE 
         ENDIF
*----------------------------------------------------------------------
*"    < 0.2. ポロイダル >
*
         IF ( OAPOR ) THEN
            DO 200 K = 2, KDIM-2
            DO 200 NML = 1, NMPMAX
*
*"        < 0.2.1 ポロイダル拡散 >
*
               CALL RESET ( WDTOR , NMDIM*KMAX )
               CALL RESET ( WDPOR , NMDIM*KMAX )
               CALL RESET ( WDDPOR, NMDIM*KMAX )
               CALL RESET ( WDT   , NMDIM*KMAX )
*
               WDDPOR( NMPLST(NML), K ) = 1.0
*
               CALL LANDYN      !" 力学 線型項
     O                ( WTTOR , WTPOR , WTT   ,
     M                  WDTOR , WDPOR , WDT   , WDDPOR ,  
     C                  0.0   , 0.0   , 1.0   , 
     C                  NMO   , FLAPLA, 
     C                  ARAD  , DRAD                    )
*
               DO 201 KJ = 2, KDIM-2
               DO 201 NMLJ = 1, NMPMAX
                  PRNDMX( NMKMTX( NMLJ,KJ,2 ), NMKMTX( NML,K,2 )  )
     &                    = WTPOR( NMPLST(NMLJ), KJ )
 201           CONTINUE 
*
*"        < 0.2.2 ポロイダル回転/温度傾度依存 >
*
               CALL RESET ( WDTOR , NMDIM*KMAX )
               CALL RESET ( WDPOR , NMDIM*KMAX )
               CALL RESET ( WDDPOR, NMDIM*KMAX )
               CALL RESET ( WDT   , NMDIM*KMAX )
*
               WDDPOR( NMPLST(NML), K ) = 1.0
*
               CALL LANDYN      !" 力学 線型項(回転, 温度傾度あり)
     O                ( WTTOR , WTPOR , WTT   ,
     M                  WDTOR , WDPOR , WDT   , WDDPOR ,  
     C                  0.0   , TAUS  , 1.0   , 
     C                  NMO   , FLAPLA, 
     C                  ARAD  , DRAD                    )
*
               IF ( OATOR ) THEN
                  DO 210 KJ = 1, KDIM-1
                  DO 210 NMLJ = 1, NMTMAX
                     RATAMX( NMKMTX( NMLJ,KJ,1 ), NMKMTX( NML,K,2 )  )
     &                    = WTTOR( NMTLST(NMLJ), KJ )
 210              CONTINUE 
               ENDIF
*
               IF ( OAPOR ) THEN
                  DO 220 KJ = 2, KDIM-2
                  DO 220 NMLJ = 1, NMPMAX
                     RATAMX( NMKMTX( NMLJ,KJ,2 ), NMKMTX( NML,K,2 ) )
     &                 = WTPOR( NMPLST(NMLJ), KJ )
     &                 - PRNDMX( NMKMTX(NMLJ,KJ,2), NMKMTX(NML,K,2) )
 220              CONTINUE 
               ENDIF
*
               IF ( OAT ) THEN
                  DO 230 KJ = 1, KDIM-1
                  DO 230 NMLJ = 1, NMQMAX
                     PRNDMX( NMKMTX( NMLJ,KJ,3 ), NMKMTX( NML,K,2 )  )
     &                    = WTT( NMQLST(NMLJ), KJ )
 230              CONTINUE 
               ENDIF
 200        CONTINUE 
         ENDIF
*----------------------------------------------------------------------
*"    < 0.3 温度 >
*
         IF ( OAT ) THEN
            DO 300 K = 1, KDIM-1
            DO 300 NML = 1, NMQMAX
*
*"        < 0.3.1 温度拡散/基本場温度傾度移流 >
*
               CALL RESET ( WDTOR , NMDIM*KMAX )
               CALL RESET ( WDPOR , NMDIM*KMAX )
               CALL RESET ( WDDPOR, NMDIM*KMAX )
               CALL RESET ( WDT   , NMDIM*KMAX )
*
               WDT( NMQLST(NML), K ) = 1.0
*
               CALL LANDYN      !" 力学 線型項
     O                ( WTTOR , WTPOR , WTT   ,
     M                  WDTOR , WDPOR , WDT   , WDDPOR ,  
     C                  1.0   , 0.0   , 1.0   , 
     C                  NMO   , FLAPLA, 
     C                  ARAD  , DRAD                    )
*
               DO 301 KJ = 1, KDIM-1
               DO 301 NMLJ = 1, NMQMAX
                  PRNDMX( NMKMTX( NMLJ,KJ,3 ), NMKMTX( NML,K,3 ) )
     &                    = WTT( NMQLST(NMLJ), KJ )
 301           CONTINUE 
*
               IF ( OAPOR ) THEN
                  DO 320 KJ = 2, KDIM-2
                  DO 320 NMLJ = 1, NMPMAX
                     RATAMX( NMKMTX( NMLJ,KJ,2 ), NMKMTX( NML,K,3 ) )
     &                    = WTPOR( NMPLST(NMLJ), KJ )
 320              CONTINUE 
               ENDIF
*
 300        CONTINUE 
         ENDIF
*
      ENDIF                 !"  <-------------------------- OMTRIX 終り
*======================================================================
*" << 1. 係数行列の設定 >>
*
*----------------------------------------------------------------------
*"     < 1. トロイダル >
*
      IF ( OATOR ) THEN
         DO 1000 K = 1, KDIM-1
         DO 1000 NML = 1, NMTMAX
*
            IF ( OATOR ) THEN
               DO 1100 KJ = 1, KDIM-1
               DO 1100 NMLJ = 1, NMTMAX
                  COFMTX( NMKMTX( NMLJ,KJ,1 ), NMKMTX( NML,K,1 ) )
     &              =  PRND 
     &                  * PRNDMX( NMKMTX(NMLJ,KJ,1), NMKMTX(NML,K,1) )
     &               + PRND * TAU / TAUS 
     &                  * RATAMX( NMKMTX(NMLJ,KJ,1), NMKMTX(NML,K,1) )
 1100          CONTINUE 
            ENDIF
*
            IF ( OAPOR ) THEN
               DO 1200 KJ = 2, KDIM-2
               DO 1200 NMLJ = 1, NMPMAX
                  COFMTX( NMKMTX( NMLJ,KJ,2 ), NMKMTX( NML,K,1 ) )
     &              = PRND * TAU / TAUS 
     &                  * RATAMX( NMKMTX(NMLJ,KJ,2), NMKMTX(NML,K,1) )
 1200          CONTINUE 
            ENDIF
*
            IF ( OAT ) THEN
               DO 1300 KJ = 1, KDIM-1
               DO 1300 NMLJ = 1, NMQMAX
                  COFMTX( NMKMTX( NMLJ,KJ,3 ), NMKMTX( NML,K,1 ) ) = 0. 
 1300          CONTINUE 
            ENDIF
 1000    CONTINUE 
      ENDIF
*----------------------------------------------------------------------
*"   < 2. ポロイダル >
*
      IF ( OAPOR ) THEN
         DO 2000 K = 2, KDIM-2
         DO 2000 NML = 1, NMPMAX
*
            IF ( OATOR ) THEN
               DO 2100 KJ = 1, KDIM-1
               DO 2100 NMLJ = 1, NMTMAX
                  COFMTX(  NMKMTX( NMLJ,KJ,1 ), NMKMTX( NML,K,2 )  )
     &              = PRND * TAU / TAUS
     &                  * RATAMX( NMKMTX(NMLJ,KJ,1), NMKMTX(NML,K,2) )
 2100          CONTINUE 
            ENDIF
*
            IF ( OAPOR ) THEN
               DO 2200 KJ = 2, KDIM-2
               DO 2200 NMLJ = 1, NMPMAX
                  COFMTX(  NMKMTX( NMLJ,KJ,2 ), NMKMTX( NML,K,2 )  )
     &              =  PRND 
     &                  * PRNDMX( NMKMTX(NMLJ,KJ,2), NMKMTX(NML,K,2) )
     &               + PRND * TAU / TAUS
     &                  * RATAMX( NMKMTX(NMLJ,KJ,2), NMKMTX(NML,K,2) )
 2200          CONTINUE 
            ENDIF
*
            IF ( OAT ) THEN
               DO 2300 KJ = 1, KDIM-1
               DO 2300 NMLJ = 1, NMQMAX
                  COFMTX( NMKMTX( NMLJ,KJ,3 ), NMKMTX( NML,K,2 ) )
     &              = PRNDMX( NMKMTX(NMLJ,KJ,3), NMKMTX(NML,K,2) )
 2300          CONTINUE 
            ENDIF
 2000    CONTINUE 
      ENDIF
*----------------------------------------------------------------------
*"     < 3. 温度 >
*
      IF ( OAT ) THEN
         DO 3000 K = 1, KDIM-1
         DO 3000 NML = 1, NMQMAX
*
            IF ( OATOR ) THEN
               DO 3100 KJ = 1, KDIM-1
               DO 3100 NMLJ = 1, NMTMAX
                  COFMTX( NMKMTX( NMLJ,KJ,1 ), NMKMTX( NML,K,3 ) ) = 0. 
 3100          CONTINUE 
            ENDIF
*
            IF ( OAPOR ) THEN
               DO 3200 KJ = 2, KDIM-2
               DO 3200 NMLJ = 1, NMPMAX
                  COFMTX( NMKMTX( NMLJ,KJ,2 ), NMKMTX( NML,K,3 ) )
     &              = PRND * RAYL
     &                 * RATAMX( NMKMTX(NMLJ,KJ,2), NMKMTX(NML,K,3) )
 3200          CONTINUE 
            ENDIF
*
            IF ( OAT ) THEN
               DO 3300 KJ = 1, KDIM-1
               DO 3300 NMLJ = 1, NMQMAX
                  COFMTX( NMKMTX( NMLJ,KJ,3 ), NMKMTX( NML,K,3 ) )
     &              =  PRNDMX( NMKMTX(NMLJ,KJ,3), NMKMTX(NML,K,3) )
 3300          CONTINUE 
            ENDIF
 3000    CONTINUE 
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE LANDYN          !" 力学解析用線型力学
     O           ( WTTOR , WTPOR , WTT   ,
     M             WDTOR , WDPOR , WDT   , WDDPOR ,  
     C             RAYL  , TAU   , PRND  , 
     C             NMO   , FLAPLA, 
     C             ARAD  , DRAD                    )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#endif
*
*  [MODIFY]
      REAL       WTTOR ( NMDIM, 0:KDIM      )  !" トロイダル Ψ
      REAL       WTPOR ( NMDIM, 0:KDIM      )  !" ポロイダル Φ
      REAL       WTT   ( NMDIM, 0:KDIM      )  !" 温度  Ｔ
*
*  [INPUT]
      REAL       WDTOR ( NMDIM , 0:KDIM     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , 0:KDIM     )  !" ポロイダル Φ
      REAL       WDDPOR( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
      REAL       WDT   ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
*
      REAL       RAYL                        !" レイリー数
      REAL       TAU                         !" √テイラー数
      REAL       PRND                        !" プランドル数
*
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
*
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
      REAL       DELT                        !" 時間刻みΔt
c$$$      DATA       DELT  / 0.00001 /
      DATA       DELT  / 1 /
*
      INTEGER    ITRUNG                        !" Runge-Kutta 制御変数
      INTEGER    IORUNG                        !" Runge-Kutta 次数
      PARAMETER( IORUNG = 1 )
*
      CALL BNDRY
     M     ( WDTOR , WDPOR , WDT  , WDDPOR, 
     C       FLAPLA,
     C       ARAD  , DRAD                   )
*
*" << Runge-Kutta Loop >>
*
      ITRUNG = 0
*
 8888 CONTINUE                    !" <---  Runge-Kutta Loop : ここから
        ITRUNG = ITRUNG + 1 
         CALL DLOPER
     M           ( WDDPOR, 
     I             WDPOR , 
     F             'POS ',
     C             FLAPLA,
     C             ARAD  , DRAD    )
*
         CALL RESET( WTTOR, NMDIM*KMAX )
         CALL RESET( WTPOR, NMDIM*KMAX )
         CALL RESET( WTT  , NMDIM*KMAX )
*
         CALL LNDYN             !" 力学 線型項
     M           ( WTTOR , WTPOR , WTT   ,
     I             WDTOR , WDPOR , WDT   , WDDPOR ,  
     C             RAYL  , TAU   , PRND  , 
     C             NMO   , FLAPLA, 
     C             ARAD  , DRAD                    )
*
*" << INTGR : 時間積分 >>
*
         CALL TINTGR            !" スペクトル時間積分
     M         ( WDTOR , WDPOR , WDDPOR, WDT   , 
     I           WTTOR , WTPOR , WTT   , 
     I           DELT  , ITRUNG, 
     C           FLAPLA, 
     C           ARAD  , DRAD             )
*
      IF ( ITRUNG . LT. IORUNG ) GOTO 8888  !" <--- Runge-Kutta Loop 終わり
*
      RETURN
      END
***********************************************************************
      SUBROUTINE MKMLST
     O            ( NMKMTX , NMKMAX ,
     O              NMTLST , NMTMAX , 
     O              NMPLST , NMPMAX , 
     O              NMQLST , NMQMAX , 
     F              OATOR  , OAPOR  , OAT   , 
     I              NMO    , MFIX              )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZLDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zldim.F"                    !" 係数行列の大きさ
#endif
*
*"  [OUTPUT]
      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 の最大添え字
*
*"  [INPUT]
      LOGICAL    OATOR , OAPOR , OAT         !" 解析対象物理量
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      INTEGER    MFIX                        !" 設定東西波数
*
*"  [INTERNAL WORK]
*
      INTEGER    LTSTR                       !" トロイダル南北波数(最小)
      INTEGER    LTEND                       !" トロイダル南北波数(最大)
      INTEGER    LTDIV                       !" トロイダル南北波数(間隔)
*
      INTEGER    LPSTR                       !" ポロイダル南北波数(最小)
      INTEGER    LPEND                       !" ポロイダル南北波数(最大)
      INTEGER    LPDIV                       !" ポロイダル南北波数(間隔)
*
      INTEGER    LQSTR                       !" 温度南北波数(最小)
      INTEGER    LQEND                       !" 温度南北波数(最大)
      INTEGER    LQDIV                       !" 温度南北波数(間隔)
*
      INTEGER    L, IP, NM, NMK, K
*
      LOGICAL    OFIRST
      DATA       OFIRST  / .TRUE. /
*
*"  < 0. 計算する南北波数を設定 >
*
      IF( OFIRST )THEN
         OFIRST = .FALSE.
         CALL NWCOMP                   
     O         (  LTSTR , LTEND , LTDIV ,
     O            LPSTR , LPEND , LPDIV ,
     O            LQSTR , LQEND , LQDIV   )
      ENDIF
*
*"  << 1. 東西波数 MFIX のリスト作成 >>
*
*"       < 1.1  トロイダル成分 >
*
      NM = 0
      IF ( MFIX .EQ. 0 )THEN
         IP = 1
         DO 110 L = LTSTR, LTEND, LTDIV
            IF ( NMO( IP, MFIX, L ) .NE. 0 )THEN
               NM = NM + 1
               NMTLST( NM ) = NMO( IP, MFIX, L )
            ENDIF
 110     CONTINUE 
      ELSE
         DO 120 L = LTSTR, LTEND, LTDIV
            DO 120 IP = 1, 2
               IF ( NMO( IP, MFIX, L ) .NE. 0 )THEN
                  NM = NM + 1
                  NMTLST( NM ) = NMO( IP, MFIX, L )
               ENDIF
 120     CONTINUE 
      ENDIF
*
      NMTMAX = NM
*
*"       < 1.2  ポロイダル成分 >
*
      NM = 0
      IF ( MFIX .EQ. 0 )THEN
         IP = 1
         DO 130 L = LPSTR, LPEND, LPDIV
            IF ( NMO( IP, MFIX, L ) .NE. 0 )THEN
               NM = NM + 1
               NMPLST( NM ) = NMO( IP, MFIX, L )
            ENDIF
 130     CONTINUE 
      ELSE
         DO 140 L = LPSTR, LPEND, LPDIV
            DO 140 IP = 1, 2
               IF ( NMO( IP, MFIX, L ) .NE. 0 )THEN
                  NM = NM + 1
                  NMPLST( NM ) = NMO( IP, MFIX, L )
               ENDIF
 140     CONTINUE 
      ENDIF
*
      NMPMAX = NM
*
*"       < 1.3  温度成分 >
*
      NM = 0
      IF ( MFIX .EQ. 0 )THEN
         IP = 1
         DO 150 L = LQSTR, LQEND, LQDIV
            IF ( NMO( IP, MFIX, L ) .NE. 0 )THEN
               NM = NM + 1
               NMQLST( NM ) = NMO( IP, MFIX, L )
            ENDIF
 150     CONTINUE 
      ELSE
         DO 160 L = LQSTR, LQEND, LQDIV
            DO 160 IP = 1, 2
               IF ( NMO( IP, MFIX, L ) .NE. 0 )THEN
                  NM = NM + 1
                  NMQLST( NM ) = NMO( IP, MFIX, L )
               ENDIF
 160     CONTINUE 
      ENDIF
*
      NMQMAX = NM
*
*"  << 2. 係数行列の番号づけ >>
*
      NMK = 0
*
*"       < 2.1 トロイダルポテンシャル >
*
      IF( OATOR ) THEN
         IP = 1
         DO 210 K = 1, KDIM-1
            DO 210 NM = 1, NMTMAX
               NMK = NMK + 1
               NMKMTX( NM, K , IP ) = NMK
 210     CONTINUE 
      ENDIF
*
*"       < 2.2 ポロイダルポテンシャル >
*
      IF( OAPOR ) THEN
         IP = 2
         DO 220 K = 2, KDIM-2
            DO 220 NM = 1, NMPMAX
               NMK = NMK + 1
               NMKMTX( NM, K , IP ) = NMK
 220     CONTINUE
      ENDIF 
*
*"       < 2.3 温度 >
*
      IF( OAT ) THEN
         IP = 3
         DO 230 K = 1, KDIM-1
            DO 230 NM = 1, NMQMAX
               NMK = NMK + 1
               NMKMTX( NM, K , IP ) = NMK
 230     CONTINUE 
      ENDIF
*
      NMKMAX = NMK
*
*"  << 3. 行列の大きさチェック >>
*
      IF ( NMKMAX .GT. MATDMX ) THEN
         CALL MSGDMP( 'E', 'MKMLST', 
     &                 'DIMENSION OF MATRIX EXCEEDS THE LIMIT' )
      ENDIF
*
      RETURN
      END
**********************************************************************
      SUBROUTINE NWCOMP                   !" 南北波数範囲 の読み込み
     O            ( LTSTR , LTEND , LTDIV ,
     O              LPSTR , LPEND , LPDIV ,
     O              LQSTR , LQEND , LQDIV   )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*"  [OUTPUT]
      INTEGER    LPSTR                  !" ポロイダル南北波数(最小)
      INTEGER    LPEND                  !" ポロイダル南北波数(最大)
      INTEGER    LPDIV                  !" ポロイダル南北波数(間隔)
*
      INTEGER    LTSTR                  !" トロイダル南北波数(最小)
      INTEGER    LTEND                  !" トロイダル南北波数(最大)
      INTEGER    LTDIV                  !" トロイダル南北波数(間隔)
*
      INTEGER    LQSTR                       !" 温度南北波数(最小)
      INTEGER    LQEND                       !" 温度南北波数(最大)
      INTEGER    LQDIV                       !" 温度南北波数(間隔)
*
*"  [INTERNAL ONCE]
      LOGICAL    OLALL                  !" 南北波数計算スイッチ
      INTEGER    LTMIN , LTMAX , LTSTP  !" トロイダル南北波数(最大,最小,間隔)
      INTEGER    LPMIN , LPMAX , LPSTP  !" ポロイダル南北波数(最大,最小,間隔)
      INTEGER    LQMIN , LQMAX , LQSTP  !" 温度南北波数(最大,最小,間隔)
*
      NAMELIST   /NMLF/  OLALL, 
     &                   LTMIN  , LTMAX , LTSTP ,
     &                   LPMIN  , LPMAX , LPSTP ,
     &                   LQMIN  , LQMAX , LQSTP 
*
*   [INTERNAL WORK]
      INTEGER    IFPAR, JFPAR
*
*"     < 1. 南北波数範囲 >
*
      OLALL = .TRUE.
      LTMIN = 0
      LTMAX = 0 
      LTSTP = 0
      LPMIN = 0
      LPMAX = 0 
      LPSTP = 0
      LQMIN = 0
      LQMAX = 0 
      LQSTP = 0
*
C      WRITE  ( JFPAR, NMLF )
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMLF , END=1190 )
 1190 WRITE  ( JFPAR, NMLF )
*
      IF( OLALL )THEN        !" 全南北波数を計算
         LTSTR = 0
         LTEND = LMAX
         LTDIV = 1
         LPSTR = 0
         LPEND = LMAX
         LPDIV = 1
         LQSTR = 0
         LQEND = LMAX
         LQDIV = 1
         CALL MSGDMP( 'W', 'NWCOMP', 'CALCULATE ALL N-W WAVES' )
*
      ELSE IF(      ( LTSTP .LE. 0 ) 
     &         .OR. ( LPSTP .LE. 0 ) 
     &         .OR. ( LQSTP .LE. 0 )     )THEN
*
         CALL MSGDMP( 'E', 'NWCOMP', 'PARAMETER IS INVALID' )
*
      ELSE IF(      ( LTMIN .LT. 0 ) .OR. ( LTMAX .GT. LMAX ) 
     &         .OR. ( LPMIN .LT. 0 ) .OR. ( LPMAX .GT. LMAX )  
     &         .OR. ( LQMIN .LT. 0 ) .OR. ( LQMAX .GT. LMAX )  )THEN 
*
         CALL MSGDMP( 'E', 'NWCOMP', 'LMIN/MAX ARE BEYOND THE RANGE' )
*
      ELSE                               !" 南北波数を制限する
         LTSTR = LTMIN
         LTEND = LTMAX
         LTDIV = LTSTP
         LPSTR = LPMIN
         LPEND = LPMAX
         LPDIV = LPSTP
         LQSTR = LQMIN
         LQEND = LQMAX
         LQDIV = LQSTP
*
         CALL MSGDMP( 'W', 'NWCOMP', 'N-W WAVES ARE RESTRICTED' )
      ENDIF
*
      RETURN
      END
