* PACKAGE LSTEDY  !" 線形力学定常解(非線形項に対する応答)
*"                  ver 1.02                       93/09/27 takepiro
**********************************************************************
      SUBROUTINE LSTEDY                      !"  線形力学定常解
     M             ( GAU   , GAV   , GAW   , GAT   ,
     M               GATOR , GAPOR , 
     I               MFIX  ,
     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
*
*   [MODIFY]
      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>  格子点データ
*"           入力-非線形項を計算, 出力-非線形項に対する応答
*
*   [INPUT]
      INTEGER    MFIX                        !" 定常解の東西波数
*
      REAL       ALON  ( IDIM )              !" 経度
      REAL       DLON  ( IDIM )              !" 経度荷重
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
      REAL       ARAD  ( 0:KDIM )            !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )            !" Δｒ(整数)
*"        : 座標格子
*
*   [INTERNAL ONCE]
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       EDEL  ( NMDIM  )            !" ζ，Ｄ→Ｕ，Ｖ
*
      REAL       UVFACT ( IDIM, JDIM )       !" u→U のファクター
*
      REAL       RAYL                        !" レイリー数
      REAL       TAU                         !" √テイラー数
      REAL       PRND                        !" プランドル数
*
      SAVE       NMO   , FLAPLA, EDEL  ,
     &           UVFACT
*
      LOGICAL    OSETC , OFIRST
      DATA       OSETC , OFIRST  / .FALSE., .TRUE.  /
      SAVE       OSETC , OFIRST
*
      LOGICAL    ONONLN                      !" 非線形力学スイッチ
*
      LOGICAL    OATOR , OAPOR , OAT         !" 解析対象物理量
      SAVE       OATOR , OAPOR , OAT
*
*   [INTERNAL WORK]
      REAL       GTUA  ( IDIM, JDIM, 0:KDIM )  !" 東西運動量変化項ＵＡ
      REAL       GTVA  ( IDIM, JDIM, 0:KDIM )  !" 南北運動量変化項ＶＡ
      REAL       GTWA  ( IDIM, JDIM, 0:KDIM )  !" 鉛直運動量変化項ＶＡ
      REAL       GTH   ( IDIM, JDIM, 0:KDIM )  !" 温度時間変化項  Ｈ
*"       : 格子点変化項を足し込む
*
      REAL       GTUT  ( IDIM, JDIM, 0:KDIM )  !" 温度東西移流項  ＵＴ
      REAL       GTVT  ( IDIM, JDIM, 0:KDIM )  !" 温度南北移流項  ＶＴ
      REAL       GTWT  ( IDIM, JDIM, 0:KDIM )  !" 温度鉛直移流項  ＷＴ
*"       : 時間変化項
*
      REAL       WTTOR ( NMDIM, 0:KDIM      )  !" トロイダル Ψ
      REAL       WTPOR ( NMDIM, 0:KDIM      )  !" ポロイダル Φ
      REAL       WTT   ( NMDIM, 0:KDIM      )  !" 温度  Ｔ
*"       : スペクトル時間変化項
*
      REAL       COFMTX( MATDMX , MATDMX )   !" 係数行列
      REAL       VSTEDY( MATDMX )            !" 非線形項 / 定常解
*
      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    IERR                        !" サブルーチンエラーコード
      REAL       EPS                         !" 特異性の判定基準
      REAL       EPS0                        !" 特異性の判定基準
      DATA       EPS0 / 3.52D-15 /
*"       : 固有値計算用行列
*
      INTEGER    IUNIT                       !" 標準出力装置番号
      INTEGER    I , J
