anchor.diff:        ELSE IF (p%TMD_DOF_MODE == DOFMode_Omni) THEN
                  ! Compute the first time derivatives of the continuous states of Omnidirectional TMD mode by sm 2015-0904
            do i_pt=1,p%NumMeshPts
  -            dxdt%tmd_x(2,i_pt) =  ( omega_P_N(2)**2 + omega_P_N(3)**2 - K(1) / p%M_XY) * x%tmd_x(1,i_pt) &
  -                           - ( m%C_ctrl( 1)/p%M_XY ) * x%tmd_x(2,i_pt)                                &
  -                           - ( m%C_Brake(1)/p%M_XY ) * x%tmd_x(2,i_pt)                                &
  -                           + B_X + 1 / p%M_XY * ( m%F_fr(1) )                                      &
  -                           - ( omega_P_N(1)*omega_P_N(2) - alpha_P_N(3) ) * x%tmd_x(3,i_pt)           &
  -                          +2 * omega_P_N(3) * x%tmd_x(4,i_pt)
  -            dxdt%tmd_x(4,i_pt) =  ( omega_P_N(1)**2 + omega_P_N(3)**2 - K(2) / p%M_XY) * x%tmd_x(3,i_pt) &
  -                           - ( m%C_ctrl( 2)/p%M_XY ) * x%tmd_x(4,i_pt)                                &
  -                           - ( m%C_Brake(2)/p%M_XY ) * x%tmd_x(4,i_pt)                                &
  -                           + B_Y + 1 / p%M_XY * ( m%F_fr(2) )                                      &
  -                           - ( omega_P_N(1)*omega_P_N(2) + alpha_P_N(3) ) * x%tmd_x(1,i_pt)           &
  -                          -2 * omega_P_N(3) * x%tmd_x(2,i_pt)
  +            dxdt%tmd_x(2,i_pt) =  ( omega_P(2,i_pt)**2 + omega_P(3,i_pt)**2 - K(1) / p%M_XY) * x%tmd_x(1,i_pt) &
  +                                - ( m%C_ctrl( 1)/p%M_XY ) * x%tmd_x(2,i_pt)                                    &
  +                                - ( m%C_Brake(1)/p%M_XY ) * x%tmd_x(2,i_pt)                                    &
  +                                +  B_X(i_pt) + 1/p%M_XY * ( m%F_fr(1) )                                        &
  +                                - ( omega_P(1,i_pt)*omega_P(2,i_pt) - alpha_P(3,i_pt) ) * x%tmd_x(3,i_pt)      &
  +                               +2 * omega_P(3,i_pt) * x%tmd_x(4,i_pt)
  +            dxdt%tmd_x(4,i_pt) =  ( omega_P(1,i_pt)**2 + omega_P(3,i_pt)**2 - K(2) / p%M_XY) * x%tmd_x(3,i_pt) &
  +                                - ( m%C_ctrl( 2)/p%M_XY ) * x%tmd_x(4,i_pt)                                    &
  +                                - ( m%C_Brake(2)/p%M_XY ) * x%tmd_x(4,i_pt)                                    &
  +                                +  B_Y(i_pt) + 1/p%M_XY * ( m%F_fr(2) )                                        &
  +                                - ( omega_P(1,i_pt)*omega_P(2,i_pt) + alpha_P(3,i_pt) ) * x%tmd_x(1,i_pt)      &
  +                               -2 * omega_P(3,i_pt) * x%tmd_x(2,i_pt)
            enddo
   
         ELSE IF (p%TMD_DOF_MODE == DOFMode_TLCD) THEN !MEG & SP
            ! Compute the first time derivatives of the continuous states of TLCD mode
            do i_pt=1,p%NumMeshPts
  -            dxdt%tmd_x(2,i_pt) = (2*p%rho_FA*p%area_FA*x%tmd_x(1,i_pt)*rddot_N_N(3)  &
  -                              +p%rho_FA*p%area_FA*p%B_FA*alpha_P_N(2)*((p%L_FA-p%B_FA)/2) &
  -                              -p%rho_FA*p%area_FA*p%B_FA*omega_P_N(1)*omega_P_N(3)*((p%L_FA-p%B_FA)/2)   &
  -                            +2*p%rho_FA*p%area_FA*omega_P_N(1)*omega_P_N(1)*x%tmd_x(1,i_pt)*(p%L_FA-p%B_FA)   &
  -                            +2*p%rho_FA*p%area_FA*omega_P_N(2)*omega_P_N(2)*x%tmd_x(1,i_pt)*(p%L_FA-p%B_FA)   &
  -                            +2*p%rho_FA*p%area_FA*x%tmd_x(1,i_pt)*a_G_N(3)   &
  -                              -p%rho_FA*p%area_FA*p%B_FA*rddot_N_N(1)   &
  -                              +p%rho_FA*p%area_FA*p%B_FA*a_G_N(1) &
  -                           -.5*p%rho_FA*p%area_FA*p%headLossCoeff_FA*p%area_ratio_FA*p%area_ratio_FA*x%tmd_x(2,i_pt)*ABS(x%tmd_x(2,i_pt)))/(p%rho_FA*p%area_FA*(p%L_FA-p%B_FA+p%area_ratio_FA*p%B_FA))
  -            dxdt%tmd_x(4,i_pt) = (2*p%rho_SS*p%area_SS*x%tmd_x(3,i_pt)*rddot_N_N(3)  &
  -                              +p%rho_SS*p%area_SS*p%B_SS*alpha_P_N(1)*((p%L_SS-p%B_SS)/2) &
  -                              -p%rho_SS*p%area_SS*p%B_SS*omega_P_N(2)*omega_P_N(3)*((p%L_SS-p%B_SS)/2)   &
  -                            +2*p%rho_SS*p%area_SS*x%tmd_x(3,i_pt)*omega_P_N(1)*omega_P_N(1)*(p%L_SS-p%B_SS)   &
  -                            +2*p%rho_SS*p%area_SS*x%tmd_x(3,i_pt)*omega_P_N(2)*omega_P_N(2)*(p%L_SS-p%B_SS)   &
  -                            +2*p%rho_SS*p%area_SS*x%tmd_x(3,i_pt)*a_G_N(3)-p%rho_SS*p%area_SS*p%B_SS*rddot_N_N(2)   &
  -                              +p%rho_SS*p%area_SS*p%B_SS*a_G_N(2)    &
  -                           -.5*p%rho_SS*p%area_SS*p%headLossCoeff_SS*p%area_ratio_SS*p%area_ratio_SS*x%tmd_x(4,i_pt)*ABS(x%tmd_x(4,i_pt)))/(p%rho_SS*p%area_SS*(p%L_SS-p%B_SS+p%area_ratio_SS*p%B_SS))
  +            dxdt%tmd_x(2,i_pt) = (2*p%rho_FA*p%area_FA*x%tmd_x(1,i_pt)*rddot_P(3,i_pt)                                       &
  +                                   +p%rho_FA*p%area_FA*p%B_FA*alpha_P(2,i_pt)*((p%L_FA-p%B_FA)/2)                            &
  +                                   -p%rho_FA*p%area_FA*p%B_FA*omega_P(1,i_pt)*omega_P(3,i_pt)*((p%L_FA-p%B_FA)/2)            &
  +                                 +2*p%rho_FA*p%area_FA*omega_P(1,i_pt)*omega_P(1,i_pt)*x%tmd_x(1,i_pt)*(p%L_FA-p%B_FA)       &
  +                                 +2*p%rho_FA*p%area_FA*omega_P(2,i_pt)*omega_P(2,i_pt)*x%tmd_x(1,i_pt)*(p%L_FA-p%B_FA)       &
  +                                 +2*p%rho_FA*p%area_FA*x%tmd_x(1,i_pt)*a_G(3,i_pt)                                           &
  +                                   -p%rho_FA*p%area_FA*p%B_FA*rddot_P(1,i_pt)                                                &
  +                                   +p%rho_FA*p%area_FA*p%B_FA*a_G(1,i_pt)                                                    &
  +                                -.5*p%rho_FA*p%area_FA*p%headLossCoeff_FA*p%area_ratio_FA*p%area_ratio_FA*x%tmd_x(2,i_pt)    &
  +                                       *ABS(x%tmd_x(2,i_pt)))/(p%rho_FA*p%area_FA*(p%L_FA-p%B_FA+p%area_ratio_FA*p%B_FA))
  +            dxdt%tmd_x(4,i_pt) = (2*p%rho_SS*p%area_SS*x%tmd_x(3,i_pt)*rddot_P(3,i_pt)                                       &
  +                                   +p%rho_SS*p%area_SS*p%B_SS*alpha_P(1,i_pt)*((p%L_SS-p%B_SS)/2)                            &
  +                                   -p%rho_SS*p%area_SS*p%B_SS*omega_P(2,i_pt)*omega_P(3,i_pt)*((p%L_SS-p%B_SS)/2)            &
  +                                 +2*p%rho_SS*p%area_SS*x%tmd_x(3,i_pt)*omega_P(1,i_pt)*omega_P(1,i_pt)*(p%L_SS-p%B_SS)       &
  +                                 +2*p%rho_SS*p%area_SS*x%tmd_x(3,i_pt)*omega_P(2,i_pt)*omega_P(2,i_pt)*(p%L_SS-p%B_SS)       &
  +                                 +2*p%rho_SS*p%area_SS*x%tmd_x(3,i_pt)*a_G(3,i_pt)-p%rho_SS*p%area_SS*p%B_SS*rddot_P(2,i_pt) &
  +                                   +p%rho_SS*p%area_SS*p%B_SS*a_G(2,i_pt)                                                    &
  +                                -.5*p%rho_SS*p%area_SS*p%headLossCoeff_SS*p%area_ratio_SS*p%area_ratio_SS*x%tmd_x(4,i_pt)    &
  +                                       *ABS(x%tmd_x(4,i_pt)))/(p%rho_SS*p%area_SS*(p%L_SS-p%B_SS+p%area_ratio_SS*p%B_SS))
            enddo
   
         END IF
  
