*" PACKAGE GTSTSPCT 時空間パワースペクトル、位相
*" [HIS]  93/11/25 堀之内 第一版（未完成部分多々あり）
*"        93/12/08 堀之内 今のところ"空間"軸(SDIM)が x方向であると仮定
*"                        して軸ファイルがつくられるようになっている。
*"                        オプション -TIMSQ は正常に働かない。
*"                     [未確認]：SDIM≠1 のときの動作
***********************************************************************
      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 = 4.5E6 )   !" 読み込む全データ数の最大値
      PARAMETER (IJKDIM= 1E5  )   !" 1レコードのデータ数の最大値
**      PARAMETER (NDIM4 = 7536700 )   !" 読み込む全データ数の最大値
**      PARAMETER (IJKDIM= 94300  )   !" 1レコードのデータ数の最大値
*
      REAL       GDAT4( NDIM4 )
*
      CHARACTER  HHEAD ( NDC )*(NCC)
      CHARACTER  HHEADS( NDC )*(NCC)
      CHARACTER  HHEADW( NDC )*(NCC)
      REAL       GDATA ( IJKDIM )
*
      CHARACTER  HFILE( 1 ) *(NFILN)
      DATA       HFILE / '$GTTMPDIR/gtool.out' /
      DATA       IFILE / 50 /
      DATA       JFILE / 60 /
      DATA       JFILEP/ 61 /
*
      CHARACTER  OUT( 1 )    *(NFILN)
      DATA       OUT / '$GTTMPDIR/gtool.out' /
      CHARACTER  OUTP( 1 )    *(NFILN)
      DATA       OUTP/'$GTTMPDIR/gtoolph.out' /
      LOGICAL    APND
      DATA       APND   / .FALSE. /
*
      CHARACTER  ITEM   *(NCC)
      CHARACTER  UNIT   *(NCC)
      CHARACTER  TITLE  *(NCC*2)
      CHARACTER  DSET   *(NCC)
      DATA       ITEM, UNIT, TITLE, DSET /4*' '/
      LOGICAL    GRESET
      DATA       GRESET / .TRUE. /
      LOGICAL    HELP
      DATA       HELP   / .FALSE. /
      LOGICAL    NOPH              !" T -> 位相を出力しない
      DATA       NOPH   / .FALSE. /
      INTEGER    SDIM       !" Fourier変換される空間次元は何次元目か
      DATA       SDIM   / 1 /
      LOGICAL    TIMSQ   !" T -> 入力データは既に時系列化されている
*"    但し、時系列データは、１レコードの最終次元が時間軸であると仮定する
      DATA       TIMSQ / .FALSE. /
      INTEGER    SEQ  !" 元が空間３次元データの場合、できた4次元データを
*                     !" 第 SEQ 次元以外の3次元データの列として出力する
      DATA       SEQ   / 4 /  !" デフォルトは、周波数毎のデータの列
      LOGICAL    RSTRD  !"  T -> 空間フーリエの前にリニアトレンドを引く
      DATA       RSTRD  / .FALSE. /
      LOGICAL    RTTRD  !"  T -> 時間フーリエの前にリニアトレンドを引く
      DATA       RTTRD  / .FALSE. /
      LOGICAL    LOG  !" T-> スペクトルのスケーリングタイプを log にする
      DATA       LOG    / .FALSE. /
      LOGICAL    NOAX   !"  T-> 軸ファイルを書かない
      DATA       NOAX   / .FALSE. /
      CHARACTER  UTIM   *4   !" 共通パラメター 'UTIM': 変えない方が無難
      DATA       UTIM   / 'DAY ' /
*
      NAMELIST  /OPTION/ OUT,  OUTP, APND,
     &                   ITEM, UNIT,  TITLE, DSET, GRESET,
     &                   HELP, HFILE, NOPH,  SDIM, TIMSQ,
     &                   SEQ,  RSTRD, RTTRD, LOG,  NOAX, UTIM
*
      CHARACTER  EDIT   *(NCC)
      CHARACTER  ETTL   *(NCC)
      DATA       EDIT   / 'STPWSP' /
      DATA       ETTL   / 's-t power spectr' /
