C-----------------------------------------------------------------------
C File: smooth.for
C Date: 09/12/1997
C-----------------------------------------------------------------------
C Version: 96B
C Enhancement: Datafield-size: 101 x 101 x 101
C-----------------------------------------------------------------------
C Tibor F. Nagy, Department of Physics, MSU
C E-mail: nagy_t@pa.msu.edu
C-----------------------------------------------------------------------
C To smooth the sharp edges and corners in the AVS/Explorer
C volume-field for better visualization.
C
C STEPS: -reading the field in
C        -smoothing
C        -rescaling 
C        -writing the field out
C-----------------------------------------------------------------------
      PROGRAM   SMOOTH
C-----------------------------------------------------------------------
      IMPLICIT  NONE
      REAL*4    W(2,0:102,0:102,0:102) ! Main field
      CHARACTER FILNAM*30              ! Filename
      INTEGER   LENG                   ! Length of the filename
      INTEGER   NX,NY,NZ               ! Dimensions of the field
      REAL*4    RMIN,RMAX    ! Min/Max value of the real physical field
      INTEGER   NSD                    ! Number of smoothing
      REAL*4    BND                    ! Boundary condition
      INTEGER   IS,ID                  ! Source and Destination index
C-----------------------------------------------------------------------
C Read the name of the inputfile
      WRITE(6,*) '--- Inputfile ---'
      CALL INQNAM (FILNAM,LENG)
C Set the source and the destination pointers
      IS = 1
      ID = 2
C Read the number of smoothing and the boundary condition
      CALL NSDBND (W,NSD,BND)
C Read the field into the source-array
      CALL RDFILE(FILNAM,LENG,W,NX,NY,NZ,IS,RMIN,RMAX)
C Smoothing the field
      CALL SMTH (W,NX,NY,NZ,NSD,IS,ID)
C Rescaling of the field
      CALL RESCALE (W,NX,NY,NZ,ID)
C Read the name of the outputfile
      WRITE(6,*) '--- Outputfile ---'
      CALL INQNAM (FILNAM,LENG)
C Write the smoothed destination-array to the ".dat"-file
      CALL WRFILE (FILNAM,LENG,W,NX,NY,NZ,ID,RMIN,RMAX)
      STOP
      END
C End of program SMOOTH
C-----------------------------------------------------------------------
      SUBROUTINE INQNAM (FILNAM,LENG)
C-----------------------------------------------------------------------
C Read a string and compute the length of the string
C-----------------------------------------------------------------------
      IMPLICIT   NONE
      CHARACTER  FILNAM*30
      INTEGER    LENG
C-----------------------------------------------------------------------
   1  WRITE(6,'(1X,A,$)') 'Enter filename. (no extension) '
      READ (5,'(A30)') FILNAM
      LENG = 30
      DO WHILE (FILNAM(LENG:LENG).EQ.' ')
        LENG = LENG - 1
      ENDDO
      IF (LENG.LT.1) GOTO 1
      RETURN
      END
C End of subroutine INQNAM
C-----------------------------------------------------------------------
      SUBROUTINE NSDBND (W,NSD,BND)
C-----------------------------------------------------------------------
C Routine the read controll-parameters
C-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL*4   W(2,0:102,0:102,0:102)
      INTEGER  NSD
      REAL*4   BND
C Local:
      INTEGER I,J,K,L
C-----------------------------------------------------------------------
      WRITE(6,'(1X,A,$)') 'Number of smoothing? (1,2,3,...) '
      READ(5,*) NSD
      WRITE(6,'(1X,A,$)') 'Boundary condition? (0-255) '
      READ(5,*) BND
C Set the boundary condition for the whole field
      DO I=1,2
        DO J=0,102
          DO K=0,102
            DO L=0,102
              W(I,J,K,L) = BND
            ENDDO
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END
C End of subroutine NSDBND
C-----------------------------------------------------------------------
      SUBROUTINE RDFILE(FILNAM,LENG,W,NX,NY,NZ,INDX,RMIN,RMAX)
C-----------------------------------------------------------------------
C Routine to read an AVS/Explorer ".dat"-file
      IMPLICIT   NONE
      CHARACTER  FILNAM*30
      INTEGER    LENG
      INTEGER    INDX
      REAL*4     W(2,0:102,0:102,0:102)
      INTEGER    NX,NY,NZ
      REAL*4     RMIN,RMAX
