* PACKAGE DLNDYN  !" 力学 線型項
*
*"  [HIS] 92/12/05(takepiro) トロイダルポロイダル interaction 微分表現
*"        93/02/15(takepiro)
*"        93/06/06(takepiro)
*
***********************************************************************
      SUBROUTINE LNDYN         !" 力学 線型項
     M         ( WTTOR , WTPOR , WTT   , 
     I           WDTOR , WDPOR , WDT   , WDDPOR,
     C           RAYL  , TAU   , PRND  , 
     C           NMO   , FLAPLA, 
     C           ARAD  , DRAD                   )
*
*"
*
#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       WDT   ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
      REAL       WDDPOR( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
*
      REAL       RAYL                          !" レイリー数
      REAL       TAU                           !" √テイラー数
      REAL       PRND                          !" プランドル数
*
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
*
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
      INTEGER    NM, K
*
*"         < 1. トロイダル Ψ >
*
      CALL GLOPER
     M           ( WTTOR , 
     I             WDTOR ,
     F             'ADD ',
     C             TAU   , PRND   , 
     C             NMO   , FLAPLA ,
     C             ARAD  , DRAD     )
*
      CALL QOPER 
     M           ( WTTOR , 
     I             WDPOR , 
     F             'SUB' ,
     C             TAU   , PRND   , 
     C             NMO   , FLAPLA , 
     C             ARAD  , DRAD     )
*
      DO 1000 K = 0, KDIM
         DO 1000 NM = 1, NMDIM
            IF( FLAPLA(NM) .EQ. 0. )THEN
               WTTOR( NM, K ) = 0.   
            ENDIF
 1000 CONTINUE 
*
*"         < 2. ポロイダル Φ >
*
      CALL GLOPER
     M           ( WTPOR , 
     I             WDDPOR, 
     F             'ADD ', 
     C             TAU   , PRND  ,
     C             NMO   , FLAPLA, 
     C             ARAD  , DRAD    )
*
      CALL QOPER 
     M           ( WTPOR , 
     I             WDTOR , 
     F             'ADD' ,
     C             TAU   , PRND   , 
     C             NMO   , FLAPLA ,
     C             ARAD  , DRAD            )
*
      DO 2000 K = 0, KDIM
         DO 2000 NM = 1, NMDIM
            IF( FLAPLA(NM) .NE. 0. )THEN
               WTPOR ( NM, K ) = WTPOR ( NM,K )
     &                         - RAYL * PRND * WDT ( NM,K )
            ELSE
               WTPOR( NM, K ) = 0.   
            ENDIF
 2000 CONTINUE 
*
*"         < 3. 温度場 Ｔ >
*
      CALL DLOPER
     M           ( WTT   , 
     I             WDT   , 
     F             'ADD ', 
     C             FLAPLA,
     C             ARAD  , DRAD    )
*
      DO 3000 K = 0, KDIM
         DO 3000 NM = 1, NMDIM
            WTT ( NM, K ) = WTT ( NM, K ) 
     &                   + ( -FLAPLA(NM) ) * WDPOR( NM, K )
 3000 CONTINUE 
*
      RETURN
      END
***********************************************************************
      SUBROUTINE GLOPER     !" 演算子 Pr * G_lm ( 力学:拡散+β )の計算 
     M           ( WODATA, 
     I             WDATA , 
     F             FUNC  , 
     C             TAU   , PRND  ,
     C             NMO   , FLAPLA,
     C             ARAD  , DRAD    )
*
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [MODIFY] 
      REAL       WODATA( NMDIM, 0:KDIM  )    !" 
*
*   [INPUT]
      REAL       WDATA ( NMDIM, 0:KDIM  )      !" 
      REAL       TAU                           !" √テイラー数
      REAL       PRND                          !" プランドル数
      CHARACTER  FUNC*4                        !" 動作スイッチ
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX )   !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )              !" ラプラシアンの係数
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
      REAL       WGLDAT( NMDIM, 0:KDIM  )      !" 
      INTEGER    M, L, LEND, K, NM
      REAL       AN
*
*"     < 1. ｉｍ をかける >
*
      CALL RESET ( WGLDAT , NMDIM*KMAX )
