anchor.diff:  
   !Impose boundary conditions for a moving circle with u velocity
   !in the air frame of reference
  -SUBROUTINE BCondVel(N, XNODE, U, CHI)
  +subroutine BCondVel(N, XNODE, U, CHI)
         IMPLICIT NONE
  -      INTEGER, INTENT(IN) :: N
  -      REAL(KIND=8), DIMENSION(N,2), INTENT(IN) :: XNODE
  -      REAL(KIND=8), DIMENSION(2), INTENT(IN) :: U
  -      REAL(KIND=8), DIMENSION(N), INTENT(OUT) :: CHI
  -      REAL(KIND=8), DIMENSION(2) :: T
  -      INTEGER :: I, NOE
  -      DO I = 1, N
  +      integer, intent(IN) :: N
  +      real(kind=8), dimension(N,2), intent(IN) :: XNODE
  +      real(kind=8), dimension(2), intent(IN) :: U
  +      real(kind=8), dimension(N), intent(OUT) :: CHI
  +      real(kind=8), dimension(2) :: T
  +      integer :: I, NOE
  +      do I = 1, N
          T = XNODE(NOE(N,I,1),:2) - XNODE(NOE(N,I,0),:2)
          T = T/SQRT(T(1)**2 + T(2)**2)
          CHI(I) = U(1)*T(2) - U(2)*T(1)
  -      END DO
  -END SUBROUTINE
  +      end do
  +end subroutine
   
   
   !We assume that U Horizontal is constant
   subroutine BCondLap(Nelem, Xnode, Ampl, ChiLap)
         IMPLICIT NONE
  -      INTEGER, INTENT(IN) :: Nelem
  -      REAL(KIND=8), DIMENSION(Nelem,2), INTENT(IN) :: Xnode
  -      REAL(KIND=8), INTENT(IN) :: Ampl
  -      complex(kind=8), dimension(Nelem), INTENT(OUT) :: ChiLap
  +      integer, intent(IN) :: Nelem
  +      real(kind=8), dimension(Nelem,2), intent(IN) :: Xnode
  +      real(kind=8), intent(IN) :: Ampl
  +      complex(kind=8), dimension(Nelem), intent(OUT) :: ChiLap
         real(kind=8), dimension(Nelem,2) :: n
         integer :: I
         call normals(Nelem, Xnode, n)
  
positive.code: 
  
  !Solve the (NxN) linear system (0.5*I - C)*phi = C*chi
  SUBROUTINE SOLVEPHI(N, B, C, PHI, CHI)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: N
        REAL(KIND=8), DIMENSION(N,N), INTENT(IN) :: B, C
        REAL(KIND=8), DIMENSION(N), INTENT(OUT) :: PHI
        REAL(KIND=8), DIMENSION(N), INTENT(IN) :: CHI
        INTEGER :: I, INFO
        REAL(KIND=8), DIMENSION(N,N) :: A
        REAL(KIND=8), DIMENSION(N) :: RHS, IPIV
        A = -C
        DO I = 1, N
         A(I,I) = A(I,I) + 0.5
        END DO
        RHS = MATMUL(B,CHI)
        CALL DGESV(N, 1, A, N, IPIV, RHS, N, INFO)
        IF (INFO .NE. 0) THEN
          WRITE(*,*) "ERROR IN LINEAR SYSTEM"
        END IF
        PHI = RHS
  END SUBROUTINE
  
  
  SUBROUTINE SolvePhiTime(N, B, C, NWake, D, NTime, ChiTime, PhiTime, DPhiW)
        IMPLICIT NONE
        integer, intent(IN) :: N, NWake, NTime
        real(kind=8), dimension(N,N), intent(IN) :: B, C
  
negative.code:       do i = 1, N
         A(i,i) = A(i,i) + dcmplx(0.5,0.)
        end do
        RHS = matmul(dcmplx(B),ChiLap)
        call ZGESV(N, 1, A, N, IPIV, RHS, N, INFO)
        if (INFO .NE. 0) then
          write(*,*) "ERROR IN LINEAR SYSTEM"
        end if
        PhiLap = RHS
  end subroutine
  
  
  !Calculate the potential in the field
  !  PHISRF   On surface                          (in)
  !  CHISRF   On surface                          (in)
  !  B, C     In the field                        (in)
  !  PHifLD   In the field                        (out)
  subroutine CalcPhiFld(N, PhiSrf, ChiSrf, NX, B, C, PHiFld)
        IMPLICIT NONE
        integer, intent(IN) :: N, NX
        real(kind=8), dimension(N), intent(IN) :: PhiSrf, ChiSrf
        real(kind=8), dimension(NX,N) :: B, C
        real(kind=8), dimension(NX), intent(OUT) :: PhiFld
        PhiFld = matmul(B, ChiSrf) + matmul(C, PhiSrf)
  end subroutine
  