subroutine max3_set (f, fp)
*
*  smooth a function by replacing the function value by an averge of
*  three function values.
*
      include 'cparam.inc'
      integer :: direction
      real, DIMENSION(mx,my,mz):: f, fp, gp
!HPF$ DISTRIBUTE (*,BLOCK,BLOCK) :: f, fp

        fp(2:mx-1,:,:) = amax1(f(2:mx-1,:,:),
     &                         f(1:mx-2,:,:),f(3:mx,:,:))
        fp(1,:,:)      = amax1(f(1,:,:),
     &                         f(mx,:,:),f(2,:,:))
        fp(mx,:,:)     = amax1(f(mx,:,:),
     &                        f(mx-1,:,:),f(1,:,:))
        gp(:,2:my-1,:) = amax1(fp(:,2:my-1,:),
     &                         fp(:,1:my-2,:),fp(:,3:my,:))
        gp(:,1,:)      = amax1(fp(:,1,:),
     &                         fp(:,my,:),fp(:,2,:))
        gp(:,my,:)     = amax1(fp(:,my,:),
     &                         fp(:,my-1,:),fp(:,1,:))
        fp(:,:,2:mz-1) = amax1(gp(:,:,2:mz-1),
     &                         gp(:,:,1:mz-2),gp(:,:,3:mz))
        fp(:,:,1)      = amax1(gp(:,:,1),
     &                         gp(:,:,mz),gp(:,:,2))
        fp(:,:,mz)     = amax1(gp(:,:,mz),
     &                         gp(:,:,mz-1),gp(:,:,1))

      nflop=nflop+9
      end
      subroutine smooth3_set (f, fp)
*
*  smooth a function by replacing the function value by an averge of
*  three function values.
*
      include 'cparam.inc'
      integer :: direction
      real, DIMENSION(mx,my,mz):: f, fp, gp
!HPF$ DISTRIBUTE (*,BLOCK,BLOCK) :: f, fp

        fp(2:mx-1,:,:) = 0.5*f(2:mx-1,:,:)+
     &                   .25*(f(1:mx-2,:,:)+f(3:mx,:,:))
        fp(1,:,:)      = 0.5*f(1,:,:)+
     &                   .25*(f(mx,:,:)+f(2,:,:))
        fp(mx,:,:)     = 0.5*f(mx,:,:)+
     &                  .25*(f(mx-1,:,:)+f(1,:,:))
        gp(:,2:my-1,:) = 0.5*fp(:,2:my-1,:)+
     &                   .25*(fp(:,1:my-2,:)+fp(:,3:my,:))
        gp(:,1,:)      = 0.5*fp(:,1,:)+
     &                   .25*(fp(:,my,:)+fp(:,2,:))
        gp(:,my,:)     = 0.5*fp(:,my,:)+
     &                   .25*(fp(:,my-1,:)+fp(:,1,:))
        fp(:,:,2:mz-1) = 0.5*gp(:,:,2:mz-1)+
     &                   .25*(gp(:,:,1:mz-2)+gp(:,:,3:mz))
        fp(:,:,1)      = 0.5*gp(:,:,1)+
     &                   .25*(gp(:,:,mz)+gp(:,:,2))
        fp(:,:,mz)     = 0.5*gp(:,:,mz)+
     &                   .25*(gp(:,:,mz-1)+gp(:,:,1))

      nflop=nflop+12
      end