*
*"  << 0.SETCON : 定数設定 >>
*
*   [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    DPARM                     !" パラメター設定
     O             ( RAYL   , TAU  , PRND  , 
     O               ONONLN                   )
*
         IF( ONONLN )THEN
            CALL MSGDMP('W','LSTDY',
     &                  'N-LIN. DYNAMICS SWITCH WAS NEGLECTED')
         ENDIF
*
         CALL HSPARM
     I          ( 'RAYLEIGH NUMBER', RAYL   )
         CALL HSPARM
     I          ( 'TAYLOR NUMBER'  , TAU**2 )
         CALL HSPARM
     I          ( 'PRANDTL NUMBER' , PRND   )
*
c$$$         CALL DHSTRG                    !" 出力設定
*
         CALL LSTPHS                         !" 解析対象物理量の設定
     O         ( OATOR  , OAPOR  , OAT  )
      ENDIF
*
*" << 1. 係数行列の計算 : 線型力学 >>
*
      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' )
*
*" << 2. 非線形項の計算 >>
*
      CALL CLCSTR( 'W2VECT' )
*
      CALL GRDDYN               !"  非線型力学項
     O         ( GTUA  , GTVA  , GTWA  , 
     O           GTUT  , GTVT  , GTWT  , GTH   ,
     I           GAU   , GAV   , GAW   , GAT   ,
     C           UVFACT, FLAPLA, 
     C           ALAT  , DLAT  , ARAD  , DRAD    )
*
      CALL GRDFRC               !" 力学外部強制項(格子点)
     M         ( GTUA  , GTVA  , GTWA  , 
     M           GTUT  , GTVT  , GTWT  , GTH   ,
     C           ALON  , DLON  , ALAT  , DLAT  ,
     C           ARAD  , DRAD                    )
*
      CALL TENG2W               !" 時間変化項をスペクトルに
     O         ( WTTOR , WTPOR , WTT   , 
     I           GTUA  , GTVA  , GTWA  , 
     I           GTUT  , GTVT  , GTWT  , GTH , 
     C           FLAPLA,
     C           ARAD  , DRAD                  )
*
      CALL W2VECT               !" 非線形項をベクトルに
     O         ( VSTEDY , 
     I           WTTOR  , WTPOR , WTT   , 
     I           NMKMTX , NMKMAX , 
     I           NMTLST , NMTMAX , 
     I           NMPLST , NMPMAX , 
     I           NMQLST , NMQMAX , 
     F           OATOR  , OAPOR  , OAT    )
*
      DO 2000 I = 1, NMKMAX
         VSTEDY( I ) = - VSTEDY( I )
 2000 CONTINUE 
*
      CALL CLCSTR( 'W2VECT' )
*
*" << 3. 定常解計算 >>
*
      CALL CLCSTR( 'LGLU' )
*
      EPS = EPS0 * PRND         !" 特異性の基準値を設定
*
      CALL LEIGLU               !" 連立 1 次方程式を解く
c$$$      CALL LSVDSL               !" 特異連立 1 次方程式を解く
     M       ( COFMTX ,
     M         VSTEDY ,
     O         IERR   , 
     I         EPS    ,
     I         NMKMAX , .TRUE.   )
*
      IF ( IERR .EQ. 0 )THEN
         CALL LEIVEC
     O       ( GAU    , GAV    , GAW    , GAT    ,
     O         GATOR  , GAPOR  , 
     I         VSTEDY , VSTEDY ,
     I         NMKMTX , 
     I         NMTLST , NMTMAX , 
     I         NMPLST , NMPMAX , 
     I         NMQLST , NMQMAX , 
     F         OATOR  , OAPOR  , OAT    ,
     C         NMO    , FLAPLA , UVFACT , 
     C         ARAD   , DRAD                       )
      ENDIF
      CALL CLCSTR( 'LGLU' )
*
*
*" << OUTPUT : 出力 >>
*
      CALL HISTIN ( GTUA  , 'UA' )  !" 出力領域へ
      CALL HISTIN ( GTVA  , 'VA' )
      CALL HISTIN ( GTWA  , 'WA' )
      CALL HISTIN ( GTUT  , 'UT' )  !" 出力領域へ
      CALL HISTIN ( GTVT  , 'VT' )
      CALL HISTIN ( GTWT  , 'WT' )
      CALL HISTIN ( GTH   , 'H'  )
*
      RETURN
      END
**********************************************************************
      SUBROUTINE  W2VECT               !" 非線形項をベクトルに
     O                ( EIGVEC , 
     I                  WDTOR  , WDDPOR , WDT   , 
     I                  NMKMTX , NMKMAX , 
     I                  NMTLST , NMTMAX , 
     I                  NMPLST , NMPMAX , 
     I                  NMQLST , NMQMAX , 
     F                  OATOR  , OAPOR  , OAT    )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZLDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zldim.F"                    !" 係数行列の大きさ
#endif
*
*"  [OUTPUT]
      REAL       EIGVEC( MATDMX )            !" 固有ベクトル
*
*"  [INPUT]
      REAL       WDTOR ( NMDIM , 0:KDIM     )  !" トロイダル Ψ
      REAL       WDDPOR( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
      REAL       WDT   ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
*
      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 の最大添え字
*
      LOGICAL    OATOR , OAPOR , OAT         !" 解析対象物理量
*
      INTEGER    K, NML
*
*"     < 1. 固有ベクトルの計算(スペクトル) >
*
      IF ( OATOR ) THEN
         DO 1100 K = 1, KDIM-1
            DO 1100 NML = 1, NMTMAX
               EIGVEC( NMKMTX(NML,K,1) ) 
     &               = WDTOR( NMTLST(NML), K ) 
 1100    CONTINUE 
      ENDIF
*
      IF ( OAPOR )THEN
         DO 1200 K = 2, KDIM-2
            DO 1200 NML = 1, NMPMAX
               EIGVEC( NMKMTX(NML,K,2) ) 
     &               = WDDPOR( NMPLST(NML), K ) 
 1200    CONTINUE 
      ENDIF
*
      IF ( OAT ) THEN
         DO 1300 K = 1, KDIM-1
            DO 1300 NML = 1, NMQMAX
               EIGVEC( NMKMTX(NML,K,3) )
     &               = WDT( NMQLST(NML), K ) 
 1300    CONTINUE 
      ENDIF
*
      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