*
      DO 1000 M = 0 , MMAX, MINT
         LEND = MIN( LMAX, NMAX-M )
         DO 1000 L = 0 , LEND
               AN = L + M
               IF(  AN .NE. 0. )THEN
                  DO 1010 K = 0, KDIM
                     WGLDAT( NMO(1,M,L),K ) 
     &                    =  - REAL( M ) * TAU / ( AN*(AN+1) ) 
     &                          *  WDATA( NMO(2,M,L),K )
                     WGLDAT( NMO(2,M,L),K )
     &                    =    REAL( M ) * TAU / ( AN*(AN+1) ) 
     &                          * WDATA( NMO(1,M,L),K )
 1010             CONTINUE 
               ENDIF
 1000 CONTINUE 
*
*"     < 2. D_l をかける >
*
      CALL DLOPER
     M           ( WGLDAT, 
     I             WDATA , 
     F             'ADD ',
     C             FLAPLA,
     C             ARAD  , DRAD    )
*
*"     < 3. 出力 >
*
      IF ( FUNC(1:1) .EQ. 'A' )THEN
         DO 2100 K = 0, KDIM
            DO 2100 NM = 1, NMDIM
               WODATA( NM,K ) = WODATA( NM,K ) + PRND * WGLDAT( NM,K )
 2100    CONTINUE 
*
      ELSE IF ( FUNC(1:1) .EQ. 'S' )THEN         
         DO 2200 K = 0, KDIM
            DO 2200 NM = 1, NMDIM
               WODATA( NM,K ) = WODATA( NM,K ) - PRND * WGLDAT( NM,K )
 2200    CONTINUE 
*
      ELSE IF ( FUNC(1:1) .EQ. 'N' )THEN         
         DO 2300 K = 0, KDIM
            DO 2300 NM = 1, NMDIM
               WODATA( NM,K ) = - PRND * WGLDAT( NM,K )
 2300    CONTINUE 
*
      ELSE 
         DO 2400 K = 0, KDIM
            DO 2400 NM = 1, NMDIM
               WODATA( NM,K ) = PRND * WGLDAT( NM,K )
 2400    CONTINUE 
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE DLOPER     !" 演算子 D_l の計算
     M           ( WODATA, 
     I             WDATA , 
     F             FUNC  ,
     C             FLAPLA,
     C             ARAD  , DRAD    )
*
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [MODIFY] 
      REAL       WODATA( NMDIM, 0:KDIM  )      !" 
*
*   [INPUT]
      REAL       WDATA ( NMDIM, 0:KDIM  )      !" 
      CHARACTER  FUNC*4                        !" 動作スイッチ
      REAL       FLAPLA( NMDIM  )              !" ラプラシアンの係数 -n(n+1)
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/ WDLDAT, WDAT1, WDAT2, CDRDAT
      REAL       WDLDAT( NMDIM, 0:KDIM  )      !"  
      REAL       WDAT2 ( NMDIM, 0:KDIM  )      !" 
      REAL       WDAT1 ( NMDIM, 0:KDIM  )      !" 
      REAL       CDRDAT( NMDIM, 0:KDIM  )      !" 
*
      INTEGER    NM, K
*
*"  < 1. 演算子 D_l の計算 >
*
      CALL W2C 
     M         ( CDRDAT,
     I           WDATA ,
     I           'DR'  , 'POS' , KMAX )
*
      CALL C2W
     O         ( WDAT2 ,
     I           CDRDAT,
     I           'DR'  , 'POS' , KMAX  )
*
      CALL C2W
     O         ( WDAT1 ,
     I           CDRDAT,
     I           ' '  , 'POS' , KMAX  )
*
      DO 1000 K = 0, KDIM
         DO 1000 NM = 1, NMDIM
            WDLDAT( NM,K ) = WDAT2( NM, K ) 
     &                 + 2.0/ ARAD( K ) * WDAT1( NM, K ) 
     &                 - ( -FLAPLA(NM) ) / ARAD(K)**2 * WDATA( NM,K )
 1000 CONTINUE 
*
*"  < 2. 出力 >
*
      IF ( FUNC(1:1) .EQ. 'A' )THEN
         DO 2100 K = 0, KDIM
            DO 2100 NM = 1, NMDIM
               WODATA( NM,K ) = WODATA( NM,K ) + WDLDAT( NM,K )
 2100    CONTINUE 
*
      ELSE IF ( FUNC(1:1) .EQ. 'S' )THEN         
         DO 2200 K = 0, KDIM
            DO 2200 NM = 1, NMDIM
               WODATA( NM,K ) = WODATA( NM,K ) - WDLDAT( NM,K )
 2200    CONTINUE 
*
      ELSE IF ( FUNC(1:1) .EQ. 'N' )THEN         
         DO 2300 K = 0, KDIM
            DO 2300 NM = 1, NMDIM
               WODATA( NM,K ) = - WDLDAT( NM,K )
 2300    CONTINUE 
