anchor.diff:        do l = 1,np+1
            
          if ((abs(spts(j)-1.).LE.eps).AND.(l.EQ.np+1)) then
  -         dPhis(l,j) = float((np+1)*np)/4.
  +         dPhis(l,j) = REAL((np+1)*np, BDKi)/4.0_BDKi
          elseif ((abs(spts(j)+1.).LE.eps).AND.(l.EQ.1)) then
  -         dPhis(l,j) = -float((np+1)*np)/4.
  +         dPhis(l,j) = -REAL((np+1)*np, BDKi)/4.0_BDKi
          elseif (abs(spts(j)-npts(l)).LE.eps) then
  -         dPhis(l,j) = 0.
  +         dPhis(l,j) = 0.0_BDKi
          else
  -         dPhis(l,j) = 0.
  -         den = 1.          
  +         dPhis(l,j) = 0.0_BDKi
  +         den = 1.0_BDKi          
            do i = 1,np+1
              if (i.NE.l) then
                den = den*(npts(l)-npts(i))
              endif
  -           dnum = 1.
  +           dnum = 1.0_BDKi
              do k = 1,np+1
                if ((k.NE.l).AND.(k.NE.i).AND.(i.NE.l)) then
                  dnum = dnum*(spts(j)-npts(k))
                elseif (i.EQ.l) then
  -               dnum = 0.
  +               dnum = 0.0_BDKi
                endif
              enddo
              dPhis(l,j) = dPhis(l,j) + dnum
  
positive.code:       do l = 1,np+1
           
         if(abs(spts(j)-npts(l)).LE.eps) then
           Ps(l,j) = 1.
         else
           dnum = 1.
           den = 1.
           do k = 1,np+1
             if (k.NE.l) then
               den = den*(npts(l) - npts(k))
               dnum = dnum*(spts(j) - npts(k))
             endif
           enddo
  
negative.code: ! This subroutine computes Global mass matrix and force vector to 
  ! calculate the forces along the beam
  !----------------------------------------------------------------------------------------
     REAL(BDKi),         INTENT(IN   ):: uuN(:)       ! Displacement of Mass 1: m
     REAL(BDKi),         INTENT(IN   ):: vvN(:)       ! Velocity of Mass 1: m/s
     REAL(BDKi),         INTENT(IN   ):: aaN(:)       ! Velocity of Mass 1: m/s
     REAL(BDKi),         INTENT(IN   ):: Stif0(:,:,:) ! Element stiffness matrix
     REAL(BDKi),         INTENT(IN   ):: Mass0(:,:,:) ! Element stiffness matrix
     TYPE(BD_InputType), INTENT(IN   ):: u            ! Inputs at t
     INTEGER(IntKi),     INTENT(IN   ):: damp_flag    ! Number of Gauss points
     REAL(BDKi),         INTENT(IN   ):: beta(:)
     INTEGER(IntKi),     INTENT(IN   ):: elem_total   ! Total number of elements
     INTEGER(IntKi),     INTENT(IN   ):: node_elem    ! Node per element
     INTEGER(IntKi),     INTENT(IN   ):: dof_node     ! Degrees of freedom per node
     INTEGER(IntKi),     INTENT(IN   ):: ngp          ! Number of Gauss points
     INTEGER(IntKi),     INTENT(IN   ):: quadrature   ! Number of Gauss points
     REAL(BDKi),         INTENT(IN   ):: gw(:)
     REAL(BDKi),         INTENT(IN   ):: hhx(:,:)
     REAL(BDKi),         INTENT(IN   ):: hpx(:,:)
     REAL(BDKi),         INTENT(IN   ):: Jaco(:,:)
     REAL(BDKi),         INTENT(IN   ):: uu0(:,:)
     REAL(BDKi),         INTENT(IN   ):: E10(:,:)
     REAL(BDKi),         INTENT(  OUT):: RHS(:)       ! Right hand side of the equation Ax=B
     INTEGER(IntKi),     INTENT(  OUT):: ErrStat      ! Error status of the operation
     CHARACTER(*),       INTENT(  OUT):: ErrMsg       ! Error message if ErrStat /= ErrID_None
  
     REAL(BDKi),           ALLOCATABLE:: Nuuu(:)      ! Nodal displacement of Mass 1 for each element
     REAL(BDKi),           ALLOCATABLE:: Nrrr(:)      ! Nodal rotation parameters for displacement of Mass 1
     REAL(BDKi),           ALLOCATABLE:: Nvvv(:)      ! Nodal velocity of Mass 1: m/s for each element
     REAL(BDKi),           ALLOCATABLE:: Naaa(:)      ! Nodal velocity of Mass 1: m/s for each element
     REAL(BDKi),           ALLOCATABLE:: EStif0_GL(:,:,:)  ! Nodal material properties for each element
     REAL(BDKi),           ALLOCATABLE:: EMass0_GL(:,:,:)  ! Nodal material properties for each element
     REAL(BDKi),           ALLOCATABLE:: DistrLoad_GL(:,:) ! Nodal material properties for each element
     REAL(BDKi),           ALLOCATABLE:: elf(:)       ! Total element force (Fc, Fd, Fb)
     INTEGER(IntKi)                   :: dof_elem     ! Degree of freedom per node
     INTEGER(IntKi)                   :: rot_elem     ! Rotational degrees of freedom
     INTEGER(IntKi)                   :: nelem        ! number of elements
     INTEGER(IntKi)                   :: j            ! Index counter
     INTEGER(IntKi)                   :: temp_id      ! Index counter
     INTEGER(IntKi)                   :: ErrStat2     ! Temporary Error status
     CHARACTER(ErrMsgLen)             :: ErrMsg2      ! Temporary Error message
     CHARACTER(*), PARAMETER          :: RoutineName = 'BD_GenerateDynamicElementForce'
  
     ErrStat    = ErrID_None
  