Code source de l'adjoint du filtre récursif 1D

!
!----------------------------------------------------------------------
!
! A D J O I N T   R E C U R S I V E    F I L T E R
!
!----------------------------------------------------------------------
!
   SUBROUTINE recursive_filter_adj_1d_r(A, alpha, pass)
!
! 1D adjoint 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))
!
      B(:) = 0.
      C(:) = 0. 
!
! Backward langevin
!
      C(1:N) = C(1:N) + A(1:N)
      A(1:N) = 0.
      DO i = 1, N - 1
        B(i) = B(i) + ( 1 - alpha ) * C(i)
        C(i+1) = C(i+1) + alpha * C(i)
        C(i) = 0.
      ENDDO
      SELECT CASE (pass)
      CASE (1)
         B(N) = B(N) + 1 / ( 1 + alpha ) * C(N)
      CASE (2)
         B(N)   = B(N)   + ( 1 - alpha ) / ( 1 - alpha**2 )**2 * C(N)
         B(N-1) = B(N-1) - alpha**3 * ( 1 - alpha ) / ( 1 - alpha**2 )**2 * C(N)
      CASE (3)
         B(N)   = B(N)   - ( 1 - alpha ) / ( 1 - alpha**2 )**3 * C(N)
         B(N-1) = B(N-1) + alpha**3 * ( 3 - alpha**2 ) * ( 1 - alpha ) / ( 1 - alpha**2 )**3 * C(N)
         B(N-2) = B(N-2) - alpha**4 * ( 1 - alpha ) / ( 1 - alpha**2 )**3 * C(N)
      CASE default
         B(N)   = B(N)   - ( 1 - alpha ) / ( 1 - alpha**2 )**3 * C(N)
         B(N-1) = B(N-1) + alpha**3 * ( 3 - alpha**2 ) * ( 1 - alpha ) / ( 1 - alpha**2 )**3 * C(N)
         B(N-2) = B(N-2) - alpha**4 * ( 1 - alpha ) / ( 1 - alpha**2 )**3 * C(N)
      END SELECT
!
! Forward Langevin
!
      DO i = N, 2, -1
        A(i) = A(i) + ( 1 - alpha ) * B(i)
        B(i-1) = B(i-1) + alpha * B(i)
        B(i) = 0.
      ENDDO
      SELECT CASE (pass)
      CASE (1)
         A(1) = A(1) + ( 1 - alpha ) * B(1)
      CASE (2)
         A(1) = A(1) + 1 / ( 1 + alpha ) * B(1)
      CASE (3)
         A(1) = A(1) + ( 1 - alpha ) / ( 1 - alpha**2 )**2 * B(1)
         A(2) = A(2) - alpha**3 * ( 1 - alpha ) / ( 1 - alpha**2 )**2 * B(1)
      CASE default
         A(1) = A(1) + ( 1 - alpha ) / ( 1 - alpha**2 )**2 * B(1)
         A(2) = A(2) - alpha**3 * ( 1 - alpha ) / ( 1 - alpha**2 )**2 * B(1)
      END SELECT
!
      DEALLOCATE (B)
      DEALLOCATE (C)
!
   END SUBROUTINE recursive_filter_adj_1d_r