positive.code:             END IF
           enddo
  
           B_X = - rddot_N_N(1) + a_G_N(1) + 1 / p%M_XY * ( m%F_ext(1) + m%F_stop(1) - m%F_table(1)*(m%F_k_x) )
           B_Y = - rddot_N_N(2) + a_G_N(2) + 1 / p%M_XY * ( m%F_ext(2) + m%F_stop(2) - m%F_table(2)*(m%F_k_y) )
  
        ENDIF
        IF (p%TMD_On_Blade) THEN
           do i_pt=1,p%NumMeshPts
           ! Compute inputs
              B_X_B(i_pt) = - rddot_B_B(1,i_pt) + a_G_B(1,i_pt) + 1 / p%M_X * ( m%F_ext(1) + m%F_stop(1) - m%F_table(1) )
              B_Y_B(i_pt) = - rddot_B_B(2,i_pt) + a_G_B(2,i_pt) + 1 / p%M_Y * ( m%F_ext(2) + m%F_stop(2) - m%F_table(2) )
           enddo
        END IF
  
  
negative.code: 
           ! tmd external forces of dependent degrees:
           do i_pt=1,p%NumMeshPts
              F_x_tmdY_P(i_pt) = - p%M_Y * (a_G(1,i_pt) - rddot_P(1,i_pt) + (alpha_P(3,i_pt) - omega_P(1,i_pt)*omega_P(2,i_pt))*x%tmd_x(3,i_pt) + 2*omega_P(3,i_pt)*x%tmd_x(4,i_pt))
              F_z_tmdY_P(i_pt) = - p%M_Y * (a_G(3,i_pt) - rddot_P(3,i_pt) - (alpha_P(1,i_pt) + omega_P(2,i_pt)*omega_P(3,i_pt))*x%tmd_x(3,i_pt) - 2*omega_P(1,i_pt)*x%tmd_x(4,i_pt))
  
              F_y_tmdX_P(i_pt) = - p%M_X *( a_G(2,i_pt) - rddot_P(2,i_pt) - (alpha_P(3,i_pt) + omega_P(1,i_pt)*omega_P(2,i_pt))*x%tmd_x(1,i_pt) - 2*omega_P(3,i_pt)*x%tmd_x(2,i_pt))
              F_z_tmdX_P(i_pt) = - p%M_X * (a_G(3,i_pt) - rddot_P(3,i_pt) + (alpha_P(2,i_pt) - omega_P(1,i_pt)*omega_P(3,i_pt))*x%tmd_x(1,i_pt) + 2*omega_P(2,i_pt)*x%tmd_x(2,i_pt))
           enddo
  
           ! forces in local coordinates
           do i_pt=1,p%NumMeshPts
              F_P(1,i_pt) =  p%K_X * x%tmd_x(1,i_pt) + m%C_ctrl(1) * x%tmd_x(2,i_pt) + m%C_Brake(1) * x%tmd_x(2,i_pt) - m%F_stop(1) - m%F_ext(1) - m%F_fr(1) - F_x_tmdY_P(i_pt) + m%F_table(1)
              F_P(2,i_pt) =  p%K_Y * x%tmd_x(3,i_pt) + m%C_ctrl(2) * x%tmd_x(4,i_pt) + m%C_Brake(2) * x%tmd_x(4,i_pt) - m%F_stop(2) - m%F_ext(2) - m%F_fr(2) - F_y_tmdX_P(i_pt) + m%F_table(2)
              F_P(3,i_pt) = - F_z_tmdX_P(i_pt) - F_z_tmdY_P(i_pt)
           enddo
  
           ! inertial contributions from mass of TMDs and acceleration of nacelle
           ! Moments on nacelle in local coordinates
           do i_pt=1,p%NumMeshPts
              M_P(1,i_pt) =  - F_z_tmdY_P(i_pt)  * x%tmd_x(3,i_pt)
              M_P(2,i_pt) =    F_z_tmdX_P(i_pt)  * x%tmd_x(1,i_pt)
              M_P(3,i_pt) = (- F_x_tmdY_P(i_pt)) * x%tmd_x(3,i_pt) + (F_y_tmdX_P(i_pt)) * x%tmd_x(1,i_pt)
           enddo
  
           do i_pt=1,p%NumMeshPts
              ! forces in global coordinates
              y%Mesh(i_pt)%Force(:,1) =  matmul(transpose(u%Mesh(i_pt)%Orientation(:,:,1)),F_P(1:3,i_pt))
              ! moments in global coordinates
              y%Mesh(i_pt)%Moment(:,1) = matmul(transpose(u%Mesh(i_pt)%Orientation(:,:,1)),M_P(1:3,i_pt))
           enddo
  
        ELSE IF (p%TMD_DOF_MODE == DOFMode_Omni) THEN
  