anchor.diff:        CHARACTER(*),                      INTENT(  OUT)  :: ErrMsg      !< Error message if ErrStat /= ErrID_None
   
         REAL(ReKi)                                        :: F_D(6), F_DP(6), D_F_I(3), kvec(3), v(3),  vf(3), vrel(3), vmag
  -      INTEGER                                           :: I, J, K, nodeIndx
  +      INTEGER                                           :: I, J, K, nodeIndx, IntWrapIndx
         REAL(ReKi)                                        :: elementWaterState
         REAL(ReKi)                                        :: AllOuts(MaxMrsnOutputs)  ! TODO: think about adding to OtherState
         REAL(ReKi)                                        :: qdotdot(6) ,qdotdot2(3)     ! The structural acceleration of a mesh node
  
positive.code:             ! Obtain the node index because WaveVel, WaveAcc, and WaveDynP are defined in the node indexing scheme, not the markers
  
           nodeIndx = p%lumpedToNodeIndx(J)
           nodeInWater = REAL( InterpWrappedStpInt( REAL(Time, SiKi), p%WaveTime(:), p%nodeInWater(:,nodeIndx), m%LastIndWave, p%NStepWave + 1 ), ReKi )
              ! Determine the dynamic pressure at the marker
           m%L_FDynP(J) = InterpWrappedStpReal ( REAL(Time, SiKi), p%WaveTime(:), p%WaveDynP(:,nodeIndx), &
                                      m%LastIndWave, p%NStepWave + 1       )
           
           
           DO I=1,3
                 ! Determine the fluid acceleration and velocity at the marker
              m%L_FA(I,J) = InterpWrappedStpReal ( REAL(Time, SiKi), p%WaveTime(:), p%WaveAcc(:,nodeIndx,I), &
                                      m%LastIndWave, p%NStepWave + 1       )
                 
              m%L_FV(I,J) = InterpWrappedStpReal ( REAL(Time, SiKi), p%WaveTime(:), p%WaveVel(:,nodeIndx,I), &
                                      m%LastIndWave, p%NStepWave + 1       )
              vrel(I)     = m%L_FV(I,J) - u%LumpedMesh%TranslationVel(I,J)
           END DO
           
  
negative.code: 
        
  END SUBROUTINE Morison_UpdateStates
  !> This routine is similar to InterpWrappedStpReal, except it returns only the slope for the interpolation.
  !! By returning the slope based on Time, we don't have to calculate this for every variable (Yary) we want to interpolate.
  FUNCTION GetInterpolationSlope(Time, p, m, IntWrapIndx) RESULT( InterpSlope )
        REAL(DbKi),                        INTENT(IN   )  :: Time        !< Current simulation time in seconds
        TYPE(Morison_ParameterType),       INTENT(IN   )  :: p           !< Parameters
        TYPE(Morison_MiscVarType),         INTENT(INOUT)  :: m           !< Misc/optimization variables            
        INTEGER, OPTIONAL,                 INTENT(  OUT)  :: IntWrapIndx
  
        REAL(SiKi)                                        :: Time_SiKi
        REAL(SiKi)                                        :: TimeMod
        REAL(ReKi)                                        :: InterpSlope
  
        Time_SiKi = REAL(Time, SiKi)
        TimeMod = MOD(Time_SiKi, p%WaveTime(p%NStepWave + 1))
        IF ( TimeMod <= p%WaveTime(2) )  THEN
           m%LastIndWave = 1
        END IF
        
        IF ( TimeMod <= p%WaveTime(1) )  THEN
           m%LastIndWave = 1
           InterpSlope = 0.0_ReKi  ! returns values at m%LastIndWave
           IF(PRESENT(IntWrapIndx)) IntWrapIndx = 1
        ELSE IF ( TimeMod >= p%WaveTime(p%NStepWave + 1) )  THEN
           m%LastIndWave = p%NStepWave
           InterpSlope = 1.0_ReKi  ! returns values at p%NStepWave + 1
           IF(PRESENT(IntWrapIndx)) IntWrapIndx = p%NStepWave + 1
        ELSE
           m%LastIndWave = MAX( MIN( m%LastIndWave, p%NStepWave-1 ), 1 )
  
           DO
  
              IF ( TimeMod < p%WaveTime(m%LastIndWave) )  THEN
  
                 m%LastIndWave = m%LastIndWave - 1
  
              ELSE IF ( TimeMod >= p%WaveTime(m%LastIndWave+1) )  THEN
  
                 m%LastIndWave = m%LastIndWave + 1
  
              ELSE
                 IF(PRESENT(IntWrapIndx)) IntWrapIndx = m%LastIndWave
                 
                 InterpSlope = ( TimeMod - p%WaveTime(m%LastIndWave) )/( p%WaveTime(m%LastIndWave+1) - p%WaveTime(m%LastIndWave) )
                 RETURN ! stop checking DO loop
              END IF
  
           END DO
     
        END IF
        
  END FUNCTION GetInterpolationSlope
  !> Use in conjunction with GetInterpolationSlope, to replace InterpWrappedStpReal here.
  FUNCTION InterpolateWithSlope(InterpSlope, Ind, YAry)
        REAL(ReKi), INTENT(IN)                            :: InterpSlope
        INTEGER(IntKi), INTENT(IN )                       :: Ind           !< Misc/optimization variables
        REAL(ReKi), INTENT(IN)                            :: YAry(:)
        REAL(ReKi)                                        :: InterpolateWithSlope
  
        InterpolateWithSlope = ( YAry(Ind+1) - YAry(Ind) )*InterpSlope + YAry(Ind)
  
  END FUNCTION InterpolateWithSlope
  !----------------------------------------------------------------------------------------------------------------------------------
  !> Routine for computing outputs, used in both loose and tight coupling.
  SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg )   
  