!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE SCALE_BCKGRD              ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!-----------------------------------------------------------------------
!

SUBROUTINE recurfilt_3d(nx,ny,nz,pgrd,ipass_filt,hradius,nradius_z) 2,3
!cc
!cc
  IMPLICIT none
  INTEGER :: i,j,k,n
  INTEGER :: nx,ny,nz,ipass_filt,nradius_z
     REAL :: hradius
  REAL :: alpha, alpha_z, ee
  REAL :: temp,temp2

  REAL :: pgrd(nx,ny,nz)
  REAL, DIMENSION(:), allocatable :: temx
  REAL, DIMENSION(:), allocatable :: temy
  REAL, DIMENSION(:), allocatable :: temz
!
  REAL,   DIMENSION (:), allocatable :: varb
!
  allocate ( temx(nx) )
  allocate ( temy(ny) )
  allocate ( temz(nz) )
!
!
  IF( hradius == 0 .and. nradius_z == 0 ) return 
!
  ee    = REAL(ipass_filt) / (hradius*hradius)
  alpha = 1+ee-SQRT( ee*(ee+2.) )
!
  IF( nradius_z /= 0 ) THEN
    ee     = REAL (ipass_filt) / REAL (nradius_z*nradius_z)
    alpha_z = 1 + ee - SQRT (ee*(ee+2.))
  ENDIF
!
!
  DO n = 1, ipass_filt/2
!
!
    allocate ( varb(nx) )
    DO k = 1, nz
      DO j = 1, ny
!
        if(n==1) varb(1) = (1-alpha) * pgrd(1,j,k)
        if(n==2) varb(1) = (1-alpha)/(1-alpha*alpha) * pgrd(1,j,k)
        if(n==3) then
          temp = (1-alpha)/((1-alpha*alpha)*(1-alpha*alpha))
          temp2 =alpha*alpha*alpha
          varb(1) = temp * (pgrd(1,j,k)-temp2*pgrd(1,j,k))
        ENDIF
        if(n>=4) then
          temp2 =alpha*alpha*alpha
          temp = (1-alpha)/(1-3*alpha*alpha+3*temp2*alpha-temp2*temp2)
          varb(1) = temp * (pgrd(1,j,k)-3*temp2*pgrd(2,j,k)+              &
                 temp2*alpha*alpha*pgrd(2,j,k)+temp2*alpha*pgrd(3,j,k))
        ENDIF
!
        DO i = 2, nx, 1
          varb(i) = alpha*varb(i-1) + (1.-alpha)*pgrd(i,j,k)
        END DO
!
        if(n==0) pgrd(nx,j,k) = (1-alpha) * varb(nx)
        if(n==1) pgrd(nx,j,k) = (1-alpha)/(1-alpha*alpha) * varb(nx)
        if(n==2) then
          temp = (1-alpha)/((1-alpha*alpha)*(1-alpha*alpha))
          temp2 =alpha*alpha*alpha
          pgrd(nx,j,k) = temp * (varb(nx)-temp2*varb(nx-1))
        ENDIF
        if(n>=3) then
          temp2 =alpha*alpha*alpha
          temp = (1-alpha)/(1-3*alpha*alpha+3*temp2*alpha-temp2*temp2)
          pgrd(nx,j,k) = temp * (varb(nx)-3*temp2*varb(nx-1)+            &
               temp2*alpha*alpha*varb(nx-1)+temp2*alpha*varb(nx-2))
        ENDIF
!
!
        DO i = nx-1, 1, -1
          pgrd(i,j,k) = alpha*pgrd(i+1,j,k) + (1.-alpha)*varb(i)
        END DO
!
      END DO
    END DO
    deallocate (varb)
!
!
    allocate ( varb(ny) )
    DO k = 1, nz
      DO i = 1, nx
