w3fm07.f

Go to the documentation of this file.

C> @file
C> @brief Nine-point smoother for rectangular grids.
C> @author P. Chase @date 1975-04-01
 
C> Smooths data on a rectangular grid using a nine-point
C> smoothing operator.
C>
C> Program history log:
C> P. Chase 1975-04-01
C> Ralph Jones 1984-07-01 Change to ibm vs fortran
C> Ralph Jones 1991-04-24 Change to cray cft77 fortran
C>
C> @param[in] FIN - Real size(ncol*nrow) array of data to be smoothed
C> @param[in] CWORK - Real size(2*ncol*(nrow+2)) work array
C> @param[in] GAMMA - Complex smoothing parameter. The imaginary part must
C> be positive.
C> @param[in] NCOL - Integer number of columns in the grid
C> @param[in] NROW - Integer number of rows in the grid
C> @param[out] FOUT - Real size(ncol*nrow) array of smoothed data. May
C> be the same array as 'fin' or overlap it in any fashion.
C>
C> @author P. Chase @date 1975-04-01
      SUBROUTINE w3fm07(FIN,FOUT,CWORK,GAMMA,NCOL,NROW)
C
      REAL    FIN(NCOL,NROW)
      REAL    FOUT(NCOL,NROW)
C
      COMPLEX CWORK(NCOL,*),GAMMA,GAMMX,GAMA,GAMB,GAMC
C
      gammx  = gamma
      xswtch = aimag(gammx)
      ncolm  = ncol-1
      nrowm  = nrow-1
C
C INITIALIZE WORK ARRAY.  WORK ARRAY STARTS UP TWO ROWS SO IT CAN SMOOTH
C DOWNWARD WITHOUT OVERLAP OF SMOOTHED AND UNSMOOTHED DATA
C
      DO 10 j = 1,nrow
      DO 10 i = 1,ncol
        cwork(i,j+2) = cmplx(fin(i,j),0.)
 10   CONTINUE
      IF (xswtch .EQ. 0.) GO TO 30
      DO 20 j = 1,nrow,nrowm
        jj = j+isign(1,nrowm-j)
      DO 20 i = 1,ncol,ncolm
        ii = i+isign(1,ncolm-i)
        cwork(i,j+2) = cmplx(fin(i,jj)+fin(ii,j)-fin(ii,jj),0.)
 20   CONTINUE
C
C SET SMOOTHING OPERATORS
C
 30   gama = 0.50 * gammx * (1.0 - gammx)
      gamb = 0.25 * gammx * gammx
      gamc = 0.50 * gammx
C
C SMOOTH WORK ARRAY, PUTTING SMOOTHED POINTS DOWN TWO ROWS
C
      cwork(1,1) = cwork(1,3)
      cwork(ncol,1) = cwork(ncol,3)
      DO 40 i = 2,ncolm
      cwork(i,1) = cwork(i,3)+gamc*(cwork(i-1,3)-2.*cwork(i,3)+
     & cwork(i+1,3))
 40   CONTINUE
      DO 60 j = 2,nrowm
      DO 50 i = 1,ncol,ncolm
      cwork(i,j) = cwork(i,j+2)+gamc*(cwork(i,j+1)-2.*cwork(i,j+2)+
     & cwork(i,j+3))
 50   CONTINUE
      DO 60 i = 2,ncolm
      cwork(i,j) = cwork(i,j+2)+gama*(cwork(i+1,j+2)+cwork(i-1,j+2)+
     & cwork(i,j+1)+cwork(i,j+3)-4.*cwork(i,j+2))+gamb*(cwork(i-1,j+1)+
     & cwork(i+1,j+1)+cwork(i-1,j+3)+cwork(i+1,j+3)-4.*cwork(i,j+2))
 60   CONTINUE
      cwork(1,nrow) = cwork(1,nrow+2)
      cwork(ncol,nrow) = cwork(ncol,nrow+2)
      DO 70 i = 2,ncolm
      cwork(i,nrow) = cwork(i,nrow+2)+gamc*(cwork(i-1,nrow+2)-2.*
     & cwork(i,nrow+2)+cwork(i+1,nrow+2))
 70   CONTINUE
C
C IF IMAGINARY PART OF SMOOTHING PARAMETER IS NOT POSITIVE, DONE
C
      IF (xswtch .LE. 0.) GO TO 90
C
C OTHERWISE MOVE WORK ARRAY BACK UP TWO ROWS
C
      DO 80 jj=1,nrow
      j = nrow+1-jj
      DO 80 i=1,ncol
      cwork(i,j+2) = cwork(i,j)
 80   CONTINUE
C
C SET SMOOTHING PARAMETER FOR CONJUGATE PASS AND GO DO IT
C
      gammx  = conjg(gammx)
      xswtch = aimag(gammx)
      GO TO 30
C
C DONE.  OUTPUT SMOOTH ARRAY
C
 90   DO 100 j = 1,nrow
      DO 100 i = 1,ncol
        fout(i,j) = real(cwork(i,j))
 100  CONTINUE
      RETURN
      END