C Function to get one character from a file (Unix only!)
      INTEGER    FGETC
C Local:
      CHARACTER  C*1
      INTEGER    I,J,K
      INTEGER    LIN,IER
      INTEGER    IMIN,IMAX
C Dummy variables for AVS field
      INTEGER    IDUM
      REAL*4     RDUM(13)
C-----------------------------------------------------------------------
      LIN = 8
      OPEN(LIN,FILE=FILNAM(1:LENG)//'.dat',FORM='UNFORMATTED',
     #     STATUS='OLD')
C 1st: read the dimensions of the field
      IER = FGETC(LIN,C)
      NX  = ICHAR(C)
      IER = FGETC(LIN,C)
      NY  = ICHAR(C)
      IER = FGETC(LIN,C)
      NZ  = ICHAR(C)
C 2nd: read the field itself
      DO I=1,NZ
        DO J=1,NY
          DO K=1,NX
            IER = FGETC(LIN,C)
            W(INDX,K,J,I) = REAL(ICHAR(C))
          ENDDO
        ENDDO
      ENDDO
C 3rd: read the dummy and the closing parameters
      READ(LIN) IDUM,RDUM,IMIN,RMIN,IMAX,RMAX
      CLOSE(LIN)
      WRITE(6,'(1X,3(A,I4))') 'Dimensions: ',NX,' X ',NY,' X ',NZ
      WRITE(6,'(1X,A)') 'Scaling:'
      WRITE(6,'(1X,A,I4,A,F8.3)') 'Min: ',IMIN,' = ',RMIN
      WRITE(6,'(1X,A,I4,A,F8.3)') 'Max: ',IMAX,' = ',RMAX
      RETURN
      END
C End of subroutine RDFILE
C-----------------------------------------------------------------------
      SUBROUTINE SMTH (W,NX,NY,NZ,NSD,IS,ID)
C-----------------------------------------------------------------------
C Routine to smooth the field
C-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL*4   W(2,0:102,0:102,0:102)
      INTEGER  NX,NY,NZ
      INTEGER  NSD
      INTEGER  IS,ID
C Local:
      INTEGER  NS
      INTEGER  S1,S2,S4,S8
      INTEGER  I,J,K
      INTEGER  IPUF
C-----------------------------------------------------------------------
C Set the smoothing counter to zero
      NS = 0
C Main smoothing loops
 10   DO I=1,NX
        DO J=1,NY
          DO K=1,NZ
            S8 = 8.0 * W(IS,I,J,K)
            S4 =      W(IS,(I+1),J,K)
            S4 = S4 + W(IS,I,(J+1),K)
            S4 = S4 + W(IS,I,J,(K+1))
            S4 = S4 + W(IS,(I-1),J,K)
            S4 = S4 + W(IS,I,(J-1),K)
            S4 = S4 + W(IS,I,J,(K-1))
            S4 = S4 * 4.0
            S2 =      W(IS,(I+1),(J+1),K)
            S2 = S2 + W(IS,(I+1),(J-1),K)
            S2 = S2 + W(IS,(I-1),(J+1),K)
            S2 = S2 + W(IS,(I-1),(J-1),K)
            S2 = S2 + W(IS,I,(J+1),(K+1))
            S2 = S2 + W(IS,I,(J+1),(K-1))
            S2 = S2 + W(IS,I,(J-1),(K+1))
            S2 = S2 + W(IS,I,(J-1),(K-1))
            S2 = S2 + W(IS,(I+1),J,(K+1))
            S2 = S2 + W(IS,(I+1),J,(K-1))
            S2 = S2 + W(IS,(I-1),J,(K+1))
            S2 = S2 + W(IS,(I-1),J,(K-1))
            S2 = S2 * 2.0
            S1 =      W(IS,(I+1),(J+1),(K+1))
            S1 = S1 + W(IS,(I+1),(J+1),(K-1))
            S1 = S1 + W(IS,(I+1),(J-1),(K+1))
            S1 = S1 + W(IS,(I+1),(J-1),(K-1))
            S1 = S1 + W(IS,(I-1),(J+1),(K+1))
            S1 = S1 + W(IS,(I-1),(J+1),(K-1))
            S1 = S1 + W(IS,(I-1),(J-1),(K+1))
            S1 = S1 + W(IS,(I-1),(J-1),(K-1))
            W(ID,I,J,K) = (S8+S4+S2+S1) / 64.0
          ENDDO
        ENDDO
      ENDDO
      NS = NS + 1
C Return if it is enough of smoothing
      IF (NS.EQ.NSD) THEN
C Delete the outest surface for visualization
        DO I=1,NX
          DO J=1,NY
            W(ID,I,J,1)  = 0.0
            W(ID,I,J,NZ) = 0.0
          ENDDO
        ENDDO
        DO I=1,NX
          DO J=1,NZ
            W(ID,I,1,J)  = 0.0
            W(ID,I,NY,J) = 0.0
          ENDDO
        ENDDO
        DO I=1,NY
          DO J=1,NZ
            W(ID,1,I,J)  = 0.0
            W(ID,NX,I,J) = 0.0
          ENDDO
        ENDDO
        RETURN
      ELSE
C If it was not enough then change source & destination pointers
C and go back again
        IPUF = IS
        IS   = ID
        ID   = IPUF
        GOTO 10
      ENDIF
      END
C End of subroutine SMTH
C-----------------------------------------------------------------------
      SUBROUTINE RESCALE (W,NX,NY,NZ,INDX)
C-----------------------------------------------------------------------
C Rescaling the field pointed by the index
C-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL*4   W(2,0:102,0:102,0:102)
      INTEGER  NX,NY,NZ,INDX
C Local:
      INTEGER  I,J,K
      REAL*4   FMIN,FMAX
C-----------------------------------------------------------------------
C Searching for the minimal value in the field
      FMIN = 255.0
      DO I=1,NX
        DO J=1,NY
          DO K=1,NZ
            IF (W(INDX,I,J,K).LT.FMIN) FMIN=W(INDX,I,J,K)
          ENDDO
        ENDDO
      ENDDO
C Shifting the whole field to zero
      DO I=1,NX
        DO J=1,NY
          DO K=1,NZ
            W(INDX,I,J,K) = W(INDX,I,J,K) - FMIN
          ENDDO
        ENDDO
      ENDDO
C Searching for the maximal value in the shifted field
      FMAX = 0.0
      DO I=1,NX
        DO J=1,NY
          DO K=1,NZ
            IF (W(INDX,I,J,K).GT.FMAX) FMAX = W(INDX,I,J,K)
          ENDDO
        ENDDO
      ENDDO
C Multiply the whole field with the factor of (255/FMAX) for rescaling
      DO I=1,NX
        DO J=1,NY
          DO K=1,NZ
            W(INDX,I,J,K) = W(INDX,I,J,K) * 255.0 / FMAX
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END
C End of subroutine RESCALE
C-----------------------------------------------------------------------
      SUBROUTINE WRFILE (FILNAM,LENG,W,NX,NY,NZ,INDX,RMIN,RMAX)
C-----------------------------------------------------------------------
C Routine to write the field to an AVS/Explorer ".dat"-file
C-----------------------------------------------------------------------
      IMPLICIT  NONE
      REAL*4    W(2,0:102,0:102,0:102)
      INTEGER   NX,NY,NZ,INDX
      CHARACTER FILNAM*30
      INTEGER   LENG
      REAL*4    RMIN,RMAX
C Function to put one character to a file (Unix only!)
      INTEGER   FPUTC
C Local:
      INTEGER   I,J,K
      INTEGER   LOUT
      INTEGER   IER
C Dummy-array for AVS file.
      REAL*4    ZERO(13)
      DATA      ZERO /13*0.0/
C-----------------------------------------------------------------------
      LOUT = 8
      OPEN(LOUT,FILE=FILNAM(1:LENG)//'.dat',FORM='UNFORMATTED',
     &     STATUS='UNKNOWN')
C 1st: write the dimensions of the field
      IER = FPUTC(LOUT,CHAR(NX))
      IER = FPUTC(LOUT,CHAR(NY))
      IER = FPUTC(LOUT,CHAR(NZ))
C 2nd: write the field into the file
      DO I=1,NZ
        DO J=1,NY
          DO K=1,NX
            IER = FPUTC(LOUT,CHAR(INT(W(INDX,K,J,I))))
          ENDDO
        ENDDO
      ENDDO
C 3rd: write the ending parameters
      WRITE(LOUT) 0,ZERO,0,RMIN,255,RMAX
      CLOSE(LOUT)
      RETURN
      END
C End of subroutine WRFILE
C-----------------------------------------------------------------------
C23456789012345678901234567890123456789012345678901234567890123456789012
C-----------------------------------------------------------------------