!
        if(n==1) varb(1) = (1-alpha) * pgrd(i,1,k)
        if(n==2) varb(1) = (1-alpha)/(1-alpha*alpha) * pgrd(i,1,k)
        if(n==3) then
          temp = (1-alpha)/((1-alpha*alpha)*(1-alpha*alpha))
          temp2 =alpha*alpha*alpha
          varb(1) = temp * (pgrd(i,1,k)-temp2*pgrd(i,1,k))
        ENDIF
        if(n>=4) then
          temp2 =alpha*alpha*alpha
          temp = (1-alpha)/(1-3*alpha*alpha+3*temp2*alpha-temp2*temp2)
          varb(1) = temp * (pgrd(i,1,k)-3*temp2*pgrd(i,2,k)+              &
                 temp2*alpha*alpha*pgrd(i,2,k)+temp2*alpha*pgrd(i,3,k))
        ENDIF
!
        DO j = 2, ny, 1
          varb(j) = alpha*varb(j-1) + (1.-alpha)*pgrd(i,j,k)
        END DO
!
        if(n==0) pgrd(i,ny,k) = (1-alpha) * varb(ny)
        if(n==1) pgrd(i,ny,k) = (1-alpha)/(1-alpha*alpha) * varb(ny)
        if(n==2) then
          temp = (1-alpha)/((1-alpha*alpha)*(1-alpha*alpha))
          temp2 =alpha*alpha*alpha
          pgrd(i,ny,k) = temp * (varb(ny)-temp2*varb(ny-1))
        ENDIF
        if(n>=3) then
          temp2 =alpha*alpha*alpha
          temp = (1-alpha)/(1-3*alpha*alpha+3*temp2*alpha-temp2*temp2)
          pgrd(i,ny,k) = temp * (varb(ny)-3*temp2*varb(ny-1)+            &
               temp2*alpha*alpha*varb(ny-1)+temp2*alpha*varb(ny-2))
        ENDIF
!
        DO j = ny-1, 1, -1
          pgrd(i,j,k) = alpha*pgrd(i,j+1,k) + (1.-alpha)*varb(j)
        END DO
!
      END DO
    END DO
    deallocate (varb)
!
!
    IF( nradius_z /= 0 ) THEN
      allocate ( varb(nz) )
      DO i = 1, nx
        DO j = 1, ny

        if(n==1) varb(1) = (1-alpha_z) * pgrd(i,j,1)
        if(n==2) varb(1) = (1-alpha_z)/(1-alpha_z*alpha_z) * pgrd(i,j,1)
        if(n==3) then
          temp = (1-alpha_z)/((1-alpha_z*alpha_z)*(1-alpha_z*alpha_z))
          temp2 =alpha_z*alpha_z*alpha_z
          varb(1) = temp * (pgrd(i,j,1)-temp2*pgrd(i,j,1))
        ENDIF
        if(n>=4) then
          temp2 =alpha_z*alpha_z*alpha_z
          temp = (1-alpha_z)/(1-3*alpha_z*alpha_z+3*temp2*alpha_z-temp2*temp2)
          varb(1) = temp * (pgrd(i,j,1)-3*temp2*pgrd(i,j,2)+              &
                 temp2*alpha_z*alpha_z*pgrd(i,j,2)+temp2*alpha_z*pgrd(i,j,3))
        ENDIF
!
        DO k = 2, nz, 1
          varb(k) = alpha_z*varb(k-1) + (1.-alpha_z)*pgrd(i,j,k)
        END DO
!
        if(n==0) pgrd(i,j,nz) = (1-alpha_z) * varb(nz)
        if(n==1) pgrd(i,j,nz) = (1-alpha_z)/(1-alpha_z*alpha_z) * varb(nz)
        if(n==2) then
          temp = (1-alpha_z)/((1-alpha_z*alpha_z)*(1-alpha_z*alpha_z))
          temp2 =alpha_z*alpha_z*alpha_z
          pgrd(i,j,nz) = temp * (varb(nz)-temp2*varb(nz-1))
        ENDIF
        if(n>=3) then
          temp2 =alpha_z*alpha_z*alpha_z
          temp = (1-alpha_z)/(1-3*alpha_z*alpha_z+3*temp2*alpha_z-temp2*temp2)
          pgrd(i,j,nz) = temp * (varb(nz)-3*temp2*varb(nz-1)+            &
               temp2*alpha_z*alpha_z*varb(nz-1)+temp2*alpha_z*varb(nz-2))
        ENDIF
