anchor.diff:              end if
               
            end if
  +      end if
  +   
  +      !..........
  +      ! dU^{ED}/du^{MAP}
  +      !..........
  +      
  +      if ( p_FAST%CompMooring == Module_MAP ) then 
  +
  +         ED_Start_mt = Indx_u_ED_Platform_Start(u_ED, y_FAST) &
  +                        + u_ED%PlatformPtMesh%NNodes * 3         ! 3 forces at each node (we're going to start at the moments)
  +      
  +            ! Transfer MAP loads to ED PlatformPtmesh input:
  +            ! we're mapping loads, so we also need the sibling meshes' displacements:
  +
  +         MAP_Start = y_FAST%Lin%Modules(MODULE_MAP)%Instance(1)%LinStartIndx(LIN_INPUT_COL)
  +      
  +         ! NOTE: Assumes at least one MAP Fairlead point    
  +      
  +         CALL Linearize_Point_to_Point( MAPp%y%ptFairleadLoad, u_ED%PlatformPtMesh, MeshMapData%Mooring_P_2_ED_P, ErrStat2, ErrMsg2, MAPp%Input(1)%PtFairDisplacement, y_ED%PlatformPtMesh) !MAPp%Input(1)%ptFairleadLoad and y_ED%PlatformPtMesh contain the displaced positions for load calculations
  +            CALL SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, RoutineName)
  +               
  +            ! HD is source in the mapping, so we want M_{uSm}
  +         if (allocated(MeshMapData%Mooring_P_2_ED_P%dM%m_us )) then
  +            call SetBlockMatrix( dUdu, MeshMapData%Mooring_P_2_ED_P%dM%m_us, ED_Start_mt, MAP_Start )
  +         end if
  +   
  +      end if
  +      
  +   end if    
  +END SUBROUTINE Linear_ED_InputSolve_du
  +
  +!----------------------------------------------------------------------------------------------------------------------------------
  +!> This routine forms the dU^{SD}/du^{HD} and dU^{SD}/du^{SD} and dU^{SD}/du^{MAP} blocks (SD row) of dUdu. (i.e., how do changes in the HD and SD inputs affect the SD inputs?)
  +SUBROUTINE Linear_SD_InputSolve_du( p_FAST, y_FAST, u_SD, y_SD, y_ED, HD, MAPp, MeshMapData, dUdu, ErrStat, ErrMsg )
  +!LIN-TODO: Augment this interface for HD and MAP
  +   TYPE(FAST_ParameterType),       INTENT(IN   )  :: p_FAST         !< Glue-code simulation parameters
  +   TYPE(FAST_OutputFileType),      INTENT(IN   )  :: y_FAST         !< Glue-code output parameters (for linearization)
  +   TYPE(SD_InputType),             INTENT(INOUT)  :: u_SD           !< SD Inputs at t
  +   TYPE(SD_OutputType),            INTENT(IN   )  :: y_SD           !< SubDyn outputs (need translation displacement on meshes for loads mapping)
  +   TYPE(ED_OutputType),            INTENT(IN   )  :: y_ED           !< ElastoDyn outputs
  +   TYPE(HydroDyn_Data),            INTENT(INOUT)  :: HD             !< HD data at t
  +   TYPE(MAP_Data),                 INTENT(INOUT)  :: MAPp           !< MAP data at t
  +   TYPE(FAST_ModuleMapType),       INTENT(INOUT)  :: MeshMapData    !< Data for mapping between modules
  +   REAL(R8Ki),                     INTENT(INOUT)  :: dUdu(:,:)      !< Jacobian matrix of which we are computing the dU^(SD)/du^(AD) block
  +   INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat        !< Error status
  +   CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg         !< Error message
  +   
  +      ! local variables
  +   INTEGER(IntKi)                                 :: HD_Start
  +   INTEGER(IntKi)                                 :: MAP_Start
  +   INTEGER(IntKi)                                 :: SD_Start, SD_Start_td, SD_Start_tr
  +   INTEGER(IntKi)                                 :: ErrStat2       ! temporary Error status of the operation
  +   CHARACTER(ErrMsgLen)                           :: ErrMsg2        ! temporary Error message if ErrStat /= ErrID_None
  +   
  +   CHARACTER(*), PARAMETER                        :: RoutineName = 'Linear_SD_InputSolve_du'
  +   
  +   
  +      ! Initialize error status
  +      
  +   ErrStat = ErrID_None
  +   ErrMsg = ""
  +
  +   IF ( p_FAST%CompSub == Module_SD ) THEN ! see routine U_ED_SD_HD_BD_Orca_Residual() in SolveOption1
  +   
  +   !..........
  +   ! dU^{SD}/du^{SD}
  +   !..........
  +   
  +   call Linearize_Point_to_Point( y_ED%PlatformPtMesh, u_SD%TPMesh, MeshMapData%ED_P_2_SD_TP, ErrStat2, ErrMsg2 )
  +      call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
  +         
  +   ! SD is destination in the mapping, so we want M_{tv_uD} and M_{ta_uD}
  +            
  +   SD_Start_td = y_FAST%Lin%Modules(MODULE_SD)%Instance(1)%LinStartIndx(LIN_INPUT_COL)  
  +   SD_Start_tr = SD_Start_td + u_SD%TPMesh%NNodes * 6 ! skip 2 fields (TranslationDisp and Orientation) with 3 components before translational velocity field      
  +         
  +      ! translational velocity:
  +   if (allocated(MeshMapData%ED_P_2_SD_TP%dM%tv_uD )) then             
  +      call SetBlockMatrix( dUdu, MeshMapData%ED_P_2_SD_TP%dM%tv_ud, SD_Start_tr, SD_Start_td )
      end if
  -   
  +         
  +      ! translational acceleration:
  +   SD_Start_tr = SD_Start_tr + u_SD%TPMesh%NNodes * 6 ! skip 2 fields ( TranslationVel and RotationVel)
  +   if (allocated(MeshMapData%ED_P_2_SD_TP%dM%ta_uD )) then            
  +      call SetBlockMatrix( dUdu, MeshMapData%ED_P_2_SD_TP%dM%ta_ud, SD_Start_tr, SD_Start_td )
  +   end if
  +
  +   
  +  
      !..........
  -   ! dU^{ED}/du^{MAP}
  +   ! dU^{SD}/du^{HD}
      !..........
  -   ! LIN-TODO: Implement
  -   if ( p_FAST%CompMooring == Module_MAP ) then 
  -
  -      ED_Start_mt = Indx_u_ED_Platform_Start(u_ED, y_FAST) &
  -                     + u_ED%PlatformPtMesh%NNodes * 3         ! 3 forces at each node (we're going to start at the moments)
  -      
  -         ! Transfer MAP loads to ED PlatformPtmesh input:
  -         ! we're mapping loads, so we also need the sibling meshes' displacements:
  -
  -      MAP_Start = y_FAST%Lin%Modules(MODULE_MAP)%Instance(1)%LinStartIndx(LIN_INPUT_COL)
  -      
  -      ! NOTE: Assumes at least one MAP Fairlead point    
  -      
  -      CALL Linearize_Point_to_Point( MAPp%y%ptFairleadLoad, u_ED%PlatformPtMesh, MeshMapData%Mooring_P_2_ED_P, ErrStat2, ErrMsg2, MAPp%Input(1)%PtFairDisplacement, y_ED%PlatformPtMesh) !MAPp%Input(1)%ptFairleadLoad and y_ED%PlatformPtMesh contain the displaced positions for load calculations
  -         CALL SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, RoutineName)
  +   
  +    ! we're just going to assume u_SD%LMesh is committed
  +   SD_Start = Indx_u_SD_LMesh_Start(u_SD, y_FAST) &
  +                        + u_SD%LMesh%NNodes * 3         ! 3 forces at each node (we're going to start at the moments)
  +   
  +   if ( p_FAST%CompHydro == Module_HD ) then ! HydroDyn-{ElastoDyn or SubDyn}
  +           
  +         
  +         ! Transfer HD load outputs to SD LMesh input:
  +         if ( HD%y%Morison%Mesh%Committed  ) then ! meshes for floating
  +             
  +            ! dU^{SD}/du^{HD}  
  +            
  +               ! we're mapping loads, so we also need the sibling meshes' displacements:
  +            HD_Start = Indx_u_HD_Morison_Start(HD%Input(1), y_FAST)
  +            
  +            call Linearize_Point_to_Point( HD%y%Morison%Mesh, u_SD%LMesh, MeshMapData%HD_M_P_2_SD_P, ErrStat2, ErrMsg2, HD%Input(1)%Morison%Mesh, y_SD%Y2Mesh) !HD%Input(1)%Mesh and y_ED%PlatformPtMesh contain the displaced positions for load calculations
  +               call SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, RoutineName)
                  
  -         ! HD is source in the mapping, so we want M_{uSm}
  -      if (allocated(MeshMapData%Mooring_P_2_ED_P%dM%m_us )) then
  -         call SetBlockMatrix( dUdu, MeshMapData%Mooring_P_2_ED_P%dM%m_us, ED_Start_mt, MAP_Start )
  +               ! HD is source in the mapping, so we want M_{uSm}
  +            if (allocated(MeshMapData%HD_M_P_2_SD_P%dM%m_us )) then
  +               call SetBlockMatrix( dUdu, MeshMapData%HD_M_P_2_SD_P%dM%m_us, SD_Start, HD_Start )
  +            end if
  +            
  +        
  +            
  +            
  +         end if      
  +         if ( HD%y%WAMITMesh%Committed  ) then ! meshes for floating
  +             
  +               ! Transfer HD load outputs to ED PlatformPtMesh input:
  +               ! we're mapping loads, so we also need the sibling meshes' displacements:
  +            HD_Start = Indx_u_HD_WAMIT_Start(HD%Input(1), y_FAST)
  +            
  +            call Linearize_Point_to_Point( HD%y%WAMITMesh, u_SD%LMesh, MeshMapData%HD_W_P_2_SD_P, ErrStat2, ErrMsg2, HD%Input(1)%WAMITMesh, y_SD%Y2Mesh) !HD%Input(1)%Mesh and y_ED%PlatformPtMesh contain the displaced positions for load calculations
  +               call SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, RoutineName)
  +               
  +               ! HD is source in the mapping, so we want M_{uSm}
  +            if (allocated(MeshMapData%HD_W_P_2_SD_P%dM%m_us )) then
  +               call SetBlockMatrix( dUdu, MeshMapData%HD_W_P_2_SD_P%dM%m_us, SD_Start, HD_Start )
  +            end if
  +            
  +           
  +         end if
  +   end if
  +   
  +   !..........
  +   ! dU^{SD}/du^{MAP}
  +   !..........
  +
  +      if ( p_FAST%CompMooring == Module_MAP ) then 
  +
  +            ! Transfer MAP loads to ED PlatformPtmesh input:
  +            ! we're mapping loads, so we also need the sibling meshes' displacements:
  +   
  +         MAP_Start = y_FAST%Lin%Modules(MODULE_MAP)%Instance(1)%LinStartIndx(LIN_INPUT_COL)
  +      
  +         ! NOTE: Assumes at least one MAP Fairlead point    
  +     
  +         CALL Linearize_Point_to_Point( MAPp%y%ptFairleadLoad, u_SD%LMesh, MeshMapData%Mooring_P_2_SD_P, ErrStat2, ErrMsg2, MAPp%Input(1)%PtFairDisplacement, y_SD%Y2Mesh) !MAPp%Input(1)%ptFairleadLoad and y_SD%Y2Mesh contain the displaced positions for load calculations
  +            CALL SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, RoutineName)
  +               
  +            ! SD is source in the mapping, so we want M_{uSm}
  +         if (allocated(MeshMapData%Mooring_P_2_SD_P%dM%m_us )) then
  +            call SetBlockMatrix( dUdu, MeshMapData%Mooring_P_2_SD_P%dM%m_us, SD_Start, MAP_Start )
  +         end if
  +   
         end if
  -   
  +   END IF   
  +END SUBROUTINE Linear_SD_InputSolve_du
  +
  +!----------------------------------------------------------------------------------------------------------------------------------
  +!> This routine forms the dU^{SD}/du^{HD} and dU^{SD}/du^{SD} blocks (SD row) of dUdu. (i.e., how do changes in the HD and SD inputs affect the SD inputs?)
  +SUBROUTINE Linear_SD_InputSolve_dy( p_FAST, y_FAST, u_SD, y_SD, y_ED, HD, MAPp, MeshMapData, dUdy, ErrStat, ErrMsg )
  +!LIN-TODO: Augment this interface for HD and MAP
  +   TYPE(FAST_ParameterType),       INTENT(IN   )  :: p_FAST         !< Glue-code simulation parameters
  +   TYPE(FAST_OutputFileType),      INTENT(IN   )  :: y_FAST         !< Glue-code output parameters (for linearization)
  +   TYPE(SD_InputType),             INTENT(INOUT)  :: u_SD           !< SD Inputs at t
  +   TYPE(SD_OutputType),            INTENT(IN   )  :: y_SD           !< SubDyn outputs (need translation displacement on meshes for loads mapping)
  +   TYPE(ED_OutputType),            INTENT(IN   )  :: y_ED           !< ElastoDyn outputs
  +   TYPE(HydroDyn_Data),            INTENT(INOUT)  :: HD             !< HD data at t
  +   TYPE(MAP_Data),                 INTENT(INOUT)  :: MAPp           !< MAP data at t
  +   TYPE(FAST_ModuleMapType),       INTENT(INOUT)  :: MeshMapData    !< Data for mapping between modules
  +   REAL(R8Ki),                     INTENT(INOUT)  :: dUdy(:,:)      !< Jacobian matrix of which we are computing the dU^(SD)/dy^(SD) block
  +   INTEGER(IntKi),                 INTENT(  OUT)  :: ErrStat        !< Error status
  +   CHARACTER(*),                   INTENT(  OUT)  :: ErrMsg         !< Error message
  +   
  +      ! local variables
  +   INTEGER(IntKi)                                 :: SD_Start, SD_Out_Start, HD_Start, HD_Out_Start, ED_Out_Start, MAP_Out_Start
  +   INTEGER(IntKi)                                 :: MAP_Start
  +   INTEGER(IntKi)                                 :: ErrStat2       ! temporary Error status of the operation
  +   CHARACTER(ErrMsgLen)                           :: ErrMsg2        ! temporary Error message if ErrStat /= ErrID_None
  +   
  +   CHARACTER(*), PARAMETER                        :: RoutineName = 'Linear_SD_InputSolve_du'
  +   
  +   
  +      ! Initialize error status
  +      
  +   ErrStat = ErrID_None
  +   ErrMsg = ""
  +   if ( p_FAST%CompSub /= Module_SD ) return  !LIN-TODO: Verify that it is ok to simply return.  What is dUdy in that case
  +   
  +   !..........
  +   ! dU^{SD}/dy^{ED}
  +   !..........
  +  
  +   !!! ! This linearization was done in forming dUdu (see Linear_SD_InputSolve_du()), so we don't need to re-calculate these matrices 
  +   !!! ! while forming dUdy, too.
  +   !!!call Linearize_Point_to_Line2( y_ED%PlatformPtMesh, u_SD%TPMesh, MeshMapData%ED_P_2_SD_TP, ErrStat2, ErrMsg2 )
  +      
  +   SD_Start     = Indx_u_SD_TPMesh_Start(u_SD, y_FAST)  ! start of u_SD%MTPMesh%TranslationDisp field     
  +   ED_Out_Start = Indx_y_ED_Platform_Start(y_ED, y_FAST) ! start of y_ED%PlatformPtMesh%TranslationDisp field
  +   call Assemble_dUdy_Motions(y_ED%PlatformPtMesh, u_SD%TPMesh, MeshMapData%ED_P_2_SD_TP, SD_Start, ED_Out_Start, dUdy, .false.)
  +   
  +   !..........
  +   ! dU^{SD}/dy^{HD}
  +   !..........
  +   ! HD
  +   ! parts of dU^{SD}/dy^{HD} and dU^{SD}/dy^{SD}:
  +   if ( p_FAST%CompHydro == Module_HD ) then ! HydroDyn-SubDyn
  +      SD_Out_Start = Indx_y_SD_Y2Mesh_Start(y_SD, y_FAST) ! start of y_SD%Y2Mesh%TranslationDisp field
  +         ! we're just going to assume u_SD%LMesh is committed
  +      if ( HD%y%Morison%Mesh%Committed  ) then ! meshes for floating
  +         !!! ! This linearization was done in forming dUdu (see Linear_ED_InputSolve_du()), so we don't need to re-calculate these matrices 
  +         !!! ! while forming dUdy, too.
  +         ! call Linearize_Point_to_Point( HD%y%Morison, u_ED%PlatformPtMesh, MeshMapData%HD_M_P_2_ED_P, ErrStat2, ErrMsg2, HD%Input(1)%Morison, y_ED%PlatformPtMesh) !HD%Input(1)%Morison and y_ED%PlatformPtMesh contain the displaced positions for load calculations
  +         HD_Out_Start = Indx_y_HD_Morison_Start(HD%y, y_FAST)
  +         SD_Start     = Indx_u_SD_LMesh_Start(u_SD, y_FAST) ! start of u_SD%LMesh%Force field
  +         call Assemble_dUdy_Loads(HD%y%Morison%Mesh, u_SD%LMesh, MeshMapData%HD_M_P_2_SD_P, SD_Start, HD_Out_Start, dUdy)
  +         
  +            ! SD translation displacement-to-SD moment transfer (dU^{SD}/dy^{SD}):
  +         SD_Start = Indx_u_SD_LMesh_Start(u_SD, y_FAST) + u_SD%LMesh%NNodes*3   ! start of u_SD%LMesh%Moment field (skip the SD forces) 
  +         call SetBlockMatrix( dUdy, MeshMapData%HD_M_P_2_SD_P%dM%m_uD, SD_Start, SD_Out_Start )
  +! maybe this should be SumBlockMatrix with future changes to linearized modules???            
  +      end if      
  +      if ( HD%y%WAMITMesh%Committed  ) then ! meshes for floating
  +         !!! ! This linearization was done in forming dUdu (see Linear_ED_InputSolve_du()), so we don't need to re-calculate these matrices 
  +         !!! ! while forming dUdy, too.
  +         ! call Linearize_Point_to_Point( HD%y%WAMITMesh, u_ED%PlatformPtMesh, MeshMapData%HD_W_P_2_ED_P, ErrStat2, ErrMsg2, HD%Input(1)%WAMITMesh, y_ED%PlatformPtMesh) !HD%Input(1)%WAMITMesh and y_ED%PlatformPtMesh contain the displaced positions for load calculations
  +         HD_Out_Start = Indx_y_HD_WAMIT_Start(HD%y, y_FAST)
  +         SD_Start     = Indx_u_SD_LMesh_Start(u_SD, y_FAST) ! start of u_SD%LMesh%Force field
  +         call Assemble_dUdy_Loads(HD%y%WAMITMesh, u_SD%LMesh, MeshMapData%HD_W_P_2_SD_P, SD_Start, HD_Out_Start, dUdy)
  +         
  +            ! SD translation displacement-to-SD moment transfer (dU^{SD}/dy^{SD}):
  +         SD_Start = Indx_u_SD_LMesh_Start(u_SD, y_FAST) + u_SD%LMesh%NNodes*3   ! start of u_SD%LMesh%Moment field (skip the SD forces)  
  +         call SumBlockMatrix( dUdy, MeshMapData%HD_W_P_2_SD_P%dM%m_uD, SD_Start, SD_Out_Start )
  +! maybe this should be SumBlockMatrix with future changes to linearized modules???            
  +      end if
      end if
  -      
  -END SUBROUTINE Linear_ED_InputSolve_du
  -
  +   
  +   !..........
  +   ! dU^{SD}/dy^{MAP}
  +   !..........
  +   if ( p_FAST%CompMooring == Module_MAP ) then
  +      if ( MAPp%y%ptFairleadLoad%Committed  ) then ! meshes for floating
  +         !!! ! This linearization was done in forming dUdu (see Linear_SD_InputSolve_du()), so we don't need to re-calculate these matrices 
  +         !!! ! while forming dUdy, too.
  +         !       CALL Linearize_Point_to_Point( MAPp%y%ptFairleadLoad, u_SD%LMesh, MeshMapData%Mooring_P_2_SD_P, ErrStat2, ErrMsg2, MAPp%Input(1)%PtFairDisplacement, y_SD%Y2Mesh) !MAPp%Input(1)%ptFairleadLoad and y_ED%Y2Mesh contain the displaced positions for load calculations
  +         MAP_Out_Start = y_FAST%Lin%Modules(MODULE_MAP)%Instance(1)%LinStartIndx(LIN_OUTPUT_COL)
  +         SD_Start      = Indx_u_SD_LMesh_Start(u_SD, y_FAST) ! start of u_SD%LMesh%TranslationDisp field
  +         call Assemble_dUdy_Loads(MAPp%y%ptFairLeadLoad, u_SD%LMesh, MeshMapData%Mooring_P_2_SD_P, SD_Start, MAP_Out_Start, dUdy)
  +      
  +         ! SD translation displacement-to-SD moment transfer (dU^{SD}/dy^{SD}):
  +         SD_Start = Indx_u_SD_LMesh_Start(u_SD, y_FAST) + u_SD%LMesh%NNodes*3   ! start of u_ED%LMesh%Moment field (skip the SD forces)
  +         SD_Out_Start = Indx_y_SD_Y2Mesh_Start(y_SD, y_FAST) ! start of y_SD%Y2Mesh%TranslationDisp field
  +         call SumBlockMatrix( dUdy, MeshMapData%Mooring_P_2_SD_P%dM%m_uD, SD_Start, SD_Out_Start )
  +      end if     
  +   end if
  +END SUBROUTINE Linear_SD_InputSolve_dy  
  +   
   
   !----------------------------------------------------------------------------------------------------------------------------------
   !> This routine forms the dU^{BD}/du^{BD} and dU^{BD}/du^{AD} blocks (BD row) of dUdu. (i.e., how do changes in the AD and BD inputs 
  
positive.code:          call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
     end if 
  
     IF (p_FAST%CompSub == Module_SD) THEN
         write(*,*)'>>> FAST_LIN: Linear_SD_InputSolve_dy, TODO'
     ELSE IF (p_FAST%CompSub == Module_ExtPtfm) THEN
         CALL WrScr('>>> FAST_LIN: Linear_ExtPtfm_InputSolve_du, TODO')
     ENDIF
  
        !............
        ! \f$ \frac{\partial U_\Lambda^{HD}}{\partial u^{HD}} \end{bmatrix} = \f$ (dUdu block row 6=HD)
        !............
     IF (p_FAST%CompHydro == MODULE_HD) THEN 
        call Linear_HD_InputSolve_du( p_FAST, y_FAST, HD%Input(1), ED%y, MeshMapData, dUdu, ErrStat2, ErrMsg2 )
           call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
     end if 
  
     ! LIN-TODO: Update the doc lines below to include HD and MAP
     !.....................................
     ! dUdy
     !> \f$ \frac{\partial U_\Lambda}{\partial y} =  
  
negative.code:       end if
     END IF
              
     !LIN-TODO: Do we need this for SD and Map?  GJH          
     !!!! SubDyn: copy final predictions to actual states
     !!!IF ( p_FAST%CompSub == Module_SD ) THEN
     !!!ELSE IF ( p_FAST%CompSub == Module_ExtPtfm ) THEN
  