* PACKAGE UCDIFF    !" 差分の計算
*
*"  [HIS] 95/08/14 (takepiro) 下位ルーチン切りだし, スキーム要チェック
***********************************************************************
      SUBROUTINE ARAJAC               !" ARAKAWA ヤコビアンの計算
     O           ( AJAC ,
     I             PSI  , G   ,
     C             DX   , DZ  ,
     I             ISX  , IFX , ISZ , IFZ ,
     D             NX   , NZ    )
*
*"      << ARAKAWA schemeを用いてヤコビアン 
*"           J(PSI,G) = dp/dx dg/dz - dp/dz dg/dx を計算する >>
* [OUTPUT]
      REAL AJAC( NX, NZ )       !" J(PSI,G)
*
* [INPUT]
      REAL PSI ( NX, NZ )       !" 流線関数
      REAL G   ( NX, NZ )       !" 物理量の2次元場
*
      REAL DX                   !" X方向の格子点間隔
      REAL DZ                   !" Z方向の格子点間隔
*
      INTEGER ISX               !" X 方向にヤコビアンを計算し始める点
      INTEGER IFX               !" X 方向にヤコビアンを計算し終える点
      INTEGER ISZ               !" Z 方向にヤコビアンを計算し始める点
      INTEGER IFZ               !" Z 方向にヤコビアンを計算し終える点
*
      INTEGER NX                !" X 方向の配列の次元
      INTEGER NZ                !" Z 方向の配列の次元
*
* [INTERNAL WORK]
      REAL RJPP
      REAL RJMP
      REAL RJPM
*
      INTEGER IX
      INTEGER IZ
*
      DO 1020 IZ = ISZ, IFZ
         DO 1010 IX =  ISX, IFX
          RJPP =     ( PSI(IX+1,IZ)-PSI(IX-1,IZ) )/(2.0D0*DX)
     &               * ( G(IX,IZ+1)-G(IX,IZ-1) )/(2.0D0*DZ)
     &           -   ( PSI(IX,IZ+1)-PSI(IX,IZ-1) )/(2.0D0*DZ)
     &             * ( G(IX+1,IZ)-G(IX-1,IZ) )/(2.0D0*DX)
          RJMP =     (   G(IX,IZ+1)*( PSI(IX+1,IZ+1)-PSI(IX-1,IZ+1) )
     &                 - G(IX,IZ-1)*( PSI(IX+1,IZ-1)-PSI(IX-1,IZ-1) ))
     &             / (4.0D0*DX*DZ)
     &           -   (   G(IX+1,IZ)*( PSI(IX+1,IZ+1)-PSI(IX+1,IZ-1) )
     &                 - G(IX-1,IZ)*( PSI(IX-1,IZ+1)-PSI(IX-1,IZ-1) ))
     &             / (4.0D0*DX*DZ)
          RJPM =     (   PSI(IX+1,IZ)*( G(IX+1,IZ+1)-G(IX+1,IZ-1) )
     &                 - PSI(IX-1,IZ)*( G(IX-1,IZ+1)-G(IX-1,IZ-1) ))
     &             / (4.0D0*DX*DZ)
     &           -   (   PSI(IX,IZ+1)*( G(IX+1,IZ+1)-G(IX-1,IZ+1) )
     &                 - PSI(IX,IZ-1)*( G(IX+1,IZ-1)-G(IX-1,IZ-1) ))
     &             / (4.0D0*DX*DZ)
          AJAC( IX,IZ ) = ( RJPP + RJMP + RJPM )/3.0D0
 1010   CONTINUE
 1020 CONTINUE
*
      RETURN
      END
***********************************************************************
      SUBROUTINE CLAPLA             !" ラプラシアンの差分計算を行う
     O           ( ALAP,
     I             F   ,
     C             DX  , DZ  , 
     I             ISX , IFX , ISZ , IFZ ,
     D             NX  , NZ    )
*
* [OUTPUT]
      REAL    ALAP( NX, NZ  )   !" 物理量のラプラシアン
*
* [INPUT]
      REAL    F( NX, NZ )       !" 物理量の2次元場
*
      REAL    DX                !" X 方向の格子点間隔
      REAL    DZ                !" Z 方向の格子点間隔
*
      INTEGER ISX               !" X 方向にラプラシアンを計算し始める点
      INTEGER IFX               !" X 方向にラプラシアンを計算し終える点
      INTEGER ISZ               !" Z 方向にラプラシアンを計算し始める点
      INTEGER IFZ               !" Z 方向にラプラシアンを計算し終える点
*
      INTEGER NX                !" X 方向の配列の次元
      INTEGER NZ                !" Z 方向の配列の次元
*
* [INTERNAL WORK]
      INTEGER IX                !" ループカウンタ
      INTEGER IZ                !" ループカウンタ
*
*"        < 1. ラプラシアンの計算 >
*
      DO 1020 IZ = ISZ, IFZ
        DO 1010 IX = ISX, IFX
          ALAP(IX,IZ) =     ( F(IX-1,IZ)-2.0D0*F(IX,IZ)+F(IX+1,IZ) )
     &                    / DX**2
     &                  +   ( F(IX,IZ-1)-2.0D0*F(IX,IZ)+F(IX,IZ+1) )
     &                    / DZ**2
 1010   CONTINUE
 1020 CONTINUE
*
      RETURN
      END
