* PACKAGE DZ2PSI    !" 渦度から流線関数を求める(境界条件 psi=const.)
*
*"  [HIS] 95/08/09 (takepiro) for psi=0 
*"        95/08/29 (takepiro) for psi=const
***********************************************************************
      SUBROUTINE DZ2PSI     !" 渦度から流線関数を求める(境界条件 psi=const)
     O           ( PSI   , 
     I             ZETA  ,
     I             PSIBTM, PSITOP, 
     I             DX    , DZ      )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [OUTPUT]
      REAL       PSI  ( 0:NX+1 , 0:NZ+1 )  !" 流線関数
*
* [INPUT]
      REAL       ZETA ( 0:NX+1 , 0:NZ+1 )  !" 渦度
      REAL       PSIBTM         !" 下端境界の流線関数値
      REAL       PSITOP         !" 下端境界の流線関数値
      REAL       DX             !" X方向の格子点間隔
      REAL       DZ             !" Z方向の格子点間隔
*
* [INTERNAL SAVE]
      REAL      ALU  ( NX , NZ, -1:1  ) !" 係数行列(3重対角)
      REAL      WTBL ( 2  , NX ) !" 三角関数表
      LOGICAL   OFIRST
*
* [INTERNAL WORK]
      REAL      Z    ( NX , NZ ) !" 作業配列
      INTEGER 	IX, IZ
      INTEGER 	IFAC
*
      DATA    OFIRST / .TRUE. /
*
      SAVE    OFIRST, ALU, WTBL
*
      IF ( OFIRST ) THEN
         OFIRST = .FALSE.
         CALL DZ2PST     !" ポアソン方程式の係数行列を設定する
     O        (  ALU  , 
     I           DX   , DZ  )
*
         CALL RFFTIM            !" FFT 3角関数表作成
     I        ( NX    ,
     O          WTBL  , IFAC  )
*
      ENDIF
*
*"        < 1. ZETAをフーリエ変換する. >
*
      CALL COPY( PSI , ZETA , (NX+2)*(NZ+2) )
*
      DO 1000 IX = 0, NX+1      !" 境界条件の組みこみ
         PSI( IX, 1 ) = PSI( IX, 1 ) - PSIBTM/(DZ)**2
         PSI( IX,NZ ) = PSI( IX,NZ ) - PSITOP/(DZ)**2
 1000 CONTINUE 
*
      CALL FFT99Y                          !" FFT を呼ぶ
     M         ( PSI   ,
     M           Z     ,
     C           WTBL  , 
     C           2     , NX  , NX+2  , 
     C           2     , NZ  , NZ+2  , 
     C           0                      )
*
*"        < 2. PSIのフーリエ成分を求める. >
*
      CALL LUSOL3   !" ＬＵ分解による解の計算 [ ３重対角行列 ]
     M      ( Z    ,
     I        ALU  ,
     D        1    , NX , NZ )
*
*"       < 3. 逆フーリエ変換してPSIを求める. >
*
      CALL FFT99Y     !" FFT を呼ぶ
     M         ( PSI   ,
     M           Z     ,
     C           WTBL  , 
     C           2     , NX  , NX+2  , 
     C           2     , NZ  , NZ+2  , 
     C           1                      )
*
*"       < 4. 境界条件の適用 >
*
      DO 4000 IX = 0, NX+1
	PSI( IX, 0 )    = PSIBTM
	PSI( IX, NZ+1 ) = PSITOP
 4000 CONTINUE
*
      DO 4100 IZ = 0, NZ+1
         PSI(    0, IZ ) = PSI( NX, IZ ) !" 周期的境界条件
         PSI( NX+1, IZ ) = PSI(  1, IZ ) !" 周期的境界条件
 4100 CONTINUE
*
      RETURN
      END
************************************************************************
      SUBROUTINE DZ2PST     !" ポアソン方程式の係数行列を設定する
     O	         (  ALU  , 
     I              DX   , DZ  )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [output]
      REAL       ALU  ( NX , NZ, -1:1  )   !" 係数行列(3重対角)
*
      REAL       DX                        !" X方向の格子点間隔
      REAL       DZ                        !" Z方向の格子点間隔
*
* [INTERNAL WORK]
      INTEGER IWN
      INTEGER IZ
      REAL    PI		           !" 円周率
*
      SAVE
*
      PI=4.0*ATAN(1.0D0)
*
*        < 1. 係数行列の設定 >
*
*"          < 1.1 対角成分の設定 >
*
      DO 1000 IZ=1,NZ
          ALU(1,IZ,0) = - 2./(DZ**2)
          ALU(2,IZ,0) = - ( 2 * PI* (NX/2)/(NX*DX) )**2
     &                  - 2./(DZ**2)
 1000 CONTINUE

      DO 1100 IWN=4,NX,2
        DO 1100 IZ=1,NZ
          ALU(IWN,IZ,0) = - ( 2 * PI* (IWN/2-1)/(NX*DX) )**2
     &                 - 2./(DZ**2)
          ALU(IWN-1,IZ,0) =  ALU(IWN,IZ,0) 
 1100 CONTINUE
*
*"          < 1.2 非対角成分の設定 >
*"		境界条件が組みこまれていることに注意
*
      DO 1200 IWN=1,NX
        DO 1200 IZ=1,NZ
          ALU(IWN, IZ,  1) = 1./(DZ**2)
          ALU(IWN, IZ, -1) = 1./(DZ**2)
 1200 CONTINUE
*
*"        < 2. LU分解の実行 >
*
      CALL LUMAK3        !" 行列のＬＵ分解 [ ３重対角行列 ]
     M     ( ALU ,
     D       NX  , NZ )
*
      RETURN
      END

