*" PACKAGE GTSTSPCT 時空間フーリエフィルター
*" FFT の都合で時系列の長さは 2,3,5 で割りきれる偶数でなければならない。
*" [HIS] 95/12/14 堀之内武  (gtspspct を利用)
***********************************************************************
      PROGRAM GTSTSP
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)
#else
#include         "gzsize.F"
#endif
*
*" [PARAM] 必要な領域は約 (2*NDIM4+3*IJKDIM)*4 バイト. 確保可能か注意
      PARAMETER (NDIM4 = 30E5)   !" 読み込む全データ数の最大値
      PARAMETER (IJKDIM= 15E4)   !" 1レコードのデータ数の最大値
**      PARAMETER (NDIM4 = 7536700 )   !" 読み込む全データ数の最大値
**      PARAMETER (IJKDIM= 94300  )   !" 1レコードのデータ数の最大値
      PARAMETER (NTMAX = 1E4)    !" 時系列の個数の最大値
*
      REAL       GDAT4( NDIM4 )
      CHARACTER  HHEAD( NDC )*(NCC)
      INTEGER    ITIME(NTMAX)
      CHARACTER  HDATE(NTMAX)  *(NCC)
      INTEGER    IDNUM(NTMAX)
*
      CHARACTER  HFILE( 1 ) *(NFILN)
      DATA       HFILE / '$GTTMPDIR/gtool.out' /
      CHARACTER  OUT( 1 )    *(NFILN)
      DATA       OUT / '$GTTMPDIR/gtool.out' /
      LOGICAL    APND
      DATA       APND   / .FALSE. /
*
      LOGICAL    EWD  !" T/F: pass only eastward/westward components
      DATA       EWD/.TRUE./
      INTEGER    MMIN,MMAX !"pass wavemumber M s.t. (0<) MIN<=M<=MAX
      DATA       MMIN,MMAX /1,-999/ !" negative MAX doesn't allow high cut
      INTEGER    FMIN,FMAX !"pass freq F(ineger) s.t. (0<) FMIN<=F<=FMAX
      DATA       FMIN,FMAX /1,-999/ !" negative MAX doesn't allow high cut
*
      INTEGER    M4(4),F4(4) !" 波数振動数の範囲を4辺形で設定する.
*"               (M4(i),F4(i)),i=1,..,4 は反時計回りになっていることを
*"               仮定する. すべての値に負でない値を与えたときのみ有効.
*"               MMIN,MMAX,FMIN,FMAXは有効のまま(だから最大で8角形切断).
      DATA       M4,F4/ 4*-999,4*-999/
*
      CHARACTER  ITEM   *(NCC)
      CHARACTER  UNIT   *(NCC)
      CHARACTER  TITLE  *(NCC*2)
      CHARACTER  DSET   *(NCC)
      DATA       ITEM, UNIT, TITLE, DSET /4*' '/
      CHARACTER  EDIT   *(NCC)
      CHARACTER  ETTL   *(NCC)
      DATA       EDIT   / ' ' /
      DATA       ETTL   / ' ' /
      LOGICAL    GRESET
      DATA       GRESET / .FALSE. /
      LOGICAL    HELP
      DATA       HELP   / .FALSE. /
      INTEGER    SDIM       !" Fourier変換される空間次元は何次元目か
      DATA       SDIM   / 1 /
**      LOGICAL    RSTRD  !"  T -> 空間フーリエの前にリニアトレンドを引く
**      DATA       RSTRD  / .FALSE. /
**      LOGICAL    RTTRD  !"  T -> 時間フーリエの前にリニアトレンドを引く
**      DATA       RTTRD  / .FALSE. /
*
      NAMELIST  /OPTION/ HFILE, OUT, APND,  SDIM, EWD, MMIN, MMAX,
     &                   FMIN, FMAX, M4,  F4,
     &                   GRESET, ITEM, UNIT,  TITLE, DSET, EDIT, ETTL,
     &                   HELP
**     &                   ,RSTRD, RTTRD
*
      DATA       IFILE / 50 /
      DATA       JFILE / 60 /
