* PACKAGE DZ2PSI    !" 渦度から流線関数を求める(境界条件 psi=const.)
*
*"  [HIS] 96/03/14 (takepiro) TZ2PSI より移植
*
***********************************************************************
*"        <<< △ζ = Ｒａ ＤＴ／Ｄｘ を解く >>>
**********************************************************************
      SUBROUTINE DTMP2Z          !"  温度から渦度の計算(P inf 近似)
     O             ( ZETA  , 
     I               T     , 
     C               RAYL  , 
     C               DX    , DZ    , DT  ,
     W               RDTDX    )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [OUTPUT]
      REAL ZETA ( 0:NX+1 , 0:NZ+1 )  !" 渦度(現STEPの値)
*
* [INPUT]
      REAL T    ( 0:NX+1 , 0:NZ+1 )  !" 温度(現STEPの値)
*
      REAL RAYL                      !" レイリー数
*
      REAL DX                        !" X 軸上の格子点間隔
      REAL DZ                        !" Z 軸上の格子点間隔
      REAL DT                        !" 時間間隔
*
* [WORK]
      REAL RDTDX ( 0:NX+1 , 0:NZ+1 ) !" 浮力項
*
*  [INTERNAL WORK]
      INTEGER IX                     !" ループカウンタ
      INTEGER IZ                     !" ループカウンタ
*
*"        < 1. 浮力項の計算 >
*
      DO 1020 IZ = 1 , NZ
        DO 1010 IX = 1 , NX
           RDTDX( IX,IZ ) = RAYL*( T(IX+1,IZ)-T(IX-1,IZ) )/(2.*DX) 
 1010   CONTINUE
 1020 CONTINUE
*
*"        < 2. ポアソン方程式を解く >
*
      CALL DTDX2Z         !" 温度傾度から渦度を求める(境界条件 zeta=0)
     O     ( ZETA , 
     I       RDTDX,
     I       DX    , DZ   )
*
      RETURN
      END
**********************************************************************
      SUBROUTINE DTDX2Z     !" 温度傾度から渦度を求める(境界条件 zeta=0)
     O           ( ZETA , 
     I             RDTDX,
     I             DX    , DZ      )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [OUTPUT]
      REAL       ZETA ( 0:NX+1 , 0:NZ+1 )  !" 渦度
*
* [INPUT]
      REAL       RDTDX( 0:NX+1 , 0:NZ+1 ) !" 浮力項
      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. RDTDX をフーリエ変換する. >
*
      CALL COPY( ZETA , RDTDX , (NX+2)*(NZ+2) )
*
      CALL FFT99Y                          !" FFT を呼ぶ
     M         ( ZETA  ,
     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         ( ZETA  ,
     M           Z     ,
     C           WTBL  , 
     C           2     , NX  , NX+2  , 
     C           2     , NZ  , NZ+2  , 
     C           1                      )
*
*"       < 4. 境界条件の適用 >
*
      DO 4000 IX = 0, NX+1
	ZETA( IX, 0 )    = 0.0
	ZETA( IX, NZ+1 ) = 0.0
 4000 CONTINUE
*
      DO 4100 IZ = 0, NZ+1
         ZETA(    0, IZ ) = ZETA( NX, IZ ) !" 周期的境界条件
         ZETA( NX+1, IZ ) = ZETA(  1, IZ ) !" 周期的境界条件
 4100 CONTINUE
*
      RETURN
      END