Code source du filtre récursif 1D

!
!----------------------------------------------------------------------
!
! R E C U R S I V E    F I L T E R
!
! from Iterative analysis using covariance functions and filters by
!      Andrew Lorenc (Q. J. R. Meteorol. Soc. (1992), 118, pp. 569-591)
!      Appendix A pp 587-589 
!
!----------------------------------------------------------------------
!
   SUBROUTINE recursive_filter_1d_r(A, alpha, pass)
!
! 1D recursive filter for reals
!
! A is the initial value
! B is the value after filtering for i = 1 to n
! C (in our case A) is the value after one pass of the filter in each 
!   direction
!
! alpha is the filter coefficient
! N is the number of time step
!
      IMPLICIT NONE
!
      REAL, DIMENSION(:), INTENT(inout) :: A
      REAL, INTENT(in)                  :: alpha
      INTEGER, INTENT(in)               :: pass
!
      INTEGER                           :: N, i
      REAL, DIMENSION(:), ALLOCATABLE   :: B, C 
!
! Look for time dimension
!
      N = SIZE(A,1)
!
! Allocate array
!
      ALLOCATE (B(N))
      ALLOCATE (C(N))
!
!
! Forward Langevin
!
      B(:) = 0.0
      SELECT CASE (pass)
      CASE (1)
         B(1) = ( 1 - alpha ) * A(1)
      CASE (2)
         B(1) = 1 / ( 1 + alpha ) * A(1)
         B(1) = ( 1 - alpha ) / ( 1 - alpha**2 ) *  A(1)
      CASE (3)
         B(1) = ( 1 - alpha ) / ( 1 - alpha**2 )**2 * ( A(1) - alpha**3 * A(2) )
      CASE default
         B(1) = ( 1 - alpha ) / ( 1 - alpha**2 )**2 * ( A(1) - alpha**3 * A(2) )
      END SELECT
      DO i = 2, N
         B(i) = alpha * B(i-1) + ( 1 - alpha ) * A(i) 
      ENDDO
!
! Backward langevin
!
      SELECT CASE (pass)
      CASE (1)
         C(N) = 1 / ( 1 + alpha ) * B(N)
      CASE (2)
         C(N) = ( 1 - alpha ) / ( 1 - alpha**2 )**2 * ( B(N) - alpha**3 * B(N-1) )
      CASE (3)
         C(N) = ( 1 - alpha ) / ( 1 - alpha**2 )**3 * &
                ( -1 * B(N) + ( 3 - alpha**2 ) * alpha**3 * B(N-1) - alpha**4 * B(N-2) )
      CASE default
         C(N) = ( 1 - alpha ) / ( 1 - alpha**2 )**3 * &
                ( -1 * B(N) + ( 3 - alpha**2 ) * alpha**3 * B(N-1) - alpha**4 * B(N-2) )
      END SELECT
      DO i = N-1, 1, -1
         C(i) = alpha * C(i+1) + ( 1 - alpha ) * B(i)
      ENDDO
!
      A(1:N) = C(1:N)
!
!
      DEALLOCATE (B)
      DEALLOCATE (C)
!
   END SUBROUTINE recursive_filter_1d_r