anchor.diff:     ! these models are valid for only small tower deflections; check for potential division-by-zero errors:   
      call CheckTwrInfl( u, ErrStat2, ErrMsg2 ); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ); if (ErrStat >= AbortErrLev) return
   
  -   MaxDiam = maxval(p%TwrDiam(:))
  -
      !$OMP PARALLEL default(shared)
  -   !$OMP do private(i,Pos,theta_tower_trans,W_tower,xbar,ybar,zbar,TwrCd,TwrClrnc,denom,u_TwrPotent,v_TwrPotent,u_TwrShadow,v) schedule(runtime)
  +   !$OMP do private(i,Pos,r_TowerBlade,theta_tower_trans,W_tower,xbar,ybar,zbar,TwrCd,TwrClrnc,TwrDiam,found,denom,u_TwrPotent,v_TwrPotent,u_TwrShadow,v) schedule(runtime)
      do i = 1, size(Positions,2)
         Pos=Positions(1:3,i)
            
  -      ! Find nearest line2 element or node of the tower 
  +      ! Find nearest line2 element or node of the tower  (see getLocalTowerProps)
         ! values are found for the deflected tower, returning theta_tower, W_tower, xbar, ybar, zbar, and TowerCd:
  -      call getLocalTowerPropsNoCheck(p, u, Pos, theta_tower_trans, W_tower, xbar, ybar, zbar, TwrCd, TwrClrnc)
  -
  -      if ( TwrClrnc>20*MaxDiam) then
  +      ! option 1: nearest line2 element
  +      call TwrInfl_NearestLine2Element(p, u, Pos, r_TowerBlade, theta_tower_trans, W_tower, xbar, ybar, zbar, TwrCd, TwrDiam, found)
  +      if ( .not. found) then 
  +         ! option 2: nearest node
  +         call TwrInfl_NearestPoint(p, u, Pos, r_TowerBlade, theta_tower_trans, W_tower, xbar, ybar, zbar, TwrCd, TwrDiam)
  +      end if
  +      TwrClrnc = TwoNorm(r_TowerBlade) - 0.5_ReKi*TwrDiam
  +
  +      if ( TwrClrnc>20*TwrDiam) then
            ! Far away, we skip the computation and keep undisturbed inflow 
  -      elseif ( TwrClrnc<=0.0_ReKi) then
  -         ! Inside the tower, will happen for vortex elements, we keep undisturbed inflow
  +      elseif ( TwrClrnc<=0.01_ReKi*TwrDiam) then
  +         ! Inside the tower, or very close, (will happen for vortex elements) we keep undisturbed inflow
  +         ! We don't want to reach the stagnation points
         else
            ! calculate tower influence:
            if ( abs(zbar) < 1.0_ReKi .and. p%TwrPotent /= TwrPotent_none ) then
  
positive.code:             u_TwrShadow = 0.0_ReKi
           end if
                       
           ! NOTE: we should introduce a safety check to avoid zero or negative velocities, have a minimal velocity
           v(1) = (u_TwrPotent + u_TwrShadow)*W_tower
           v(2) = v_TwrPotent*W_tower
           v(3) = 0.0_ReKi
  
           
           Inflow(1:3,i) = Inflow(1:3,i) + matmul( theta_tower_trans, v ) 
        endif ! Check if point far away or in tower
     enddo ! loop on points
     !$OMP END DO 
     !$OMP END PARALLEL
  contains 
     !> This is a simpler version of getLocalTowerProps
     SUBROUTINE getLocalTowerPropsNoCheck(p, u, BladeNodePosition, theta_tower_trans, W_tower, xbar, ybar, zbar, TwrCd, TwrClrnc)
        TYPE(AD_InputType),           INTENT(IN   )  :: u                       !< Inputs at Time t
        TYPE(AD_ParameterType),       INTENT(IN   )  :: p                       !< Parameters
        REAL(ReKi)                   ,INTENT(IN   )  :: BladeNodePosition(3)    !< local blade node position
        REAL(ReKi)                   ,INTENT(  OUT)  :: theta_tower_trans(3,3)  !< transpose of local tower orientation expressed as a DCM
        REAL(ReKi)                   ,INTENT(  OUT)  :: W_tower                 !< local relative wind speed normal to the tower
        REAL(ReKi)                   ,INTENT(  OUT)  :: xbar                    !< local x^ component of r_TowerBlade normalized by tower radius
        REAL(ReKi)                   ,INTENT(  OUT)  :: ybar                    !< local y^ component of r_TowerBlade normalized by tower radius
        REAL(ReKi)                   ,INTENT(  OUT)  :: zbar                    !< local z^ component of r_TowerBlade normalized by tower radius
        REAL(ReKi)                   ,INTENT(  OUT)  :: TwrCd                   !< local tower drag coefficient
        REAL(ReKi)                   ,INTENT(  OUT)  :: TwrClrnc                !< tower clearance for potential output 
        ! local variables
        real(ReKi)                                   :: r_TowerBlade(3)         ! distance vector from tower to blade
        real(ReKi)                                   :: TwrDiam                 ! local tower diameter  
        logical                                      :: found   
        character(*), parameter                      :: RoutineName = 'getLocalTowerPropsNoCheck'
        ! option 1: nearest line2 element
        call TwrInfl_NearestLine2Element(p, u, BladeNodePosition, r_TowerBlade, theta_tower_trans, W_tower, xbar, ybar, zbar, TwrCd, TwrDiam, found)
        if ( .not. found) then 
           ! option 2: nearest node
           call TwrInfl_NearestPoint(p, u, BladeNodePosition, r_TowerBlade, theta_tower_trans, W_tower, xbar, ybar, zbar, TwrCd, TwrDiam)
        end if
        TwrClrnc = TwoNorm(r_TowerBlade) - 0.5_ReKi*TwrDiam
     END SUBROUTINE getLocalTowerPropsNoCheck
  END SUBROUTINE TwrInflArray
  !----------------------------------------------------------------------------------------------------------------------------------
  !> This routine returns the tower constants necessary to compute the tower influence. 
  
negative.code:    real(ReKi)                                   :: denom                   ! denominator
     real(ReKi)                                   :: v(3)                    ! temp vector
     integer(IntKi)                               :: i                       ! loop counters for points
     real(ReKi)                                   :: TwrClrnc                ! local tower clearance
     real(ReKi)                                   :: r_TowerBlade(3)         ! distance vector from tower to blade
     real(ReKi)                                   :: TwrDiam                 ! local tower diameter  
     logical                                      :: found   
     integer(intKi)                               :: ErrStat2
     character(ErrMsgLen)                         :: ErrMsg2
     character(*), parameter                      :: RoutineName = 'TwrInflArray'
  