diff --git a/docs/source/dev/glue-code/ModVar.svg b/docs/source/dev/glue-code/ModVar.svg
deleted file mode 100644
index 8271b93f32..0000000000
--- a/docs/source/dev/glue-code/ModVar.svg
+++ /dev/null
@@ -1,558 +0,0 @@
-
-
\ No newline at end of file
diff --git a/docs/source/user/general.rst b/docs/source/user/general.rst
index e68f7e2cba..347e9c097d 100644
--- a/docs/source/user/general.rst
+++ b/docs/source/user/general.rst
@@ -47,6 +47,9 @@ of publications, presentations, and past studies on the subject.
Glue Code and Mesh Mapping
**************************
+For current documentation on the glue code structure, module variable API,
+solver, and linearization see :ref:`glue-code`.
+
- `FAST Modular Wind Turbine CAE Tool: Nonmatching Spatial and Temporal Meshes `_
- `FAST Modular Framework for Wind Turbine Simulation: New Algorithms and Numerical Examples `_
- :download:`Predictor-Corrector Approach <../../OtherSupporting/ProposedPCApproach_Rev4.docx>`
diff --git a/docs/source/user/glue-code/ModVar.svg b/docs/source/user/glue-code/ModVar.svg
new file mode 100644
index 0000000000..46620dc081
--- /dev/null
+++ b/docs/source/user/glue-code/ModVar.svg
@@ -0,0 +1,572 @@
+
+
+
+
+
diff --git a/docs/source/user/glue-code/index.rst b/docs/source/user/glue-code/index.rst
new file mode 100644
index 0000000000..cebb573088
--- /dev/null
+++ b/docs/source/user/glue-code/index.rst
@@ -0,0 +1,18 @@
+.. _glue-code:
+
+Glue Code
+=========
+
+The OpenFAST *glue code* is the layer of software that initializes each physics
+module, manages the flow of data between them, orchestrates the time-stepping
+loop, and—optionally—linearizes the assembled system. This section documents
+the glue code from a user and module-developer perspective.
+
+.. toctree::
+ :maxdepth: 2
+
+ overview
+ modvar
+ modglue
+ solver
+ linearization
diff --git a/docs/source/user/glue-code/linearization.rst b/docs/source/user/glue-code/linearization.rst
new file mode 100644
index 0000000000..a2b3e413cc
--- /dev/null
+++ b/docs/source/user/glue-code/linearization.rst
@@ -0,0 +1,269 @@
+.. _glue-code-linearization:
+
+Linearization
+=============
+
+OpenFAST can linearise the full multi-physics system about a periodic (or
+static) operating point to produce continuous-time, first-order state-space
+matrices of the form
+
+.. math::
+
+ \dot{\mathbf{x}} &= A\,\mathbf{x} + B\,\mathbf{u} \\
+ \mathbf{y} &= C\,\mathbf{x} + D\,\mathbf{u}
+
+together with the coupling matrices *dUdu* (input-to-input feed-through) and
+*dUdy* (output-to-input coupling). The linearization engine lives in
+``modules/openfast-library/src/FAST_ModGlue.f90``.
+
+.. contents::
+ :local:
+ :depth: 2
+
+User inputs for linearization
+------------------------------
+
+The following parameters appear in the main OpenFAST input file (``*.fst``)
+under the **Linearization** section.
+
+.. list-table::
+ :header-rows: 1
+ :widths: 25 12 63
+
+ * - Parameter
+ - Type
+ - Description
+ * - ``Linearize``
+ - logical
+ - Master switch. Set to ``True`` to enable all linearization
+ functionality. When ``False`` all other linearization parameters are
+ ignored.
+ * - ``CalcSteady``
+ - logical
+ - When ``True``, OpenFAST first runs the simulation forward until the
+ outputs at each target azimuth converge from one rotor revolution to the
+ next (steady-state trimming), then performs linearization at each
+ azimuth. When ``False``, linearization is performed at user-specified
+ absolute simulation times (``LinTimes``).
+ * - ``TrimCase``
+ - integer
+ - Controller degree of freedom trimmed during ``CalcSteady`` to achieve
+ periodic steady state.
+
+ * ``1`` – yaw
+ * ``2`` – generator torque
+ * ``3`` – collective blade pitch
+ * - ``TrimTol``
+ - real
+ - RMS convergence tolerance on normalised output error across one
+ rotor revolution. Trimming stops when the error falls below this
+ value. Typical value: ``1.0e-5``.
+ * - ``TrimGain``
+ - real
+ - Proportional gain used by the built-in trim controller.
+ Units are rad/(rad/s) for yaw/pitch cases and N·m/(rad/s) for
+ the torque case.
+ * - ``Twr_Kdmp``
+ - real
+ - Artificial tower damping coefficient (N/(m/s)) added during the
+ ``CalcSteady`` run to help damp transients and reach steady state
+ faster. Set to 0 to disable.
+ * - ``Bld_Kdmp``
+ - real
+ - Artificial blade damping coefficient (N/(m/s)) during ``CalcSteady``.
+ * - ``NLinTimes``
+ - integer
+ - Number of linearization time points per rotor revolution (or number of
+ equally spaced absolute time instants when ``CalcSteady=False``).
+ Must be ≥ 1. For a periodic model at least 12 azimuths are typically
+ needed to resolve the per-revolution variation.
+ * - ``LinTimes``
+ - real array
+ - Absolute simulation times (seconds) at which to linearise when
+ ``CalcSteady=False``. Length must equal ``NLinTimes``. Ignored when
+ ``CalcSteady=True``.
+ * - ``LinInputs``
+ - integer
+ - Controls which input variables appear in the **B** and **D** matrices.
+
+ * ``0`` (``LIN_NONE``) – no inputs; produces state matrix only.
+ * ``1`` (``LIN_STANDARD``) – inputs flagged ``VF_Linearize`` by the
+ module (default set by each module's ``InitVars``).
+ * ``2`` (``LIN_ALL``) – all module inputs including debug ones.
+ * - ``LinOutputs``
+ - integer
+ - Controls which output variables appear in the **C** and **D** matrices.
+
+ * ``0`` (``LIN_NONE``) – no outputs.
+ * ``1`` (``LIN_STANDARD``) – ``WriteOutput`` channels only
+ (``VF_WriteOut`` flag).
+ * ``2`` (``LIN_ALL``) – all module outputs.
+ * - ``LinOutJac``
+ - logical
+ - When ``True`` (requires ``LinInputs=LinOutputs=2``), the full module
+ Jacobian matrices are written to the linearization output file for
+ debugging.
+ * - ``LinOutMod``
+ - logical
+ - When ``True``, per-module ``.lin`` files are written in addition to the
+ full-system file.
+
+Module support for linearization
+----------------------------------
+
+Modules that appear in the linearization variable ordering (set in
+``ModGlue_Init``) are:
+
+InflowWind → SeaState → ServoDyn → ElastoDyn → BeamDyn → AeroDyn →
+HydroDyn → SubDyn → MAP++ → MoorDyn
+
+A module that is not in this ordered list causes a fatal error if
+``Linearize=True``.
+
+Variable selection
+------------------
+
+During ``ModGlue_Init``, the ``VF_Linearize`` flag is applied to variables
+according to the ``LinInputs`` and ``LinOutputs`` settings:
+
+* **States (x)**: the ``VF_Linearize`` flag is always set on all continuous
+ state variables of every participating module.
+* **Inputs (u)**:
+
+ * ``LIN_NONE`` → flag cleared on all input variables.
+ * ``LIN_STANDARD`` → keeps whatever ``VF_Linearize`` flag was set in the
+ module's ``InitVars``; module developers choose the *standard* input set.
+ * ``LIN_ALL`` → flag set on all input variables.
+ * Variables with ``VF_NoLin`` always have ``VF_Linearize`` cleared,
+ regardless of the above setting.
+
+* **Outputs (y)**:
+
+ * ``LIN_NONE`` → flag cleared on all output variables.
+ * ``LIN_STANDARD`` → flag set only on outputs that also carry ``VF_WriteOut``.
+ * ``LIN_ALL`` → flag set on all output variables.
+ * Variables with ``VF_NoLin`` are always excluded.
+
+The combined variable set is assembled into a ``ModGlueType`` named ``Lin``
+via ``ModGlue_CombineModules``.
+
+Steady-state trimming (``CalcSteady``)
+---------------------------------------
+
+When ``CalcSteady=True``, ``ModGlue_CalcSteady`` is called at each time step
+to detect periodicity:
+
+1. The module outputs tagged ``VF_Linearize`` (excluding ``VF_WriteOut``) are
+ collected into a buffer indexed by azimuth angle.
+2. After each complete revolution the outputs at each of the ``NLinTimes``
+ azimuth targets are compared against the previous revolution via the
+ normalised RMS error:
+
+ .. math::
+
+ \varepsilon = \sqrt{\frac{1}{N} \sum_{i=1}^{N}
+ \left(\frac{y_i^{\rm current} - y_i^{\rm previous}}{r_i}\right)^2}
+
+ where :math:`r_i = \max(y_{i,\rm max} - y_{i,\rm min},\, 0.01)` is the
+ output range from the current revolution (with a floor to avoid division
+ by near-zero).
+
+3. When :math:`\varepsilon < \texttt{TrimTol}`, ``FoundSteady=True`` and
+ linearization at all ``NLinTimes`` azimuths proceeds automatically.
+
+4. If the simulation reaches within approximately two revolutions of ``TMax``
+ without converging, a warning is issued and linearization is forced.
+
+The azimuth interpolation between buffer samples uses the extrapolation
+routines from ``MV_ExtrapInterp`` (supports constant, linear, and quadratic
+schemes depending on the number of available samples).
+
+linearization at an operating point
+-------------------------------------
+
+``ModGlue_Linearize_OP`` assembles the full-system matrices at a single
+operating point (time / azimuth):
+
+1. **Module Jacobians**: for each module,
+ ``FAST_JacobianPInput`` and ``FAST_JacobianPContState`` are called to
+ compute the per-module sub-matrices *dYdu*, *dXdu*, *dYdx*, *dXdx* by
+ central-difference finite differentiation. The perturbation magnitudes are
+ taken from each variable's ``Perturb`` field (see :ref:`glue-code-modvar`).
+
+2. **Operating point extraction**: ``FAST_GetOP`` packs the current states,
+ inputs, and outputs into the linearization arrays
+ (``ModGlue%Lin%x``, ``%u``, ``%y``).
+
+3. **Coupling matrices**: the input-output coupling matrices *dUdu* and *dUdy*
+ are assembled from the mesh-mapping Jacobians to account for the fact that
+ some module inputs are functions of other modules' outputs.
+
+4. **Full-system assembly**: the per-module sub-matrices are placed into the
+ combined glue-level matrices using the ``iGlu`` index ranges stored in
+ each ``ModVarType``.
+
+5. **Output**: ``ModGlue_CalcWriteLinearMatrices`` writes the ``.lin`` file
+ containing:
+
+ * Operating point values (**x_op**, **u_op**, **y_op**)
+ * linearization channel names (from ``LinNames``)
+ * Derivative order indicators (``VF_DerivOrder1``, ``VF_DerivOrder2``)
+ * Rotating-frame flags (``VF_RotFrame``)
+ * Full-system matrices **A**, **B**, **C**, **D**, **dUdu**, **dUdy**
+ * Per-module matrices (if ``LinOutMod=True``)
+ * Full Jacobians (if ``LinOutJac=True``)
+
+Output file format
+-------------------
+
+Each linearization call produces a file named
+``..lin`` where *N* is the linearization index (1 … ``NLinTimes``).
+The file is a plain-text ASCII file that can be read by the
+`openfast_io `_ Python library or the
+`pyFAST `_ post-processing tools.
+
+Key fields in the file header:
+
+* ``Rotor_Speed`` – rotor speed at linearization time (RPM)
+* ``Azimuth`` – blade-1 azimuth at linearization time (deg)
+
+Variable naming conventions
+----------------------------
+
+In linearization output files each channel label follows the pattern:
+
+`` [, component [, node [, unit]]]``
+
+Examples:
+
+* ``ED BlPitch1, rad`` – ElastoDyn individual blade-1 pitch state
+* ``AD B1N001Fx force, node 1, N`` – AeroDyn blade 1 node 1 X-force input
+* ``BD_1 B1TipTDxr translation displacement, node 10, m`` – BeamDyn instance 1
+
+Module developers should ensure that the ``Name`` argument to ``MV_AddVar`` /
+``MV_AddMeshVar`` and the entries in ``LinNames`` follow this convention for
+consistency with post-processing tools.
+
+Module developer responsibilities
+-----------------------------------
+
+To participate in linearization a module must:
+
+1. Call ``MV_AddVar`` / ``MV_AddMeshVar`` with appropriate ``VF_Linearize``
+ flags and supply ``LinNames`` for all variables that may appear in the
+ standard linearization set.
+
+2. Implement ``_JacobianPInput`` and ``_JacobianPContState``
+ subroutines (or supply analytical Jacobians through the registry). The
+ glue code calls these via the ``FAST_JacobianPInput`` /
+ ``FAST_JacobianPContState`` wrappers in ``FAST_Funcs.f90``.
+
+3. Implement ``_GetOP`` (via the registry) to extract the operating-point
+ values of states, inputs, and outputs.
+
+4. Mark variables that should **not** participate in linearization with
+ ``VF_NoLin``.
+
+5. Mark variables in the rotating reference frame with ``VF_RotFrame`` so that
+ multi-blade coordinate (MBC) transformations applied by post-processing
+ tools are aware of these variables.
diff --git a/docs/source/user/glue-code/modglue.rst b/docs/source/user/glue-code/modglue.rst
new file mode 100644
index 0000000000..c07a39d5ec
--- /dev/null
+++ b/docs/source/user/glue-code/modglue.rst
@@ -0,0 +1,265 @@
+.. _glue-code-modglue:
+
+ModGlue – Combining Modules into Global Arrays
+===============================================
+
+``ModGlue_CombineModules`` is the single routine that transforms the
+per-module variable descriptions registered by ``MV_AddVar`` /
+``MV_AddMeshVar`` / ``MV_AddModule`` (see :ref:`glue-code-modvar`) into a
+*monolithic* view of the coupled system. The output of the routine is a
+``ModGlueType`` structure whose ``Vars`` member holds globally indexed
+arrays for every state, input, and output variable that belongs to a
+specified set of modules. These global arrays underpin both the solver
+Jacobian construction (:ref:`glue-code-solver`) and the linearization
+procedure (:ref:`glue-code-linearization`).
+
+.. _fig-modvar-types:
+
+.. figure:: ModVar.svg
+ :alt: ModVar data structure hierarchy
+ :align: center
+ :width: 100%
+
+ Complete data-structure hierarchy used by the glue code. Types shown on
+ the left (``DatLoc``, ``Field ID``, ``Variable Flag``) are referenced from
+ each ``ModVarType``. ``ModDataType`` (centre, per module) is collected into
+ the top-level ``ModGlueType`` (right) by ``ModGlue_CombineModules``.
+
+.. contents::
+ :local:
+ :depth: 2
+
+The ``ModGlueType`` structure
+------------------------------
+
+``ModGlueType`` is defined in ``modules/openfast-library/src/Glue_Types.f90``:
+
+.. code-block:: fortran
+
+ type ModGlueType
+ character(ChanLen) :: Name ! Label (e.g. 'Solver', 'Lin')
+ type(ModDataType), allocatable :: ModData(:) ! Per-module view
+ type(ModVarsType) :: Vars ! Combined variable arrays
+ type(ModLinType) :: Lin ! Linearization matrices
+ type(VarMapType), allocatable :: VarMaps(:) ! Relevant mappings
+ end type ModGlueType
+
+Key members:
+
+.. list-table::
+ :header-rows: 1
+ :widths: 20 80
+
+ * - Member
+ - Description
+ * - ``ModData(:)``
+ - One entry per module included in this glue instance. Each entry is a
+ filtered copy of the original ``ModDataType``, containing only the
+ variable subset selected by ``FlagFilter``. Index ordering matches the
+ ``iModAry`` argument passed to ``ModGlue_CombineModules``.
+ * - ``Vars%x(:)``
+ - Concatenated array of ``ModVarType`` descriptors for all *continuous-state*
+ variables across the selected modules.
+ * - ``Vars%u(:)``
+ - Concatenated array of ``ModVarType`` descriptors for all *input* variables.
+ * - ``Vars%y(:)``
+ - Concatenated array of ``ModVarType`` descriptors for all *output* variables.
+ * - ``Vars%Nx / %Nu / %Ny``
+ - Total number of scalar values in each group (sum of ``Var%Num`` across
+ all variables in ``Vars%x / %u / %y`` respectively). These are the
+ row/column dimensions of the global data arrays and Jacobian matrices.
+ * - ``Lin``
+ - Holds linearization operating-point arrays (``x``, ``dx``, ``u``, ``y``)
+ and the full-system matrices (``dXdx``, ``dXdu``, ``dYdx``, ``dYdu``,
+ ``dUdu``, ``dUdy``). Only allocated when ``Linearize=.true.``.
+ * - ``VarMaps(:)``
+ - Filtered subset of the global ``Mappings`` array containing only the
+ mappings whose source **and** destination modules both appear in
+ ``iModAry``. Used during Jacobian finite-differencing to account for
+ output-to-input coupling.
+
+The ``iLoc`` index range
+-------------------------
+
+After ``ModGlue_CombineModules`` returns, each variable in
+``ModGlue%Vars%x / %u / %y`` carries an ``iLoc(1:2)`` range that locates it
+inside the *glue-level* data vectors. Specifically, for a variable at position
+*k* in ``ModGlue%Vars%x``:
+
+* ``iLoc(1)`` – index of its first scalar value in a length-``Vars%Nx`` array.
+* ``iLoc(2)`` – index of its last scalar value.
+
+This is the *glue-local* index after filtration; the corresponding
+per-module position is still available through the variable's ``iGlu`` range
+(set earlier by ``FAST_SolverInit → CalcVarGlobalIndices``).
+The ``iLoc / iGlu`` separation means the same variable descriptor can live
+simultaneously in the per-module ``ModData`` view, the solver ``m%Mod``
+view, and the linearization ``m%ModGlue`` view with consistent, non-overlapping
+index ranges in each context.
+
+What ``ModGlue_CombineModules`` does
+--------------------------------------
+
+The subroutine signature is:
+
+.. code-block:: fortran
+
+ subroutine ModGlue_CombineModules(ModGlue, ModDataAry, Mappings, &
+ iModAry, FlagFilter, Linearize, &
+ ErrStat, ErrMsg, Name)
+
+.. list-table::
+ :header-rows: 1
+ :widths: 22 12 66
+
+ * - Argument
+ - Intent
+ - Description
+ * - ``ModGlue``
+ - ``out``
+ - The ``ModGlueType`` structure to populate.
+ * - ``ModDataAry(:)``
+ - ``in``
+ - Full array of per-module data registered by the glue code.
+ * - ``Mappings(:)``
+ - ``in``
+ - Full array of mesh and variable mappings.
+ * - ``iModAry(:)``
+ - ``in``
+ - Ordered list of indices into ``ModDataAry`` specifying *which* modules
+ to include and in *what order*. The order determines the row/column
+ layout of the global vectors and Jacobians.
+ * - ``FlagFilter``
+ - ``in``
+ - Bitmask of ``VF_*`` flags. Only variables that have **at least one** of
+ these flags set (i.e. ``MV_HasFlagsAny(Var, FlagFilter)`` is true) are
+ copied into the glue ``Vars`` arrays.
+ * - ``Linearize``
+ - ``in``
+ - When ``.true.``, allocates the ``Lin`` operating-point arrays and the
+ full-system Jacobian matrices.
+ * - ``Name``
+ - ``in`` (optional)
+ - Human-readable label stored in ``ModGlue%Name`` (e.g. ``'Solver'``,
+ ``'Lin'``).
+
+The four main steps performed internally are:
+
+1. **Count and allocate**. Iterate over each module in ``iModAry`` and count
+ how many variable descriptors (and how many total scalar values) pass the
+ ``FlagFilter`` test for each of the *x*, *u*, *y* groups. Allocate
+ ``ModGlue%Vars%x / %u / %y`` to exactly those sizes.
+
+2. **Copy and re-index**. For each module, copy the filtered
+ ``ModVarType`` descriptors into ``ModGlue%Vars`` and assign contiguous
+ ``iLoc`` ranges so the glue-level index is consecutive across all
+ modules. The per-module ``ModGlue%ModData(i)%Vars`` sub-arrays receive
+ the same descriptors with the same ``iLoc`` so that
+ scatter/gather routines can operate directly on global vectors.
+ Linear name prefixes (e.g. ``"ED "``, ``"BD_1 "``) are prepended to
+ ``LinNames`` at this step to produce globally unique channel labels.
+
+3. **Filter mappings**. Iterate over the full ``Mappings`` array and retain
+ only those where both ``iModSrc`` and ``iModDst`` appear in ``iModAry``.
+ Re-index the retained entries against the local ``ModData`` position (not
+ the global ``ModDataAry`` position) and store them in
+ ``ModGlue%VarMaps``.
+
+4. **Allocate linearization storage** (only when ``Linearize=.true.``).
+ Allocate the ``Lin`` operating-point vectors (``x``, ``dx``, ``u``, ``y``)
+ and all six Jacobian matrices (``dXdx``, ``dXdu``, ``dYdx``, ``dYdu``,
+ ``dUdu``, ``dUdy``) dimensioned by ``Vars%Nx`` and ``Vars%Nu / Ny``.
+
+Where ``ModGlue_CombineModules`` is called
+------------------------------------------
+
+The routine is called twice during OpenFAST initialisation, producing two
+distinct ``ModGlueType`` instances with different variable selections:
+
+``m%Mod`` – the Solver glue module (in ``FAST_SolverInit``)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. code-block:: fortran
+
+ iMod = [p%iModTC, p%iModOpt1] ! TC + Option-1 indices
+
+ call ModGlue_CombineModules(m%Mod, GlueModData, GlueModMaps, iMod, &
+ VF_Solve, .true., ErrStat2, ErrMsg2, &
+ Name='Solver')
+
+* **Modules included**: tight-coupling (TC) modules plus Option-1 modules.
+* **Variable filter**: ``VF_Solve``. This flag is set by
+ ``FAST_SolverInit → SetVarSolveFlags`` on every variable that must appear
+ in the Newton iteration: TC continuous states, motion/load mesh
+ inputs/outputs involved in inter-module couplings, and any
+ variable-to-variable mapped inputs/outputs (see :ref:`glue-code-solver`
+ for the full selection criteria).
+* **Linearize**: ``.true.`` — the Jacobian matrices are allocated here
+ because the solver's Newton linear system reuses the same storage as the
+ operating-point linearization.
+* **Result**: ``m%Mod%Vars%Nx`` equals the number of TC displacement/velocity
+ scalars (``p%NumQ``); ``m%Mod%Vars%Nu`` covers all TC and Option-1 inputs
+ flagged ``VF_Solve``. The Jacobian dimension ``p%NumJ = p%NumQ + p%NumU``
+ follows directly.
+
+``m_FAST%ModGlue`` – the Linearization glue module (in ``ModGlue_Init``)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. code-block:: fortran
+
+ LinFlags = VF_Linearize + VF_Mapping + VF_Mesh
+
+ call ModGlue_CombineModules(m%ModGlue, m%ModData, m%Mappings, &
+ p%Lin%iMod, LinFlags, &
+ p_FAST%Linearize, ErrStat2, ErrMsg2, &
+ Name="Lin")
+
+* **Modules included**: all modules participating in linearization, in the
+ canonical order set by ``p%Lin%iMod`` (InflowWind → SeaState → ServoDyn →
+ ElastoDyn → BeamDyn → AeroDyn → HydroDyn → SubDyn → MAP++ → MoorDyn).
+* **Variable filter**: ``VF_Linearize + VF_Mapping + VF_Mesh``. The
+ ``VF_Linearize`` flag is applied per variable according to the
+ ``LinInputs`` / ``LinOutputs`` settings in the input file;
+ ``VF_Mapping`` / ``VF_Mesh`` ensure that mesh-coupled variables are always
+ included so that coupling Jacobians can be assembled even when a variable
+ is not a formal linearization output.
+* **Linearize**: ``p_FAST%Linearize`` — only allocates the full full-system
+ matrices when linearization is requested.
+* **Result**: ``m%ModGlue%Vars%Nx / %Nu / %Ny`` give the dimensions of the
+ **A**, **B**, **C**, **D**, **dUdu**, **dUdy** matrices written to the
+ ``.lin`` file.
+
+How the global index enables matrix assembly
+--------------------------------------------
+
+Because every variable descriptor carries its ``iLoc`` range in the glue
+``Vars`` array, scatter and gather operations on global data vectors become
+trivial index-range assignments. For example, to pack a module's state
+vector into the global solver vector ``x_global``:
+
+.. code-block:: fortran
+
+ do iVar = 1, size(ModData%Vars%x)
+ associate (Var => ModData%Vars%x(iVar))
+ x_global(Var%iLoc(1):Var%iLoc(2)) = x_mod(Var%iGlu(1):Var%iGlu(2))
+ end associate
+ end do
+
+And the corresponding gather from the global Jacobian column into a
+per-module sub-column requires no offset arithmetic: the ``iLoc`` index
+already encodes the correct global position.
+
+Similarly, the ``VarMaps`` array stored in ``ModGlueType`` makes the
+Jacobian coupling terms self-contained. During
+``BuildJacobianTC`` / ``BuildJacobianIO`` (in ``FAST_Solver.f90``) and
+``ModGlue_Linearize_OP`` (in ``FAST_ModGlue.f90``), the loop is simply:
+
+.. code-block:: fortran
+
+ do i = 1, size(ModGlue%VarMaps)
+ ! perturb source output, evaluate destination input change
+ ! scatter result into J(iVarDst%iLoc, iVarSrc%iLoc)
+ end do
+
+No module-specific knowledge is needed at this level — the ``ModGlueType``
+instance is fully self-describing.
diff --git a/docs/source/dev/glue-code/modvar.dot b/docs/source/user/glue-code/modvar.dot
similarity index 98%
rename from docs/source/dev/glue-code/modvar.dot
rename to docs/source/user/glue-code/modvar.dot
index 4e329730f3..ba821233e9 100644
--- a/docs/source/dev/glue-code/modvar.dot
+++ b/docs/source/user/glue-code/modvar.dot
@@ -47,6 +47,7 @@ digraph ModVarDiagram {
VF_DerivOrder1
2048
VF_DerivOrder2
4096
VF_Mapping
8192
+
VF_NoLin
16384
>;];
ModVarType [shape = plain;label = <
@@ -97,7 +98,6 @@ digraph ModVarDiagram {
dXdy(:,:)
R8Ki
dUdu(:,:)
R8Ki
dUdy(:,:)
R8Ki
-
StateRotation(:,:)
R8Ki
>;];
ModDataType [shape = plain;label = <
@@ -106,7 +106,9 @@ digraph ModVarDiagram {
iMod
IntKi
ID
IntKi
Ins
IntKi
+
iRotor
IntKi
SubSteps
IntKi
+
Category
IntKi
DT
R8Ki
Vars
ModVarsType
Lin
ModLinType
diff --git a/docs/source/user/glue-code/modvar.rst b/docs/source/user/glue-code/modvar.rst
new file mode 100644
index 0000000000..fe1b18677b
--- /dev/null
+++ b/docs/source/user/glue-code/modvar.rst
@@ -0,0 +1,459 @@
+.. _glue-code-modvar:
+
+Module Variables (``ModVar``)
+=============================
+
+The ``ModVar`` module (``modules/nwtc-library/src/ModVar.f90``) provides the data
+structures and subroutines that allow each physics module to declare its
+continuous states, inputs, and outputs in a way that the glue code can
+manipulate generically—without knowing the internals of any particular module.
+
+The complete type hierarchy—from ``DatLoc`` and ``ModVarType`` through
+``ModVarsType`` and ``ModDataType`` up to the top-level ``ModGlueType``—is
+illustrated in :numref:`fig-modvar-types`.
+
+.. contents::
+ :local:
+ :depth: 2
+
+Data structures
+---------------
+
+``DatLoc``
+~~~~~~~~~~
+
+A ``DatLoc`` (data location) is a small structure used to uniquely identify
+where a particular variable lives inside a module's derived-type hierarchy:
+
+.. code-block:: fortran
+
+ TYPE :: DatLoc
+ INTEGER :: Num ! Data identification number (from _Types.f90 file)
+ INTEGER :: i1 ! First index
+ INTEGER :: i2 ! Second index
+ INTEGER :: i3 ! Third index
+ INTEGER :: i4
+ INTEGER :: i5
+ END TYPE DatLoc
+
+A ``DatLoc`` value is created once per variable by the registry-generated
+``DatLoc()`` constructor calls inside each module's ``InitVars`` subroutine and
+passed to ``MV_AddVar`` / ``MV_AddMeshVar``. The glue code stores this value
+inside each ``ModVarType`` and uses ``MV_EqualDL`` to match variables across
+source and destination modules when setting up mesh mappings.
+
+``ModVarType``
+~~~~~~~~~~~~~~
+
+Describes a single variable (or group of variables for a mesh field):
+
+.. code-block:: fortran
+
+ TYPE :: ModVarType
+ INTEGER :: Field ! Field type (FieldForce, FieldTransDisp, ...)
+ INTEGER :: Nodes ! Number of nodes (mesh variables only)
+ INTEGER :: Num ! Total number of scalar values
+ INTEGER :: Flags ! Bit-mask of VF_* flags
+ INTEGER :: DerivOrder ! 0=disp/orientation, 1=velocity, 2=acceleration
+ INTEGER :: iLoc(2) ! [start, end] in module-local array
+ INTEGER :: iGlu(2) ! [start, end] in glue-level array
+ INTEGER :: iq(2) ! [start, end] in solver state (q) array
+ INTEGER :: iLB, iUB ! User array bounds (for array-valued scalars)
+ INTEGER :: j, k, m, n ! Additional user-defined indices
+ REAL(R8Ki) :: Perturb ! Perturbation size for Jacobian finite differences
+ TYPE(DatLoc) :: DL ! Data location
+ CHARACTER :: Name ! Human-readable variable name
+ CHARACTER(:), ALLOCATABLE :: LinNames(:) ! Per-value linearization labels
+ END TYPE ModVarType
+
+``ModVarsType``
+~~~~~~~~~~~~~~~
+
+Holds all variables for one module, partitioned into three groups:
+
+.. code-block:: fortran
+
+ TYPE :: ModVarsType
+ INTEGER :: Nx, Nu, Ny
+ TYPE(ModVarType), ALLOCATABLE :: x(:) ! Continuous states
+ TYPE(ModVarType), ALLOCATABLE :: u(:) ! Inputs
+ TYPE(ModVarType), ALLOCATABLE :: y(:) ! Outputs
+ END TYPE ModVarsType
+
+``ModDataType``
+~~~~~~~~~~~~~~~
+
+Top-level container that the glue code holds for each module instance:
+
+.. code-block:: fortran
+
+ TYPE :: ModDataType
+ CHARACTER :: Abbr ! Module abbreviation ("ED", "BD", ...)
+ INTEGER :: iMod ! Index in glue module array
+ INTEGER :: ID ! Module_ED, Module_BD, ...
+ INTEGER :: Ins ! Instance number
+ INTEGER :: iRotor ! Rotor index (0 = all rotors)
+ INTEGER :: SubSteps ! Module sub-steps per solver step
+ INTEGER :: Category ! Bit-mask of MC_* coupling flags
+ REAL(R8Ki) :: DT ! Module time step
+ TYPE(ModVarsType) :: Vars
+ TYPE(ModLinType) :: Lin
+ END TYPE ModDataType
+
+Variable flags (``VF_*``)
+-------------------------
+
+Flags are combined via ``IOR`` and tested with ``MV_HasFlagsAll`` /
+``MV_HasFlagsAny``.
+
+.. list-table::
+ :header-rows: 1
+ :widths: 25 10 65
+
+ * - Flag
+ - Value
+ - Meaning
+ * - ``VF_None``
+ - 0
+ - No flags set; used as a wildcard that matches any variable.
+ * - ``VF_Mesh``
+ - 1
+ - Variable is a mesh field (set automatically by ``MV_AddMeshVar``).
+ * - ``VF_Line``
+ - 2
+ - Mesh is a *line* mesh (loads per unit length); linearization labels get
+ ``/m`` suffix.
+ * - ``VF_RotFrame``
+ - 4
+ - Variable lives in the rotating reference frame.
+ * - ``VF_Linearize``
+ - 8
+ - Variable is included in the full-system linearization.
+ * - ``VF_ExtLin``
+ - 16
+ - Variable is included in extended linearization output.
+ * - ``VF_SmallAngle``
+ - 32
+ - Use small-angle approximation when computing orientation differences for
+ linearization.
+ * - ``VF_2PI``
+ - 64
+ - Scalar angle with range [0, 2π] (e.g. generator azimuth).
+ * - ``VF_WriteOut``
+ - 256
+ - Output variable associated with a ``WriteOutput`` channel.
+ * - ``VF_Solve``
+ - 512
+ - Variable participates in the tight-coupling Jacobian or input-output
+ convergence solve. Set automatically by ``FAST_SolverInit``;
+ module developers should not set this flag manually.
+ * - ``VF_AeroMap``
+ - 1024
+ - Variable used in aeromap computation.
+ * - ``VF_Mapping``
+ - 8192
+ - Variable participates in a module-to-module transfer mapping.
+ * - ``VF_NoLin``
+ - 16384
+ - Explicitly excludes a variable from both linearization and the solver
+ (overrides ``VF_Linearize`` and ``VF_Solve``).
+
+Field types (``Field*``)
+------------------------
+
+Used in the ``Field`` member of ``ModVarType`` and in ``MV_AddMeshVar``'s
+``Fields`` argument:
+
+.. list-table::
+ :header-rows: 1
+ :widths: 30 70
+
+ * - Constant
+ - Meaning
+ * - ``FieldForce``
+ - Nodal force (3 components per node, N)
+ * - ``FieldMoment``
+ - Nodal moment (3 components per node, N·m)
+ * - ``FieldTransDisp``
+ - Translational displacement (m)
+ * - ``FieldOrientation``
+ - Orientation, stored internally as unit-quaternion parameters (rad)
+ * - ``FieldTransVel``
+ - Translational velocity (m/s)
+ * - ``FieldAngularVel``
+ - Angular velocity (rad/s)
+ * - ``FieldTransAcc``
+ - Translational acceleration (m/s²)
+ * - ``FieldAngularAcc``
+ - Angular acceleration (rad/s²)
+ * - ``FieldAngularDisp``
+ - Angular displacement (rad)
+ * - ``FieldScalar``
+ - Generic scalar values
+
+Convenience arrays defined in ``ModVar.f90``:
+
+* ``LoadFields`` = ``[FieldForce, FieldMoment]``
+* ``TransFields`` = ``[FieldTransDisp, FieldTransVel, FieldTransAcc]``
+* ``AngularFields``= ``[FieldOrientation, FieldAngularVel, FieldAngularAcc, FieldAngularDisp]``
+* ``MotionFields`` = all translational and angular motion fields
+
+Adding module variables
+-----------------------
+
+Each module that participates in the glue code must implement an
+``InitVars`` (or equivalent) subroutine that populates a ``ModVarsType``
+structure by calling the ``MV_Add*`` subroutines documented below. This
+subroutine is called during the module's ``_Init`` routine and the resulting
+``Vars`` is passed immediately to ``MV_AddModule``.
+
+``MV_AddVar``
+~~~~~~~~~~~~~
+
+Adds a single (possibly multi-element) scalar variable to a variable array.
+
+.. code-block:: fortran
+
+ subroutine MV_AddVar(VarAry, Name, Field, DL, &
+ Num, iAry, jAry, kAry, &
+ Flags, DerivOrder, Perturb, LinNames, Active)
+
+.. list-table::
+ :header-rows: 1
+ :widths: 20 10 70
+
+ * - Argument
+ - Intent
+ - Description
+ * - ``VarAry``
+ - ``INOUT``
+ - Allocatable array of ``ModVarType``; the new variable is appended.
+ * - ``Name``
+ - ``IN``
+ - Human-readable name used in debug output and linearization labels.
+ * - ``Field``
+ - ``IN``
+ - Field type constant (``FieldScalar``, ``FieldTransDisp``, etc).
+ * - ``DL``
+ - ``IN``
+ - ``DatLoc`` identifying where the data live in the module's derived type.
+ * - ``Num``
+ - ``IN`` (optional)
+ - Number of scalar values. Defaults to 1. If 0, the call is a no-op.
+ * - ``iAry``
+ - ``IN`` (optional)
+ - Starting lower-bound index if the data are stored in an array.
+ * - ``jAry``, ``kAry``
+ - ``IN`` (optional)
+ - Second and third array indices (for 2-D or 3-D arrays).
+ * - ``Flags``
+ - ``IN`` (optional)
+ - Initial ``VF_*`` flag bit-mask. Defaults to ``VF_None``.
+ * - ``DerivOrder``
+ - ``IN`` (optional)
+ - Override the automatically-inferred derivative order (0, 1, or 2).
+ * - ``Perturb``
+ - ``IN`` (optional)
+ - Finite-difference perturbation magnitude. A good default is roughly
+ 1 % of the expected variable magnitude. For mesh fields the default
+ computed inside ``ModVarType_Init`` is used when this is omitted.
+ * - ``LinNames``
+ - ``IN`` (optional)
+ - Array of per-value linearization channel labels (length = ``Num``).
+ **Required** for non-mesh, non-scalar variables.
+ * - ``Active``
+ - ``IN`` (optional)
+ - Set to ``.false.`` to conditionally skip adding the variable.
+
+**Example** – registering the generator torque input of ServoDyn:
+
+.. code-block:: fortran
+
+ call MV_AddVar(Vars%u, 'Generator torque command', FieldScalar, &
+ DatLoc(SrvD_u_GenTrq), &
+ Flags=VF_Linearize, &
+ Perturb=1.0e3_R8Ki, & ! N*m
+ LinNames=['SrvD GenTrq, N*m'])
+
+``MV_AddMeshVar``
+~~~~~~~~~~~~~~~~~
+
+Adds all requested mesh fields for a single ``MeshType`` to a variable array.
+It is a convenience wrapper around ``MV_AddVar`` that iterates over the
+``Fields`` argument and skips fields that are absent from the committed mesh.
+
+.. code-block:: fortran
+
+ subroutine MV_AddMeshVar(VarAry, Name, Fields, DL, Mesh, &
+ Flags, Perturbs, Active, iVar)
+
+.. list-table::
+ :header-rows: 1
+ :widths: 20 10 70
+
+ * - Argument
+ - Intent
+ - Description
+ * - ``VarAry``
+ - ``INOUT``
+ - Allocatable array to append the new variable(s) to.
+ * - ``Name``
+ - ``IN``
+ - Base name for all fields on this mesh.
+ * - ``Fields``
+ - ``IN``
+ - Integer array of field-type constants; use the ``LoadFields``,
+ ``MotionFields``, etc. convenience parameters where appropriate.
+ * - ``DL``
+ - ``IN``
+ - ``DatLoc`` for this mesh within the module's data type.
+ * - ``Mesh``
+ - ``INOUT``
+ - The committed ``MeshType``. Its ``ID`` field is set to identify it in
+ mesh-mapping operations. The subroutine returns without adding anything
+ if the mesh has not been committed.
+ * - ``Flags``
+ - ``IN`` (optional)
+ - Extra ``VF_*`` flags added on top of ``VF_Mesh``.
+ * - ``Perturbs``
+ - ``IN`` (optional)
+ - Array of perturbation values, one per entry in ``Fields``.
+ * - ``Active``
+ - ``IN`` (optional)
+ - Conditionally disable the entire mesh variable registration.
+ * - ``iVar``
+ - ``OUT`` (optional)
+ - Returns the assigned mesh ``ID`` so the caller can store it for later
+ field look-ups.
+
+**Example** – registering ElastoDyn's blade-root output motion mesh:
+
+.. code-block:: fortran
+
+ call MV_AddMeshVar(Vars%y, 'BladeRootMotion', MotionFields, &
+ DatLoc(ED_y_BladeRootMotion, i), & ! blade i
+ p%BladeRootMotion(i), &
+ Flags=VF_Linearize + VF_RotFrame)
+
+``MV_AddModule``
+~~~~~~~~~~~~~~~~
+
+After a module's ``InitVars`` subroutine is complete, the caller registers the
+module with the glue code using ``MV_AddModule``.
+
+.. code-block:: fortran
+
+ subroutine MV_AddModule(ModDataAry, ModID, ModAbbr, Instance, &
+ ModDT, SolverDT, Vars, Linearize, &
+ ErrStat, ErrMsg, iRotor)
+
+.. list-table::
+ :header-rows: 1
+ :widths: 20 10 70
+
+ * - Argument
+ - Intent
+ - Description
+ * - ``ModDataAry``
+ - ``INOUT``
+ - Allocatable array of ``ModDataType``; the new entry is appended.
+ * - ``ModID``
+ - ``IN``
+ - Module identifier constant (``Module_ED``, ``Module_BD``, etc.).
+ * - ``ModAbbr``
+ - ``IN``
+ - Short abbreviation string used in output labels (``"ED"``, ``"BD"``).
+ * - ``Instance``
+ - ``IN``
+ - Instance number (1-based). Most modules have a single instance.
+ * - ``ModDT``
+ - ``IN``
+ - Module time step (seconds). Must be an exact integer divisor of
+ ``SolverDT``.
+ * - ``SolverDT``
+ - ``IN``
+ - Solver (global) time step.
+ * - ``Vars``
+ - ``IN``
+ - Populated ``ModVarsType`` from the module's ``InitVars`` call.
+ * - ``Linearize``
+ - ``IN``
+ - Whether linearization is enabled. When ``.false.``, ``LinNames``
+ arrays are deallocated to save memory.
+ * - ``ErrStat`` / ``ErrMsg``
+ - ``OUT``
+ - Error status and message.
+ * - ``iRotor``
+ - ``IN`` (optional)
+ - Rotor number for multi-rotor turbines (0 = all rotors).
+
+**Sub-stepping logic**: if ``ModDT < SolverDT``, ``MV_AddModule`` calculates
+``ModData%SubSteps = NINT(SolverDT/ModDT)`` and validates that the module DT
+divides the solver DT exactly. An error is returned if ``ModDT > SolverDT``.
+
+**Typical call sequence inside FAST_Subs.f90**:
+
+.. code-block:: fortran
+
+ ! Module computes its own Vars in Init
+ call ED_Init(InitInp, u, p, ..., InitOut, ErrStat, ErrMsg)
+ ! Register with glue code
+ call MV_AddModule(m%ModData, Module_ED, 'ED', 1, p%DT, p_FAST%DT, &
+ InitOut%Vars, p_FAST%Linearize, ErrStat, ErrMsg, iRotor=1)
+
+``MV_InitVarsJac``
+~~~~~~~~~~~~~~~~~~
+
+Called inside each module's ``InitVars`` after all ``MV_AddVar`` /
+``MV_AddMeshVar`` calls are complete. It assigns the module-local ``iLoc``
+index ranges to each variable and allocates the ``ModJacType`` working arrays
+used during Jacobian calculations.
+
+.. code-block:: fortran
+
+ subroutine MV_InitVarsJac(Vars, Jac, Linearize, ErrStat, ErrMsg)
+
+Perturbation values
+-------------------
+
+Every variable carries a ``Perturb`` value used for central-difference
+finite-differencing when building module-level Jacobians. The ``Perturb``
+argument to ``MV_AddVar`` / ``MV_AddMeshVar`` should be chosen so that the
+resulting output change is large enough to distinguish from numerical noise
+but small enough to stay in the linear regime. Typical values:
+
+* Translational displacement: ``1.0e-4`` m
+* Rotational (orientation): ``2.0e-5`` rad
+* Translational velocity: ``1.0e-3`` m/s
+* Angular velocity: ``2.0e-4`` rad/s
+* Translational acceleration: ``1.0e-2`` m/s²
+* Force: ``1.0e1`` N
+* Moment: ``1.0e1`` N·m
+* Generic scalar: context-dependent
+
+The ``UJacSclFact`` input parameter (see :ref:`glue-code-solver-inputs`) is a
+global conditioning factor that the solver applies to load variables in the
+Jacobian to improve matrix conditioning when force/moment magnitudes are very
+different from state magnitudes.
+
+Orientation representation
+--------------------------
+
+Orientations are **not** stored or manipulated as direction cosine matrices
+(DCMs) inside the glue-variable arrays. Instead, a compact three-component
+unit-quaternion parameterization is used:
+
+.. math::
+
+ \mathbf{q}_p = [q_1, q_2, q_3] \quad \text{where } q_0 = \sqrt{1 - q_1^2 - q_2^2 - q_3^2}
+
+This parameterization avoids the redundancy in a full DCM and enables
+straightforward finite-differencing via quaternion composition
+(``quat_compose``). Conversion utilities exported from ``ModVar`` include
+``dcm_to_quat``, ``quat_to_dcm``, ``quat_compose``, ``quat_inv``,
+``quat_to_rvec``, ``rvec_to_quat``, ``wm_to_quat``, and ``quat_to_wm``.
+
+When computing orientation differences for Jacobian rows, ``MV_ComputeDiff``
+computes the relative rotation between the negative and positive-perturbation
+quaternions and converts it to a rotation vector (small-angle approximation
+or full Rodrigues formula depending on ``VF_SmallAngle``).
diff --git a/docs/source/user/glue-code/overview.rst b/docs/source/user/glue-code/overview.rst
new file mode 100644
index 0000000000..e775961555
--- /dev/null
+++ b/docs/source/user/glue-code/overview.rst
@@ -0,0 +1,103 @@
+.. _glue-code-overview:
+
+Glue Code Overview
+==================
+
+The OpenFAST glue code is the software layer that connects individual physics modules
+(ElastoDyn, AeroDyn, HydroDyn, ServoDyn, etc.) into a coupled simulation. It lives
+primarily in ``modules/openfast-library/src/`` and relies on the ``ModVar`` module in
+``modules/nwtc-library/src/ModVar.f90`` for an abstract description of every
+variable exchanged between modules.
+
+High-level responsibilities include:
+
+* Initialising each module and registering its variables with the glue code
+* Managing multi-rate sub-stepping (modules whose time step is a divisor of the
+ global time step)
+* Mapping outputs of one module to inputs of another (motion meshes, load meshes,
+ and scalar variables)
+* Running the time-stepping loop under either loose coupling or tight
+ generalized-alpha coupling
+* Linearising the assembled system to produce state-space matrices
+
+Source files
+------------
+
+.. list-table::
+ :header-rows: 1
+ :widths: 35 65
+
+ * - File
+ - Purpose
+ * - ``FAST_Subs.f90``
+ - Top-level initialisation: reads input file, calls each module's ``_Init``,
+ and calls ``MV_AddModule`` to register every module with the glue code.
+ * - ``FAST_ModGlue.f90``
+ - Combines per-module variable descriptions into a monolithic ``ModGlueType``
+ structure via ``ModGlue_CombineModules``; performs linearization
+ (``ModGlue_Linearize_OP``) and steady-state trimming
+ (``ModGlue_CalcSteady``).
+ * - ``FAST_Solver.f90``
+ - Implements the generalized-alpha tight-coupling solver
+ (``FAST_SolverStep``), the input-output convergence loop
+ (``FAST_CalcOutputsAndSolveForInputs``), and Jacobian assembly.
+ * - ``FAST_Mapping.f90``
+ - Mesh-to-mesh and variable-to-variable transfer mappings.
+ * - ``FAST_Funcs.f90``
+ - Wrappers around module-level ``CalcOutput``, ``UpdateStates``,
+ ``CalcContStateDeriv``, ``GetOperatingPoint``, and ``SetOperatingPoint``
+ that dispatch to the correct module instance.
+ * - ``ModVar.f90`` (nwtc-library)
+ - The ``ModVar`` module: data structures (``ModVarType``, ``ModVarsType``,
+ ``ModDataType``, ``DatLoc``) and all ``MV_*`` subroutines.
+
+Module coupling categories
+--------------------------
+
+Each module is assigned to exactly one coupling category during initialisation in
+``FAST_SolverInit``:
+
+.. list-table::
+ :header-rows: 1
+ :widths: 20 15 65
+
+ * - Category
+ - Flag
+ - Description
+ * - Tight Coupling (TC)
+ - ``MC_Tight``
+ - States and accelerations are solved simultaneously via the generalized-alpha
+ Newton iteration. ElastoDyn, BeamDyn, and SubDyn are tight-coupling modules
+ when ``ModCoupling`` ≥ 2.
+ * - Option 1
+ - ``MC_Option1``
+ - Modules whose inputs depend on TC outputs and are converged in the same
+ Newton loop (e.g. HydroDyn, MoorDyn, ServoDyn with structural controllers).
+ * - Option 2
+ - ``MC_Option2``
+ - Loosely coupled modules that are called once per step before the convergence
+ loop (InflowWind, SeaState, AeroDyn, etc.).
+ * - Post
+ - ``MC_Post``
+ - Modules whose input solve is deferred until after the convergence loop
+ (ServoDyn, ExternalInflow).
+
+Time-stepping loop (overview)
+------------------------------
+
+Each call to ``FAST_SolverStep`` follows this sequence:
+
+1. **Correction iterations** (outer loop) – at most ``p%MaxConvIter`` iterations.
+2. **Option 2** – input solve + state update + ``CalcOutput`` for loosely coupled modules.
+3. **Option 1** – input solve + state update for semi-implicit modules.
+4. **TC input solve** – gather inputs for the tight-coupling modules.
+5. **Convergence iterations** (inner loop) – Newton-Raphson updates of TC
+ states and inputs until the update norm falls below ``ConvTol`` or the
+ iteration limit is reached.
+6. **Post-solve input solves** – ServoDyn, ExternalInflow.
+
+Module registration and variable ordering are described in detail in
+:ref:`glue-code-modvar`. How the per-module variables are assembled into
+global arrays and Jacobian matrices is covered in :ref:`glue-code-modglue`.
+The solver algorithm and Jacobian construction are covered in
+:ref:`glue-code-solver`.
diff --git a/docs/source/user/glue-code/solver.rst b/docs/source/user/glue-code/solver.rst
new file mode 100644
index 0000000000..62bf2a6256
--- /dev/null
+++ b/docs/source/user/glue-code/solver.rst
@@ -0,0 +1,338 @@
+.. _glue-code-solver:
+
+Solver
+======
+
+The OpenFAST tight-coupling solver is implemented in
+``modules/openfast-library/src/FAST_Solver.f90``. It integrates the continuous
+states and resolves the input-output coupling between modules using a
+generalized-alpha scheme with Newton-Raphson convergence iterations.
+
+.. contents::
+ :local:
+ :depth: 2
+
+.. _glue-code-solver-inputs:
+
+User input parameters
+---------------------
+
+All solver parameters are set in the main OpenFAST input file
+(``*.fst``) under the **Feature Switches and Flags** and
+**Tight-Coupling / Solver** sections.
+
+.. list-table::
+ :header-rows: 1
+ :widths: 22 12 66
+
+ * - Parameter
+ - Type
+ - Description
+ * - ``DT``
+ - real
+ - Global (solver) time step in seconds. All module time steps must be
+ equal to or an integer sub-divisor of ``DT``.
+ * - ``ModCoupling``
+ - integer
+ - Coupling method.
+
+ * ``1`` – Loose coupling: structural modules (ED/BD/SD) are treated as
+ Option 1 and do **not** participate in the tight Newton loop.
+ * ``2`` – Tight coupling with fixed Jacobian updates (``DT_UJac``
+ controls update frequency).
+ * ``3`` – Tight coupling with adaptive Jacobian updates (the Jacobian is
+ rebuilt whenever the Newton loop fails to converge within the
+ iteration budget).
+ * - ``RhoInf``
+ - real
+ - Numerical damping parameter ρ∞ for the generalized-alpha integrator.
+ Range [0, 1]; 1 = no numerical damping (second-order accurate), 0 =
+ maximum damping (first-order accurate). Typical value: **0.9**.
+ Reducing ``RhoInf`` below 1 damps high-frequency numerical noise at the
+ cost of slightly reduced accuracy.
+ * - ``MaxConvIter``
+ - integer
+ - Maximum number of Newton convergence iterations per time step before the
+ solver declares convergence failure. Typical value: **20**.
+ With ``ModCoupling=2`` or ``1``, a fatal error is issued on failure;
+ with ``ModCoupling=3`` the Jacobian is rebuilt first and the step is
+ retried before a warning is emitted.
+ * - ``ConvTol``
+ - real
+ - Convergence tolerance. The iteration stops when the average
+ `L2`-norm of the Newton update vector falls below this value.
+ Typical value: ``1.0e-4``. Tighter tolerances increase
+ computational cost but may be needed for stiff problems.
+ * - ``DT_UJac``
+ - real
+ - Time interval (seconds) between Jacobian rebuilds when
+ ``ModCoupling=2``.
+
+ * If ``DT_UJac < DT``: the Jacobian is rebuilt at a fraction of the
+ convergence-iteration budget.
+ * If ``DT_UJac ≥ DT``: the Jacobian is rebuilt every
+ ``CEILING(DT_UJac/DT)`` time steps.
+ * Setting ``DT_UJac`` very large (e.g. ``9999``) freezes the Jacobian
+ for the entire simulation; useful for profiling or when the system
+ is nearly linear and the Jacobian is expensive.
+ * - ``UJacSclFact``
+ - real
+ - Conditioning scale factor applied to load rows and columns of the
+ Jacobian. Force and moment variables are divided by this factor before
+ the linear solve and multiplied back afterwards, equalising the magnitude
+ of load entries relative to displacement/velocity entries. Typical
+ value: **1.0e5** for offshore systems; may need adjustment for very
+ large or very small turbines.
+ * - ``CompElast``
+ - integer
+ - Select the structural dynamics module: ``1`` = ElastoDyn,
+ ``2`` = BeamDyn (blades only, ElastoDyn still handles the tower/platform),
+ ``3`` = Simplified ElastoDyn. The chosen modules become TC members when
+ ``ModCoupling ≥ 2``.
+ * - ``CompSub``
+ - integer
+ - Sub-structural module: ``0`` = none, ``1`` = SubDyn, ``2`` = ExtPtfm,
+ ``3`` = SlD (SoilDyn). SubDyn joins the TC set when
+ ``ModCoupling ≥ 2``.
+ * - ``CompHydro``
+ - integer
+ - ``0`` = none, ``1`` = HydroDyn. HydroDyn is always Option 1.
+ * - ``CompMooring``
+ - integer
+ - ``0`` = none, ``1`` = MAP++, ``2`` = FEAMooring, ``3`` = MoorDyn,
+ ``4`` = OrcaFlex. Mooring modules are always Option 1.
+ * - ``CompAero``
+ - integer
+ - Aerodynamics module: ``0`` = none, ``1`` = AeroDisk, ``2`` = AeroDyn.
+ AeroDyn is Option 2 for land-based turbines and Option 1 for MHK.
+ * - ``CompServo``
+ - integer
+ - Controller module: ``0`` = none, ``1`` = ServoDyn. ServoDyn is
+ Post-solve by default but becomes Option 1 when structural controllers
+ (tower, blade, nacelle StC) are active.
+
+Generalized-alpha integration
+------------------------------
+
+The tight-coupling solver integrates second-order ODEs of the form
+
+.. math::
+
+ \mathbf{M}\,\ddot{\mathbf{q}} + \mathbf{f}(\mathbf{q}, \dot{\mathbf{q}}, t) = 0
+
+using the **generalized-alpha method** (Chung & Hulbert, 1993). Given the
+spectral radius ρ∞ specified by ``RhoInf``, the method parameters are:
+
+.. math::
+
+ \alpha_m &= \frac{2\rho_\infty - 1}{\rho_\infty + 1} \\
+ \alpha_f &= \frac{\rho_\infty}{\rho_\infty + 1} \\
+ \gamma &= \tfrac{1}{2} - \alpha_m + \alpha_f \\
+ \beta &= \tfrac{1}{4}(1 - \alpha_m + \alpha_f)^2
+
+Two derived coefficients used throughout the convergence loop are:
+
+.. math::
+
+ \beta' &= h^2 \beta \frac{1 - \alpha_f}{1 - \alpha_m} \\
+ \gamma' &= h \gamma \frac{1 - \alpha_f}{1 - \alpha_m}
+
+where *h* = ``DT``.
+
+**State vector layout** – the solver maintains a per-module *generalized
+coordinate* (q) vector with four columns:
+
+.. list-table::
+ :header-rows: 1
+ :widths: 15 85
+
+ * - Column
+ - Meaning
+ * - ``q``
+ - Displacement / orientation states (``DerivOrder = 0``)
+ * - ``v``
+ - Velocity states (``DerivOrder = 1``)
+ * - ``vd``
+ - Acceleration (physical, from module ``CalcContStateDeriv``)
+ * - ``a``
+ - Algorithmic acceleration (generalized-alpha internal variable)
+
+State prediction at the start of each step:
+
+.. math::
+
+ q_{n+1}^{\rm pred} &= q_n + h v_n + h^2[(\tfrac{1}{2} - \beta)a_n + \beta\, a_{n+1}] \\
+ v_{n+1}^{\rm pred} &= v_n + h[(1-\gamma)a_n + \gamma\, a_{n+1}]
+
+Module ordering
+---------------
+
+During ``FAST_SolverInit`` each module is categorised based on ``ModCoupling``
+and its own physics type, and assigned to one of the ordered index arrays
+in the ``Glue_TCParam`` structure:
+
+.. list-table::
+ :header-rows: 1
+ :widths: 20 80
+
+ * - Array
+ - Modules (in order)
+ * - ``iModTC``
+ - ElastoDyn, BeamDyn, SubDyn (when ``ModCoupling ≥ 2``)
+ * - ``iModOpt1``
+ - ServoDyn (when StC active), SED, AD (MHK), ExtPtfm, HydroDyn, OrcaFlex,
+ MoorDyn; ED/BD/SD also appear here when ``ModCoupling = 1``
+ * - ``iModOpt2``
+ - ServoDyn, SED, ED, BD, SD, InflowWind, SeaState, AeroDyn (land),
+ AeroDisk, ExtLoads, MAP++, FEAMooring, IceDyn, IceFloe, SoilDyn
+ * - ``iModPost``
+ - ServoDyn, ExternalInflow
+ * - ``iModInit``
+ - SED, ED, BD, SD, InflowWind, ExtLoads (Step 0 initialisation only)
+
+Jacobian construction
+---------------------
+
+Two separate Jacobians are assembled:
+
+1. **TC/Option-1 Jacobian** (``BuildJacobianTC``) — for the main time-stepping
+ convergence loop.
+2. **IO Jacobian** (``BuildJacobianIO``) — for the initial and linearization
+ input-output solve.
+
+Variable selection (``VF_Solve`` flag)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+During ``FAST_SolverInit → SetVarSolveFlags``, the ``VF_Solve`` flag is set on
+the variables that must appear in the Jacobian:
+
+* **Continuous states** of all TC modules (automatically).
+* **Motion mesh** inputs/outputs of TC-to-TC mappings (all fields).
+* **Motion mesh** input accelerations of TC-to-Option1 or
+ Option1-to-TC mappings.
+* **Load mesh** inputs and outputs involved in any TC/Option1 mapping.
+* **Load mesh** displacement outputs of the destination module when the mapping
+ carries moments (needed for moment-arm Jacobian terms).
+* **Variable-to-variable** mapped inputs/outputs of TC/Option1 modules.
+* Any variable with ``VF_NoLin`` is excluded from ``VF_Solve``.
+
+Jacobian structure (TC Jacobian)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The assembled TC Jacobian **J** has size ``NumJ × NumJ``, where:
+
+.. math::
+
+ N_J = \underbrace{N_Q}_{\text{TC states}} +
+ \underbrace{N_{U_T}}_{\text{TC inputs}} +
+ \underbrace{N_{U_1}}_{\text{Option-1 inputs}}
+
+The columns and rows are partitioned as:
+
+.. math::
+
+ \mathbf{J} = \begin{bmatrix}
+ J_{11} & J_{12} \\
+ J_{21} & J_{22}
+ \end{bmatrix}
+
+where
+
+* **J₁₁** (``NumQ × NumQ``) — derivative of the acceleration residual with
+ respect to TC displacement/velocity states (formed from the module
+ ``dXdx`` sub-Jacobians plus the generalized-alpha tangent).
+* **J₁₂** (``NumQ × NumU_T``) — derivative of the acceleration residual with
+ respect to TC inputs (from ``dXdu``).
+* **J₂₁** (``NumU_T × NumQ``) — derivative of the input residual with
+ respect to TC states (from ``dUdx = dUdy · dydx``).
+* **J₂₂** (``NumU × NumU``) — derivative of the input residual with
+ respect to inputs, including load conditioning rows/columns.
+
+The right-hand side (XB) contains the residuals:
+
+* **State residual** (rows ``iJX``): difference between the predicted
+ velocity derivative and the module-computed accelerations.
+* **Input residual** (rows ``iJU``): difference between the inputs computed
+ from mesh mappings (``FAST_InputSolve``) and the current iterate.
+
+The loads portion (rows ``iJL``) is pre-divided by ``UJacSclFact`` before the
+factorisation to improve conditioning.
+
+Jacobian update strategy
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+``ModCoupling = 2`` (fixed updates)
+ The Jacobian is rebuilt if either of these counters reaches zero:
+
+ * ``UJacStepsRemain`` — steps remaining; initialised to
+ ``CEILING(DT_UJac/DT)`` each time the Jacobian is rebuilt.
+ * ``UJacIterRemain`` — iteration budget; initialised to
+ ``CEILING(DT_UJac/DT · MaxConvIter)`` when ``DT_UJac < DT``.
+
+ On convergence failure the solver returns a fatal error immediately.
+
+``ModCoupling = 3`` (adaptive updates)
+ The Jacobian is rebuilt the first time the convergence loop fails. If the
+ step still does not converge after the forced rebuild, a non-fatal warning
+ is issued and the simulation proceeds.
+
+Per-module Jacobian contributions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The module-level Jacobian sub-matrices are computed by finite differencing
+inside ``BuildJacobianTC`` and ``BuildJacobianIO`` using the ``MV_Perturb`` /
+``MV_ComputeDiff`` / ``MV_ComputeCentralDiff`` utilities from ``ModVar``. For
+each variable flagged ``VF_Solve``:
+
+1. Apply a positive perturbation of magnitude ``Var%Perturb`` to the working
+ state/input array.
+2. Call ``FAST_CalcOutput`` (or ``FAST_GetContStateDeriv``).
+3. Apply an equal negative perturbation.
+4. Call again.
+5. Compute the central difference: ``(y_plus - y_minus) / (2·Perturb)``.
+
+For orientation variables (``FieldOrientation``), perturbations are applied by
+quaternion composition rather than direct addition (``MV_Perturb``), and
+differences are extracted as rotation vectors (``MV_ComputeDiff``).
+
+Linear solve
+~~~~~~~~~~~~
+
+The LU factorisation of **J** is computed with ``LAPACK_getrf`` and the system
+is solved with ``LAPACK_getrs`` (packed in ``NWTC_LAPACK``). The same
+factored matrix is reused across convergence iterations until the update
+strategy triggers a rebuild.
+
+Convergence check
+~~~~~~~~~~~~~~~~~
+
+After each Newton step the convergence error is the average `L2`-norm of the
+update vector:
+
+.. math::
+
+ e = \frac{\|\Delta \mathbf{z}\|_2}{N_J}
+
+where :math:`\Delta \mathbf{z}` combines state and input updates. The loop
+exits if ``e < ConvTol`` (``ErrID_None``) or the iteration count reaches
+``MaxConvIter`` (``ErrID_Fatal`` or ``ErrID_Warn`` depending on
+``ModCoupling``).
+
+Output channels from the solver
+--------------------------------
+
+Three output channels are written to ``DriverWriteOutput`` each step and
+appear in the output file when enabled:
+
+.. list-table::
+ :header-rows: 1
+ :widths: 20 80
+
+ * - Index
+ - Content
+ * - 1
+ - Total convergence iterations in the step (``TotalIter``)
+ * - 2
+ - Final convergence error (``ConvError``)
+ * - 3
+ - Number of Jacobian rebuilds in the step (``NumUJac``)
diff --git a/docs/source/user/index.rst b/docs/source/user/index.rst
index f06e595df3..eb331f8113 100644
--- a/docs/source/user/index.rst
+++ b/docs/source/user/index.rst
@@ -14,6 +14,7 @@ This section contains documentation for the OpenFAST module-coupling environment
:maxdepth: 1
General considerations
+ Glue Code
AeroDyn
OLAF
Aeroacoustics
diff --git a/modules/nwtc-library/src/NWTC_Library_Types.f90 b/modules/nwtc-library/src/NWTC_Library_Types.f90
index 72489ece86..93fa73a811 100644
--- a/modules/nwtc-library/src/NWTC_Library_Types.f90
+++ b/modules/nwtc-library/src/NWTC_Library_Types.f90
@@ -62,7 +62,7 @@ MODULE NWTC_Library_Types
INTEGER(IntKi), PUBLIC, PARAMETER :: VF_DerivOrder1 = 2048 ! Variable is derivative order 1 in linearization file [-]
INTEGER(IntKi), PUBLIC, PARAMETER :: VF_DerivOrder2 = 4096 ! Variable is derivative order 2 in linearization file [-]
INTEGER(IntKi), PUBLIC, PARAMETER :: VF_Mapping = 8192 ! Variable is used in a module-to-module transfer mapping [-]
- INTEGER(IntKi), PUBLIC, PARAMETER :: VF_NoLin = 16384 ! Variable is used in a module-to-module transfer mapping [-]
+ INTEGER(IntKi), PUBLIC, PARAMETER :: VF_NoLin = 16384 ! Flag to exclude variable from linearization or solver [-]
INTEGER(IntKi), PUBLIC, PARAMETER :: MC_None = 0 ! [-]
INTEGER(IntKi), PUBLIC, PARAMETER :: MC_Tight = 1 ! [-]
INTEGER(IntKi), PUBLIC, PARAMETER :: MC_Option1 = 2 ! [-]
@@ -121,7 +121,7 @@ MODULE NWTC_Library_Types
! =======================
! ========= DatLoc =======
TYPE, PUBLIC :: DatLoc
- INTEGER(IntKi) :: Num = 0 !< Mesh number in module [-]
+ INTEGER(IntKi) :: Num = 0 !< Field number in module input, output, or continuous state derived data types [-]
INTEGER(IntKi) :: i1 = 0 !< Index 1 [-]
INTEGER(IntKi) :: i2 = 0 !< Index 2 [-]
INTEGER(IntKi) :: i3 = 0 !< Index 3 [-]
diff --git a/modules/nwtc-library/src/Registry_NWTC_Library.txt b/modules/nwtc-library/src/Registry_NWTC_Library.txt
index 48e3f779db..e63a1498a6 100644
--- a/modules/nwtc-library/src/Registry_NWTC_Library.txt
+++ b/modules/nwtc-library/src/Registry_NWTC_Library.txt
@@ -72,7 +72,7 @@ param ^ - IntKi VF_AeroMap - 1024 -
param ^ - IntKi VF_DerivOrder1 - 2048 - "Variable is derivative order 1 in linearization file" -
param ^ - IntKi VF_DerivOrder2 - 4096 - "Variable is derivative order 2 in linearization file" -
param ^ - IntKi VF_Mapping - 8192 - "Variable is used in a module-to-module transfer mapping" -
-param ^ - IntKi VF_NoLin - 16384 - "Variable is used in a module-to-module transfer mapping" -
+param ^ - IntKi VF_NoLin - 16384 - "Flag to exclude variable from linearization or solver" -
param ^ - IntKi MC_None - 0 - "" -
param ^ - IntKi MC_Tight - 1 - "" -
@@ -80,7 +80,7 @@ param ^ - IntKi MC_Option1 - 2 -
param ^ - IntKi MC_Option2 - 4 - "" -
param ^ - IntKi MC_Post - 8 - "" -
-typedef ^ DatLoc IntKi Num - 0 - "Mesh number in module"
+typedef ^ DatLoc IntKi Num - 0 - "Field number in module input, output, or continuous state derived data types"
typedef ^ ^ IntKi i1 - 0 - "Index 1"
typedef ^ ^ IntKi i2 - 0 - "Index 2"
typedef ^ ^ IntKi i3 - 0 - "Index 3"
diff --git a/modules/nwtc-library/src/Registry_NWTC_Library_base.txt b/modules/nwtc-library/src/Registry_NWTC_Library_base.txt
index 7f0449cbeb..ce6839eea3 100644
--- a/modules/nwtc-library/src/Registry_NWTC_Library_base.txt
+++ b/modules/nwtc-library/src/Registry_NWTC_Library_base.txt
@@ -72,7 +72,7 @@ param ^ - IntKi VF_AeroMap - 1024 -
param ^ - IntKi VF_DerivOrder1 - 2048 - "Variable is derivative order 1 in linearization file" -
param ^ - IntKi VF_DerivOrder2 - 4096 - "Variable is derivative order 2 in linearization file" -
param ^ - IntKi VF_Mapping - 8192 - "Variable is used in a module-to-module transfer mapping" -
-param ^ - IntKi VF_NoLin - 16384 - "Variable is used in a module-to-module transfer mapping" -
+param ^ - IntKi VF_NoLin - 16384 - "Flag to exclude variable from linearization or solver" -
param ^ - IntKi MC_None - 0 - "" -
param ^ - IntKi MC_Tight - 1 - "" -
@@ -80,7 +80,7 @@ param ^ - IntKi MC_Option1 - 2 -
param ^ - IntKi MC_Option2 - 4 - "" -
param ^ - IntKi MC_Post - 8 - "" -
-typedef ^ DatLoc IntKi Num - 0 - "Mesh number in module"
+typedef ^ DatLoc IntKi Num - 0 - "Field number in module input, output, or continuous state derived data types"
typedef ^ ^ IntKi i1 - 0 - "Index 1"
typedef ^ ^ IntKi i2 - 0 - "Index 2"
typedef ^ ^ IntKi i3 - 0 - "Index 3"