anchor.diff:  END FUNCTION
   
   
  -!The integral of G(XI,XJ) from XJM to XJP, nodes of the J element
  -REAL FUNCTION BIJ(XI, XJM, XJP)
  -      IMPLICIT NONE
  -      REAL, DIMENSION(2) :: XI, XJM, XJP, EX, VXP, VXM
  -      REAL :: XM, YM, XP, YP, RP, RM
  -      REAL, PARAMETER :: PI = 4.*ATAN(1.)
  -      EX = XJP - XJM
  -      EX = EX/SQRT(EX(1)**2 + EX(2)**2)
  -      VXP = (XJP - XI)
  -      VXM = (XJM - XI)
  -      XP = VXP(1)*EX(1) + VXP(2)*EX(2)
  -      XM = VXM(1)*EX(1) + VXM(2)*EX(2)
  -!     YP and YM should be equals
  -      YP =-VXP(1)*EX(2) + VXP(2)*EX(1)
  -      YM =-VXM(1)*EX(2) + VXM(2)*EX(1)
  -      RP = SQRT(XP**2 + YP**2)
  -      RM = SQRT(XM**2 + YM**2)
  -      BIJ = (XP*LOG(RP) - XP + YP*ATAN2(XP,YP) - XM*LOG(RM) + XM - YM*ATAN2(XM,YM))/(2.*PI)
  -      RETURN
  -END FUNCTION
  -
  -
  -!The integral of dG/dn(XS,XJ) from XJM to XJP, nodes of the J element
  -REAL FUNCTION CIJ(XS, XJM, XJP)
  -      IMPLICIT NONE
  -      REAL, DIMENSION(2) :: XS, XJM, XJP, EX, VXP, VXM
  -      REAL :: XM, XP, YS
  +!BIJ is the integral of G(XI,XJ) from XJM to XJP, nodes of the J element
  +!CIJ is the integral of dG/dn(XI,XJ) from XJM to XJP, nodes of the J element
  +SUBROUTINE BCIJ(XI, XJM, XJP, BIJ, CIJ)
  +      IMPLICIT NONE
  +      REAL, DIMENSION(2), INTENT(IN) :: XI, XJM, XJP
  +      REAL, INTENT(OUT) :: BIJ, CIJ
  +      REAL, DIMENSION(2) :: EX, VXP, VXM
  +      REAL :: XM, XP, YS, RM, RP
         REAL, PARAMETER :: PI = 4.*ATAN(1.)
         EX = XJP - XJM
         EX = EX/SQRT(EX(1)**2 + EX(2)**2)
   !     Vector X_+
  -      VXP = (XS - XJP)
  +      VXP = (XI - XJP)
   !     Vector X_-
  -      VXM = (XS - XJM)
  +      VXM = (XI - XJM)
         XP = -VXP(1)*EX(1) - VXP(2)*EX(2)
         XM = -VXM(1)*EX(1) - VXM(2)*EX(2)
         YS = -VXP(1)*EX(2) + VXP(2)*EX(1)
  -      CIJ = (ATAN2(XP,YS) - ATAN2(XM,YS))/(2.*PI)
  -      RETURN
  -END FUNCTION
  +      RP = SQRT(XP**2 + YS**2)
  +      RM = SQRT(XM**2 + YS**2)
  +      BIJ = (XP*LOG(RP) - XP + YS*ATAN2(XP,YS) - XM*LOG(RM) + XM - YS*ATAN2(XM,YS))/(2.*PI)
  +      CIJ = (ATAN2(XP,YS) - ATAN2(XM,YS))/(2.*PI)      
  +END SUBROUTINE
   
   
   !Calculate the matrix B and C
  
positive.code: SUBROUTINE SRFMATBC(N, XNODE, B, C)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: N
        REAL, DIMENSION(N,2), INTENT(IN) :: XNODE
        REAL, DIMENSION(N,N), INTENT(OUT) :: B, C
        REAL, DIMENSION(2) :: XI
        REAL :: CIJ, BIJ
        INTEGER :: I, J, NODOFEL
        DO I = 1, N
         !Midpoint of I element
         XI = (XNODE(NODOFEL(N,I,1),:) + XNODE(NODOFEL(N,I,2),:))/2.
         C(I,I) = 0.
         DO J = 1, N
          B(I,J) = BIJ(XI, XNODE(NODOFEL(N,J,1),:), XNODE(NODOFEL(N,J,2),:))
          IF (I .NE. J) THEN
            C(I,J) = CIJ(XI, XNODE(NODOFEL(N,J,1),:), XNODE(NODOFEL(N,J,2),:))
          END IF
         END DO
        END DO
  
negative.code: SUBROUTINE BCONDVEL(N, XNODE, U, CHI)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: N
        REAL, DIMENSION(N,3), INTENT(IN) :: XNODE
        REAL, DIMENSION(2), INTENT(IN) :: U
        REAL, DIMENSION(N), INTENT(OUT) :: CHI
        REAL, DIMENSION(2) :: T
        INTEGER :: I, NODOFEL
        DO I = 1, N
         T = XNODE(NODOFEL(N,I,2),:2) - XNODE(NODOFEL(N,I,1),:2)
         T = T/SQRT(T(1)**2 + T(2)**2)
         CHI(I) = U(1)*T(2) - U(2)*T(1)
        END DO
  