!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE RECURSIVE_FILTER            ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  recursive filter in xy two direction. according to paper
!  published by Lorenc and Purser.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!
!  Jidong Gao, CAPS, July, 2000
!
!-----------------------------------------------------------------------
!

SUBROUTINE  recurfilt_2d(nx,ny,pgrd,ipass_filt,radius)
!
!
  INTEGER :: nx, ny, ipass_filt, radius
  DIMENSION pgrd(nx,ny)
  REAL, DIMENSION (:), allocatable :: temx
  REAL, DIMENSION (:), allocatable :: temy
!
!
  allocate (temx(nx))
  allocate (temy(ny))
!
!
  ee    = REAL(ipass_filt) / REAL(radius* radius)
  alpha = 1 + ee - SQRT(ee*(ee+2.))
!
!
  DO n = 1, ipass_filt
!
    DO j = 1, ny
!
      DO i = 1, nx
        temx (i) = pgrd(i,j)
      END DO
!
      CALL recurfilt_1d( temx,nx,alpha )
!
      DO i = 1, nx
        pgrd(i,j) = temx(i)
      END DO
!
    END DO
!
!
    DO i = 1, nx
!
      DO j = 1, ny
        temy (j) = pgrd(i,j)
      END DO
!
      CALL recurfilt_1d( temy,ny,alpha )
!
      DO j = 1, ny
        pgrd(i,j) = temy (j)
      END DO
!
    END DO
!
  END DO
!
!
  deallocate (temx)
  deallocate (temy)
!
  RETURN
END SUBROUTINE  recurfilt_2d
!
!


SUBROUTINE  arecurfilt_2d(nx,ny,pgrd,ipass_filt,radius)
!
!
  INTEGER :: nx, ny, ipass_filt
  DIMENSION pgrd(nx,ny)
  REAL, DIMENSION (:), allocatable :: temx
  REAL, DIMENSION (:), allocatable :: temy
!
!
  allocate (temx(nx))
  allocate (temy(ny))
!
!
  ee    = REAL(ipass_filt) / REAL( radius*radius )
  alpha = 1 + ee - SQRT(ee*(ee+2.))
!
!
  DO i = 1, nx
    temx (i) = 0.
  END DO
!
  DO j = 1, ny
    temy (j) = 0.
  END DO
!
!
  DO n = 1, ipass_filt
!
!
    DO i =  nx, 1, -1
      DO j = ny, 1, -1
        temy(j)   = temy(j) + pgrd(i,j)
        pgrd(i,j) = 0.
      END DO
!
      CALL arecurfilt_1d ( temy,ny,alpha )
!
      DO j = ny, 1, -1
        pgrd(i,j) = pgrd(i,j)+ temy(j)
        temy(j)   = 0.
      END DO
!
    END DO
!
!
    DO j = ny, 1, -1
!
      DO i = nx, 1, -1
        temx(i)   = temx(i) + pgrd(i,j)
        pgrd(i,j) = 0.
      END DO
!
      CALL arecurfilt_1d( temx,nx,alpha )
!
      DO i = nx, 1, -1
        pgrd(i,j) = pgrd(i,j) + temx(i)
        temx(i)   = 0.
      END DO
!
    END DO
!
  END DO
!
!
  deallocate (temx)
  deallocate (temy)
!
!
  RETURN
END SUBROUTINE  arecurfilt_2d