* PACKAGE GTXPWSPCT x 軸パワースペクトル                 A.Numaguti
*"                                    bug fix : 97/02/04 S.Takehiro
*********************************************************************
      PROGRAM GTXPSP
*
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (GZSIZE)
#else
#include         "gtsinc.F"
#include         "gzsize.F"
#endif
      COMMON     /GMWORK/ MWORK
      REAL       MWORK  ( IJKDIM )
*
      CHARACTER  HHEAD ( NDC )*(NCC)
      REAL       GDATA ( IJKDIM )
*
      CHARACTER  HFILE( 1 ) *(NFILN)
      DATA       HFILE / '$GTTMPDIR/gtool.out' /
      DATA       IFILE / 50 /
      DATA       JFILE / 60 /
*
      CHARACTER  OUT    *(NFILN)
      DATA       OUT    / '$GTTMPDIR/gtool.out' /
      LOGICAL    APND
      DATA       APND   / .FALSE. /
      CHARACTER  ITEM   *(NCC)
      CHARACTER  UNIT   *(NCC)
      CHARACTER  TITLE  *(NCC*2)
      CHARACTER  DSET   *(NCC)
      CHARACTER  EDIT   *(NCC)
      CHARACTER  ETTL   *(NCC)
      DATA       ITEM, UNIT, TITLE, DSET, EDIT, ETTL /6*' '/
      LOGICAL    GRESET
      DATA       GRESET / .TRUE. /
      LOGICAL    HELP
      DATA       HELP   / .FALSE. /
*
      EXTERNAL   XPWSPC
*     
      NAMELIST  /OPTION/ OUT, APND,
     &                   ITEM, UNIT, TITLE, DSET, EDIT, ETTL, GRESET,
     &                   HELP, HFILE
*
      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 GTSIZE ( HHEADP, IJKDIM )
      CALL GMSIZE ( IJKDIM  )
*
      CALL GURNTF ( HFILE( 1 ), OUT  , '$GTTMPDIR/gtool.in' )
*
      CALL GFROPN ( IFILE , HFILE( 1 ) )
      CALL GFOOPN ( JFILE ,  OUT , APND )
*
      CALL GUNENV( OUT,'.',.FALSE. )
      IL=LENC(OUT)
      WRITE (6,*) 'output='//OUT(1:IL)
*
 1100 CONTINUE
         CALL   GFREAD
     O        ( HHEAD , GDATA , IEOD  ,
     I          IFILE , 1               )
*
C         IF ( MAX(IEOD,IEODP) .EQ.0 ) THEN
         IF ( IEOD.EQ.0 ) THEN
            CALL GMCAL1
     I         ( XPWSPC,
     M           HHEAD , GDATA ,
     I           EDIT  , ETTL  )
*
            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
*
            CALL  GFWRIT
     I                 ( HHEAD , GDATA ,
     I                   JFILE , 1     , 0       )
*
      GOTO 1100
         ENDIF
*
      STOP
      END
********************************************************************
      SUBROUTINE XPWSPC
     I         ( HDT   , DT    ,
     O           HPW   , PW    ,
     D           IDIM  , JDIM  , KDIM  )
*
      CHARACTER  HDT ( * )*(*)
      REAL       DT  ( IDIM, JDIM, KDIM )
      CHARACTER  HPW ( * )*(*)
      REAL       PW  ( IDIM/2-1, JDIM, KDIM )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)
#else
#include         "gzsize.F"
#endif
      PARAMETER  ( NW=512 )
      REAL       WN  ( NW )
      REAL       DWN ( NW )
      CHARACTER  HWN ( NDC )*(NCC)
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
*
      IF ( IDIM .GE. NW*2 ) THEN
         CALL MSGDMP( 'E','XPWSPC','WORK AREA TOO SMALL' )
      ENDIF
*
      CALL       GTSIZE( HWN, NW )
*
      DO 1100 K = 1, KDIM
         DO 1100 J = 1, JDIM
            CALL RFFTI( IDIM, PW )
            CALL RFFTF(IDIM,DT(1,J,K),PW )
 1100 CONTINUE
*
      DO 2100 K = 1, KDIM
         DO 2100 J = 1, JDIM
            DO 2100 I = 1, IDIM/2-1
               PW( I,J,K ) = ( DT(I*2,J,K) **2 + DT(I*2+1,J,K) **2 )
     &                     / REAL(IDIM) **2
 2100 CONTINUE
*
      CALL GHCSET( HPW, 'AITM1', '@WVN01' )
      CALL GHPSET( HPW, 'AEND1', IDIM/2-1 )
*
*"       < 波数空間軸の生成 >
*
      IF( OFIRST ) THEN
         OFIRST = .FALSE.
*
         IF ( IDIM/2-1 .GT. NW ) THEN
            CALL MSGDMP('E','XPWSPC','WORK SPACE NOT ENOUGH')
         ENDIF
         DO 3100 I = 1, IDIM/2-1
            WN ( I ) = REAL(I)
            DWN( I ) = 1.
 3100    CONTINUE
*
         CALL    GHPGET ( HDT, 'MISS', VMISS )
         CALL    GHPACA
     O         ( HWN,
     I           'IAXLOC', '@WVN01',
     I           'wave number', ' ',
     I           IDIM/2-1,
     I           'UR4' , VMISS ,
     I           VMISS , VMISS , VMISS , VMISS , 1 )
*
         CALL    GUSIAX
     I         ( HWN , WN , DWN ,
     I           1                )
*
         CALL  GUWIAX
     I         ( HPW )
      ENDIF
*
      RETURN
      END