*
      ELSE 
         DO 2400 K = 0, KDIM
            DO 2400 NM = 1, NMDIM
               WODATA( NM,K ) = WDLDAT( NM,K )
 2400    CONTINUE 
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE QOPER     !" 演算子 Pr * τ * Ｑ/n(n+1) の計算
     M           ( WODATA, 
     I             WDATA , 
     F             FUNC  ,
     C             TAU   , PRND   ,
     C             NMO   , FLAPLA ,
     C             ARAD  , DRAD           )
*
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [MODIFY] 
      REAL       WODATA( NMDIM, 0:KDIM  )      !" 
*
*   [INPUT]
      REAL       WDATA ( NMDIM, 0:KDIM  )      !" 
*
      CHARACTER  FUNC*4                        !" 動作スイッチ
*
      REAL       TAU                           !" √テイラー数
      REAL       PRND                          !" プランドル数
*
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX )   !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )              !" ラプラシアンの係数
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
c$$$      COMMON    /COMWRK/ WKGRD, WKGL2, WORK
      REAL       WKGRD ( NMDIM, 0:KDIM  )      !" 
      REAL       WKGL2 ( NMDIM, 0:KDIM  )      !" 
      REAL       WORK  ( NMDIM, 0:KDIM  )      !" 
*
      INTEGER    NM, K
*
*"    < 1. ｋ・▽ の計算 > 
*
      CALL KGRAD 
     O     ( WKGRD , 
     I       WDATA , 
     C       NMO   ,
     C       ARAD  , DRAD    )
*
*"    < 2. ｋ・▽ L_2 の計算 > 
*
      DO 2000 K = 0, KDIM
         DO 2000 NM = 1, NMDIM
            WORK( NM,K ) = - FLAPLA( NM ) * WDATA( NM,K )
 2000 CONTINUE 
*
      CALL KGRAD 
     O     ( WKGL2 , 
     I       WORK  , 
     C       NMO   ,
     C       ARAD  , DRAD    )
*
*"    < 3. Ｑ の計算 > 
*
      DO 3000 K = 0, KDIM
         DO 3000 NM = 1, NMDIM
            IF( FLAPLA( NM ) .NE. 0 ) THEN
               WORK( NM,K ) = TAU / ( -FLAPLA( NM ) )
     &                   * ( WKGRD( NM,K ) 
     &                       - 0.5 * ( - FLAPLA( NM )* WKGRD( NM,K ) 
     &                                 + WKGL2( NM,K ) ) )
            ELSE
               WORK( NM,K ) = 0.
            ENDIF
 3000 CONTINUE 
*
*"    < 4. 出力 >
*
      IF ( FUNC(1:1) .EQ. 'A' )THEN
         DO 3100 K = 0, KDIM
            DO 3100 NM = 1, NMDIM
               WODATA( NM,K ) = WODATA( NM,K ) + PRND * WORK( NM,K )
 3100    CONTINUE 
*
      ELSE IF ( FUNC(1:1) .EQ. 'S' )THEN         
         DO 3200 K = 0, KDIM
            DO 3200 NM = 1, NMDIM
               WODATA( NM,K ) = WODATA( NM,K ) - PRND * WORK( NM,K )
 3200    CONTINUE 
*
      ELSE IF ( FUNC(1:1) .EQ. 'N' )THEN         
         DO 3300 K = 0, KDIM
            DO 3300 NM = 1, NMDIM
               WODATA( NM,K ) = - PRND * WORK( NM,K )
 3300    CONTINUE 
*
      ELSE 
         DO 3400 K = 0, KDIM
            DO 3400 NM = 1, NMDIM
               WODATA( NM,K ) = PRND * WORK( NM,K )
 3400    CONTINUE 
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE KGRAD     !" 演算子 ｋ・▽ の計算
     O           ( WODATA, 
     I             WDATA , 
     C             NMO   ,
     C             ARAD  , DRAD    )
*
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [MODIFY] 
      REAL       WODATA( NMDIM, 0:KDIM  )      !" 
*
*   [INPUT]
      REAL       WDATA ( NMDIM, 0:KDIM  )      !" 
*
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX )   !" スペクトルの添字順番
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/ WDAT1, WDAT2, WORK, GWORK
      REAL       WDAT2 ( NMDIM, 0:KDIM  )      !" 
      REAL       WDAT1 ( NMDIM, 0:KDIM  )      !" 
      REAL       WORK  ( NMDIM, 0:KDIM  )      !" 
      REAL       GWORK ( IDIM*JDIM, 0:KDIM  )  !" 
