anchor.diff:          !dq = matmul(jac, -res)
           dq = direct_solve(jac, -res, TOT_NDOF)
   
  -        update_norm  =  maxval(abs(dq));
  -        res_norm     =  maxval(abs(res));
  +        !-------------------------------------------------------------!
  +        ! Calcualte residual tolratances                    
  +        !-------------------------------------------------------------!
  +
  +        update_norm  =  maxval(abs(dq)); res_norm  =  maxval(abs(res));
   
           !-------------------------------------------------------------!
           ! Update the state variables and then update the body
           !-------------------------------------------------------------!
  +
           q            = get_updated_q(q, dq)
           q_dot        = get_updated_q_dot(q_dot, dq)
           ! q_double_dot = get_updated_q_double_dot(q_double_dot, dq) 
  
positive.code:        &  '             KE', '          PE','               TE', &
         &'       Niter',' FCNT'
  
    ! initialize the time tracker
    time = start_time
  
    time_march: do while ( time .le. end_time)  ! loop for time-marching
  
       time = time + dT ! update the time
  
       ! Newton iterations
       newton_cnt = 0
  
       newton: do k = 1, MAX_NEWTON_ITER
  
negative.code:   use utils,only:deg2rad
    use types, only: body, vector
    use solver_utils, only: get_updated_q, get_updated_q_dot,&
         & get_updated_q_double_dot, get_approx_q, get_approx_q_dot
    use bodies, only: create_body, set_state
    use residual, only: get_residual
    use jacobian, only: get_jacobian
  