C*********************************************************************
C*********************************************************************
C
      SUBROUTINE CLPRES
     M   (  DUTOTL, DWTOTL,
     I      CONVM,
     O      PFUNC, PPAI,
     W      PDIV ,  PDIVM, PFUNCM, ALPHA, GAMMA,
     M      VMODE, TMODE, FMODE, TRIGS, IFAC,
     W      UUFT , ANFT ,
     G         DX,
     C      CPDRY,
     B      FDZM, TPOT0,
     D      NX, NZ, NXB, NZB                         )
C
C*********************************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C*********************************************************************
C
      DIMENSION
     \   DUTOTL( -2:NXB, -2:NZB ),  DWTOTL( -2:NXB, -2:NZB )
C
      DIMENSION
     \   PFUNC ( -2:NXB, -2:NZB ), PPAI ( -2:NXB, -2:NZB )
C
      DIMENSION
     \   PDIV  ( -2:NXB, -2:NZB ), PDIVM ( -2:NXB, -2:NZB ),
     \   GAMMA ( -2:NXB, -2:NZB ), CONVM ( -2:NXB, -2:NZB ),
     \   PFUNCM( -2:NXB, -2:NZB ), ALPHA ( -2:NXB, -2:NZB )
C
      DIMENSION
     \   VMODE ( -2:NZB, NZ  ),  TMODE( -2:NZB, NZ  )
C
      DIMENSION
     \   FMODE  ( -2:NXB, NZ )

c     dimension
c    \   UUFT ( -2:NXB, NZ ),
c    \   ANFT ( -2:NXB, NZ )

      dimension
     \   UUFT ( NZ, NX ),  ANFT ( NZ, NX )
C
      DIMENSION
     \   TRIGS( * ), IFAC(15)
C
      DIMENSION
     \   FDZM ( -2:NZB ), TPOT0 ( -2:NZB )
C
C=====================================================================
C
C     WRITE(6,*)'BUOY'
C     CALL ADUMP( BUOY, NXB, NZB, 1, 1 )
C
      DO 10 IX = 1, NX
         GAMMA(IX,0) = FDZM(0) * DWTOTL(IX,0)
   10 CONTINUE
C
      DO 11 IZ = 1, NZ
      DO 11 IX = 1, NX
         GAMMA(IX,IZ) =    GAMMA(IX,IZ-1)
     \                    + FDZM(IZ) * DWTOTL(IX,IZ)
   11 CONTINUE
C
      CALL BOUND( GAMMA, NX, NZ, NXB, NZB, 1, 1, 0, 1 )
C     WRITE(6,*)'GAMMA'
C     CALL ADUMP( GAMMA, NXB, NZB, 1, 1 )
C
      DO 12 IZ = 0, NZ-1
      DO 12 IX = 1, NX
         ALPHA(IX,IZ) = - ( GAMMA(IX,IZ) - GAMMA(IX-1,IZ) )
     \           + DUTOTL(IX,IZ)
   12 CONTINUE
