* PACKAGE DBNDR  !" 力学 境界条件(チェビシェフ空間での計算)
*
*"  [HIS] 92/11/19(takepiro)
*"        92/12/22(takepiro) トロイダルをチェビシェフ空間で計算
*"        92/12/30(takepiro) rigid 条件追加
*
***********************************************************************
      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
      INTEGER    IFPAR, JFPAR
*
      LOGICAL    RIGIDB                      !" 下端力学境界条件スイッチ
      LOGICAL    RIGIDT                      !" 上端力学境界条件スイッチ
      LOGICAL    TFIXB                       !" 下端熱境界条件スイッチ
      LOGICAL    TFIXT                       !" 上端熱境界条件スイッチ
      DATA       RIGIDB, RIGIDT / .FALSE., .FALSE./
      DATA       TFIXB , TFIXT  / .TRUE. , .TRUE. /
*
      NAMELIST   / NMBNDR / RIGIDB, RIGIDT, TFIXB, TFIXT
*
      LOGICAL    ORIGDT                      !" 上端力学境界条件スイッチ
      LOGICAL    ORIGDB                      !" 下端力学境界条件スイッチ
      LOGICAL    OTFIXT                      !" 上端熱境界条件スイッチ
      LOGICAL    OTFIXB                      !" 下端熱境界条件スイッチ
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
*
      SAVE
*
      IF( OFIRST )THEN
         OFIRST = .FALSE.
         WRITE ( 6, * ) ' BOUNDARY CONDITIONS : DATE=92/12/22' 
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMBNDR, END=1190 )
 1190    WRITE  ( JFPAR, NMBNDR  )
*
         ORIGDB = RIGIDB
         ORIGDT = RIGIDT
         OTFIXB = TFIXB
         OTFIXT = TFIXT
*
         CALL SETMTX
      ENDIF
*
*"     < 1. トロイダルポテンシャル  d(Ψ/r)/dr = 0 >
*
      CALL BNDTOR
     M     (   WDTOR , 
     F         ORIGDB, ORIGDT, 
     C         FLAPLA,
     C         ARAD  , DRAD  )
*
*"     < 2. ポロイダルポテンシャル  Φ = 0,  d^2Φ/dr^2 = 0 >
*
      CALL DLINVS
     O     (   WDPOR , 
     I         WDDPOR, 
     F         ORIGDB, ORIGDT,
     C         FLAPLA,
     C         ARAD  , DRAD  )
*
*"     < 3. 温度  T = 0 >
*
      DO 1300 NM = 1 , NMDIM
         IF( OTFIXT )THEN
            WDT( NM,  0  ) = 0.
         ELSE
            WDT( NM,  0  ) = WDT( NM, 1 )
         ENDIF
         IF( OTFIXB )THEN
            WDT( NM,KDIM ) = 0.
         ELSE
            WDT( NM,KDIM ) = WDT( NM, KDIM-1 )
         ENDIF
 1300 CONTINUE 
*
      RETURN
      END
************************************************************************
      SUBROUTINE BNDTOR     !" チェビシェフ空間での逆ラプラシアン計算
     O             (   WDTOR , 
     F                 ORIGDB, ORIGDT, 
     C                 FLAPLA,
     C                 ARAD  , DRAD  )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       WDTOR ( NMDIM , KMAX  )  !" トロイダル Ψ
*
*   [INPUT] 
      LOGICAL    ORIGDB                      !" 下端力学境界条件スイッチ
      LOGICAL    ORIGDT                      !" 上端力学境界条件スイッチ
*
      REAL       FLAPLA( NMDIM  )         !" ラプラシアンの係数
*
      REAL       ARAD  ( KMAX )           !" ｒレベル(整数)
      REAL       DRAD  ( KMAX )           !" Δｒ(整数)
*
*   [INTERNAL WORK] 
      REAL       TMATRX( NMDIM, KMAX, KMAX )   !" トロイダルの行列
*                                              !" 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 TORMTX                   !" D_l 演算子の行列設定
     O        (   TMATRX,
     F            ORIGDB, ORIGDT, 
     C            FLAPLA,
     C            ARAD  , DRAD  )
*
         CALL  LUMAKE                  !" 行列のＬＵ分解
     M         ( TMATRX,
     O           KP    ,
     D           NMDIM , KMAX )
      ENDIF
*
*"   < 1. D_l の逆演算 >
*
      DO 100 K = 1, KMAX
         DO 100 NM = 1, NMDIM
            WORK( NM, K ) = WDTOR( NM, K )
 100  CONTINUE 
