anchor.diff:              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. 
  
positive.code:    real(ReKi)                                   :: denom                   ! denominator
     real(ReKi)                                   :: v(3)                    ! temp vector
     integer(IntKi)                               :: i                       ! loop counters for points
     real(ReKi)                                   :: MaxDiam                 ! Maximum tower diameter
     real(ReKi)                                   :: TwrClrnc                ! local tower clearance
     integer(intKi)                               :: ErrStat2
     character(ErrMsgLen)                         :: ErrMsg2
     character(*), parameter                      :: RoutineName = 'TwrInflArray'
  
negative.code:    ! 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
  
     !$OMP PARALLEL default(shared)
     !$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  (see getLocalTowerProps)
        ! values are found for the deflected tower, returning theta_tower, W_tower, xbar, ybar, zbar, and TowerCd:
        ! 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.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
  