C
C     WRITE(6,*)'ALPHA'
C     CALL ADUMP(ALPHA, NXB, NZB, 1, 1 )
      CALL BOUND( ALPHA, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
      DO 21 IZ = 0, NZ-1
      DO 21 IX = 1, NX
         PDIV(IX,IZ) =    CONVM(IX,IZ)
     \          + 1.125D0 * DX
     \              * ( ALPHA(IX+1,IZ) - ALPHA(IX,IZ) )
     \          - 1.0D0 / 24.0D0 * DX
     \              * ( ALPHA(IX+2,IZ) - ALPHA(IX-1,IZ) )
   21 CONTINUE
C
C  *** MODE SEPARATION ***
C
      DO 31 IM = 1, NZ
C
         DO 301 IX = 1, NX
            PDIVM(IX,IM) = 0.0
  301    CONTINUE
C
         DO 302 IZ = 0, NZ-1
         DO 302 IX = 1, NX
            PDIVM(IX,IM) =    PDIVM(IX,IM)
     \                      + TMODE(IZ,IM) * PDIV(IX,IZ)
  302    CONTINUE
C
   31 CONTINUE
C
      DO 32 IZ = 0, NZ-1
         DO 322 IM = 1, NZ
         DO 322 IX = 1, NX
            PDIV(IX,IZ) =    PDIV(IX,IZ)
     \                      - VMODE(IZ,IM) * PDIVM(IX,IM)
  322    CONTINUE
   32 CONTINUE
C
C     WRITE(6,*)'PDIVM '
C     CALL ADUMP( PDIVM, NXB, NZB, 1, 1 )
C
C     WRITE(6,*)'PDIV DECOMPOSITION ERROR'
C     CALL ADUMP( PDIV, NXB, NZB, 1, 1 )
C                                        下のループはデバッグ用
C     DO 361 IM = 1, NZ
C     DO 361 IX = 1, NX
C        PDIV(IX,IZ) = PDIVM(IX,IZ)
C 361 CONTINUE
C
C     CALL BOUND ( PDIV, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C     ここから　ＦＦＴ　を使う．
C
      DO 51 IM = 1, NZ
      DO 51 IX = 1, NX
c        UUFT(IX,IM) = PDIVM(IX,IM)
         UUFT(IM,IX) = PDIVM(IX,IM)
   51 CONTINUE
C
C     WRITE(6,*) ' IFAC IN CLPRES ', IFAC
c     CALL RFFTFR ( NX, NZ, 1, NX+5, 4, 1, NX+5, 4,
c    \              UUFT, TRIGS, IFAC, ANFT )

      call vrfftf ( nz, nx, uuft, anft, nz, trigs )

C
C     WRITE(6,*)'UUFT  '
C     DO 5001 IM = 1, NZ
C5001 WRITE(6,*) ( UUFT(KX,IM), KX = 1, NX )
C
C     WRITE(6,*)'FMODE '
C     DO 5002 IM = 1, NZ
C5002 WRITE(6,*) ( FMODE(KX,IM), KX = 1, NX )
C
      DO 55 IM = 1, NZ
      DO 55 IK = 1, NX
c        ANFT(IK,IM) = FMODE(IK,IM) * UUFT(IK,IM)
        ANFT(IM,IK) = FMODE(IK,IM) * UUFT(IM,IK)
   55 CONTINUE
C
C     WRITE(6,*)'ANFT  '
C     DO 5004 IM = 1, NZ
C5004 WRITE(6,*) ( ANFT(KX,IM), KX = 1, NX )
C
C
c     CALL RFFTBR ( NX, NZ, 1, NX+5, 4, 1, NX+5, 4,
c    \              ANFT, TRIGS, IFAC, UUFT )

      call vrfftb ( nz, nx, anft, uuft, nz, trigs )

C
      DO 61 IM = 1, NZ
      DO 61 IX = 1, NX
         PFUNCM(IX,IM) = ANFT(IM,IX)
   61 CONTINUE
C
      CALL BOUND ( PFUNCM, NX, NZ, NXB, NZB, 1, 1, 0, 1)
C     WRITE(6,*)'PFUNCM'
C     CALL ADUMP( PFUNCM, NXB, NZB, 1, 1 )
C
C     DO 661 IM = 1, NZ-1
C     DO 661 IX = 1, NX
C        PDIVM(IX,IM) =
C    \     (  - 54.0D0 * PFUNCM(IX,IM)
C    \        + 28.0D0 * ( PFUNCM(IX-1,IM) + PFUNCM(IX+1,IM) )
C    \        -        ( PFUNCM(IX-2,IM) + PFUNCM(IX+2,IM) ) ) /24.0D0
C    \     + PFUNCM(IX,IM) / FMODE(1,IM)
C 661 CONTINUE
C     DO 662 IX = 1, NX
C        IM = NZ
C        PDIVM(IX,IM) =
C    \     (  - 54.0D0 * PFUNCM(IX,IM)
C    \        + 28.0D0 * ( PFUNCM(IX-1,IM) + PFUNCM(IX+1,IM) )
C    \        -        ( PFUNCM(IX-2,IM) + PFUNCM(IX+2,IM) ) ) /24.0D0
C 662 CONTINUE
C
C     WRITE(6,*)'PDIVM RECONSTRUCTED FROM PFUNC'
C     CALL ADUMP( PDIVM , NXB, NZB, 1, 1 )
C
      DO 71 IZ = 0, NZ
         DO 701 IX = 1, NX
            PFUNC(IX,IZ) = 0.0
  701    CONTINUE
         DO 702 IM = 1, NZ
         DO 702 IX = 1, NX
            PFUNC(IX,IZ) = PFUNC(IX,IZ)
     \         +  VMODE(IZ,IM) * PFUNCM(IX,IM)
  702    CONTINUE
   71 CONTINUE
C
      CALL BOUND ( PFUNC, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C---------------------------------------------------------------------
C     Perturbation non-dimensional pressure function.
C---------------------------------------------------------------------
C
      DO 81 IZ = 0, NZ-1
      DO 81 IX = 1, NX
         PPAI (IX,IZ) =
     \      PFUNC(IX,IZ) + GAMMA(IX,IZ) * DX
   81 CONTINUE
C
      CALL BOUND ( PPAI , NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
      DO 91 IX =  -2, NXB
         PPAI (IX,-1) =    PPAI (IX,0)
     \                   - FDZM(0) * DX * DWTOTL(IX,0)
         PPAI (IX,NZ) =    PPAI (IX,NZ-1)
     \                   + FDZM(NZ)* DX * DWTOTL(IX,NZ)
   91 CONTINUE
C
      DO 92 IX =  -2, NXB
         PPAI (IX,-2) =    PPAI (IX,-1)
     \    - FDZM(-1)* DX * DWTOTL(IX,-1)
     \    - FDZM( 1)* DX * DWTOTL(IX, 1)
     \            +  ( PPAI (IX,1) - PPAI (IX,0) )
         PPAI (IX,NZ+1) =    PPAI (IX,NZ  )
     \    + FDZM(NZ+1) * DX * DWTOTL(IX,NZ+1)
     \    + FDZM(NZ-1) * DX * DWTOTL(IX,NZ-1)
     \            -  ( PPAI (IX,NZ-1) - PPAI (IX,NZ-2) )
   92 CONTINUE
C
      DO 101 IZ = -2, NZB
      DO 101 IX = -2, NXB
         PPAI(IX,IZ) = PPAI (IX,IZ) / CPDRY / TPOT0(IZ)
  101 CONTINUE
C
C--------------------------------------------------------------------
C     Time derivative of U and W.
C--------------------------------------------------------------------
C
      DO 201 IZ = 0, NZ-1
      DO 201 IX = 1, NX
         DUTOTL(IX,IZ) =
     \       + ALPHA(IX,IZ)
     \       - ( PFUNC(IX,IZ) - PFUNC(IX-1,IZ) ) / DX
  201 CONTINUE
C
      DO 202 IZ = 1, NZ-1
      DO 202 IX = 1, NX
         DWTOTL(IX,IZ) =
     \       - ( PFUNC(IX,IZ) - PFUNC(IX,IZ-1) ) / DX / FDZM(IZ)
  202 CONTINUE
C
      DO 203 IX = 1, NX
         DWTOTL(IX,0 ) = 0.0D0
         DWTOTL(IX,NZ) = 0.0D0
  203 CONTINUE
C
      CALL BOUND ( DUTOTL, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( DWTOTL, NX, NZ, NXB, NZB, 1,-1, 0, 1 )
C
C=====================================================================
C
      RETURN
      END
