anchor.diff:     real(ReKi)              :: xb              !< 
      real(ReKi)              :: yb              !< 
      real(ReKi)              :: zb              !< 
  -   real(ReKi)              :: visc_param2     !< Square of viscous param
  +   real(ReKi)              :: SmoothParam2     !< Square of viscous param
      real(ReKi)              :: exp_value       !< 
      real(ReKi),parameter    :: fourpi_inv =  0.25_ReKi / ACOS(-1.0_Reki )
   
  -   !$OMP PARALLEL default(shared)
  -   !$OMP do private(&
  -   !$OMP& icp,is,Uind,P1,P2,DP1,DP2,norm2_r0,Gam,visc_param,&
  -   !$OMP& crossprod,denom,denominator,h2,h,Kv,norm_a,norm_b,norm2_crossprod,xa,ya,za,xb,yb,zb,visc_param2,exp_value&
  -   !$OMP& ) schedule(runtime)
  +   !OMP PARALLEL default(shared)
  +   !OMP do private(&
  +   !OMP& icp,is,Uind,P1,P2,DP1,DP2,norm2_r0,Gam,SmoothParam,&
  +   !OMP& crossprod,denom,denominator,h2,h,Kv,norm_a,norm_b,norm2_crossprod,xa,ya,za,xb,yb,zb,SmoothParam2,exp_value&
  +   !OMP& ) schedule(runtime)
      ! loop on CPs 
      do icp=iCPStart,iCPEnd
         do is=iSegStart,iSegEnd ! loop on selected segments 
  
positive.code:             crossprod(2) = za*xb-xa*zb
              crossprod(3) = xa*yb-ya*xb
              norm2_crossprod=crossprod(1)**2+ crossprod(2)**2+ crossprod(3)**2
              !print "(A,6F21.16)","xa xb:",xa,ya,za,xb,yb,zb
              ! check for singularity  
              ! TODO this check is viscous model dependent(e.g. crossprod)... have to go towards different funcitons soon..
              ! to avoid singularity, save a bit of computation and use limit values
              if (denominator<PRECISION_UI .or. norm2_crossprod<PRECISION_UI) then
                 !--- Exactly on the Singularity, velocity is zero
                 !Uind(1:3)=0.0_ReKi
  
negative.code:          Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
        end do ! Loop on segments
     enddo ! Loop on control points
     !OMP END DO 
     !OMP END PARALLEL
  end subroutine  
  
  
  !> Velocity induced by one vortex quad on nCPs Control Points
  subroutine ui_quad_n1(CPs, nCPs, P1, P2, P3, P4, Gamm, SmoothModel, SmoothParam, Uind)
     ! Arguments declarations 
     integer, intent(in)                        :: nCPs        !< 
     real(ReKi), dimension(3,nCPs), intent(in)  :: CPs         !< 
     real(ReKi), dimension(3), intent(in)       :: P1,P2,P3,P4 !< Coordinates of vortex quadrilateral
     real(ReKi),                    intent(in)  :: Gamm
     integer(IntKi) , intent(in)                :: SmoothModel !< Most likely should be 0
     real(ReKi),                    intent(in)  :: SmoothParam !< 
     real(ReKi), dimension(3,nCPs), intent(out) :: Uind        !< no side effects!!!
     ! Variable declarations 
     real(ReKi), dimension(3) :: CP           !< 
     real(ReKi), dimension(3) :: Uindtmp        !< 
     real(ReKi), dimension(3) :: DP1          !< 
     real(ReKi), dimension(3) :: DP2          !< 
     integer                :: icp
     ! 
     !OMP PARALLEL DEFAULT(SHARED)
     !OMP DO PRIVATE(icp,CP,Uindtmp,DP1,DP2) schedule(runtime)
     do icp=1,nCPs
        CP(1:3)=CPs(1:3,icp)
        ! 1-2 segment
        DP1=CP-P1; DP2=CP-P2; 
        call ui_seg_11 ( DP1, DP2, Gamm, SmoothModel, SmoothParam, Uindtmp)
        Uind(1:3,icp) = Uind(1:3,icp)+Uindtmp(1:3)
        ! 3-4 segment
        DP1=CP-P3; DP2=CP-P4; 
        call ui_seg_11 ( DP1, DP2, Gamm, SmoothModel, SmoothParam, Uindtmp)
        Uind(1:3,icp) = Uind(1:3,icp)+Uindtmp(1:3)
        ! 2-3 segment
        DP1=CP-P2; DP2=CP-P3; 
        call ui_seg_11 ( DP1, DP2, Gamm, SmoothModel, SmoothParam, Uindtmp)
        Uind(1:3,icp) = Uind(1:3,icp)+Uindtmp(1:3)
        ! 4-1 segment
        DP1=CP-P4; DP2=CP-P1; 
        call ui_seg_11 ( DP1, DP2, Gamm, SmoothModel, SmoothParam, Uindtmp)
        Uind(1:3,icp) = Uind(1:3,icp)+Uindtmp(1:3)
     end do  ! loop on CPs
     !OMP END DO 
     !OMP END PARALLEL
  end subroutine  ui_quad_n1
  
  end module FVW_BiotSavart
  