*
      CHARACTER  HHEADS( NDC )*(NCC)
*
      CHARACTER  EDITP  *(NCC)
      CHARACTER  ETTLP  *(NCC)
      DATA       EDITP  / 'STPH' /
      DATA       ETTLP  / 's-t phase' /
*
      CHARACTER  AITMF  *(NCC)  !" 時間軸ファイル名
      DATA       AITMF  / '@FREQ' /  !" その頭の部分
      CHARACTER  AITMS  *(NCC)  !" 空間軸ファイル名
      DATA       AITMS  / '@WVNWE' /  !" その頭の部分
      LOGICAL    LFIRST
      DATA       LFIRST  /.TRUE./
      LOGICAL    LINIT
      DATA       LINIT  /.TRUE./
      CHARACTER  HAITM1*(NCC),HAITM2*(NCC),HAITM3*(NCC)
      INTEGER    IXDIM,IYDIM,IZDIM,IWDIM !" ４次元データの各要素数
      INTEGER    IRDIM,ITDIM     !" 変換に使われる空間、時間の要素数
      INTEGER    IQDIM,ISDIM !" 時系列は(IQ,IR,IS,IT)DIM の４次元として扱われる
*
*  [WORK]
      REAL       WX(IJKDIM),WT(IJKDIM)      !" work(consts) for FFT
      REAL       WY(NDIM4)                  !" work for FFT
      REAL       SP(NDIM4/2),PH(NDIM4/2)    !" power spectrum, phase
      EQUIVALENCE (WY(1),SP)
      EQUIVALENCE (WY(NDIM4/2+1),PH)
*
      INTEGER    IWT(5),IWX(5)
*
      CALL OPTARG ( 91, 'OPTION', 'HFILE', NOPT, NFILE )
      READ (91,OPTION,IOSTAT=IOS)
      CLOSE(91)
      IF ( HELP ) THEN
         WRITE(6,OPTION)
         WRITE(6,*) 'Filter for Eastward or Westward moving component'
         WRITE(6,*) '  on the basis of space-time Fourier transf.'
         WRITE(6,*) 'Input data has to be a time sequence.'
         WRITE(6,*) 'Options:'
         WRITE(6,*) ' sdim: assume sdim-th axis to be the "space" axis'
         WRITE(6,*) ' ewd: T/F pass eastward/westward moving component'
         WRITE(6,*) ' mmax,mmin: select wavenumber M s.t. min<=M<=max'
         WRITE(6,*) 
     +          ' fmax,fmin: select freq F(integer) s.t. fmin<=F<=fmax'
         WRITE(6,*) 'If all of M4(4),F4(4) are .le. 0,'
         WRITE(6,*) '  select (M,F)s within a quadrilateral ',
     +              '(M4(i),F4(i)),i=1,..,4 (assume anti-clockwise)'
         WRITE(6,*) 
         STOP
      ELSE IF ( IOS.NE.0 ) THEN
         WRITE(6,OPTION)
         STOP
      ENDIF
*
      CALL GTOPEN
      CALL GTSIZE ( HHEAD, IJKDIM )
*
      CALL GURNTF ( HFILE( 1 ), OUT(1)  , '$GTTMPDIR/gtool.in' )
*
      CALL GFROPN ( IFILE , HFILE( 1 ) )
      CALL GFOOPN ( JFILE ,  OUT(1) , APND )
      CALL GUNENV( OUT(1),'.',.FALSE. )
      IL=LENC(OUT(1))
      WRITE (6,*) 'output='//OUT(1)(1:IL)
*
*  < read data >
*
      IANUM = 0
      NGDAT=0
 1100 CONTINUE
        IF ( (IANUM+1)*NGDAT .GT. NDIM4 )
     &         CALL MSGDMP('E',' MAIN ','"NDIM4" INSUFFICIENT')
        CALL   GFREAD
     O        ( HHEAD , GDAT4(1+IANUM*NGDAT) , IEOD  ,
     I          IFILE , 1               )