*
      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' /  !" その頭の部分
      INTEGER    ITIMD   !" 入力の何次元目が時間軸か。
      DATA       ITIMD  / 4 / 
      INTEGER    ISN
      DATA       ISN    / 1 /
      CHARACTER  HFTIME  *(NCC)
      CHARACTER  HFTIM1  *(NCC)
      LOGICAL    LFIRST
      DATA       LFIRST  /.TRUE./
      LOGICAL    LINIT
      DATA       LINIT  /.TRUE./
      CHARACTER  HAITM1*(NCC),HAITM2*(NCC),HAITM3*(NCC)
      CHARACTER  HAXSEQ*(NCC)
      CHARACTER  HAXSDM*(NCC)
      INTEGER    IXDIM,IYDIM,IZDIM,IWDIM !" ４次元データの各要素数
      INTEGER    IRDIM,ITDIM     !" 変換に使われる空間、時間の要素数
      INTEGER    IQDIM,ISDIM !" 時系列は(IQ,IR,IS,IT)DIM の４次元として扱われる
      INTEGER    ITIMDO    !"   出力の何次元目が時間軸か
      LOGICAL    LSEQ   !"  T -> 入力が空間３次元データの(時系)列である。
      DATA       LSEQ  /.FALSE./
*
*  [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 ( IOS.NE.0 .OR. HELP ) THEN
         WRITE(6,OPTION)
         STOP
      ENDIF
*
      CALL GTOPEN
      CALL GTSIZE ( HHEAD, IJKDIM )
*
      CALL GURNTF ( HFILE( 1 ), OUT(1)  , '$GTTMPDIR/gtool.in' )
      CALL GURNTF ( HFILE( 1 ), OUTP(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)
      IF(.NOT.NOPH) THEN
        CALL GFOOPN ( JFILEP,  OUTP(1) , APND )
        CALL GUNENV( OUTP(1),'.',.FALSE. )
        IL=LENC(OUTP(1))
        WRITE (6,*) 'output='//OUTP(1)(1:IL)
      ENDIF
*
      CALL GTCSET ( 'UTIM', UTIM )
*
*  < read data >
*
      IANUM = 0
 1100 CONTINUE
        CALL   GHCOPY( HHEAD,HHEADS)
        CALL   GFREAD
     O        ( HHEAD , GDATA , IEOD  ,
     I          IFILE , 1               )
*
        IF ( IEOD .EQ.0 ) THEN
          IANUM = IANUM + 1
          IF (LFIRST) THEN
            LFIRST = .FALSE.
            CALL GHCGET( HHEAD, 'TIME' , HFTIM1 )
            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  )
            IF (TIMSQ) THEN
              IF  ( HAITM2 .EQ. ' ' )
     &            CALL MSGDMP( 'E',' MAIN ','ZERO SPACE DIMENSION' )
              IF  ( HAITM3 .EQ. ' ' ) THEN
                ITIMD = 2
                ITIMDO= 2
                ITDIM = IYDIM
              ELSE
                ITIMD = 3
                ITIMDO= 3
                ITDIM = IZDIM
              ENDIF
              IF (.NOT.(SDIM.LT.ITIMDO))
     &            CALL MSGDMP( 'E',' MAIN ','INVALID SPACE DIM.NUM.' )
            ELSE
              IF(HAITM2.EQ.' ') THEN
                ITIMDO = 2
              ELSE IF(HAITM3.EQ.' ') THEN
                ITIMDO = 3
              ELSE 
                LSEQ = .TRUE.
                IF (SEQ.EQ.4) THEN
                  ITIMDO = 4
                ELSE
                  ITIMDO = 3
                ENDIF
              ENDIF
            ENDIF
          ENDIF
          IF ( IANUM*IXDIM*IYDIM*IZDIM .GT. NDIM4 )
     &         CALL MSGDMP('E',' MAIN ','TOO SMALL "NDIM4"')
          CALL COPY(GDAT4(1 + (IANUM-1)*IXDIM*IYDIM*IZDIM ),
     &              GDATA, IXDIM*IYDIM*IZDIM )
          IF (MOD(IANUM,2).EQ.0) CALL GHCGET( HHEAD, 'TIME' , HFTIME )
      GOTO 1100
        ENDIF
        CALL   GHCOPY( HHEADS,HHEAD)
        CALL   GHCOPY( HHEADS,HHEADW)
*
*  < set dimension >
*
      IWDIM = IANUM
      IF(.NOT.TIMSQ) THEN
        ITDIM = IWDIM
      ENDIF
*
      ITDIM = 2*(ITDIM/2)      !" 偶数に限る.
      IF ( .NOT.TIMSQ  ) CALL GHCSET(HHEAD, 'TIME', HFTIME)
*
      IF (SDIM.EQ.1) THEN
        IQDIM = 1
        IRDIM = IXDIM
        IF (ITIMD.EQ.2) THEN
          ISDIM = 1
        ELSE IF (ITIMD.EQ.3) THEN
          ISDIM = IYDIM
        ELSE
          ISDIM = IYDIM*IZDIM
        ENDIF
      ELSE IF (SDIM.EQ.2) THEN
        IQDIM = IXDIM
        IRDIM = IYDIM
        IF (ITIMD.EQ.3) THEN
          ISDIM = 1
        ELSE
          ISDIM = IZDIM
        ENDIF
      ELSE IF (SDIM.EQ.3) THEN
        IQDIM = IXDIM*IYDIM
        IRDIM = IZDIM
        ISDIM = 1
      ELSE
        CALL MSGDMP('E',' MAIN ','INVALID SPACE DIM.NUMBER' )
      ENDIF
*
*  < Fourier trans,  spectrum,phase > 
*
      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.
*
      CALL  FST2SP  !" space-time fourier --> space-time power spectrum
     O   ( SP    ,
     I     GDAT4 , IQDIM, IRDIM, ISDIM, ITDIM )
*
      IF (.NOT.NOPH) THEN
        IF ( IQDIM*(IRDIM+1)*ISDIM*ITDIM/2 .GT. NDIM4/2 )
     &       CALL MSGDMP('E',' MAIN ','TOO Small "NDIM4"')    
        CALL  FST2PH !" space-time fourier --> space-time fourier phase
     O   ( PH    ,
     I     GDAT4 , IQDIM, IRDIM, ISDIM, ITDIM )
      ENDIF
*
*  < header set >
*
      IF ( ITEM .NE. ' ' ) THEN
        CALL GHCSET( HHEAD , 'ITEM', ITEM )
      ENDIF
      IF ( UNIT .NE. ' ' ) THEN
        CALL GHCSET( HHEAD , 'UNIT', UNIT )
      ELSE
        CALL GHCGET( HHEAD , 'UNIT', UNIT )
        IL = LENZ(UNIT)
        IF (IL.NE.0) THEN
          UNIT = '('//UNIT(1:IL)//')|2"'
          CALL GHCSET( HHEAD , 'UNIT', UNIT )
        ENDIF
      ENDIF

      IF ( TITLE .NE. ' ' ) THEN
        CALL GHCSTS( HHEAD , 'TITL', TITLE )
      ELSE
        CALL GHCGTS( HHEAD , 'TITL', TITLE )
        IL = LENC(TITLE)
        TITLE = TITLE(1:IL)//' s-t pw.sp'
        CALL GHCSTS( HHEAD , 'TITL', TITLE )
      ENDIF
      IF ( DSET .NE. ' ' ) THEN
        CALL GHCSET( HHEAD , 'DSET', DSET )
      ENDIF
      IF ( GRESET ) THEN
        CALL GHRSGP( HHEAD  )
      ENDIF
*
      CALL GHEADD( HHEAD , EDIT , ETTL )
*
      IL = LENC(AITMS)
      CALL CHVAL( '(I8)', REAL(IRDIM/2) , AITMS(IL+1:NCC) )
*
      IF (SDIM.EQ.1) THEN
        CALL GHCGET( HHEAD, 'AITM1', HAXSDM) 
        IF (HAXSDM(1:4).EQ.'CSIG') THEN
          AITMS(5:6) = 'VT'
          HAXSDM(1:4) = 'CLNS'
          CALL GHCSET( HHEADW, 'AITM1', HAXSDM) 
        ENDIF
        CALL GHCSET( HHEAD, 'AITM1', AITMS) 
        CALL GHPSET( HHEAD, 'ASTR1', 1    ) 
        CALL GHPSET( HHEAD, 'AEND1', IRDIM+1) 
      ELSE IF (SDIM.EQ.2) THEN
        CALL GHCGET( HHEAD, 'AITM2', HAXSDM) 
        IF (HAXSDM(1:4).EQ.'CSIG') THEN
          AITMS(5:6) = 'VT'
          HAXSDM(1:4) = 'CLNS'
          CALL GHCSET( HHEADW, 'AITM2', HAXSDM) 
        ENDIF
        CALL GHCSET( HHEAD, 'AITM2', AITMS) 
        CALL GHPSET( HHEAD, 'ASTR2', 1    ) 
        CALL GHPSET( HHEAD, 'AEND2', IRDIM+1) 
      ELSE IF (SDIM.EQ.3) THEN
        CALL GHCGET( HHEAD, 'AITM2', HAXSDM) 
        IF (HAXSDM(1:4).EQ.'CSIG') THEN
          AITMS(5:6) = 'VT'
          HAXSDM(1:4) = 'CLNS'
          CALL GHCSET( HHEADW, 'AITM3', HAXSDM) 
        ENDIF
        CALL GHCSET( HHEAD, 'AITM3', AITMS) 
        CALL GHPSET( HHEAD, 'ASTR3', 1    ) 
        CALL GHPSET( HHEAD, 'AEND3', IRDIM+1) 
      ENDIF
*
      IL = LENC(AITMF)
      CALL CHVAL( '(I8)', REAL(ITDIM/2) , AITMF(IL+1:NCC) )
      CALL FFTTIM
     I           ( TIMSQ,  ITDIM, 
     I             HHEAD,  HFTIM1,
     I             ITIMDO,
     O             FTIM )
      IL = LENC(AITMF)
      AITMF(IL+1:IL+1) = '-'
      CALL CHVAL( '(I8)', REAL(FTIM) , AITMF(IL+2:NCC) )
*
      CALL HEDAXF
     M           ( HHEAD,
     O             HAXSEQ,
     I             AITMF, ITDIM, ITIMDO, SEQ ,LSEQ )
*
      CALL   GHCOPY( HHEAD,HHEADS)
      IF ( TITLE .NE. ' ' ) THEN
        CALL GHCSTS( HHEADS , 'TITL', TITLE )
      ELSE
        CALL GHCGTS( HHEADS , 'TITL', TITLE )
        IL = LENC(TITLE)
        TITLE = TITLE(1:IL)//' s-t phase'
        CALL GHCSTS( HHEADS , 'TITL', TITLE )
      ENDIF
      CALL GHEADD( HHEADS , EDITP , ETTLP )
      CALL GHCSET( HHEADS ,'UNIT', 'RAD')
*
      IF (LOG) THEN
        CALL GHPSET( HHEAD , 'STYP' , 2 )
      ENDIF
*
*  < output >
*
      IF (LSEQ) THEN
        IXDIMO = IXDIM
        IYDIMO = IYDIM
        IZDIMO = IZDIM
        IF(SDIM.EQ.1) IXDIMO = IXDIM + 1
        IF(SDIM.EQ.2) IYDIMO = IYDIM + 1
        IF(SDIM.EQ.3) IZDIMO = IZDIM + 1
        IF(SEQ.EQ.1) NSEQ=IXDIMO
        IF(SEQ.EQ.2) NSEQ=IYDIMO
        IF(SEQ.EQ.3) NSEQ=IZDIMO
        IF(SEQ.EQ.4) NSEQ=IWDIM/2
        DO 100 ISEQ = 1,NSEQ
          CALL COPY43(GDATA,SP,SEQ,ISEQ,IXDIMO,IYDIMO,IZDIMO,IWDIM/2)
          IF (SEQ.EQ.4) THEN
            CALL EADDPR
     M           ( HHEAD ,
     I             ISEQ, FTIM, AITMF )
          ELSE IF (SEQ.EQ.SDIM) THEN
            CALL EADDSX
     M           ( HHEAD,
     I             ISEQ, IRDIM  , AITMS )
          ELSE
            CALL EADDSQ
     M           ( HHEAD,
     I             HAXSEQ, SEQ, ISEQ ) 
          ENDIF
          CALL GFWRIT
     I         ( HHEAD , GDATA ,
     I         JFILE , 1     , 0      )
          IF (.NOT.NOPH) THEN
            CALL COPY43(GDATA,PH,SEQ,ISEQ,IXDIMO,IYDIMO,IZDIMO,IWDIM/2)
            IF (SEQ.EQ.4) THEN
              CALL EADDPR
     M             ( HHEADS ,
     I               ISEQ, FTIM, AITMF )
            ELSE IF (SEQ.EQ.SDIM) THEN
              CALL EADDSX
     M             ( HHEADS,
     I               ISEQ, IRDIM  , AITMS )
            ELSE
              CALL EADDSQ
     M             ( HHEADS,
     I               HAXSEQ, SEQ, ISEQ ) 
            ENDIF
            CALL GFWRIT
     I           ( HHEADS , GDATA ,
     I           JFILEP, 1     , 0      )
          ENDIF
  100   CONTINUE
      ELSE
        CALL COPY( GDAT4,SP,IQDIM*(IRDIM+1)*ISDIM*(ITDIM/2) )
        CALL GFWRIT
     I       ( HHEAD , GDAT4 ,
     I       JFILE , 1     , 0      )
        IF (.NOT.NOPH) THEN
          CALL COPY( GDAT4,PH,IQDIM*(IRDIM+1)*ISDIM*(ITDIM/2) )
          CALL GFWRIT
     I         ( HHEADS , GDAT4 ,
     I         JFILEP, 1     , 0      )
        ENDIF
      ENDIF
*
*  < make axis file >
*
      IF (.NOT.NOAX) THEN
        CALL JIKFRQ
     I      ( HHEAD,  AITMF, FTIM,
     I        ITDIM/2 )
        IF (HAXSDM(1:4).EQ.'GLON') THEN
          CALL JIKSPC
     I         ( HHEAD,  AITMS,
     I           IRDIM,  SDIM )
        ELSE IF(HAXSDM(1:4).EQ.'CLNS') THEN
          
          CALL JIKVRT
     I         ( HHEAD, HHEADW,  HAXSDM,
     I           IRDIM,  SDIM )
        ENDIF
      ENDIF
*
      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
*********************************************************************
*  copy 4dim array --> 3dim array
      SUBROUTINE COPY43(X,Y,KD,ID,N1,N2,N3,N4)
      REAL    X(*),Y(*)
      IF (KD.EQ.1) THEN
        CALL COP431(X,Y,ID,N1,N2,N3,N4)
      ELSE IF (KD.EQ.2) THEN
        CALL COP432(X,Y,ID,N1,N2,N3,N4)
      ELSE IF (KD.EQ.3) THEN
        CALL COP433(X,Y,ID,N1,N2,N3,N4)
      ELSE IF (KD.EQ.4) THEN
        CALL COP434(X,Y,ID,N1,N2,N3,N4)
      ELSE
        CALL MSGDMP('E','COPY43','INVALID "KD"')
      ENDIF
      RETURN
      END
*********************************************************************
*  copy 4dim array --> 3dim array
      SUBROUTINE COP431(X,Y,ID,N1,N2,N3,N4)
      REAL    X(N2,N3,N4),Y(N1,N2,N3,N4)
      DO 10 I4 = 1,N4
      DO 10 I3 = 1,N3
      DO 10 I2 = 1,N2
        X(I2,I3,I4) = Y(ID,I2,I3,I4)
   10 CONTINUE
      RETURN
      END
*********************************************************************
*  copy 4dim array --> 3dim array
      SUBROUTINE COP432(X,Y,ID,N1,N2,N3,N4)
      REAL    X(N1,N3,N4),Y(N1,N2,N3,N4)
      DO 10 I4 = 1,N4
      DO 10 I3 = 1,N3
      DO 10 I1 = 1,N1
        X(I1,I3,I4) = Y(I1,ID,I3,I4)
   10 CONTINUE
      RETURN
      END
*********************************************************************
*  copy 4dim array --> 3dim array
      SUBROUTINE COP433(X,Y,ID,N1,N2,N3,N4)
      REAL    X(N1,N2,N4),Y(N1,N2,N3,N4)
      DO 10 I4 = 1,N4
      DO 10 I2 = 1,N2
      DO 10 I1 = 1,N1
        X(I1,I2,I4) = Y(I1,I2,ID,I4)
   10 CONTINUE
      RETURN
      END
*********************************************************************
*  copy 4dim array --> 3dim array
      SUBROUTINE COP434(X,Y,ID,N1,N2,N3,N4)
      REAL    X(N1,N2,N3),Y(N1,N2,N3,N4)
      DO 10 I3 = 1,N3
      DO 10 I2 = 1,N2
      DO 10 I1 = 1,N1
        X(I1,I2,I3) = Y(I1,I2,I3,ID)
   10 CONTINUE
      RETURN
      END
*********************************************************************
      SUBROUTINE JIKFRQ
     I      ( HHEAD,  AITMF,  FTIM,
     I        NF )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#endif
      CHARACTER  HHEAD ( NDC )*(NCC)
      CHARACTER  AITMF   *(NCC)
      CHARACTER  HUTIM   *4  !"共通パラメタ UTIM (default:'DAY')
      CHARACTER  HUNIT   *6
      PARAMETER  ( NAX=1024 )
      REAL       AX  ( NAX )
      REAL       DAX ( NAX )
      CHARACTER  HAX (NDC)*(NCC)
*
      CALL GTCGET ( 'UTIM', HUTIM )
      IL = LENC(HUTIM)
      HUNIT = '1/'//HUTIM(1:IL)
      CALL GHPGET ( HHEAD, 'MISS', VMISS )
*
      DF = 1. / FTIM
      DO 10 I = 1,NF
        AX(I) = DF*I
        DAX(I) = 1.
   10 CONTINUE
*
      CALL GHCSET( HAX, 'AITM1', AITMF )
      CALL GHPSET( HAX, 'AEND1', NF )
*
      CALL    GHPACA
     O         ( HAX,
     I           'IAXLOC', AITMF,
     I           'frequency', HUNIT ,
     I           NF,
     I           'UR4' , VMISS ,
     I           VMISS , VMISS , VMISS , VMISS , 1 )
*
      CALL    GUSIAX
     I         ( HAX , AX , DAX ,
     I           1                )
*
      CALL  GUWIAX
     I         ( HAX )
*
      RETURN
      END
*********************************************************************
      SUBROUTINE JIKSPC
     I      ( HHEAD,  AITMS,
     I        NS , IAX )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#endif
      CHARACTER  HHEAD ( NDC )*(NCC)
      CHARACTER  AITMS   *(NCC)
      CHARACTER  HAITM   *(NCC)
      PARAMETER  ( NAX=1024 )
      REAL       AX  ( NAX )
      REAL       DAX ( NAX )
      CHARACTER  HAX (NDC)*(NCC)
*
      CALL GHPGET ( HHEAD, 'MISS', VMISS )
      CALL GHPSET( HAX, 'AEND1', NS+1 )
*
      IF (IAX.EQ.1) THEN
        CALL GHCGET( HHEAD, 'AITM1', HAITM )
      ELSE IF (IAX.EQ.2) THEN
        CALL GHCGET( HHEAD, 'AITM2', HAITM )
      ELSE IF (IAX.EQ.3) THEN
        CALL GHCGET( HHEAD, 'AITM3', HAITM )
      ENDIF
*
        DO 10 I = -NS/2,NS/2
          AX(I+NS/2+1) = 1.*I
          DAX(I+NS/2+1) = 1.
   10   CONTINUE
*
        CALL    GHPACA
     O         ( HAX,
     I           'IAXLOC', AITMS,
     I           'wave number', ' ',
     I           NS+1,
     I           'UR4' , VMISS ,
     I           VMISS , VMISS , VMISS , VMISS , 1 )
*
        CALL    GUSIAX
     I         ( HAX , AX , DAX ,
     I           1                )
*
        CALL  GUWIAX
     I         ( HAX )
*
      RETURN
      END
*********************************************************************
      SUBROUTINE JIKVRT
     I      ( HHEAD,  HHEADW, HAXSDM,
     I        NS , IAX )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#endif
      CHARACTER  HHEAD ( NDC )*(NCC)
      CHARACTER  HAITM   *(NCC)
      CHARACTER  HAXSDM   *(NCC)
      PARAMETER  ( NAX=1024 )
      REAL       WGT ( NAX )
      REAL       AX  ( NAX )
      REAL       DAX ( NAX )
      CHARACTER  HAX (NDC)*(NCC)
      CHARACTER  HHEADO( NDC )*(NCC)
      LOGICAL    LREQ
*
      CALL GHPGET ( HHEAD, 'MISS', VMISS )
      CALL GHPSET( HAX, 'AEND1', NS+1 )
*
      IF (IAX.EQ.1) THEN
        CALL GHCGET( HHEAD, 'AITM1', HAITM )
      ELSE IF (IAX.EQ.2) THEN
        CALL GHCGET( HHEAD, 'AITM2', HAITM )
      ELSE IF (IAX.EQ.3) THEN
        CALL GHCGET( HHEAD, 'AITM3', HAITM )
      ENDIF
*
      CALL GTSIZE ( HHEADO  ,  NAX  )
      CALL GUQAXV
     I     ( HHEADW, IAX, 'WGT' ,
     O       HHEADO,WGT, IEOD  )
*
      W = WGT(1)
      S = WGT(1)
      DO 5 I = 2,NS
        W2 = WGT(I)
        IF(.NOT.LREQ(W,W2)) 
     &       CALL MSGDMP('W','JIKVRT','AXIS IN NOT EQUALLY DIVDED')
        W = W2
        S = S + WGT(I)
    5 CONTINUE
*
      DS = -1 / S
*
      DO 10 I = -NS/2,NS/2
        AX(I+NS/2+1) = DS*I
        DAX(I+NS/2+1) = 1.
   10 CONTINUE
*
      CALL    GHPACA
     O         ( HAX,
     I           'IAXLOC', HAITM,
     I           'vertical wn', ' ',
     I           NS+1,
     I           'UR4' , VMISS ,
     I           VMISS , VMISS , VMISS , VMISS , 1 )
*
      CALL    GUSIAX
     I         ( HAX , AX , DAX ,
     I           1                )
*
      CALL  GUWIAX
     I         ( HAX )
*
      RETURN
      END
***********************************************************************
* GET TIME FOR TIME FFT
*
      SUBROUTINE FFTTIM
     I           ( LTIMSQ, ITDIM, 
     I             HHEAD,  HFTIM1,
     I             ITIMD,
     O             FTIM )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#endif
*  [INPUT]
      LOGICAL    LTIMSQ   !"  T -> 入力データは既に時系列化されている。
      INTEGER    ITDIM    !" 変換に使われる時間の要素数
      CHARACTER  HHEAD ( NDC )*(NCC) !"          (for .NOT.LTMSQ)
      CHARACTER  HFTIM1  *(NCC)  !" 始めの時刻  (for .NOT.LTMSQ)
      INTEGER    ITIMD    !" 時間軸は何次元目か (for LTIMSQ)
*  [OUTPUT]
      REAL       FTIM     !" フーリエ変換に使われた時間 (単位 HUTIM)
*  [WORK]
      CHARACTER  HUTIM   *4  !"共通パラメタ UTIM (default:'DAY') となる
      CHARACTER  HUTIMF  *4  !" (for LTIMSQ)
      CHARACTER  HAW     *5  !" (for LTIMSQ)
      CHARACTER  HHEADW( NDC )*(NCC)
      PARAMETER  ( MAXAX = 1024 )
      REAL       AXIS(MAXAX)
*
      IF (.NOT.LTIMSQ) THEN
        CALL GUQTIM ( HHEAD , RTF , HUTIM ) !" 'DAY'?単位で時刻を読む
        CALL GHCOPY( HHEAD,HHEADW)
        CALL GHCSET(HHEADW, 'TIME', HFTIM1)
        CALL GUQTIM ( HHEADW , RTI , HUTIM ) !" 'DAY'?単位で時刻を読む
        FTIM = ( RTF - RTI ) * ITDIM / (ITDIM-1)
      ELSE
        IF (ITDIM.GT.MAXAX)
     &      CALL MSGDMP('E','FFTTIM','TOO LARGE TIME SEQUENCE')
        CALL GUQAXV
     I       ( HHEAD, ITIMD, 'LOC' ,
     O         HHEADW,AXIS , IEOD )
        HAW = 'ASTR#'
        WRITE(HAW(5:5),'(I1)') ITIMD
        CALL GHPGET(HHEAD,HAW,IASTR)
        RTI = AXIS(IASTR)
        RTF = AXIS(IASTR+ITDIM)
        FTIM = ( RTF - RTI ) * (ITDIM+1) / ITDIM
        CALL GHCGET ( HHEADW, 'UNIT',  HUTIMF )
        CALL GTCGET ( 'UTIM', HUTIM )
        IF ( HUTIMF(1:1) .NE. HUTIM(1:1)  ) THEN
          IT = FTIM
          CALL GUCT2R
     I            ( IT , HUTIMF, HUTIM , 1      ,
     O              FTIM                                )
        ENDIF
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE HEDAXF
     M           ( HHEAD,
     O             HAXSEQ,
     I             AITMF, ITDIM, ITIMDO, ISEQ ,LSEQ )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#endif
*  [MODIFY]
      CHARACTER  HHEAD ( NDC )*(NCC)
*  [OUTPUT]
      CHARACTER  HAXSEQ *(NCC)
*  [INPUT]
      CHARACTER  AITMF  *(NCC)  !" 時間軸ファイル名
      INTEGER    ITIMDO    !"   出力の何次元目が時間軸か
      INTEGER    ISEQ  !" 元が空間３次元データの場合、できた4次元データを
*                     !" 第 ISEQ 次元以外の3次元データの列として出力する。
      LOGICAL    LSEQ   !"  T -> 入力が空間３次元データの(時系)列である。
*  [WORK]
      CHARACTER*(NCC)  CAITM,CASTR,CAEND
      CHARACTER*(NCC)  HAITM,HASTR,HAEND
*
      CAITM = 'AITM#'
      CASTR = 'ASTR#'
      CAEND = 'AEND#'
*
      IF (ISEQ.NE.4) THEN
        WRITE(CAITM(5:5),'(I1)') ISEQ
        CALL GHCGET( HHEAD, CAITM, HAXSEQ )
      ENDIF
*
      IF (LSEQ .AND. ISEQ.LT.3) THEN
        DO 10 KS = ISEQ+1,3
          WRITE(CAITM(5:5),'(I1)') KS
          WRITE(CASTR(5:5),'(I1)') KS
          WRITE(CAEND(5:5),'(I1)') KS
          CALL GHCGET( HHEAD, CAITM, HAITM ) 
          CALL GHPGET( HHEAD, CASTR, HASTR ) 
          CALL GHPGET( HHEAD, CAEND, HAEND ) 
          WRITE(CAITM(5:5),'(I1)') KS-1
          WRITE(CASTR(5:5),'(I1)') KS-1
          WRITE(CAEND(5:5),'(I1)') KS-1
          CALL GHCSET( HHEAD, CAITM, HAITM ) 
          CALL GHPSET( HHEAD, CASTR, HASTR ) 
          CALL GHPSET( HHEAD, CAEND, HAEND ) 
   10   CONTINUE
      ENDIF
*
      WRITE(CAITM(5:5),'(I1)') ITIMDO
      WRITE(CASTR(5:5),'(I1)') ITIMDO
      WRITE(CAEND(5:5),'(I1)') ITIMDO
      IF (ITIMDO.NE.4) THEN
        CALL GHCSET( HHEAD, CAITM, AITMF  ) 
        CALL GHPSET( HHEAD, CASTR, 1      ) 
        CALL GHPSET( HHEAD, CAEND, ITDIM/2) 
      ENDIF
*
      RETURN
      END
*********************************************************************
*"    編集記述子追加 (周期)
      SUBROUTINE EADDPR
     M      ( HHEAD,
     I        IFRQ, FTIM, AITMF )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#endif
      CHARACTER  HHEAD ( NDC )*(NCC)
*
      INTEGER    IFRQ     !" 何番目か
      REAL       FTIM     !" フーリエ変換に使われた時間 (単位 'UTIM')
      CHARACTER  AITMF  *(NCC)  !" 時間軸ファイル名
*
      CHARACTER  EDIT   *(NCC)
      CHARACTER  ETTL   *(NCC)
      CHARACTER  HUTIM   *4  !"共通パラメタ UTIM (default:'DAY') となる
      CHARACTER  HNUM   *(NCC)
      REAL       PERIOD
      INTEGER    IENUM
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE       IENUM , OFIRST
*
      DF = 1. / FTIM
      FRQ = DF*IFRQ
      PERIOD = 1/FRQ
*
      EDIT = AITMF
*
      CALL GTCGET ( 'UTIM', HUTIM )
      CALL CHVAL( 'B' , PERIOD , HNUM )
      IL = LENC(HNUM)
      ETTL = 'PERIOD'//HNUM(1:IL)//HUTIM
*
      IF (OFIRST) THEN
        CALL GHQENM ( HHEAD, IENUM )
        OFIRST = .FALSE.
      ENDIF
      CALL GHESET( HHEAD , EDIT , ETTL , IENUM+1 )
*
      RETURN
      END
*********************************************************************
*"    編集記述子追加 (zonal w.n.)
      SUBROUTINE EADDSX
     M      ( HHEAD,
     I        ISX, IRDIM  , AITMS )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#endif
      CHARACTER  HHEAD ( NDC )*(NCC)
*
      INTEGER    ISX      !" 何番目か
      INTEGER    IRDIM    !" フーリエ変換に使われたデータ数
      CHARACTER  AITMS  *(NCC)  !" 軸ファイル名
*
      CHARACTER  EDIT   *(NCC)
      CHARACTER  ETTL   *(NCC)
      CHARACTER  HNUM   *(NCC)
      REAL       WN
      INTEGER    IENUM
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE       IENUM , OFIRST
*
      WN = -IRDIM/2 + (ISX-1)
*
      EDIT = AITMS
*
      CALL CHVAL( '(I8)' , WN , HNUM )
      ETTL = 'WVN'//HNUM
*
      IF (OFIRST) THEN
        CALL GHQENM ( HHEAD, IENUM )
        OFIRST = .FALSE.
      ENDIF
      CALL GHESET( HHEAD , EDIT , ETTL , IENUM+1 )
*
      RETURN
      END
***********************************************************************
*"    編集記述子追加
      SUBROUTINE EADDSQ
     M     ( HHEAD,
     I       HAXSEQ, KD,  IS )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)                !" NCC, NDC, NED
#else
#include         "gzsize.F"              !" NCC, NDC, NED
#endif
      CHARACTER  HHEAD ( NDC )*(NCC)
*
      CHARACTER  HAXSEQ       *(NCC)
      INTEGER    KD      !" 何次元目か
      INTEGER    IS      !" 何番目か
*
      CHARACTER  HEDIT       *(NCC)        !" 編集略記号
      CHARACTER  HETTL       *(NCC)        !" 編集タイトル
      CHARACTER  HED         *2
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      INTEGER    IENUM
      SAVE       IENUM , OFIRST
*
      IF(KD.EQ.1) THEN
        HED = 'XS'
      ELSE IF(KD.EQ.2) THEN
        HED = 'YS'
      ELSE IF(KD.EQ.3) THEN
        HED = 'ZS'
      ELSE
        HED = 'SS'
      ENDIF
*
      CALL GUEAXN
     I               ( HAXSEQ, HED    , '-select',
     I                 IS    , IS     ,
     O                 HEDIT , HETTL           )
      IF (OFIRST) THEN
        CALL GHQENM ( HHEAD, IENUM )
        OFIRST = .FALSE.
      ENDIF
      CALL GHESET( HHEAD , HEDIT , HETTL , IENUM+1 )
*
      RETURN
      END