*
*"  < 境界条件に対応する行を 0 に設定 >
*
      DO 1000 NM = 1, NMDIM
            WORK( NM , 1 ) = 0.0
            WORK( NM , KMAX   ) = 0.0
 1000 CONTINUE 
*
*"  < ＬＵ分解による解の計算 >
*
      CALL LUSOLV
     M         ( WORK  ,
     I           TMATRX, KP    , 
     D              1  , NMDIM , KMAX )
*
      CALL C2W
     O         ( WDTOR ,
     I           WORK  ,
     I           ' '   , 'POS ' , KMAX )
*
      RETURN
      END
************************************************************************
      SUBROUTINE DLINVS     !" チェビシェフ空間での逆ラプラシアン計算
     O             (   WDPOR , 
     I                 WDDPOR, 
     F                 ORIGDB, ORIGDT, 
     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 Φ
*
      LOGICAL    ORIGDB                      !" 下端力学境界条件スイッチ
      LOGICAL    ORIGDT                      !" 上端力学境界条件スイッチ
*
      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 PORMTX                   !" D_l 演算子の行列設定
     O        (   DLMTRX,
     F            ORIGDB, ORIGDT, 
     C            FLAPLA,
     C            ARAD  , DRAD  )
*
         CALL  LUMAKE                  !" 行列のＬＵ分解
     M         ( DLMTRX,
     O           KP    ,
     D           NMDIM , KMAX )
      ENDIF
*
*"   < 1. D_l の逆演算 >
*
      DO 100 K = 1, KMAX
         DO 100 NM = 1, NMDIM
            WORK( NM, K ) = WDDPOR( NM, K )
 100  CONTINUE 
*
*"  < 境界条件に対応する行を 0 に設定 >
*
      DO 1000 NM = 1, NMDIM
            WORK( NM , 1 ) = 0.0
            WORK( NM , 2 ) = 0.0
            WORK( NM , KMAX-1 ) = 0.0
            WORK( NM , KMAX   ) = 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 SETMTX    
*
#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] 
      LOGICAL    ORIGDB                      !" 下端力学境界条件スイッチ
      LOGICAL    ORIGDT                      !" 上端力学境界条件スイッチ
*
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
*
      REAL       ARAD  ( KMAX )              !" ｒレベル(整数)
      REAL       DRAD  ( KMAX )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
*
      REAL       T      ( KMAX , KMAX )       !" チェビシェフ多項式 Tj(xi)
      REAL       DTDR   ( KMAX , KMAX )       !" 1 階微分 dTj/dr(xi)
      REAL       D2TDR2 ( KMAX , KMAX )       !" 2 階微分 d^2Tj/dr^2(xi)
*
      REAL       CDATA  ( NMDIM, KMAX )      !" チェビシェフ作業領域
      REAL       CDDATA ( NMDIM, KMAX )      !" チェビシェフ作業領域
      REAL       WDATA  ( NMDIM, KMAX )      !" スペクトル作業領域
*
      INTEGER    NM, KI, KJ
*
      SAVE       T , DTDR , D2TDR2
*
*" << 1. 微分 etc.の行列表現  >>
*
      DO 1000 KJ = 1, KMAX
         CALL RESET ( CDATA, NMDIM*KMAX )
         CDATA( 1, KJ ) = 1.
*
*"     < 1.0  0 階微分 >
*
         CALL C2W
     M         ( WDATA ,
     I           CDATA ,
     I           ' ' , 'POS' , KMAX )
*
         DO 1010 KI = 1, KMAX
            T( KI, KJ ) = WDATA( 1, KI )
 1010    CONTINUE 
*
*"     < 1.1  1 階微分 >
*
         CALL CDERIV
     O         ( CDDATA,
     I           CDATA ,
     I           KMAX    )
*
         CALL C2W
     M         ( WDATA ,
     I           CDDATA,
     I           ' ' , 'POS' , KMAX )
*
         DO 1100 KI = 1, KMAX
            DTDR( KI, KJ ) = WDATA( 1, KI )
 1100    CONTINUE 
*
*"     < 1.2  2 階微分 >
*
         CALL C2W
     M         ( WDATA ,
     I           CDDATA,
     I           'DR' , 'POS' , KMAX )
*
         DO 1200 KI = 1, KMAX
            D2TDR2( KI, KJ ) = WDATA( 1, KI )
 1200    CONTINUE 
*
 1000 CONTINUE 
      RETURN
**********************************************************************
      ENTRY TORMTX      !" トロイダルポテンシャルの境界条件行列設定
     O         (   DLMTRX,
     F             ORIGDB, ORIGDT, 
     C             FLAPLA,
     C             ARAD  , DRAD  )
*
*"  << 2. 行列の設定 >>
*
      CALL RESET ( DLMTRX, NMDIM*KMAX*KMAX )
*
*"    < 2.1. トロイダルポテンシャルの行列表現設定 >
*
      DO 2100 KJ = 1, KMAX
         DO 2100 KI = 1, KMAX
            DO 2100 NM = 1, NMDIM
               DLMTRX( NM,KI,KJ )  =  T( KI,KJ )
 2100 CONTINUE 
*
*"    < 2.2 境界条件 >
*
*"        < 上端境界 : rigid or stress free >
*
      IF( ORIGDT )THEN
         DO 2210 KJ = 1, KMAX
            DO 2210 NM = 1, NMDIM
               DLMTRX( NM, 1, KJ ) = T( 1, KJ )
 2210    CONTINUE 
      ELSE
         DO 2220 KJ = 1, KMAX
            DO 2220 NM = 1, NMDIM
               DLMTRX( NM, 1, KJ ) 
     &             =   DTDR( 1, KJ ) / ARAD( 1 ) 
     &               - T( 1, KJ ) / ARAD( 1 )**2 
 2220    CONTINUE 
      ENDIF
*
*"        < 下端境界 : rigid or stress free >
*
      IF( ORIGDB )THEN
         DO 2230 KJ = 1, KMAX
            DO 2230 NM = 1, NMDIM
               DLMTRX( NM, KMAX, KJ ) = T( KMAX, KJ )
 2230    CONTINUE 
      ELSE
         DO 2240 KJ = 1, KMAX
            DO 2240 NM = 1, NMDIM
               DLMTRX( NM, KMAX, KJ ) 
     &             =   DTDR( KMAX, KJ ) / ARAD( KMAX ) 
     &               - T( KMAX, KJ ) / ARAD( KMAX )**2 
 2240    CONTINUE 
      ENDIF
*
      RETURN
**********************************************************************
      ENTRY PORMTX      !" ポロイダルポテンシャルの境界条件行列設定
     O         (   DLMTRX, 
     F             ORIGDB, ORIGDT, 
     C             FLAPLA,
     C             ARAD  , DRAD    )
*
*"  < D_l 演算子の行列設定(チェビシェフ空間) >
*
*"  << 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 ) 
     &              =    D2TDR2( KI,KJ ) 
     &                 + 2 * DTDR( KI,KJ ) / ARAD( KI )
     &                 + FLAPLA( NM ) * T( KI,KJ ) / ARAD( KI )**2 
 3100 CONTINUE 
*
*"    < 3.2 境界条件 >
*
*"        < 上下共通 : no penetration condition >
*
      DO 3200 KJ = 1, KMAX
            DO 3200 NM = 1, NMDIM
               DLMTRX( NM, 1, KJ ) = T( 1, KJ ) 
               DLMTRX( NM, KMAX, KJ ) = T( KMAX, KJ )
 3200 CONTINUE 
*
*"        < 上端境界 : rigid or stress free >
*
      IF( ORIGDT )THEN
         DO 3210 KJ = 1, KMAX
            DO 3210 NM = 1, NMDIM
               DLMTRX( NM, 2, KJ ) = DTDR( 1, KJ )
 3210    CONTINUE 
      ELSE
         DO 3220 KJ = 1, KMAX
            DO 3220 NM = 1, NMDIM
               DLMTRX( NM, 2, KJ ) = D2TDR2( 1, KJ )
 3220    CONTINUE 
      ENDIF
*
*"        < 下端境界 : rigid or stress free >
*
      IF( ORIGDB )THEN
         DO 3230 KJ = 1, KMAX
            DO 3230 NM = 1, NMDIM
               DLMTRX( NM, KMAX-1, KJ ) = DTDR( KMAX, KJ )
 3230    CONTINUE 
      ELSE
         DO 3240 KJ = 1, KMAX
            DO 3240 NM = 1, NMDIM
               DLMTRX( NM, KMAX-1, KJ ) = D2TDR2( KMAX, KJ )
 3240    CONTINUE 
      ENDIF
*
      RETURN
      END
