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:                case ( idSegSmoothRankine ) !!Rankine - t<=>rc
                    ! orthogonal distance r1xr2/r0 
                    h2          = norm2_crossprod/ norm2_r0
                    visc_param2 = visc_param**2
                    if (h2<visc_param2) then 
                       Kv=h2/visc_param2
                    else
                       Kv=1.0_ReKi 
                    end if 
                 case ( idSegSmoothLambOseen ) !!Lamb-Oseen - vsic_param<=>rc 
                    ! orthogonal distance r1xr2/r0 
                    h2          = norm2_crossprod/ norm2_r0
                    visc_param2 = visc_param**2
                    exp_value = -1.25643_ReKi*(h2)/(visc_param2)
                    if(exp_value<MIN_EXP_VALUE) then
                       Kv = 1.0_ReKi
                    else
                       Kv = 1.0_ReKi-exp(exp_value)
                    endif
                 case ( idSegSmoothVatistas ) !!Vatistas n=2 - visc_param<=>rc 
                    ! orthogonal distance r1xr2/r0 
                    h2          = norm2_crossprod/ norm2_r0
                    visc_param2 = visc_param**2
                    ! h = (norm_a+norm_b)/2; 
                    Kv = h2/sqrt(visc_param2**2+h2**2)
                 case ( idSegSmoothCompact ) !!Cut-off radius as function of length: delta^2*norm(r0)^2  
                    Kv=1.0_ReKi
                    denominator=denominator+visc_param*visc_param*norm2_r0
                 case ( idSegSmoothCompactDim ) !!Cut-off radius dimension fixed: (delta l_0)^2  
                    Kv=1.0_ReKi
                    denominator=denominator+visc_param*visc_param
                 case ( 33 ) !!Vatistas n=2 - visc_paramt<=>rc  
                    ! See Hoydonck 2012
                    h           = (norm_a+norm_b)/2.0_ReKi
                    visc_param2 = visc_param**2
                    if(h<2*sqrt(norm2_r0)) then
                       ! use Vatistas n=2, normal one
                       h2 = norm2_crossprod/norm2_r0
                    end if
                    Kv = (h2)/sqrt(visc_param2**2+h2**2)
                 case default
                    Kv=1.0_ReKi
                 end select 
  
negative.code:    real(ReKi), dimension(3) :: DP1, DP2       !<
     real(ReKi)               :: norm2_r0       !<
     real(ReKi)               :: Gam            !<
     real(ReKi)               :: SmoothParam1    !< 
     ! Variables declaration 
     real(ReKi),dimension(3) :: crossprod       !< 
     real(ReKi)              :: denom           !< 
  