* PACKAGE TDZ2PSI    !" テスト : 渦度から流線関数を求める(境界条件 psi=const.)
*
*"  [HIS] 95/07/11 (takepiro) 
*"        95/08/28 (takepiro) 
*
***********************************************************************
      PROGRAM TDZ2PS
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [OUTPUT]
      REAL       PSI  ( 0:NX+1 , 0:NZ+1 )  !" 流線関数
*
* [INPUT]
      REAL       ZETA ( 0:NX+1 , 0:NZ+1 )  !" 渦度
*
      REAL       DX                        !" X方向の格子点間隔
      REAL       DZ                        !" Z方向の格子点間隔
*
* [INTERNAL WORK]
      INTEGER IX
      INTEGER IZ
      REAL    PSIBTM, PSITOP
      REAL    PI, XX, ZZ		           !" 円周率
*
      PI =  ATAN2( 0.0, -1.0 )
      DX=2*PI/REAL(NX)
      DZ=1/REAL(NZ+1)
*
*"   < 1. 境界条件の設定 >
*
      PSIBTM = 0.0
      PSITOP = 1.0
*
*"   < 0. 渦度分布の設定 >
*
      DO 100 IZ = 0, NZ+1
        DO 100 IX = 0, NX+1
c$$$          XX = (IX-1)*DX
c$$$          ZZ = IZ*DZ
          ZETA( IX, IZ ) = 0.0
c$$$          ZETA( IX, IZ ) = COS( XX ) * SIN( PI*ZZ )
c$$$          ZETA( IX, IZ ) = SIN( PI*ZZ )
c$$$           IF ( IZ.EQ.NZ/2) THEN
c$$$              ZETA( IX, IZ ) = COS( 2 *  PI * (IX-1) / REAL(NX) )
c$$$           ENDIF
  100 CONTINUE
*
      WRITE ( 6, * ) 'Vorticity'
      DO 110 IZ=0,NZ+1
	   WRITE( 6, '(8F9.3)' ) (ZETA(IX,IZ), IX=1,NX)
  110 CONTINUE
*
      CALL ZTOPSI
     O     ( PSI   , 
     I       ZETA  , 
     I       PSIBTM, PSITOP,
     D       DX    , DZ      )
*
*"   < 3. 出力 >
*
      WRITE ( 6, * ) 'Stream Function'
      DO 300 IZ=0,NZ+1 
	   WRITE( 6, '(8F9.3)' ) (PSI(IX,IZ), IX=1,NX)
  300 CONTINUE
*
      STOP
      END

