anchor.diff:  
      t = rho*A*L;
   
  -   M(1:12,1:12) = 0
  +   M(1:12,1:12) = 0.0_ReKi
   
      M( 1,  1) = 13._ReKi/35._ReKi * t
      M( 2,  2) = 13._ReKi/35._ReKi * t
  
positive.code:    K(8,8)= EE
     K(9,9)= EAL0
  
     DC = 0
     DC( 1: 3,  1: 3) = DirCos
     DC( 4: 6,  4: 6) = DirCos
     DC( 7: 9,  7: 9) = DirCos
  
negative.code:    REAL(ReKi) :: TempCoeff
     REAL(ReKi) :: w            ! weight per unit length
     
     F = 0.0_ReKi      ! initialize whole array to zero, then set the non-zero portions
     w = rho*A*g       ! weight per unit length
     
     ! lumped forces on both nodes (z component only):
     F(3) = -0.5_ReKi*L*w 
     F(9) = F(3)
            
     ! lumped moments on node 1 (x and y components only):
     ! bjj: note that RRD wants factor of 1/12 because of boundary conditions. Our MeshMapping routines use factor of 1/6 (assuming generic/different boundary  
     !      conditions), so we may have some inconsistent behavior. JMJ suggests using line2 elements for SubDyn's input/output meshes to improve the situation.
     TempCoeff = L*L*w/12.0_ReKi ! let's not calculate this twice  
     F(4) = -TempCoeff * DirCos(2,3) ! = -L*w*Dy/12._ReKi   !bjj: DirCos(2,3) = Dy/L
     F(5) =  TempCoeff * DirCos(1,3) ! =  L*w*Dx/12._ReKi   !bjj: DirCos(1,3) = Dx/L
  
        ! lumped moments on node 2: (note the opposite sign of node 1 moment)
     F(10) = -F(4)
  