d_invert_tridiag_block_matrix Subroutine

subroutine d_invert_tridiag_block_matrix(Nb,N,sub_,diag_,over_,Ainv)
  integer                              :: ib,i,j,Nb,N
  real(8),dimension(Nb,N,N)            :: diag_
  real(8),dimension(Nb-1,N,N)          :: sub_
  real(8),dimension(Nb-1,N,N)          :: over_
  real(8),dimension(Nb,N,N)            :: Ainv
  real(8),dimension(N,N)               :: Foo,Cleft,Cright
  real(8),dimension(Nb,N,N)            :: Dleft
  real(8),dimension(Nb,N,N)            :: Dright
  !
  !DOWNWARD:
  dleft(1,:,:) = diag_(1,:,:)
  do i=2,Nb
     foo  = dleft(i-1,:,:) ; call inv(foo)
     cleft= matmul(sub_(i-1,:,:),foo)
     dleft(i,:,:) = diag_(i,:,:) - matmul(cleft,over_(i-1,:,:))
  enddo
  !
  !BACKWARD:
  dright(Nb,:,:) = diag_(Nb,:,:)
  do i=Nb-1,1,-1
     foo   = dright(i+1,:,:) ; call inv(foo)
     cright= matmul(over_(i,:,:),foo)
     dright(i,:,:) = diag_(i,:,:) - matmul(cright,sub_(i,:,:))
  enddo
  !
  do ib=1,Nb
     Ainv(ib,:,:)    =  dleft(ib,:,:) + dright(ib,:,:) - diag_(ib,:,:)
     call inv(Ainv(ib,:,:))
  end do
end subroutine d_invert_tridiag_block_matrix