* PACKAGE GTBPF   !" band-pass filter (M.Murakami)
**************************************************************************
      PROGRAM GTBPF
*
#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 /
*
      LOGICAL    X, Y, Z
      DATA       X, Y, Z / .FALSE.,.FALSE.,.FALSE. /
      REAL       PERI( 2 )
      DATA       PERI    / 2.1, 1000. /
      CHARACTER  OUT    *(NFILN)
      DATA       OUT    / '$GTTMPDIR/gtool.out' /
      CHARACTER  ITEM   *(NCC)
      LOGICAL    APND
      DATA       APND   / .FALSE. /
      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. /
*
      NAMELIST  /OPTION/ X, Y, Z, PERI, OUT, APND,
     &                   ITEM, UNIT, TITLE, DSET, EDIT, ETTL, 
     &                   GRESET, HELP, HFILE
*
      EXTERNAL   BNDPS1
*
      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
*
      IF       ( X .AND.(.NOT.Y).AND.(.NOT.Z) ) THEN
         CALL SETBPF( 'X', PERI(1), PERI(2) )         
      ELSE IF  ( Y .AND.(.NOT.X).AND.(.NOT.Z) ) THEN
         CALL SETBPF( 'Y', PERI(1), PERI(2) )         
      ELSE IF  ( Z .AND.(.NOT.X).AND.(.NOT.Y) ) THEN
         CALL SETBPF( 'Z', PERI(1), PERI(2) )         
      ELSE
         CALL MSGDMP( 'E','GTBPF',
     &                'ONE OF -x OR -y OR -z MUST BE SPECIFIED' )
      ENDIF
*
      CALL GTOPEN
      CALL GTSIZE ( HHEAD, 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)
*
      II = 0
 1100 CONTINUE
         CALL   GFREAD
     O        ( HHEAD , GDATA , IEOD  ,
     I          IFILE , 1               )
*
         IF ( IEOD.EQ.0 ) THEN
*
                  II = II + 1
*
                  CALL GMCAL1
     I            ( BNDPS1,
     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 BNDPS1
     I         ( HHEAD , GDATA ,
     O           HHEADF, GDATAF,
     D           IXDIM , IYDIM , IZDIM )
*
      CHARACTER  HHEAD  ( * ) *(*)
      CHARACTER  HHEADF ( * ) *(*)
      REAL       GDATA  ( * )
      REAL       GDATAF ( * )
      CHARACTER  HBPFZ*1, HBPF*1
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)
#else
#include         "gzsize.F"
#endif
      CHARACTER  HEDIT   *(NCC)
      CHARACTER  HETTL1 *(NCC)
      CHARACTER  HETTL2 *(NCC)
      SAVE
*
      IF      ( ( HBPF .EQ. 'X' ).OR.( HBPF .EQ. 'x' ) ) THEN
         CALL BNDPAS
     O      ( GDATAF,
     I        GDATA , PERI1, PERI2, 1, IXDIM, IYDIM*IZDIM )
      ELSE IF ( ( HBPF .EQ. 'Y' ).OR.( HBPF .EQ. 'y' ) ) THEN
         CALL BNDPAS
     O      ( GDATAF,
     I        GDATA , PERI1, PERI2, IXDIM, IYDIM, IZDIM )
      ELSE IF ( ( HBPF .EQ. 'Z' ).OR.( HBPF .EQ. 'z' ) ) THEN
         CALL BNDPAS
     O      ( GDATAF,
     I        GDATA , PERI1, PERI2, IXDIM*IYDIM, IZDIM, 1 )
      ELSE
         CALL MSGDMP( 'E', 'BNDPS1','INVALID BPF SPECIFICATION' )
      ENDIF
*
      HEDIT   = HBPF//'BPF'
      HETTL1  = HBPF//'cut='
      WRITE ( HETTL1(5:13),'(1PE9.2)' ) PERI1
      HETTL2 = HBPF//'cut='
      WRITE ( HETTL2(5:13),'(1PE9.2)' ) PERI2
      CALL GHEADD ( HHEADF, HEDIT, HETTL1 )
      CALL GHEADD ( HHEADF, HEDIT, HETTL2 )
*
      RETURN
*=====================================================================
      ENTRY      SETBPF
     I         ( HBPFZ , PERIZ1, PERIZ2 )
*
      HBPF  = HBPFZ 
      PERI1 = PERIZ1
      PERI2 = PERIZ2
*      
      END
********************************************************************
      SUBROUTINE BNDPAS
     O      ( FD    ,
     I        D     , PERI1 , PERI2 , 
     D        IMAX  , JMAX  , KMAX   )
*
      REAL       FD  ( IMAX, JMAX, KMAX )
      REAL       D   ( IMAX, JMAX, KMAX )
      REAL       PERI1, PERI2
*
      LOGICAL    OMISS
*
      CALL GLPGET( 'RMISS', VMISS )
*
      IF ( MIN( PERI1, PERI2 ) .LE. 2. ) THEN
         CALL MSGDMP( 'E','BNDPAS','PERIOD MUST BE > 2.0' )
      ENDIF
*
      PI = ATAN( 1.) *4.
      AV1 = 2.*PI/ PERI1
      AV2 = 2.*PI/ PERI2
*
      DOM = 2.* ABS( SIN(AV1)/(1.+COS(AV1)) - SIN(AV2)/(1.+COS(AV2)) )
      OM2 = 4.* SIN(AV1)*SIN(AV2)/(1.+COS(AV1))/(1.+COS(AV2))
*
      A0  = 2.*DOM           /( 4.+2.*DOM+OM2 )
      B1  = 2.*(OM2-4.)      /( 4.+2.*DOM+OM2 )
      B2  = ( 4.-2.*DOM+OM2 )/( 4.+2.*DOM+OM2 )
*
      DO 4000 K = 1, KMAX
*
         DO 1100 I = 1, IMAX 
            FD( I,1,K ) = 0.
            FD( I,2,K ) = 0.
 1100    CONTINUE 
*
         DO 1200 J = 3, JMAX
            DO 1210 I = 1, IMAX 
               FD( I,J,K ) = A0*( D( I,J,K ) - D( I,J-2,K ) )
     &                     - B1* FD( I,J-1,K ) - B2* FD( I,J-2,K )
 1210       CONTINUE 
 1200    CONTINUE 
*
         DO 2100 J = JMAX-2, 1, -1
            DO 2110 I = 1, IMAX
               FD( I,J,K ) = A0*( D( I,J,K ) - D( I,J+2,K ) )
     &                     - B1* FD( I,J+1,K ) - B2* FD( I,J+2,K )
 2110       CONTINUE 
 2100    CONTINUE 
*
         DO 3100 I = 1, IMAX
            OMISS = .FALSE.
            DO 3110 J = 1, JMAX
               IF ( D( I,J,K ) .EQ. VMISS ) THEN
                  OMISS = .TRUE.
               ENDIF
 3110       CONTINUE 
            IF ( OMISS ) THEN
               DO 3120 J = 1, JMAX
                  FD( I,J,K ) = VMISS
 3120          CONTINUE 
            ENDIF
 3100    CONTINUE 
*
 4000 CONTINUE 
*        
      RETURN
      END