!
        DO k = nz-1, 1, -1
          pgrd(i,j,k) = alpha_z*pgrd(i,j,k+1) + (1.-alpha_z)*varb(k)
        END DO

        END DO
      END DO
      deallocate (varb)
    ENDIF
!
  END DO
!
!
  deallocate (temx)
  deallocate (temy)
  deallocate (temz)
!
!
  RETURN
END SUBROUTINE recurfilt_3d
!
!


SUBROUTINE  arecurfilt_3d(nx,ny,nz,pgrd,ipass_filt,hradius,nradius_z)
!c
!c
  IMPLICIT none
  INTEGER :: i,j,k,n
  INTEGER :: nx,ny,nz,ipass_filt,nradius_z
  REAL    :: hradius
  REAL :: alpha, alpha_z, ee
  REAL :: pgrd(nx, ny, nz)
  REAL, DIMENSION(:), allocatable :: temx
  REAL, DIMENSION(:), allocatable :: temy
  REAL, DIMENSION(:), allocatable :: temz
!
  IF( hradius == 0 .and. nradius_z == 0 ) return 
!
  allocate ( temx(nx) )
  allocate ( temy(ny) )
  allocate ( temz(nz) )
!
!
  DO i = 1, nx
    temx(i) = 0.
  END DO
!
  DO j = 1, ny
    temy(j) = 0.
  END DO
!
  DO k = 1, nz
    temz(k) = 0.
  END DO
!
  IF( nradius_z /= 0 ) THEN
    ee    = REAL(ipass_filt)/REAL(nradius_z*nradius_z)
    alpha_z = 1 + ee - SQRT (ee*(ee+2.))
  ENDIF
!
  ee    = REAL(ipass_filt) / (hradius*hradius)
  alpha = 1 + ee - SQRT (ee*(ee+2.))
!
!
  DO n = ipass_filt/2, 1 , -1
!
!
  IF( nradius_z /= 0 ) THEN
      DO i = nx, 1, -1
        DO j = ny, 1, -1
          DO k = nz, 1, -1
            temz(k) = temz(k)+pgrd(i,j,k)
            pgrd(i,j,k) = 0.
          END DO
!
          CALL arecurfilt_1d(temz,nz,alpha_z,n)
!
          DO k = nz, 1, -1
            pgrd(i,j,k) = pgrd(i,j,k)+temz(k)
            temz(k) = 0.
          END DO
        END DO
      END DO
    ENDIF
!
!
    DO k =  nz, 1, -1
      DO i =  nx, 1, -1
        DO j = ny, 1, -1
          temy(j) = temy(j)+pgrd(i,j,k)
          pgrd(i,j,k) = 0.
        END DO
!
        CALL arecurfilt_1d(temy,ny,alpha,n)
!
        DO j = ny, 1, -1
          pgrd(i,j,k) = pgrd(i,j,k)+temy(j)
          temy(j)   = 0.
        END DO
      END DO
    END DO
!
!
    DO k = nz, 1, -1
      DO j = ny, 1, -1
        DO i = nx, 1, -1
          temx(i) = temx(i)+pgrd(i,j,k)
          pgrd(i,j,k) = 0.
        END DO
!
        CALL arecurfilt_1d(temx,nx,alpha,n)
!
        DO i = nx, 1, -1
          pgrd(i,j,k) = pgrd(i,j,k)+temx(i)
          temx(i) = 0.
        END DO
      END DO
    END DO
!
  END DO
!
!
  deallocate (temx)
  deallocate (temy)
!
!
  RETURN
END SUBROUTINE  arecurfilt_3d