* 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