*
        IF ( IEOD .EQ.0 ) THEN
          IANUM = IANUM + 1
          IF (LFIRST) THEN
            LFIRST = .FALSE.
            CALL   GHCOPY( HHEAD,HHEADS)
            CALL GHCGET( HHEAD, 'AITM1', HAITM1 )
            CALL GHCGET( HHEAD, 'AITM2', HAITM2 )
            CALL GHCGET( HHEAD, 'AITM3', HAITM3 )
            IF  ( HAITM1 .EQ. ' ' ) 
     &          CALL MSGDMP( 'E',' MAIN ','ZERO SPACE DIMENSION' )
            CALL GUQSIZ
     I         ( HHEAD ,
     O           IXSTR , IXEND , IXDIM ,
     O           IYSTR , IYEND , IYDIM ,
     O           IZSTR , IZEND , IZDIM  )
            NGDAT=IXDIM*IYDIM*IZDIM
          ENDIF
          CALL GHPGET( HHEAD, 'TIME' , ITIME(IANUM) )
          CALL GHCGET( HHEAD, 'DATE' , HDATE(IANUM) )
          CALL GHPGET( HHEAD, 'DNUM' , IDNUM(IANUM) )
      GOTO 1100
        ENDIF
        CALL   GHCOPY( HHEADS,HHEAD)
*
*  < set dimension >
*
      IWDIM = IANUM
      ITDIM = IWDIM
*
      ITDIM = 2*(ITDIM/2)      !" 偶数に限る.
*
      IF (SDIM.EQ.1) THEN
        IQDIM = 1
        IRDIM = IXDIM
        ISDIM = IYDIM*IZDIM
      ELSE IF (SDIM.EQ.2) THEN
        IQDIM = IXDIM
        IRDIM = IYDIM
        ISDIM = IZDIM
      ELSE IF (SDIM.EQ.3) THEN
        IQDIM = IXDIM*IYDIM
        IRDIM = IZDIM
        ISDIM = 1
      ELSE
        CALL MSGDMP('E',' MAIN ','INVALID SPACE DIM.NUMBER' )
      ENDIF
*
*  < Fourier tranf. > 
*
      ISN = 1
*
      WRITE(6,*) IXDIM,IYDIM,IZDIM,IWDIM
      WRITE(6,*) IQDIM,IRDIM,ISDIM,ITDIM
*
**      IF (RSTRD) THEN   !" 空間リニアトレンドを引く
**        CALL SUBTRD
**     M      ( GDAT4,
**     I        IQDIM, IRDIM, IRDIM, ISDIM*ITDIM )
**      ENDIF
      CALL  AFFTRU      !" 空間 Fourier
     M      ( GDAT4, 
     W        WY,
     C        WX, IWX, 
     I        IQDIM, IRDIM, IRDIM, ISDIM*ITDIM,  ISN, LINIT )
*
**      IF (RTTRD) THEN   !" 時間リニアトレンドを引く
**        CALL SUBTRD
**     M      ( GDAT4,
**     I        IQDIM*IRDIM*ISDIM, ITDIM, ITDIM, 1 )
**      ENDIF
      CALL  AFFTRU      !" 時間 Fourier
     M      ( GDAT4, 
     W        WY,
     C        WT, IWT, 
     I        IQDIM*IRDIM*ISDIM, ITDIM, ITDIM, 1,  ISN, LINIT )
      LINIT = .FALSE.
*
*  < Filtering >
*
      CALL STFILT  !" 時空間フィルター(東進or西進のみ通す)
     M      ( GDAT4,
     I        IQDIM, IRDIM, ISDIM, ITDIM,
     I        EWD,
     M        MMIN,  MMAX,  FMIN,  FMAX,
     M        M4,    F4  )
*
*  < Fourier inverse tranf. > 
*
      ISN = -1
      CALL  AFFTRU      !" 空間逆Fourier
     M      ( GDAT4, 
     W        WY,
     C        WX, IWX, 
     I        IQDIM, IRDIM, IRDIM, ISDIM*ITDIM,  ISN, LINIT )
*
      CALL  AFFTRU      !" 時間逆Fourier
     M      ( GDAT4, 
     W        WY,
     C        WT, IWT, 
     I        IQDIM*IRDIM*ISDIM, ITDIM, ITDIM, 1,  ISN, LINIT )
