anchor.diff:     INTEGER(IntKi)                :: i, indx_last
      REAL(ReKi)                    :: smallAngles(3)
      REAL(R8Ki)                    :: orientation(3,3)
  -
  -   
  -   do i=1,M_p%NNodes
  -      indx_last  = indx_first + 2 
  -      dY(indx_first:indx_last) = M_p%TranslationDisp(:,i) - M_m%TranslationDisp(:,i)
  -      indx_first = indx_last + 1
  -   end do
  +   LOGICAL                       :: Mask(FIELDMASK_SIZE)               !< flags to determine if this field is part of the packing
  +
  +   if (present(FieldMask)) then
  +      Mask = FieldMask
  +   else
  +      Mask = .true.
  +   end if
  +
  +   
  +   if (Mask(MASKID_TRANSLATIONDISP)) then
  +      do i=1,M_p%NNodes
  +         indx_last  = indx_first + 2 
  +         dY(indx_first:indx_last) = M_p%TranslationDisp(:,i) - M_m%TranslationDisp(:,i)
  +         indx_first = indx_last + 1
  +      end do
  +   end if
  +   
  +   if (Mask(MASKID_ORIENTATION)) then
  +      do i=1,M_p%NNodes
  +         orientation = transpose(M_m%Orientation(:,:,i))
  +         orientation = matmul(orientation, M_p%Orientation(:,:,i))
  +            
  +         smallAngles = GetSmllRotAngs( orientation, ErrStat2, ErrMsg2 )
  +
  +         indx_last  = indx_first + 2 
  +         dY(indx_first:indx_last) = smallAngles
  +         indx_first = indx_last + 1
  +      end do
  +   end if
  +      
  +   if (Mask(MASKID_TRANSLATIONVEL)) then
  +      do i=1,M_p%NNodes
  +         indx_last  = indx_first + 2 
  +         dY(indx_first:indx_last) = M_p%TranslationVel(:,i) - M_m%TranslationVel(:,i)
  +         indx_first = indx_last + 1
  +      end do         
  +   end if
  +      
  +   if (Mask(MASKID_ROTATIONVEL)) then
  +      do i=1,M_p%NNodes
  +         indx_last  = indx_first + 2 
  +         dY(indx_first:indx_last) = M_p%RotationVel(:,i) - M_m%RotationVel(:,i)
  +         indx_first = indx_last + 1
  +      end do          
  +   end if
            
  -   do i=1,M_p%NNodes
  -      orientation = transpose(M_m%Orientation(:,:,i))
  -      orientation = matmul(orientation, M_p%Orientation(:,:,i))
  -            
  -      smallAngles = GetSmllRotAngs( orientation, ErrStat2, ErrMsg2 )
  -
  -      indx_last  = indx_first + 2 
  -      dY(indx_first:indx_last) = smallAngles
  -      indx_first = indx_last + 1
  -   end do
  -                  
  -   do i=1,M_p%NNodes
  -      indx_last  = indx_first + 2 
  -      dY(indx_first:indx_last) = M_p%TranslationVel(:,i) - M_m%TranslationVel(:,i)
  -      indx_first = indx_last + 1
  -   end do         
  -         
  -   do i=1,M_p%NNodes
  -      indx_last  = indx_first + 2 
  -      dY(indx_first:indx_last) = M_p%RotationVel(:,i) - M_m%RotationVel(:,i)
  -      indx_first = indx_last + 1
  -   end do          
  -   
  -   !...............................
  -   ! if we want to skip the acceleration fields, lets return early:
  -   if (present(SkipAcc)) then
  -      if (SkipAcc) return
  +   if (Mask(MASKID_TRANSLATIONACC)) then
  +      do i=1,M_p%NNodes
  +         indx_last  = indx_first + 2 
  +         dY(indx_first:indx_last) = M_p%TranslationAcc(:,i) - M_m%TranslationAcc(:,i)
  +         indx_first = indx_last + 1
  +      end do         
      end if
  -   !...............................
  -   
  -   do i=1,M_p%NNodes
  -      indx_last  = indx_first + 2 
  -      dY(indx_first:indx_last) = M_p%TranslationAcc(:,i) - M_m%TranslationAcc(:,i)
  -      indx_first = indx_last + 1
  -   end do         
  -         
  -   do i=1,M_p%NNodes
  -      indx_last  = indx_first + 2 
  -      dY(indx_first:indx_last) = M_p%RotationAcc(:,i) - M_m%RotationAcc(:,i)
  -      indx_first = indx_last + 1
  -   end do      
  +   
  +   if (Mask(MASKID_ROTATIONACC)) then
  +      do i=1,M_p%NNodes
  +         indx_last  = indx_first + 2 
  +         dY(indx_first:indx_last) = M_p%RotationAcc(:,i) - M_m%RotationAcc(:,i)
  +         indx_first = indx_last + 1
  +      end do      
  +   end if         
  +                                             
   
   END SUBROUTINE PackMotionMesh_dY
   !----------------------------------------------------------------------------------------------------------------------------------
  -!> This subroutine returns the names of the output rows/columns in the Jacobian matrices. Make sure it is consistant with
  -!! elastodyn::packmotionmesh_dy.
  -SUBROUTINE PackMotionMesh_Names(M, MeshName, Names, indx_first, SkipAcc)
  -   
  -   TYPE(MeshType)                    , INTENT(IN   ) :: M                          !< ED output mesh
  -   CHARACTER(*)                      , INTENT(IN   ) :: MeshName                   !< name of mesh 
  -   CHARACTER(LinChanLen)             , INTENT(INOUT) :: Names(:)                   !< name of column of dYdu 
  -   INTEGER(IntKi)                    , INTENT(INOUT) :: indx_first                 !< index into Names array; gives location of next array position to fill
  -   LOGICAL, OPTIONAL                 , INTENT(IN   ) :: SkipAcc                    !< flag that allows us to skip the acceleration fields
  -   
  -      ! local variables:
  -   INTEGER(IntKi)                :: i, j
  -   character(1), parameter       :: Comp(3) = (/'X','Y','Z'/)
  -
  -   
  -   do i=1,M%NNodes
  -      do j=1,3
  -         Names(indx_first) = trim(MeshName)//Comp(j)//' translation displacement, node '//trim(num2lstr(i))
  -         indx_first = indx_first + 1
  -      end do      
  -   end do
  -         
  -   do i=1,M%NNodes
  -      do j=1,3
  -         Names(indx_first) = trim(MeshName)//Comp(j)//' orientation angle, node '//trim(num2lstr(i))
  -         indx_first = indx_first + 1
  -      end do      
  -   end do
  -   
  -   do i=1,M%NNodes
  -      do j=1,3
  -         Names(indx_first) = trim(MeshName)//Comp(j)//' translation velocity, node '//trim(num2lstr(i))
  -         indx_first = indx_first + 1
  -      end do      
  -   end do
  -                  
  -   do i=1,M%NNodes
  -      do j=1,3
  -         Names(indx_first) = trim(MeshName)//Comp(j)//' rotation velocity, node '//trim(num2lstr(i))
  -         indx_first = indx_first + 1
  -      end do      
  -   end do
  +!> Routine to pack the data structures representing the operating points into arrays for linearization.
  +SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, y_op, x_op, dx_op, xd_op, z_op )
  +
  +   REAL(DbKi),                           INTENT(IN   )           :: t          !< Time in seconds at operating point
  +   TYPE(ED_InputType),                   INTENT(IN   )           :: u          !< Inputs at operating point (may change to inout if a mesh copy is required)
  +   TYPE(ED_ParameterType),               INTENT(IN   )           :: p          !< Parameters
  +   TYPE(ED_ContinuousStateType),         INTENT(IN   )           :: x          !< Continuous states at operating point
  +   TYPE(ED_DiscreteStateType),           INTENT(IN   )           :: xd         !< Discrete states at operating point
  +   TYPE(ED_ConstraintStateType),         INTENT(IN   )           :: z          !< Constraint states at operating point
  +   TYPE(ED_OtherStateType),              INTENT(IN   )           :: OtherState !< Other states at operating point
  +   TYPE(ED_OutputType),                  INTENT(IN   )           :: y          !< Output at operating point
  +   TYPE(ED_MiscVarType),                 INTENT(INOUT)           :: m          !< Misc/optimization variables
  +   INTEGER(IntKi),                       INTENT(  OUT)           :: ErrStat    !< Error status of the operation
  +   CHARACTER(*),                         INTENT(  OUT)           :: ErrMsg     !< Error message if ErrStat /= ErrID_None
  +   REAL(ReKi), ALLOCATABLE, OPTIONAL,    INTENT(INOUT)           :: u_op(:)    !< values of linearized inputs
  +   REAL(ReKi), ALLOCATABLE, OPTIONAL,    INTENT(INOUT)           :: y_op(:)    !< values of linearized outputs
  +   REAL(ReKi), ALLOCATABLE, OPTIONAL,    INTENT(INOUT)           :: x_op(:)    !< values of linearized continuous states
  +   REAL(ReKi), ALLOCATABLE, OPTIONAL,    INTENT(INOUT)           :: dx_op(:)   !< values of first time derivatives of linearized continuous states
  +   REAL(ReKi), ALLOCATABLE, OPTIONAL,    INTENT(INOUT)           :: xd_op(:)   !< values of linearized discrete states
  +   REAL(ReKi), ALLOCATABLE, OPTIONAL,    INTENT(INOUT)           :: z_op(:)    !< values of linearized constraint states
  +
  +
  +
  +   INTEGER(IntKi)                                    :: i, k, index
  +   INTEGER(IntKi)                                    :: ny
  +   INTEGER(IntKi)                                    :: ErrStat2
  +   CHARACTER(ErrMsgLen)                              :: ErrMsg2
  +   CHARACTER(*), PARAMETER                           :: RoutineName = 'ED_GetOP'
  +   TYPE(ED_ContinuousStateType)                      :: dx          !< derivative of continuous states at operating point
  +   LOGICAL                                           :: Mask(FIELDMASK_SIZE)               !< flags to determine if this field is part of the packing
  +   
  +   
  +      ! Initialize ErrStat
  +
  +   ErrStat = ErrID_None
  +   ErrMsg  = ''
  +
  +   !..................................
  +   IF ( PRESENT( u_op ) ) THEN
  +      if (.not. allocated(u_op)) then         
  +         call AllocAry(u_op, size(p%Jac_u_indx,1),'u_op',ErrStat2,ErrMsg2)
  +            call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
  +         if (ErrStat>=AbortErrLev) return
  +      end if
               
  -   !...............................
  -   ! if we want to skip the acceleration fields, lets return early:
  -   if (present(SkipAcc)) then
  -      if (SkipAcc) return
  -   end if
  -   !...............................
  -   
  -   do i=1,M%NNodes
  -      do j=1,3
  -         Names(indx_first) = trim(MeshName)//Comp(j)//' translation acceleration, node '//trim(num2lstr(i))
  -         indx_first = indx_first + 1
  -      end do      
  -   end do
  -   
  -   do i=1,M%NNodes
  -      do j=1,3
  -         Names(indx_first) = trim(MeshName)//Comp(j)//' rotation acceleration, node '//trim(num2lstr(i))
  -         indx_first = indx_first + 1
  -      end do      
  -   end do
  -
  -
  -END SUBROUTINE PackMotionMesh_Names
  -!----------------------------------------------------------------------------------------------------------------------------------
  -!> This subroutine returns the names of the output rows/columns in the Jacobian matrices.
  -SUBROUTINE PackLoadMesh_Names(M, MeshName, Names, indx_first)
  -   
  -   TYPE(MeshType)                    , INTENT(IN   ) :: M                          !< ED output mesh
  -   CHARACTER(*)                      , INTENT(IN   ) :: MeshName                   !< name of mesh 
  -   CHARACTER(LinChanLen)             , INTENT(INOUT) :: Names(:)                   !< name of column of dYdu 
  -   INTEGER(IntKi)                    , INTENT(INOUT) :: indx_first                 !< index into Names array; gives location of next array position to fill
  -   
  -      ! local variables:
  -   INTEGER(IntKi)                :: i, j
  -   character(1), parameter       :: Comp(3) = (/'X','Y','Z'/)
  -
  -   
  -   do i=1,M%NNodes
  -      do j=1,3
  -         Names(indx_first) = trim(MeshName)//Comp(j)//' force, node '//trim(num2lstr(i))
  -         indx_first = indx_first + 1
  -      end do      
  -   end do
  -                           
  -   do i=1,M%NNodes
  -      do j=1,3
  -         Names(indx_first) = trim(MeshName)//Comp(j)//' moment, node '//trim(num2lstr(i))
  -         indx_first = indx_first + 1
  -      end do      
  -   end do
  +      index = 1
  +      if (allocated(u%BladePtLoads)) then
  +         do k=1,p%NumBl
  +            call PackLoadMesh(u%BladePtLoads(k), u_op, index)   
  +         end do
  +      end if
  +      call PackLoadMesh(u%PlatformPtMesh, u_op, index)   
  +      call PackLoadMesh(u%TowerPtLoads, u_op, index)   
  +      call PackLoadMesh(u%HubPtLoad, u_op, index)   
  +      call PackLoadMesh(u%NacelleLoads, u_op, index)   
  +      
  +      do k = 1,p%NumBl ! scalars
  +         u_op(index) = u%BlPitchCom(k)
  +         index = index + 1
  +      end do
  +      u_op(index) = u%YawMom ; index = index + 1
  +      u_op(index) = u%GenTrq
  +   END IF
  +
  +   !..................................
  +   IF ( PRESENT( y_op ) ) THEN
  +      
  +      if (.not. allocated(y_op)) then 
  +            ! our operating point includes DCM (orientation) matrices, not just small angles like the perturbation matrices do
  +         ny = p%Jac_ny + y%PlatformPtMesh%NNodes * 6 & ! Jac_ny has 3 for Orientation, but we need 9 at each node
  +                       + y%TowerLn2Mesh%NNodes   * 6 & ! Jac_ny has 3 for Orientation, but we need 9 at each node 
  +                       + y%HubPtMotion%NNodes    * 6 & ! Jac_ny has 3 for Orientation, but we need 9 at each node
  +                       + y%NacelleMotion%NNodes  * 6   ! Jac_ny has 3 for Orientation, but we need 9 at each node
               
  -
  -
  -END SUBROUTINE PackLoadMesh_Names
  +         if (allocated(y%BladeLn2Mesh)) then
  +            do k=1,p%NumBl
  +               ny = ny + y%BladeLn2Mesh(k)%NNodes * 6  ! Jac_ny has 3 for Orientation, but we need 9 (at each node on each blade)
  +            end do      
  +         end if
  +         do k=1,p%NumBl
  +            ny = ny + y%BladeRootMotion(k)%NNodes * 6  ! Jac_ny has 3 for Orientation, but we need 9 at each node on each blade
  +         end do
  +                                    
  +         call AllocAry(y_op, ny,'y_op',ErrStat2,ErrMsg2)
  +            call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
  +         if (ErrStat>=AbortErrLev) return
  +      end if
  +            
  +      
  +      
  +      Mask  = .false.
  +      Mask(MASKID_TRANSLATIONDISP) = .true.
  +      Mask(MASKID_ORIENTATION) = .true.
  +      Mask(MASKID_ROTATIONVEL) = .true.
  +   
  +      index = 1
  +      if (allocated(y%BladeLn2Mesh)) then
  +         do k=1,p%NumBl
  +            call PackMotionMesh(y%BladeLn2Mesh(k), y_op, index)
  +         end do      
  +      end if
  +      call PackMotionMesh(y%PlatformPtMesh, y_op, index)
  +      call PackMotionMesh(y%TowerLn2Mesh, y_op, index)
  +      call PackMotionMesh(y%HubPtMotion, y_op, index, FieldMask=Mask)
  +      do k=1,p%NumBl
  +         call PackMotionMesh(y%BladeRootMotion(k), y_op, index)
  +      end do   
  +      call PackMotionMesh(y%NacelleMotion, y_op, index)
  +      
  +      y_op(index) = y%Yaw     ; index = index + 1    
  +      y_op(index) = y%YawRate ; index = index + 1    
  +      y_op(index) = y%HSS_Spd 
  +   
  +      do i=1,p%NumOuts
  +         y_op(i+index) = y%WriteOutput(i)
  +      end do   
  +                        
  +   END IF
  +
  +   !..................................
  +   IF ( PRESENT( x_op ) ) THEN
  +
  +      if (.not. allocated(x_op)) then                           
  +         call AllocAry(x_op, p%DOFs%NActvDOF * 2,'x_op',ErrStat2,ErrMsg2)
  +            call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
  +         if (ErrStat>=AbortErrLev) return
  +      end if
  +      
  +      do i=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs
  +         x_op(i) = x%QT( p%DOFs%PS(i) )
  +      end do
  +      do i=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs
  +         x_op(i+p%DOFs%NActvDOF) = x%QDT( p%DOFs%PS(i) )
  +      end do                                
  +      
  +   END IF
  +
  +   !..................................
  +   IF ( PRESENT( dx_op ) ) THEN
  +
  +      if (.not. allocated(dx_op)) then                           
  +         call AllocAry(dx_op, p%DOFs%NActvDOF * 2,'dx_op',ErrStat2,ErrMsg2)
  +            call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
  +         if (ErrStat>=AbortErrLev) return
  +      end if
  +      
  +      call ED_CalcContStateDeriv( t, u, p, x, xd, z, OtherState, m, dx, ErrStat2, ErrMsg2 ) 
  +         call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) 
  +         if (ErrStat>=AbortErrLev) then
  +            call ED_DestroyContState( dx, ErrStat2, ErrMsg2)
  +            return
  +         end if
  +                     
  +      do i=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs
  +         dx_op(i) = dx%QT( p%DOFs%PS(i) )
  +      end do
  +      do i=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs
  +         dx_op(i+p%DOFs%NActvDOF) = dx%QDT( p%DOFs%PS(i) )
  +      end do                                
  +      
  +      call ED_DestroyContState( dx, ErrStat2, ErrMsg2)
  +            
  +   END IF
  +
  +   !..................................
  +   IF ( PRESENT( xd_op ) ) THEN
  +   END IF
  +   
  +   !..................................
  +   IF ( PRESENT( z_op ) ) THEN
  +   END IF
  +
  +END SUBROUTINE ED_GetOP
   !----------------------------------------------------------------------------------------------------------------------------------
   
   
  