*
      INTEGER    NM, K
*
*"    < 1. μ d/dr の計算 > 
*
      CALL MLTIMU
     O     ( WDAT1 , 
     I       WDATA , 
     C       NMO     )
*
      CALL W2C 
     M         ( WORK  ,
     I           WDAT1 ,
     I           'DR'  , 'POS' , KMAX )
*
      CALL C2W
     O         ( WDAT1 ,
     I           WORK  ,
     I           ' '  , 'POS' , KMAX  )
*
*"    < 2. (1-μ^2)/r d/dμ の計算 > 
*
      CALL W2G
     O         ( GWORK  ,
     I           WDATA ,
     I           'YGRA', 'POS ', KMAX )
*
      CALL G2W
     O         ( WDAT2 ,
     I           GWORK  ,
     I           ' ', 'POS ', KMAX )
*
*"    < 3. ｋ・▽ の計算 > 
*
      DO 1000 K = 0, KDIM
         DO 1000 NM = 1, NMDIM
            WODATA( NM,K ) = WDAT1( NM,K ) + WDAT2( NM,K ) / ARAD( K )
 1000 CONTINUE 
*
      RETURN
      END
***********************************************************************
      SUBROUTINE MLTIMU     !" μW_n^m Y_n^m の計算
     O           ( WODATA, 
     I             WDATA , 
     C             NMO     ) 
*
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [MODIFY] 
      REAL       WODATA( NMDIM, 0:KDIM  )      !" 
*
*   [INPUT]
      REAL       WDATA ( NMDIM, 0:KDIM  )      !" 
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX )   !" スペクトルの添字順番
*
*   [INTERNAL WORK] 
      INTEGER    NMOL   ( NMDIM )    !" 計算するスペクトルのリスト
      INTEGER    NMLSTL ( NMDIM )    !" 
      INTEGER    NMOU   ( NMDIM )    !" 
      INTEGER    NMLSTU ( NMDIM )    !" 
*
      INTEGER    NML, NMU
*
      REAL       COEFLO( NMDIM )               !" 演算子の係数
      REAL       COEFUP( NMDIM )               !" 演算子の係数
*
      SAVE       NMOL, NMOU, NMLSTL, NMLSTU, NML, NMU, COEFLO, COEFUP
*
      INTEGER    NM, L, M, N, LEND, K, I, IE
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
*
*"     < 1. 係数, リスト作成 >
*
      IF ( OFIRST )THEN
         OFIRST = .FALSE.
         NML = 0
         NMU = 0
*
         IF ( MMAX .EQ. 0 )THEN     !" 東西平均 or not
            IE = 1
         ELSE
            IE = 2
         ENDIF
*
         DO 1000 M = 0 , MMAX, MINT
            LEND = MIN( LMAX, NMAX-M )
            DO 1000 L = 0 , LEND
               DO 1000 I = 1, IE
                  N = L + M
*
                  IF( L-1 .GE. 0 )THEN
                     NML = NML + 1
                     NMOL  ( NML ) = NMO( I,M,L )
                     NMLSTL( NML ) = NMO( I,M,L-1 )
                     COEFLO( NML ) = SQRT( REAL( L*(N+M) ) 
     &                                     / REAL( (2*N-1)*(2*N+1) ) )
                  ENDIF
*
                  IF( L+1 .LE. LEND )THEN
                     NMU = NMU + 1
                     NMOU  ( NMU ) = NMO( I,M,L )
                     NMLSTU( NMU ) = NMO( I,M,L+1 )
                     COEFUP( NMU ) = SQRT( REAL( (L+1)*(N+M+1) ) 
     &                                     / REAL( (2*N+1)*(2*N+3) ) )
                  ENDIF
 1000    CONTINUE 
      ENDIF
*
*"      < 2. スペクトル係数の計算 >
*
      CALL RESET( WODATA, NMDIM*KMAX )
*
      DO 2000 K = 0, KDIM
         DO 2000 NM = 1, NML
            WODATA( NMOL(NM), K ) = WODATA( NMOL(NM), K ) 
     &                  + COEFLO( NM ) * WDATA( NMLSTL( NM ), K )
 2000 CONTINUE 
*
      DO 2100 K = 0, KDIM
         DO 2100 NM = 1, NMU
            WODATA( NMOU(NM), K ) = WODATA( NMOU(NM), K ) 
     &                  + COEFUP( NM ) * WDATA( NMLSTU( NM ), K )
 2100 CONTINUE 
*
      RETURN
      END