*
*  < output >
*
      IF ( ITEM .NE. ' ' ) THEN
        CALL GHCSET( HHEAD , 'ITEM', ITEM )
      ENDIF
      IF ( UNIT .NE. ' ' ) THEN
        CALL GHCSET( HHEAD , 'UNIT', UNIT )
      ENDIF
      IF ( TITLE .NE. ' ' ) THEN
        CALL GHCSTS( HHEAD , 'TITL', TITLE )
      ENDIF
      IF ( DSET .NE. ' ' ) THEN
        CALL GHCSET( HHEAD , 'DSET', DSET )
      ENDIF
      IF ( GRESET ) THEN
        CALL GHRSGP( HHEAD  )
      ENDIF
*
      DO 100 ISEQ = 1,ITDIM
        CALL GHPSET( HHEAD, 'TIME' , ITIME(ISEQ) )
        CALL GHCSET( HHEAD, 'DATE' , HDATE(ISEQ) )
        CALL GHPSET( HHEAD, 'DNUM' , IDNUM(ISEQ) )
*
        CALL GFWRIT
     I         ( HHEAD , GDAT4(1+(ISEQ-1)*NGDAT) ,
     I           JFILE , 1     , 0      )
  100 CONTINUE
*
      STOP
      END
*********************************************************************
      SUBROUTINE COPY(X,Y,N)
      REAL    X(N),Y(N)
      DO 10 I=1,N
        X(I) = Y(I) 
   10 CONTINUE
      RETURN
      END
*********************************************************************
      SUBROUTINE STFILT  !" 時空間フィルター(東進or西進のみ通す)
     M      ( GDAT4,
     I        IDIM, JDIM, KDIM, LDIM,
     I        LEWD,
     M        JWMIN,JWMAX,LFMIN,LFMAX,
     M        JW4,  LF4  )
*
      REAL    GDAT4(IDIM,JDIM,KDIM,LDIM) !"JDIM:wn,LDIM:freq
      LOGICAL LEWD !" T/F: pass only eastward/westward components
      INTEGER JWMIN,JWMAX !"pass wavemumber J s.t. JWMIN<=J<=JWMAX
      INTEGER LFMIN,LFMAX !"pass freq L s.t. LFMIN<=L<=LFMAX
      INTEGER JW4(4),LF4(4) !"波数振動数の範囲を反時計回りの4辺形で設定
*
      IF(JWMIN.LE.0) THEN
        JWMIN=1         !" No low wavenum cut filtering (except for J=0)
      ENDIF
      IF(JWMAX.LT.0) THEN
        JWMAX=JDIM/2-1  !" No high wavenum cut filtering
      ENDIF
*
      IF(LFMIN.LE.0) THEN
        LFMIN=1         !" No low freq cut filtering (except for L=0)
      ENDIF
      IF(LFMAX.LT.0) THEN
        LFMAX=LDIM/2-1  !" No high freq cut filtering
      ENDIF
*
      DO 1000 L=1,LDIM
      DO 1000 I=1,IDIM  !" IDIM=1がほとんどだろう
      DO 1000 K=1,KDIM
      DO 1000 J=1,MIN(2*JWMIN,JDIM) !" WN=0,..,JWMIN-1, and JDIM/2  -> ZERO
        GDAT4(I,J,K,L)=0
 1000 CONTINUE
*
      DO 1100 L=1,LDIM
      DO 1100 I=1,IDIM
      DO 1100 K=1,KDIM
      DO 1100 J=MAX(2*JWMAX+3,3),JDIM !" WN=JWMAX+1,..,JDIM/2-1  -> ZERO
          GDAT4(I,J,K,L)=0
 1100 CONTINUE
*
      DO 1200 L=1,MIN(2*LFMIN,LDIM) !" FREQ=0,..,LFMIN-1 and LDIM/2 -> ZERO
      DO 1200 I=1,IDIM
      DO 1200 K=1,KDIM
      DO 1200 J=1,JDIM
        GDAT4(I,J,K,L)=0
 1200 CONTINUE
