anchor.diff:     
      
      !.................................
  -   ! y%WriteOutput
  +   ! y%WriteOutput (for writing columns to output file)
      !.................................
      
      call AllocAry( y%WriteOutput, p%numOuts, 'WriteOutput', errStat2, errMsg2 )
  
positive.code:          call cleanup()
           return
        end if
     x%q(:) = 0.0D0
     x%dqdt(:) = 0.0D0
     OtherState%acc(:) = 0.0D0
     OtherState%xcc(:) = 0.0D0
  
  
     ! Define system output initializations (set up mesh) here:
     CALL MeshCreate( BlankMesh        = u%HubMotion            &
                     ,IOS              = COMPONENT_INPUT        &
                     ,NNodes           = 1                      &
                     , TranslationDisp = .TRUE. &
                     , Orientation     = .TRUE. &
                     ,ErrStat         = ErrStat2               &
                     ,ErrMess         = ErrMsg2                )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL MeshPositionNode ( Mesh = u%HubMotion       &
                           , INode = 1                &
                           , Pos = InitInp%HubPos     &
                           , ErrStat   = ErrStat2     &
                           , ErrMess   = ErrMsg2      &
                           , Orient = InitInp%HubRot )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL MeshConstructElement ( Mesh = u%HubMotion            &
                               , Xelement = ELEMENT_POINT      &
                               , P1       = 1                  &
                               , ErrStat  = ErrStat2            &
                               , ErrMess  = ErrMsg2            )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL MeshCommit(u%HubMotion, ErrStat2, ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     
     
     CALL MeshCreate( BlankMesh        = u%RootMotion            &
                     ,IOS              = COMPONENT_INPUT        &
                     ,NNodes           = 1                      &
                     , TranslationDisp = .TRUE. &
                     , TranslationVel  = .TRUE. &
                     , TranslationAcc  = .TRUE. &
                     , Orientation     = .TRUE. &
                     , RotationVel     = .TRUE. &
                     , RotationAcc     = .TRUE. &
                     ,ErrStat         = ErrStat2               &
                     ,ErrMess         = ErrMsg2                )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  
     CALL MeshCreate( BlankMesh        = u%PointLoad            &
                     ,IOS              = COMPONENT_INPUT        &
                     ,NNodes           = p%node_total           &
                     ,Force            = .TRUE. &
                     ,Moment           = .TRUE. &
                     ,ErrStat         = ErrStat2               &
                     ,ErrMess         = ErrMsg2                )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  
     IF(p%quadrature .EQ. 1) THEN
         temp_int = p%ngp * p%elem_total + 2
         CALL MeshCreate( BlankMesh        = u%DistrLoad          &
                         ,IOS              = COMPONENT_INPUT      &
                         ,NNodes           = temp_int             &
                         ,Force            = .TRUE. &
                         ,Moment           = .TRUE. &
                         ,ErrStat         = ErrStat2               &
                         ,ErrMess         = ErrMsg2                )
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     ELSEIF(p%quadrature .EQ. 2) THEN
         temp_int = p%ngp
         CALL MeshCreate( BlankMesh        = u%DistrLoad          &
                         ,IOS              = COMPONENT_INPUT      &
                         ,NNodes           = temp_int             &
                         ,Force            = .TRUE. &
                         ,Moment           = .TRUE. &
                         ,ErrStat         = ErrStat2               &
                         ,ErrMess         = ErrMsg2                )
         
     ENDIF
  
  
     CALL MeshConstructElement ( Mesh = u%RootMotion            &
                               , Xelement = ELEMENT_POINT      &
                               , P1       = 1                  &
                               , ErrStat  = ErrStat2            &
                               , ErrMess  = ErrMsg2            )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  
     DO i=1,p%node_total
         CALL MeshConstructElement( Mesh     = u%PointLoad      &
                                   ,Xelement = ELEMENT_POINT    &
                                   ,P1       = i                &
                                   ,ErrStat  = ErrStat2          &
                                   ,ErrMess  = ErrMsg2           )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         
     ENDDO
  
     IF(p%quadrature .EQ. 1) THEN
         temp_int = p%ngp * p%elem_total + 2
         DO i=1,temp_int-1
             CALL MeshConstructElement( Mesh     = u%DistrLoad      &
                                       ,Xelement = ELEMENT_LINE2    &
                                       ,P1       = i                &
                                       ,P2       = i+1              &
                                       ,ErrStat  = ErrStat2          &
                                       ,ErrMess  = ErrMsg2           )
               CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ENDDO
     ELSEIF(p%quadrature .EQ. 2) THEN
         temp_int = p%ngp
         DO i=1,temp_int - 1
             CALL MeshConstructElement( Mesh     = u%DistrLoad      &
                                       ,Xelement = ELEMENT_LINE2    &
                                       ,P1       = i                &
                                       ,P2       = i+1              &
                                       ,ErrStat  = ErrStat2          &
                                       ,ErrMess  = ErrMsg2           )
               CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ENDDO
     ENDIF
  
     ! place single node at origin; position affects mapping/coupling with other modules
     temp_POS(:) = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(1:3,1))
     TmpPos(1) = temp_POS(2)
     TmpPos(2) = temp_POS(3)
     TmpPos(3) = temp_POS(1)
  
     temp_CRV(1) = p%Glb_crv(2)
     temp_CRV(2) = p%Glb_crv(3)
     temp_CRV(3) = p%Glb_crv(1)
    CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     TmpDCM(:,:) = TRANSPOSE(TmpDCM)
     CALL MeshPositionNode ( Mesh = u%RootMotion      &
                           , INode = 1                &
                           , Pos = TmpPos             &
                           , ErrStat   = ErrStat2     &
                           , ErrMess   = ErrMsg2      &
                           , Orient = TmpDCM ) !bjj: add orient=DCM ! Orientation (direction cosine matrix) of node; identity by default
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  
  
     DO i=1,p%elem_total
         DO j=1,p%node_elem
             temp_id = (j-1) * p%dof_node
             temp_POS(1:3) = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(temp_id+1:temp_id+3,i))
             TmpPos(1) = temp_POS(2)
             TmpPos(2) = temp_POS(3)
             TmpPos(3) = temp_POS(1)
             temp_CRV(:) = MATMUL(p%GlbRot,p%uuN0(temp_id+4:temp_id+6,i))
             CALL BD_CrvCompose(temp_POS,p%Glb_crv,temp_CRV,0,ErrStat2,ErrMsg2)
                CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             temp_CRV(1) = temp_POS(2)
             temp_CRV(2) = temp_POS(3)
             temp_CRV(3) = temp_POS(1)
             CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
                CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             TmpDCM(:,:) = TRANSPOSE(TmpDCM)
             temp_id = (i-1)*(p%node_elem-1)+j
             CALL MeshPositionNode ( Mesh    = u%PointLoad  &
                                    ,INode   = temp_id      &
                                    ,Pos     = TmpPos       &
                                    ,ErrStat = ErrStat2      &
                                    ,ErrMess = ErrMsg2       &
                                    , Orient = TmpDCM ) !bjj: add orient=DCM ! Orientation (direction cosine matrix) of node; identity by default
              CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ENDDO
     ENDDO
  
  
     IF(p%quadrature .EQ. 1) THEN
         DO i=1,p%ngp*p%elem_total+2
             temp_POS(1:3) = p%GlbPos(1:3) + MATMUL(p%GlbRot,temp_L2(1:3,i))
             TmpPos(1) = temp_POS(2)
             TmpPos(2) = temp_POS(3)
             TmpPos(3) = temp_POS(1)
             temp_CRV(:) = MATMUL(p%GlbRot,temp_L2(4:6,i))
             CALL BD_CrvCompose(temp_POS,p%Glb_crv,temp_CRV,0,ErrStat2,ErrMsg2)
                CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             temp_CRV(1) = temp_POS(2)
             temp_CRV(2) = temp_POS(3)
             temp_CRV(3) = temp_POS(1)
             CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
                CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             TmpDCM(:,:) = TRANSPOSE(TmpDCM)
             CALL MeshPositionNode ( Mesh    = u%DistrLoad  &
                                    ,INode   = i            &
                                    ,Pos     = TmpPos       &
                                    ,ErrStat = ErrStat2      &
                                    ,ErrMess = ErrMsg2      &
                                    , Orient = TmpDCM ) !bjj: add orient=DCM ! Orientation (direction cosine matrix) of node; identity by default
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ENDDO
     ELSEIF(p%quadrature .EQ. 2) THEN
         temp_int = p%ngp
         DO i=1,temp_int
             temp_POS(1:3) = p%GlbPos(1:3) + MATMUL(p%GlbRot,temp_L2(1:3,i))
             TmpPos(1) = temp_POS(2)
             TmpPos(2) = temp_POS(3)
             TmpPos(3) = temp_POS(1)
             temp_CRV(:) = MATMUL(p%GlbRot,temp_L2(4:6,i))
             CALL BD_CrvCompose(temp_POS,p%Glb_crv,temp_CRV,0,ErrStat2,ErrMsg2)
                CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             temp_CRV(1) = temp_POS(2)
             temp_CRV(2) = temp_POS(3)
             temp_CRV(3) = temp_POS(1)
             CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
                CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             TmpDCM(:,:) = TRANSPOSE(TmpDCM)
             CALL MeshPositionNode ( Mesh    = u%DistrLoad  &
                                    ,INode   = i            &
                                    ,Pos     = TmpPos       &
                                    ,ErrStat = ErrStat2      &
                                    ,ErrMess = ErrMsg2      &
                                    , Orient = TmpDCM ) !bjj: add orient=DCM ! Orientation (direction cosine matrix) of node; identity by default
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ENDDO
     ENDIF
  
     CALL MeshCommit ( Mesh    = u%RootMotion    &
                      ,ErrStat = ErrStat2         &
                      ,ErrMess = ErrMsg2          )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL MeshCommit ( Mesh    = u%PointLoad     &
                      ,ErrStat = ErrStat2         &
                      ,ErrMess = ErrMsg2          )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL MeshCommit ( Mesh    = u%DistrLoad     &
                      ,ErrStat = ErrStat2         &
                      ,ErrMess = ErrMsg2          )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     
  
  
        if (ErrStat >= AbortErrLev) then
           call cleanup()
           return
        end if
        
     ! Define initialization-routine input here:
     u%HubMotion%TranslationDisp(1:3,1) = 0.0_ReKi
     u%HubMotion%Orientation(1:3,1:3,1) = InitInp%HubRot
  
     u%RootMotion%TranslationDisp(1:3,1) = InitInp%RootDisp(1:3)
     u%RootMotion%Orientation(1:3,1:3,1) = InitInp%RootOri(1:3,1:3)
     u%RootMotion%TranslationVel(1:3,1)  = InitInp%RootVel(1:3)
     u%RootMotion%RotationVel(1:3,1)     = InitInp%RootVel(4:6)
     u%RootMotion%TranslationAcc(:,:)  = 0.0D0
     u%RootMotion%RotationAcc(:,:)     = 0.0D0
  
     DO i=1,u%PointLoad%ElemTable(ELEMENT_POINT)%nelem
         j = u%PointLoad%ElemTable(ELEMENT_POINT)%Elements(i)%ElemNodes(1)
         u%PointLoad%Force(:,j)  = 0.0D0
         u%PointLoad%Moment(:,j) = 0.0D0
     ENDDO
  
     DO i = 1, u%DistrLoad%ElemTable(ELEMENT_LINE2)%nelem
         j = u%DistrLoad%ElemTable(ELEMENT_LINE2)%Elements(i)%ElemNodes(1)
         k = u%DistrLoad%ElemTable(ELEMENT_LINE2)%Elements(i)%ElemNodes(2)
         u%DistrLoad%Force(:,j)  = 0.0D0
         u%DistrLoad%Force(:,k)  = 0.0D0
         u%DistrLoad%Moment(:,j) = 0.0D0
         u%DistrLoad%Moment(:,k) = 0.0D0
     ENDDO
  
     CALL BD_CopyInput(u, u_tmp, MESH_NEWCOPY, ErrStat2, ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL BD_InputGlobalLocal(p,u_tmp,ErrStat2,ErrMsg2)
  
negative.code:       
     
  end subroutine Init_y
  !-----------------------------------------------------------------------------------------------------------------------------------
  subroutine Init_u( temp_L2, InitInp, p, u, ErrStat, ErrMsg )
  ! this routine initializes the inputs, u, that are used in the BeamDyn interface for coupling in the FAST framework.
  
     real(BDKi),                   intent(in   )  :: temp_L2(:,:)      ! 
     type(BD_InitInputType),       intent(in   )  :: InitInp           ! Input data for initialization routine
     type(BD_ParameterType),       intent(in   )  :: p                 ! Parameters  ! intent(out) only because it changes p%NdIndx
     type(BD_InputType),           intent(inout)  :: u                 ! Inputs
     integer(IntKi),               intent(  out)  :: ErrStat           ! Error status of the operation
     character(*),                 intent(  out)  :: ErrMsg            ! Error message if ErrStat /= ErrID_None
     
     
     real(ReKi)                                   :: DCM(3,3)          ! must be same type as mesh orientation fields
     real(ReKi)                                   :: Pos(3)            ! must be same type as mesh position fields
  
     real(BDKi)                                   :: TmpDCM(3,3)
     real(BDKi)                                   :: temp_POS(3)
     real(BDKi)                                   :: temp_CRV(3)
     
     integer(intKi)                               :: temp_id           
     integer(intKi)                               :: i,j               ! loop counters
     integer(intKi)                               :: NNodes            ! number of nodes in mesh
     
     integer(intKi)                               :: ErrStat2          ! temporary Error status
     character(ErrMsgLen)                         :: ErrMsg2           ! temporary Error message
     character(*), parameter                      :: RoutineName = 'Init_u'
     
     ErrStat = ErrID_None
     ErrMsg  = ""
     
     !.................................
     ! u%HubMotion (from ElastoDyn for pitch actuator)
     !.................................
  
     CALL MeshCreate( BlankMesh        = u%HubMotion        &
                     ,IOS              = COMPONENT_INPUT    &
                     ,NNodes           = 1                  &
                     , TranslationDisp = .TRUE.             &
                     , Orientation     = .TRUE.             &
                     ,ErrStat          = ErrStat2           &
                     ,ErrMess          = ErrMsg2            )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        if (ErrStat>=AbortErrLev) return
        
        ! possible type conversions here:
     DCM = InitInp%HubRot 
     Pos = InitInp%HubPos
     CALL MeshPositionNode ( Mesh    = u%HubMotion          &
                           , INode   = 1                    &
                           , Pos     = Pos                  &
                           , ErrStat = ErrStat2             &
                           , ErrMess = ErrMsg2              &
                           , Orient  = InitInp%HubRot       )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        
     CALL MeshConstructElement ( Mesh = u%HubMotion         &
                               , Xelement = ELEMENT_POINT   &
                               , P1       = 1               &
                               , ErrStat  = ErrStat2        &
                               , ErrMess  = ErrMsg2         )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        
     CALL MeshCommit(u%HubMotion, ErrStat2, ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )   
     
        ! initial guesses
     u%HubMotion%TranslationDisp(1:3,1) = 0.0_ReKi
     u%HubMotion%Orientation(1:3,1:3,1) = InitInp%HubRot
        
     !.................................
     ! u%RootMotion (for coupling with ElastoDyn)
     !.................................
  
     CALL MeshCreate( BlankMesh        = u%RootMotion          &
                     ,IOS              = COMPONENT_INPUT       &
                     ,NNodes           = 1                     &
                     , TranslationDisp = .TRUE.                &
                     , TranslationVel  = .TRUE.                &
                     , TranslationAcc  = .TRUE.                &
                     , Orientation     = .TRUE.                &
                     , RotationVel     = .TRUE.                &
                     , RotationAcc     = .TRUE.                &
                     ,ErrStat         = ErrStat2               &
                     ,ErrMess         = ErrMsg2                )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        if (ErrStat>=AbortErrLev) return
              
    !! ! place single node at origin; position affects mapping/coupling with other modules
    !! temp_POS(:) = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(1:3,1))
    !! Pos(1) = temp_POS(2)
    !! Pos(2) = temp_POS(3)
    !! Pos(3) = temp_POS(1)
    !!
    !! temp_CRV(1) = p%Glb_crv(2)
    !! temp_CRV(2) = p%Glb_crv(3)
    !! temp_CRV(3) = p%Glb_crv(1)
    !!CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
    !!    CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
    !! DCM = TRANSPOSE(TmpDCM)
        
     DCM = InitInp%GlbRot 
     Pos = InitInp%GlbPos      
     CALL MeshPositionNode ( Mesh    = u%RootMotion &
                           , INode   = 1            &
                           , Pos     = Pos          &
                           , ErrStat = ErrStat2     &
                           , ErrMess = ErrMsg2      &
                           , Orient  = DCM          ) 
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
           
     CALL MeshConstructElement ( Mesh     = u%RootMotion       &
                               , Xelement = ELEMENT_POINT      &
                               , P1       = 1                  &
                               , ErrStat  = ErrStat2           &
                               , ErrMess  = ErrMsg2            )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     
     
     CALL MeshCommit ( Mesh    = u%RootMotion    &
                      ,ErrStat = ErrStat2        &
                      ,ErrMess = ErrMsg2         )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        
        ! initial guesses
     u%RootMotion%TranslationDisp(1:3,1) = InitInp%RootDisp(1:3)
     u%RootMotion%Orientation(1:3,1:3,1) = InitInp%RootOri(1:3,1:3)
     u%RootMotion%TranslationVel(1:3,1)  = InitInp%RootVel(1:3)
     u%RootMotion%RotationVel(1:3,1)     = InitInp%RootVel(4:6)
     u%RootMotion%TranslationAcc(:,:)    = 0.0_ReKi
     u%RootMotion%RotationAcc(:,:)       = 0.0_ReKi
              
     !.................................
     ! u%PointLoad (currently not used in FAST)
     !.................................
     
     CALL MeshCreate( BlankMesh       = u%PointLoad      &
                     ,IOS             = COMPONENT_INPUT  &
                     ,NNodes          = p%node_total     &
                     ,Force           = .TRUE.           &
                     ,Moment          = .TRUE.           &
                     ,ErrStat         = ErrStat2         &
                     ,ErrMess         = ErrMsg2          )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        if (ErrStat>=AbortErrLev) return
  
     DO i=1,p%elem_total
         DO j=1,p%node_elem
             temp_id = (j-1) * p%dof_node
             temp_POS = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(temp_id+1:temp_id+3,i))
             Pos(1) = temp_POS(2)
             Pos(2) = temp_POS(3)
             Pos(3) = temp_POS(1)
             
             temp_CRV = MATMUL(p%GlbRot,p%uuN0(temp_id+4:temp_id+6,i))
             CALL BD_CrvCompose(temp_POS,p%Glb_crv,temp_CRV,0,ErrStat2,ErrMsg2)
                CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             temp_CRV(1) = temp_POS(2)
             temp_CRV(2) = temp_POS(3)
             temp_CRV(3) = temp_POS(1)
             CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
                CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             DCM = TRANSPOSE(TmpDCM)
             
             temp_id = (i-1)*(p%node_elem-1)+j           
             CALL MeshPositionNode ( Mesh    = u%PointLoad  &
                                    ,INode   = temp_id      &
                                    ,Pos     = Pos          &
                                    ,ErrStat = ErrStat2     &
                                    ,ErrMess = ErrMsg2      &
                                    , Orient = DCM ) 
                 CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ENDDO
     ENDDO      
        
     DO i=1,p%node_total
         CALL MeshConstructElement( Mesh     = u%PointLoad      &
                                   ,Xelement = ELEMENT_POINT    &
                                   ,P1       = i                &
                                   ,ErrStat  = ErrStat2         &
                                   ,ErrMess  = ErrMsg2          )
           CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         
     ENDDO
        
     CALL MeshCommit ( Mesh    = u%PointLoad     &
                      ,ErrStat = ErrStat2        &
                      ,ErrMess = ErrMsg2         )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     
        ! initial guesses
     u%PointLoad%Force  = 0.0_ReKi
     u%PointLoad%Moment = 0.0_ReKi      
        
     !.................................
     ! u%DistrLoad (for coupling with AeroDyn)
     !.................................
              
     IF(p%quadrature .EQ. 1) THEN
         NNodes = p%ngp * p%elem_total + 2
     ELSEIF(p%quadrature .EQ. 2) THEN
        NNodes = p%ngp
     ELSE
        CALL SetErrStat( ErrID_Fatal, "Invalid quadrature selected.", ErrStat, ErrMsg, RoutineName )
        RETURN
     ENDIF
     
     CALL MeshCreate( BlankMesh  = u%DistrLoad      &
                     ,IOS        = COMPONENT_INPUT  &
                     ,NNodes     = NNodes           &
                     ,Force      = .TRUE.           &
                     ,Moment     = .TRUE.           &
                     ,ErrStat    = ErrStat2         &
                     ,ErrMess    = ErrMsg2          )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )   
        if (ErrStat>=AbortErrLev) return
  
     DO i=1,NNodes
        temp_POS(1:3) = p%GlbPos(1:3) + MATMUL(p%GlbRot,temp_L2(1:3,i))
        Pos(1) = temp_POS(2)
        Pos(2) = temp_POS(3)
        Pos(3) = temp_POS(1)
        
        temp_CRV(:) = MATMUL(p%GlbRot,temp_L2(4:6,i))
        CALL BD_CrvCompose(temp_POS,p%Glb_crv,temp_CRV,0,ErrStat2,ErrMsg2)
           CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        temp_CRV(1) = temp_POS(2)
        temp_CRV(2) = temp_POS(3)
        temp_CRV(3) = temp_POS(1)
        CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
           CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        DCM = TRANSPOSE(TmpDCM)
        
        CALL MeshPositionNode ( Mesh    = u%DistrLoad  &
                               ,INode   = i            &
                               ,Pos     = Pos          &
                               ,ErrStat = ErrStat2     &
                               ,ErrMess = ErrMsg2      &
                               ,Orient  = DCM ) 
           CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     ENDDO   
     
     DO i=1,NNodes-1
        CALL MeshConstructElement( Mesh      = u%DistrLoad      &
                                   ,Xelement = ELEMENT_LINE2    &
                                   ,P1       = i                &
                                   ,P2       = i+1              &
                                   ,ErrStat  = ErrStat2         &
                                   ,ErrMess  = ErrMsg2          )
           CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     ENDDO
  
     
     CALL MeshCommit ( Mesh    = u%DistrLoad     &
                      ,ErrStat = ErrStat2        &
                      ,ErrMess = ErrMsg2         )
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     
        ! initial guesses
     u%DistrLoad%Force  = 0.0_ReKi
     u%DistrLoad%Moment = 0.0_ReKi      
     
        
  end subroutine Init_u
  !-----------------------------------------------------------------------------------------------------------------------------------
  SUBROUTINE BD_End( u, p, x, xd, z, OtherState, y, ErrStat, ErrMsg )
     !
  