positive.code:          
     ! note: all other inputs are ignored
        
        !....................
        ! set allocate arrays needed in numerical differentiation and list of input names
        !....................
     call AllocAry(p%du,               17,    'p%du', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) ! 17 = number of unique values in p%Jac_u_indx(:,1)      
     call AllocAry(p%Jac_u_indx,       nu, 3, 'p%Jac_u_indx', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)   
     call AllocAry(InitOut%LinNames_u, nu,    'LinNames_u', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
     if (ErrStat >= AbortErrLev) return
                       
        
     
     
     ! fill matrix to store index to help us figure out what the ith value of the u vector really means
     ! (see elastodyn::ed_perturb_u ... these MUST match )
     ! column 1 indicates module's mesh and field
     ! column 2 indicates the first index of the acceleration/load field
     ! column 3 is the node
        
     !...............
     ! ED input mappings stored in p%Jac_u_indx:   
     !...............
              
     index = 1
     if (allocated(u%BladePtLoads)) then
        !Module/Mesh/Field: u%BladePtLoads(1)%Force  = 1;
  
negative.code:       
        ! set linearization output names:
     do i=1,p%DOFs%NActvDOF
        InitOut%LinNames_x(i) = p%DOF_Desc( p%DOFs%PS(i) )
     end do
     
     do i=1,p%DOFs%NActvDOF
        InitOut%LinNames_x(i+p%DOFs%NActvDOF) = 'First time derivative of '//trim(p%DOF_Desc( p%DOFs%PS(i) ))//'/s'
     end do      
     
  END SUBROUTINE ED_Init_Jacobian_x
  