*
      DO 1300 L=MAX(2*LFMAX+3,3),LDIM !" FREQ=LFMAX,..,LDIM/2-1 -> ZERO
      DO 1300 I=1,IDIM
      DO 1300 K=1,KDIM
      DO 1300 J=1,JDIM
        GDAT4(I,J,K,L)=0
 1300 CONTINUE
*
      IF( JW4(1).GE.0 .AND. JW4(2).GE.0 .AND.
     &    JW4(3).GE.0 .AND. JW4(4).GE.0 .AND.
     &    LF4(1).GE.1 .AND. LF4(2).GE.1 .AND.
     &    LF4(3).GE.1 .AND. LF4(4).GE.1      ) THEN !" 全部正の場合
        DO 1500 L=3,LDIM
        DO 1500 J=3,JDIM
          LW=(L-1)/2
          JW=(J-1)/2
          IF( (JW4(2)-JW4(1))*(LW-LF4(1)).LT.(LF4(2)-LF4(1))*(JW-JW4(1))
     &    .OR.(JW4(3)-JW4(2))*(LW-LF4(2)).LT.(LF4(3)-LF4(2))*(JW-JW4(2))
     &    .OR.(JW4(4)-JW4(3))*(LW-LF4(3)).LT.(LF4(4)-LF4(3))*(JW-JW4(3))
     &    .OR.(JW4(1)-JW4(4))*(LW-LF4(4)).LT.(LF4(1)-LF4(4))*(JW-JW4(4))
     &      ) THEN     !" (J,L)平面で反時計回りの4辺形の外にある条件
            DO 1600 I=1,IDIM
            DO 1600 K=1,KDIM
              GDAT4(I,J,K,L)=0
 1600       CONTINUE
          ENDIF
 1500   CONTINUE
      ENDIF
*
      IF(LEWD) THEN
        DO 2000 LW=1,LDIM/2-1 !" =FREQ; COS:L=2*LW+1 , -SIN:L=2*LN+2
        DO 2000 K=1,KDIM
        DO 2000 JW=JWMIN,JWMAX !" =WN;   COS:J=2*JW+1 , -SIN:J=2*JN+2
        DO 2000 I=1,IDIM
          GDAT4(I,2*JW+1,K,2*LW+1) = ( GDAT4(I,2*JW+1,K,2*LW+1)
     &                                +GDAT4(I,2*JW+2,K,2*LW+2) )/2
          GDAT4(I,2*JW+1,K,2*LW+2) = ( GDAT4(I,2*JW+1,K,2*LW+2)
     &                                -GDAT4(I,2*JW+2,K,2*LW+1) )/2
          GDAT4(I,2*JW+2,K,2*LW+1) = -GDAT4(I,2*JW+1,K,2*LW+2)
          GDAT4(I,2*JW+2,K,2*LW+2) =  GDAT4(I,2*JW+1,K,2*LW+1)
 2000   CONTINUE
      ELSE
        DO 3000 LW=1,LDIM/2-1 !" =FREQ; COS:L=2*LW+1 , -SIN:L=2*LN+2
        DO 3000 K=1,KDIM
        DO 3000 JW=JWMIN,JWMAX !" =WN; COS:J=2*JW+1 , -SIN:J=2*JN+2
        DO 3000 I=1,IDIM
          GDAT4(I,2*JW+1,K,2*LW+1) = ( GDAT4(I,2*JW+1,K,2*LW+1)
     &                                -GDAT4(I,2*JW+2,K,2*LW+2) )/2
          GDAT4(I,2*JW+1,K,2*LW+2) = ( GDAT4(I,2*JW+1,K,2*LW+2)
     &                                +GDAT4(I,2*JW+2,K,2*LW+1) )/2
          GDAT4(I,2*JW+2,K,2*LW+1) =  GDAT4(I,2*JW+1,K,2*LW+2)
          GDAT4(I,2*JW+2,K,2*LW+2) = -GDAT4(I,2*JW+1,K,2*LW+1)
 3000   CONTINUE
      ENDIF
*
      RETURN
      END
