anchor.diff:                   , ErrStat  = ErrStat2         &
                    , ErrMess  = ErrMsg2          )
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +      if (ErrStat>=AbortErrLev) RETURN
            
         ! initialization (not necessary)
         
  
positive.code:    character(ErrMsgLen)                         :: ErrMsg2           ! temporary Error message
     character(*), parameter                      :: RoutineName = 'Init_y'
     
     ErrStat = ErrID_None
     ErrMsg  = ""
     
     
     !.................................
     ! y%BldForce
     !.................................
     CALL MeshCopy ( SrcMesh  = u%PointLoad      &
                   , DestMesh = y%BldForce       &
  
negative.code:    ErrStat = ErrID_None
     ErrMsg  = ""
  
     CALL BD_CrvExtractCrv(u%RootMotion%Orientation(1:3,1:3,1),temp3,ErrStat2,ErrMsg2) !returns temp3
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL BD_CrvCompose(temp_rv,temp3,p%Glb_crv,2,ErrStat2,ErrMsg2)  !returns temp_rv
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL BD_CrvMatrixR(temp_rv,temp_R,ErrStat2,ErrMsg2) !returns temp_R
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        
        
     !Initialize displacements and rotations
     k = 1 !when i=1, k=1
     DO i=1,p%elem_total
        DO j=k,p%node_elem
           temp_id = (j-1)*p%dof_node
           temp_p0 = MATMUL(p%GlbRot,p%uuN0(temp_id+1:temp_id+3,i))
           temp_p0 = MATMUL(temp_R,temp_p0) - temp_p0
           temp_p0 = MATMUL(temp_p0,p%GlbRot) != transpose(MATMUL(transpose(temp_p0),p%GlbRot)) = MATMUL(TRANSPOSE(p%GlbRot),temp_p0)  [transpose of a 1-d array temp_p0 is temp_p0]
           
           temp_id = ((i-1)*(p%node_elem-1)+j-1)*p%dof_node
           x%q(temp_id+1:temp_id+3) = u%RootMotion%TranslationDisp(1:3,1) + temp_p0
        ENDDO
        k = 2 ! start j loop at k=2 for remaining elements (i>1)
     ENDDO
     
     k = 1 !when i=1, k=1
     DO i=1,p%elem_total
        DO j=k,p%node_elem
           temp_id = ((i-1)*(p%node_elem-1)+j-1)*p%dof_node
           x%q(temp_id+4:temp_id+6) = MATMUL(temp_rv,p%GlbRot) != transpose(MATMUL(TRANSPOSE(temp_rv),p%GlbRot) = MATMUL(TRANSPOSE(p%GlbRot),temp_rv) because temp_rv is 1-dimension
        ENDDO
        k = 2 ! start j loop at k=2 for remaining elements (i>1)
     ENDDO
  
     !Initialize velocities and angular velocities
     x%dqdt(:) = 0.0_BDKi
     
     ! these values don't change in the loop:
     GlbRot_TransVel     = MATMUL(p%GlbRot,u%RootMotion%TranslationVel(:,1))
     GlbRot_RotVel_tilde = BD_Tilde(MATMUL(p%GlbRot,u%RootMotion%RotationVel(:,1)))
     k=1 !when i=1, k=1
     DO i=1,p%elem_total
        DO j=k,p%node_elem
           temp_id = (j-1)*p%dof_node
           temp3 = MATMUL(p%GlbRot,p%uuN0(temp_id+1:temp_id+3,i))
           temp3 = GlbRot_TransVel + MATMUL(GlbRot_RotVel_tilde,temp3)
           
           temp_id = ((i-1)*(p%node_elem-1)+j-1)*p%dof_node
           x%dqdt(temp_id+1:temp_id+3) = MATMUL(temp3,p%GlbRot) ! = transpose(MATMUL(transpose(temp3),p%GlbRot)) = MATMUL(TRANSPOSE(p%GlbRot),temp3)  because temp3 is 1-dimension
           x%dqdt(temp_id+4:temp_id+6) = u%RootMotion%RotationVel(1:3,1)
        ENDDO
        k = 2 ! start j loop at k=2 for remaining elements (i>1)
     ENDDO
  
  END SUBROUTINE BD_CalcIC
  