anchor.diff:        IF (ErrStat >= AbortErrLev) RETURN
   
      p%DT  = Interval
  -   p%Gravity = InitInp%Gravity      ! Gravity vector pointed in negative (/0,0,-g/)
  +   p%Gravity = InitInp%Gravity      ! Gravity vector pointed in negative global Z-axis (/0,0,-g/)
   
         !............................................................................................
         ! Define initial system states here:
  
positive.code:       ErrMsg  = ""
  
  !FIXME: move this to miscvars, or states.
  !      call AllocAry(a_G_N,       3, p%NumMeshPts,'a_G_N',         ErrStat2,ErrMsg2); if (Failed()) return;
  !      call AllocAry(F_P_N,       3, p%NumMeshPts,'F_P_N',         ErrStat2,ErrMsg2); if (Failed()) return;
  !      call AllocAry(M_P_N,       3, p%NumMeshPts,'M_P_N',         ErrStat2,ErrMsg2); if (Failed()) return;
  
  !FIXME: move this to miscvars, or states.
        call AllocAry(a_G_B,       3, p%NumMeshPts,'a_G_B',         ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(F_P_B,       3, p%NumMeshPts,'F_P_B',         ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(M_P_B,       3, p%NumMeshPts,'M_P_B',         ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(r_ddot_P_B  ,3, p%NumMeshPts,'r_ddot_P_B' ,   ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(omega_B_O_B ,3, p%NumMeshPts,'omega_B_O_B',   ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(alpha_B_O_B ,3, p%NumMeshPts,'alpha_B_O_B',   ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(F_x_btmdY_P_B,  p%NumMeshPts,'F_x_btmdY_P_B', ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(F_z_btmdY_P_B,  p%NumMeshPts,'F_z_btmdY_P_B', ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(F_y_btmdX_P_B,  p%NumMeshPts,'F_y_btmdX_P_B', ErrStat2,ErrMsg2); if (Failed()) return;
        call AllocAry(F_z_btmdX_P_B,  p%NumMeshPts,'F_z_btmdX_P_B', ErrStat2,ErrMsg2); if (Failed()) return;
  
  
        IF (.not. p%TMD_On_Blade) THEN
           i_pt=1
              ! Compute nacelle and gravitational acceleration in nacelle coordinates
           a_G_N       = matmul(u%Mesh(i_pt)%Orientation(:,:,1),p%Gravity)
           r_ddot_P_N  = matmul(u%Mesh(i_pt)%Orientation(:,:,1),u%Mesh(i_pt)%TranslationAcc(:,1))
           omega_N_O_N = matmul(u%Mesh(i_pt)%Orientation(:,:,1),u%Mesh(i_pt)%RotationVel(:,1))
           alpha_N_O_N = matmul(u%Mesh(i_pt)%Orientation(:,:,1),u%Mesh(i_pt)%RotationAcc(:,1))
  
        ELSE IF (p%TMD_On_Blade) THEN
  
           do i_pt=1,p%NumMeshPts
              a_G_B(:,i_pt)       = matmul(u%Mesh(i_pt)%Orientation(:,:,1),p%Gravity)
              r_ddot_P_B(:,i_pt)  = matmul(u%Mesh(i_pt)%Orientation(:,:,1),u%Mesh(i_pt)%TranslationAcc(:,1))
              omega_B_O_B(:,i_pt) = matmul(u%Mesh(i_pt)%Orientation(:,:,1),u%Mesh(i_pt)%RotationVel(:,1))
              alpha_B_O_B(:,i_pt) = matmul(u%Mesh(i_pt)%Orientation(:,:,1),u%Mesh(i_pt)%RotationAcc(:,1))
           enddo
  
        ENDIF
  
negative.code:             END IF
           enddo
  
           do i_pt=1,p%NumMeshPts
              B_X(i_pt) = - rddot_P(1,i_pt) + a_G(1,i_pt) + 1 / p%M_XY * ( m%F_ext(1) + m%F_stop(1) - m%F_table(1)*(m%F_k_x) )
              B_Y(i_pt) = - rddot_P(2,i_pt) + a_G(2,i_pt) + 1 / p%M_XY * ( m%F_ext(2) + m%F_stop(2) - m%F_table(2)*(m%F_k_y) )
           enddo
  
        ENDIF
        IF (p%TMD_On_Blade) THEN
           do i_pt=1,p%NumMeshPts
           ! Compute inputs
              B_X(i_pt) = - rddot_P(1,i_pt) + a_G(1,i_pt) + 1 / p%M_X * ( m%F_ext(1) + m%F_stop(1) - m%F_table(1) )
              B_Y(i_pt) = - rddot_P(2,i_pt) + a_G(2,i_pt) + 1 / p%M_Y * ( m%F_ext(2) + m%F_stop(2) - m%F_table(2) )
           enddo
        END IF
  
  