* PACKAGE DBNDR  !" 力学 境界条件(チェビシェフ空間での計算)
*
*"  [HIS] 92/11/19(takepiro)
*
***********************************************************************
      SUBROUTINE BNDRY
     M       ( WDTOR , WDPOR , WDT  , WDDPOR, 
     C         FLAPLA,
     C         ARAD  , DRAD                   )
*
*"        free-slep , 温度固定境界条件
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [MODIFY] 
      REAL       WDTOR ( NMDIM , 0:KDIM     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , 0:KDIM     )  !" ポロイダル Φ
      REAL       WDDPOR( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
      REAL       WDT   ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
*
*   [INPUT] 
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
*
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
      INTEGER    NM
*
*"     < 1. トロイダルポテンシャル  d(Ψ/r)/dr = 0 >
*
      DO 1100 NM = 1 , NMDIM
         WDTOR( NM,  0  ) = ARAD(0) 
     &                       * WDTOR( NM,1 ) / ARAD(1)
         WDTOR( NM,KDIM ) = ARAD(KDIM) 
     &                       * WDTOR( NM,KDIM-1 ) / ARAD(KDIM-1)
 1100 CONTINUE 
*
*"     < 2. ポロイダルポテンシャル  Φ = 0,  d^2Φ/dr^2 = 0 >
*
      CALL DLINVS
     O     (   WDPOR , 
     I         WDDPOR, 
     C         FLAPLA,
     C         ARAD  , DRAD  )
*
*"     < 3. 温度  T = 0 >
*
      DO 1300 NM = 1 , NMDIM
         WDT( NM,  0  ) = 0.
         WDT( NM,KDIM ) = 0.
 1300 CONTINUE 
*
      RETURN
      END
************************************************************************
      SUBROUTINE DLINVS     !" チェビシェフ空間での逆ラプラシアン計算
     O             (   WDPOR , 
     I                 WDDPOR, 
     C                 FLAPLA,
     C                 ARAD  , DRAD  )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       WDPOR ( NMDIM , KMAX  )  !" ポロイダル Φ
*
*   [INPUT] 
      REAL       WDDPOR( NMDIM , KMAX  )  !" ポロイダル D_l Φ
*
      REAL       FLAPLA( NMDIM  )         !" ラプラシアンの係数
*
      REAL       ARAD  ( KMAX )           !" ｒレベル(整数)
      REAL       DRAD  ( KMAX )           !" Δｒ(整数)
*
*   [INTERNAL WORK] 
      REAL       DLMTRX( NMDIM, KMAX, KMAX )   !" D_l の行列
*                                              !" DL_lko WP_lo = WDP_lk
*                                              !" 各 NM 最後 4 行は境界条件
*                                              !" 対応する WDP_lk を 0 に設定
      INTEGER    KP  ( NMDIM, KMAX )           !" ピボット
*
      REAL       WORK  ( NMDIM, KMAX )         !" 作業用行列
      INTEGER    NM, K
      LOGICAL    OFIRST
      DATA       OFIRST  / .TRUE. /
*
      SAVE
*
*"      << 行列設定 >>
*
*   [ONCE]
      IF ( OFIRST )THEN
         OFIRST = .FALSE.
*
         CALL STDLMT                   !" D_l 演算子の行列設定
     O        (   DLMTRX,
     C            FLAPLA,
     C            ARAD  , DRAD  )
*
         CALL  LUMAKE                  !" 行列のＬＵ分解
     M         ( DLMTRX,
     O           KP    ,
     D           NMDIM , KMAX )
      ENDIF
*
*"   < 1. D_l の逆演算 >
*
      CALL W2C
     M         ( WORK  ,
     I           WDDPOR,
     I           ' '   , 'POS ' , KMAX )
*
*"  < 境界条件に対応する行を 0 に設定 >
*
      DO 1000 K = KMAX-4, KMAX
         DO 1000 NM = 1, NMDIM
            WORK( NM, K ) = 0.0
 1000 CONTINUE 
*
*"  < ＬＵ分解による解の計算 >
*
      CALL LUSOLV
     M         ( WORK  ,
     I           DLMTRX, KP    , 
     D              1  , NMDIM , KMAX )
*
      CALL C2W
     M         ( WDPOR ,
     I           WORK  ,
     I           ' '   , 'POS ' , KMAX )
*
      RETURN
      END
************************************************************************
      SUBROUTINE STDLMT    !" D_l 演算子の行列設定(チェビシェフ空間)
     O             (   DLMTRX,
     C                 FLAPLA,
     C                 ARAD  , DRAD  )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       DLMTRX( NMDIM, KMAX, KMAX ) !" D_l の行列
*"                                           !" DL_lko WP_lo = WDP_lk
*   [INPUT] 
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
*
      REAL       ARAD  ( KMAX )              !" ｒレベル(整数)
      REAL       DRAD  ( KMAX )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
*
      REAL       DR1   ( KMAX , KMAX )       !" 1 階微分の表現
      REAL       DR2   ( KMAX , KMAX )       !" 2 階微分の表現
      REAL       RM1DR1( KMAX , KMAX )       !" 1/r d/dr の表現
      REAL       RMIN2 ( KMAX , KMAX )       !" 1/r^2 の表現
*
      REAL       BNDRYC ( 4, KMAX )          !" 境界条件の表現
*
      REAL       CDATA  ( NMDIM, KMAX )      !" チェビシェフ作業領域
      REAL       CDDATA ( NMDIM, KMAX )      !" チェビシェフ作業領域
      REAL       CD2DAT ( NMDIM, KMAX )      !" チェビシェフ作業領域
      REAL       CRDDAT ( NMDIM, KMAX )      !" チェビシェフ作業領域
      REAL       CRDATA ( NMDIM, KMAX )      !" チェビシェフ作業領域
      REAL       WDATA  ( NMDIM, KMAX )      !" スペクトル作業領域
*
      REAL       CHBTOP ( KMAX )             !" 上側境界
      REAL       CHBBTM ( KMAX )             !" 下側境界
*
      INTEGER    NM, K, KI, KJ
*
*" << 1. 微分 etc.の行列表現  >>
*
      DO 1000 K = 1, KMAX
         CALL RESET ( CDATA, NMDIM*KMAX )
         CDATA( 1, K ) = 1.
*
*"     < 1.1  1 階微分 >
*
         CALL CDERIV
     O         ( CDDATA,
     I           CDATA ,
     I           KMAX    )
*
         DO 1100 KI = 1, KMAX
            DR1( KI, K ) = CDDATA( 1, KI )
 1100    CONTINUE 
*
*"     < 1.2  2 階微分 >
*
         CALL CDERIV
     O         ( CD2DAT ,
     I           CDDATA ,
     I           KMAX     )
*
         DO 1200 KI = 1, KMAX
            DR2( KI, K ) = CD2DAT( 1, KI )
 1200    CONTINUE 
*
*"     < 1.3  1/r d/dr 微分 >
*
            CALL C2W     
     M         ( WDATA ,
     I           CDDATA,
     I           ' '  , 'POS ' , KMAX )
*
            DO 1300 KI = 1, KMAX
                  WDATA ( 1,KI ) = WDATA( 1,KI ) / ARAD( KI )
 1300       CONTINUE 
*
            CALL W2C     
     M         ( CRDDAT,
     I           WDATA ,
     I           ' '  , 'POS ' , KMAX )
*
            DO 1400 KI = 1, KMAX
               RM1DR1( KI, K ) = CRDDAT( 1, KI )
 1400       CONTINUE 
*
*"     < 1.4  1/r^2 のかけ算 >
*
            CALL C2W     
     M         ( WDATA ,
     I           CDATA,
     I           ' '  , 'POS ' , KMAX )
*
            DO 1500 KI = 1, KMAX
                  WDATA ( 1,KI ) = WDATA( 1,KI ) / ( ARAD( KI )**2 )
 1500       CONTINUE 
*
            CALL W2C     
     M         ( CRDATA,
     I           WDATA ,
     I           ' '  , 'POS ' , KMAX )
*
            DO 1600 KI = 1, KMAX
               RMIN2( KI, K ) = CRDATA( 1, KI )
 1600       CONTINUE 
*
 1000 CONTINUE 
*
*"  << 2. 境界条件の表現 >>
*
      CALL RESET( BNDRYC, 4*KMAX )
*
      CHBTOP( 1 ) = 0.5 
      CHBBTM( 1 ) = 0.5 
*
      DO 2100 K = 2, KMAX
         CHBTOP( K ) = 1.
         CHBBTM( K ) = REAL( (-1)**(K-1) )
 2100 CONTINUE 
*
*"      < 2.1 境界条件 Φ = 0 >
*
      DO 2200 KJ = 1, KMAX
         BNDRYC( 1,KJ ) = CHBTOP( KJ )  
         BNDRYC( 2,KJ ) = CHBBTM( KJ )  
 2200 CONTINUE 
*
*"      < 2.2 境界条件 dr^2 d^2 Φ = 0 >
*
      DO 2300 KJ = 1, KMAX
         DO 2300 KI = 1, KMAX
            BNDRYC( 3,KJ ) = BNDRYC( 3,KJ ) 
c$$$     &                        + CHBTOP( KI ) * DR2( KI,KJ )
     &                        + CHBTOP( KI ) * DR1( KI,KJ )
            BNDRYC( 4,KJ ) = BNDRYC( 4,KJ ) 
c$$$     &                        + CHBBTM( KI ) * DR2( KI,KJ )
     &                        + CHBBTM( KI ) * DR1( KI,KJ )
 2300 CONTINUE 
*
*"  << 3. 行列の設定 >>
*
      CALL RESET ( DLMTRX, NMDIM*KMAX*KMAX )
*
*"    < 3.1. D_l 演算子の行列表現設定 d/dr^2 + 2/r d/dr - l(l+1)/r^2 >
*
      DO 3100 KJ = 1, KMAX
         DO 3100 KI = 1, KMAX
            DO 3100 NM = 1, NMDIM
               DLMTRX( NM,KI,KJ ) 
     &              =    DR2( KI,KJ ) + 2.* RM1DR1( KI, KJ ) 
     &                 + FLAPLA( NM ) * RMIN2( KI, KJ )
 3100 CONTINUE 
*
*"    < 3.2 境界条件 >
*
      DO 3200 KJ = 1, KMAX
         DO 3200 KI = 1, 4
            DO 3200 NM = 1, NMDIM
               DLMTRX( NM, KMAX-4+KI, KJ ) = BNDRYC( KI,KJ )
 3200 CONTINUE 
*
      RETURN
      END
