diff --git a/pygsti/baseobjs/statespace.py b/pygsti/baseobjs/statespace.py
index 277edb0f4..b99fb9da5 100644
--- a/pygsti/baseobjs/statespace.py
+++ b/pygsti/baseobjs/statespace.py
@@ -708,6 +708,7 @@ def qudit_udims(self):
         """Integer Hilbert (unitary operator) space dimensions of the qudits in ths quantum state space."""
         return self._qudit_udims
 
+
     @property
     def udim(self):
         """
diff --git a/pygsti/modelmembers/instruments/__init__.py b/pygsti/modelmembers/instruments/__init__.py
index 8ef2c0ed2..002fd230e 100644
--- a/pygsti/modelmembers/instruments/__init__.py
+++ b/pygsti/modelmembers/instruments/__init__.py
@@ -21,7 +21,8 @@
 
 
 def instrument_type_from_op_type(op_type):
-    """Decode an op type into an appropriate instrument type.
+    """
+    Decode an op type into an appropriate instrument type.
 
     Parameters:
     -----------
@@ -38,13 +39,13 @@ def instrument_type_from_op_type(op_type):
     # Limited set (only matching what is in convert)
     instr_conversion = {
         'auto': 'full',
-        'static unitary': 'static unitary',
-        'static clifford': 'static clifford',
+        'static unitary': 'static',
+        'static clifford': 'static',
         'static': 'static',
         'full': 'full',
         'full TP': 'full TP',
-        'full CPTP': 'full CPTP',
-        'full unitary': 'full unitary',
+        'full CPTP': 'full TP',
+        'full unitary': 'full',
     }
 
     instr_type_preferences = []
@@ -71,6 +72,8 @@ def instrument_type_from_op_type(op_type):
 
     return instr_type_preferences
 
+def create_from_instrument_effects():
+    pass
 
 def convert(instrument, to_type, basis, ideal_instrument=None, flatten_structure=False):
     """
@@ -121,7 +124,7 @@ def convert(instrument, to_type, basis, ideal_instrument=None, flatten_structure
 
             if to_type == "full TP":
                 return TPInstrument(list(instrument.items()), instrument.evotype, instrument.state_space)
-            elif to_type in ("full", "static", "static unitary"):
+            elif to_type in ("full", "static"):
                 from ..operations import convert as _op_convert
                 ideal_items = dict(ideal_instrument.items()) if (ideal_instrument is not None) else {}
                 members = [(k, _op_convert(g, to_type, basis, ideal_items.get(k, None), flatten_structure))
diff --git a/pygsti/modelmembers/operations/embeddedop.py b/pygsti/modelmembers/operations/embeddedop.py
index da014dc03..78ff2a405 100644
--- a/pygsti/modelmembers/operations/embeddedop.py
+++ b/pygsti/modelmembers/operations/embeddedop.py
@@ -20,6 +20,7 @@
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 from pygsti.baseobjs.errorgenlabel import GlobalElementaryErrorgenLabel as _GlobalElementaryErrorgenLabel, LocalElementaryErrorgenLabel as _LocalElementaryErrorgenLabel
 from pygsti import SpaceT
+
 class EmbeddedOp(_LinearOperator):
     """
     An operation containing a single lower (or equal) dimensional operation within it.
diff --git a/pygsti/models/cloudnoisemodel.py b/pygsti/models/cloudnoisemodel.py
index e6aefe8a7..761b78df4 100644
--- a/pygsti/models/cloudnoisemodel.py
+++ b/pygsti/models/cloudnoisemodel.py
@@ -18,7 +18,7 @@
 import scipy.sparse as _sps
 
 from pygsti.baseobjs import statespace as _statespace
-from pygsti.models.implicitmodel import ImplicitOpModel as _ImplicitOpModel, _init_spam_layers
+from pygsti.models.implicitmodel import ImplicitOpModel as _ImplicitOpModel
 from pygsti.models.layerrules import LayerRules as _LayerRules
 from pygsti.models.memberdict import OrderedMemberDict as _OrderedMemberDict
 from pygsti.evotypes import Evotype as _Evotype
@@ -306,7 +306,7 @@ def __init__(self, processor_spec, gatedict,
                     # used to specify which gate parameters should be amplifiable by germs for a given cloud (?)
                     # TODO CHECK THIS
 
-        _init_spam_layers(self, prep_layers, povm_layers)  # SPAM
+        self._init_spam_layers(self, prep_layers, povm_layers)  # SPAM
 
         printer.log("DONE! - created Model with nqudits=%d and op-blks=" % self.state_space.num_qudits)
         for op_blk_lbl, op_blk in self.operation_blks.items():
diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py
index 5d276d9e8..cafb33700 100644
--- a/pygsti/models/explicitmodel.py
+++ b/pygsti/models/explicitmodel.py
@@ -137,7 +137,9 @@ def flagfn(typ): return {'auto_embed': True, 'match_parent_statespace': True,
         self.preps = _OrderedMemberDict(self, default_prep_type, prep_prefix, flagfn("state"))
         self.povms = _OrderedMemberDict(self, default_povm_type, povm_prefix, flagfn("povm"))
         self.operations = _OrderedMemberDict(self, default_gate_type, gate_prefix, flagfn("operation"))
-        self.instruments = _OrderedMemberDict(self, default_instrument_type, instrument_prefix, flagfn("instrument"))
+        self.instruments = _OrderedMemberDict(self, default_instrument_type, instrument_prefix, 
+                                              {'auto_embed': False, 'match_parent_statespace': True,
+                                               'match_parent_evotype': True, 'cast_to_type': 'instrument'})
         self.factories = _OrderedMemberDict(self, default_gate_type, gate_prefix, flagfn("factory"))
         self.effects_prefix = effect_prefix
         self._default_gauge_group = None
@@ -222,6 +224,8 @@ def _embed_operation(self, op_target_labels, op_val, force=False):
             return op_val  # if gate operates on full dimension, no need to embed.
 
         return _op.EmbeddedOp(self.state_space, op_target_labels, op_val)
+    
+    #TODO: Add an equivalent method for embedding instruments automatically.
 
     @property
     def default_gauge_group(self):
@@ -1552,7 +1556,10 @@ def create_processor_spec(self, qudit_labels='auto'):
         gate_unitaries = _collections.OrderedDict()
         all_sslbls = self.state_space.sole_tensor_product_block_labels
         all_udims = [self.state_space.label_udimension(lbl) for lbl in all_sslbls]
+        sslbl_to_udim = {lbl:udim for lbl,udim in zip(all_sslbls, all_udims)}
+        all_qubits = all([udim == 2 for udim in all_udims])
         availability = {}
+        instrument_names = set(lbl.name for lbl in self.instruments)
 
         def extract_unitary(Umx, U_sslbls, extracted_sslbls):
             if extracted_sslbls is None: return Umx  # no extraction to be done
@@ -1575,7 +1582,7 @@ def extract_unitary(Umx, U_sslbls, extracted_sslbls):
                     U_extracted[ii, jj] = Umx[i, j]
             return U_extracted
 
-        def add_availability(opkey, op):
+        def add_availability_gate(opkey, op):
             if opkey == _Label(()) or opkey.is_simple:
                 if opkey == _Label(()):  # special case: turn empty tuple labels into "{idle}" gate in processor spec
                     gn = "{idle}"
@@ -1598,16 +1605,57 @@ def add_availability(opkey, op):
                             availability[gn].append(sslbls)
                     else:
                         availability[gn] = [sslbls]
+                
                 elif sslbls not in availability[gn]:
                     availability[gn].append(sslbls)
 
             else:  # a COMPOUND label with components => process each component separately
                 for component in opkey.components:
-                    add_availability(component, None)  # recursive call - the reason we need this to be a function!
+                    add_availability_gate(component, None)  # recursive call - the reason we need this to be a function!
+
+        def add_availability_inst(opkey):
+            if opkey.is_simple:
+                gn = opkey.name
+                sslbls = opkey.sslbls
+                #if sslbls is not None:
+                #    observed_sslbls.update(sslbls)
+                if gn in availability:
+                    if sslbls not in availability[gn]:
+                        availability[gn].append(sslbls)
+                else:
+                    availability[gn] = [sslbls]              
 
         #observed_sslbls = set()
         for opkey, op in self.operations.items():  # TODO: need to deal with special () idle label
-            add_availability(opkey, op)
+            add_availability_gate(opkey, op)
+
+        #process instruments:
+        nonstd_instruments = dict()
+        if self.instruments:
+            for instkey, inst in self.instruments.items():
+                add_availability_inst(instkey)
+
+            #create instrument specifiers.
+            seen_names = set()
+            from pygsti.processors.processorspec import InstrumentSpec as _InstrumentSpec
+            from pygsti.baseobjs.statespace import QubitSpace as _QubitSpace, QuditSpace as _QuditSpace
+            for instkey, inst in self.instruments.items():
+                instrument_specifiers = []
+                instrument_effect_labels = []
+                for lbl, inst_effect in inst.items():
+                    instrument_effect_labels.append(lbl)
+                    instrument_specifiers.append(inst_effect.to_dense())
+                name = instkey.name
+                if name not in seen_names:
+                    instkey.sslbls 
+                    if all_qubits:
+                        nonstd_instruments[name] = _InstrumentSpec(name, instrument_effect_labels, instrument_specifiers, 
+                                                                   _QubitSpace(len(instkey.sslbls)))
+                    else:
+                        
+                        nonstd_instruments[name] = _InstrumentSpec(name, instrument_effect_labels, instrument_specifiers, 
+                                                                   _QuditSpace(instkey.sslbls, [sslbl_to_udim[lbl] for lbl in instkey.sslbls]))
+                    seen_names.add(name)
 
         #Check that there aren't any undetermined unitaries
         unknown_unitaries = [k for k, v in gate_unitaries.items() if v is None]
@@ -1623,17 +1671,21 @@ def add_availability(opkey, op):
         if qudit_labels is None:  # special case of legacy explicit models where all gates have availability [None]
             qudit_labels = tuple(range(nqudits))
 
+        
+
         assert(len(qudit_labels) == nqudits), \
             "Length of `qudit_labels` must equal %d (not %d)!" % (nqudits, len(qudit_labels))
 
         if all([udim == 2 for udim in all_udims]):
             return _QubitProcessorSpec(nqudits, list(gate_unitaries.keys()), gate_unitaries, availability,
                                        qubit_labels=qudit_labels,
-                                       instrument_names=list(self.instruments.keys()), nonstd_instruments=self.instruments)
+                                       instrument_names=list(nonstd_instruments.keys()), 
+                                       nonstd_instruments=nonstd_instruments)
         else:
             return _QuditProcessorSpec(qudit_labels, all_udims, list(gate_unitaries.keys()), gate_unitaries,
                                        availability,
-                                       instrument_names=list(self.instruments.keys()), nonstd_instruments=self.instruments)
+                                       instrument_names=list(nonstd_instruments.keys()), 
+                                       nonstd_instruments=nonstd_instruments)
 
     def create_modelmember_graph(self):
         return _MMGraph({
diff --git a/pygsti/models/implicitmodel.py b/pygsti/models/implicitmodel.py
index 9ce554f7e..d4573e018 100644
--- a/pygsti/models/implicitmodel.py
+++ b/pygsti/models/implicitmodel.py
@@ -26,6 +26,8 @@
 from pygsti.models.layerrules import LayerRules as _LayerRules
 from pygsti.forwardsims.forwardsim import ForwardSimulator as _FSim
 
+from typing import Dict
+ordereddict = Dict
 
 class ImplicitOpModel(_mdl.OpModel):
     """
@@ -58,17 +60,12 @@ class ImplicitOpModel(_mdl.OpModel):
         represented, allowing compatibility checks with (super)operator
         objects.
     """
-    def __init__(self,
-                 state_space,
-                 layer_rules,
-                 basis="pp",
-                 simulator="auto",
-                 evotype="densitymx"):
-        self.prep_blks = _collections.OrderedDict()
-        self.povm_blks = _collections.OrderedDict()
-        self.operation_blks = _collections.OrderedDict()
-        self.instrument_blks = _collections.OrderedDict()
-        self.factories = _collections.OrderedDict()
+    def __init__(self, state_space, layer_rules, basis="pp", simulator="auto", evotype="densitymx"):
+        self.prep_blks: ordereddict = dict()
+        self.povm_blks: ordereddict = dict()
+        self.operation_blks: ordereddict = dict()
+        self.instrument_blks: ordereddict = dict()
+        self.factories: ordereddict = dict()
 
         super(ImplicitOpModel, self).__init__(state_space, basis, evotype, layer_rules, simulator)
 
@@ -370,27 +367,27 @@ def _from_nice_serialization(cls, state):
         return mdl
 
 
-def _init_spam_layers(model, prep_layers, povm_layers):
-    """ Helper function for initializing the .prep_blks and .povm_blks elements of an implicit model"""
-    # SPAM (same as for cloud noise model)
-    if prep_layers is None:
-        pass  # no prep layers
-    elif isinstance(prep_layers, dict):
-        for rhoname, layerop in prep_layers.items():
-            model.prep_blks['layers'][_Label(rhoname)] = layerop
-    elif isinstance(prep_layers, _op.LinearOperator):  # just a single layer op
-        model.prep_blks['layers'][_Label('rho0')] = prep_layers
-    else:  # assume prep_layers is an iterable of layers, e.g. isinstance(prep_layers, (list,tuple)):
-        for i, layerop in enumerate(prep_layers):
-            model.prep_blks['layers'][_Label("rho%d" % i)] = layerop
-
-    if povm_layers is None:
-        pass  # no povms
-    elif isinstance(povm_layers, _povm.POVM):  # just a single povm - must precede 'dict' test!
-        model.povm_blks['layers'][_Label('Mdefault')] = povm_layers
-    elif isinstance(povm_layers, dict):
-        for povmname, layerop in povm_layers.items():
-            model.povm_blks['layers'][_Label(povmname)] = layerop
-    else:  # assume povm_layers is an iterable of layers, e.g. isinstance(povm_layers, (list,tuple)):
-        for i, layerop in enumerate(povm_layers):
-            model.povm_blks['layers'][_Label("M%d" % i)] = layerop
+    def _init_spam_layers(self, prep_layers, povm_layers):
+        """ Helper function for initializing the .prep_blks and .povm_blks elements of an implicit model"""
+        # SPAM (same as for cloud noise model)
+        if prep_layers is None:
+            pass  # no prep layers
+        elif isinstance(prep_layers, dict):
+            for rhoname, layerop in prep_layers.items():
+                self.prep_blks['layers'][_Label(rhoname)] = layerop
+        elif isinstance(prep_layers, _op.LinearOperator):  # just a single layer op
+            self.prep_blks['layers'][_Label('rho0')] = prep_layers
+        else:  # assume prep_layers is an iterable of layers, e.g. isinstance(prep_layers, (list,tuple)):
+            for i, layerop in enumerate(prep_layers):
+                self.prep_blks['layers'][_Label("rho%d" % i)] = layerop
+
+        if povm_layers is None:
+            pass  # no povms
+        elif isinstance(povm_layers, _povm.POVM):  # just a single povm - must precede 'dict' test!
+            self.povm_blks['layers'][_Label('Mdefault')] = povm_layers
+        elif isinstance(povm_layers, dict):
+            for povmname, layerop in povm_layers.items():
+                self.povm_blks['layers'][_Label(povmname)] = layerop
+        else:  # assume povm_layers is an iterable of layers, e.g. isinstance(povm_layers, (list,tuple)):
+            for i, layerop in enumerate(povm_layers):
+                self.povm_blks['layers'][_Label("M%d" % i)] = layerop
diff --git a/pygsti/models/localnoisemodel.py b/pygsti/models/localnoisemodel.py
index 739fb8f7d..4802511ae 100644
--- a/pygsti/models/localnoisemodel.py
+++ b/pygsti/models/localnoisemodel.py
@@ -15,7 +15,7 @@
 import warnings as _warnings
 import numpy as _np
 
-from pygsti.models.implicitmodel import ImplicitOpModel as _ImplicitOpModel, _init_spam_layers
+from pygsti.models.implicitmodel import ImplicitOpModel as _ImplicitOpModel
 from pygsti.models.layerrules import LayerRules as _LayerRules
 from pygsti.models.memberdict import OrderedMemberDict as _OrderedMemberDict
 from pygsti.baseobjs import qubitgraph as _qgraph, statespace as _statespace
@@ -26,6 +26,7 @@
 from pygsti.modelmembers import operations as _op
 from pygsti.modelmembers import povms as _povm
 from pygsti.modelmembers import states as _state
+from pygsti.modelmembers import instruments as _instruments
 from pygsti.modelmembers.operations import opfactory as _opfactory
 from pygsti.modelmembers.modelmembergraph import ModelMemberGraph as _MMGraph
 from pygsti.baseobjs.basis import BuiltinBasis as _BuiltinBasis
@@ -38,6 +39,8 @@
 from pygsti.tools import listtools as _lt
 from pygsti.processors.processorspec import ProcessorSpec as _ProcessorSpec, QubitProcessorSpec as _QubitProcessorSpec
 
+from typing import Dict
+ordereddict = Dict
 
 class LocalNoiseModel(_ImplicitOpModel):
     """
@@ -54,8 +57,7 @@ class LocalNoiseModel(_ImplicitOpModel):
         processor.
 
     gatedict : dict
-        A dictionary (an `OrderedDict` if you care about insertion order) that
-        associates with gate names (e.g. `"Gx"`) :class:`LinearOperator`,
+        A dictionary that associates gate names (e.g. `"Gx"`) with :class:`LinearOperator` or
         `numpy.ndarray` objects. When the objects may act on fewer than the total
         number of qudits (determined by their dimension/shape) then they are
         repeatedly embedded into operation on the entire state space as specified
@@ -102,8 +104,7 @@ class LocalNoiseModel(_ImplicitOpModel):
           a circuit for any preparation and outcome.  High memory demand;
           best for a small number of (1 or 2) qubits.
         - "map" : op_matrix-state_vector products are repeatedly computed
-          to simulate circuits.  Slower for a small number of qubits, but
-          faster and more memory efficient for higher numbers of qubits (3+).
+          to simulate circuits.
 
     on_construction_error : {'raise','warn',ignore'}
         What to do when the conversion from a value in `gatedict` to a
@@ -142,8 +143,8 @@ class LocalNoiseModel(_ImplicitOpModel):
         this idle operation is *not* added to other layers as in `"add_global"`.
     """
 
-    def __init__(self, processor_spec, gatedict, prep_layers=None, povm_layers=None, evotype="default",
-                 simulator="auto", on_construction_error='raise',
+    def __init__(self, processor_spec, gatedict, instdict=None, prep_layers=None, povm_layers=None, 
+                 evotype="default", simulator="auto", on_construction_error='raise',
                  independent_gates=False, ensure_composed_gates=False, implicit_idle_mode="none"):
 
         qudit_labels = processor_spec.qudit_labels
@@ -158,7 +159,8 @@ def __init__(self, processor_spec, gatedict, prep_layers=None, povm_layers=None,
         # For later processing, we'll create mm_gatedict to contain each item as a ModelMember.  In local noise
         # models, these gates can be parameterized however the user desires - the LocalNoiseModel just embeds these
         # operators appropriately.
-        mm_gatedict = _collections.OrderedDict()  # ops as ModelMembers
+        mm_gatedict: ordereddict = dict()  # ops as ModelMembers
+        mm_instdict: ordereddict = dict()
 
         for key, gate in gatedict.items():
             if isinstance(gate, (_op.LinearOperator, _opfactory.OpFactory)):
@@ -167,8 +169,17 @@ def __init__(self, processor_spec, gatedict, prep_layers=None, povm_layers=None,
                 mm_gatedict[key] = _op.StaticArbitraryOp(gate, basis=None, evotype=evotype,
                                                          state_space=state_space)  # static gates by default
 
+
+        for key, inst in instdict.items():
+            if isinstance(inst, (_instruments.Instrument, _instruments.TPInstrument)):
+                mm_instdict[key] = inst
+            else:
+                if isinstance(inst, dict):
+                    mm_instdict[key] = _instruments.Instrument(inst)
+                else:
+                    raise ValueError(f'Unrecognized instdict member {inst}')
+
         self.processor_spec = processor_spec
-        idle_names = processor_spec.idle_gate_names
         global_idle_layer_label = processor_spec.global_idle_layer_label
 
         layer_rules = _SimpleCompLayerRules(qudit_labels, implicit_idle_mode, None, global_idle_layer_label)
@@ -178,23 +189,32 @@ def __init__(self, processor_spec, gatedict, prep_layers=None, povm_layers=None,
 
         flags = {'auto_embed': False, 'match_parent_statespace': False,
                  'match_parent_evotype': True, 'cast_to_type': None}
-        self.prep_blks['layers'] = _OrderedMemberDict(self, None, None, flags)
-        self.povm_blks['layers'] = _OrderedMemberDict(self, None, None, flags)
+        self.prep_blks['layers']= _OrderedMemberDict(self, None, None, flags)
+        self.povm_blks['layers']= _OrderedMemberDict(self, None, None, flags)
         self.operation_blks['gates'] = _OrderedMemberDict(self, None, None, flags)
         self.operation_blks['layers'] = _OrderedMemberDict(self, None, None, flags)
+        self.instrument_blks['gates'] = _OrderedMemberDict(self, None, None, flags)
         self.instrument_blks['layers'] = _OrderedMemberDict(self, None, None, flags)
         self.factories['gates'] = _OrderedMemberDict(self, None, None, flags)
         self.factories['layers'] = _OrderedMemberDict(self, None, None, flags)
 
-        _init_spam_layers(self, prep_layers, povm_layers)  # SPAM
+        self._init_spam_layers(self, prep_layers, povm_layers)  # SPAM
+
+        self._init_gate_layers(independent_gates, mm_gatedict, ensure_composed_gates, on_construction_error, implicit_idle_mode)
+        
+        self._init_instrument_layers(independent_gates, mm_instdict, on_construction_error)
+
+        self._clean_paramvec()
 
+    def _init_gate_layers(self, independent_gates, mm_gatedict, ensure_composed_gates, on_construction_error, implicit_idle_mode):
+        """Helper method for initialization of gate layers."""
         for gateName in self.processor_spec.gate_names:
             # process gate names (no sslbls, e.g. "Gx", not "Gx:0") - we'll check for the
             # latter when we process the corresponding gate name's availability
             gate_unitary = self.processor_spec.gate_unitaries[gateName]
             resolved_avail = self.processor_spec.resolved_availability(gateName)
 
-            gate_is_idle = gateName in idle_names
+            gate_is_idle = gateName in self.processor_spec.idle_gate_names
             gate_is_factory = callable(gate_unitary)
 
             if not independent_gates:  # then get our "template" gate ready
@@ -213,12 +233,12 @@ def __init__(self, processor_spec, gatedict, prep_layers=None, povm_layers=None,
                     else:
                         self.operation_blks['gates'][_Lbl(gateName)] = gate
 
-                    if gate_is_idle and gate.state_space.num_qudits == 1 and global_idle_layer_label is None:
+                    if gate_is_idle and gate.state_space.num_qudits == 1 and self.processor_spec.global_idle_layer_label is None:
                         # then attempt to turn this 1Q idle into a global idle (for implied idle layers)
-                        global_idle = _op.ComposedOp([_op.EmbeddedOp(state_space, (qlbl,), gate)
-                                                      for qlbl in qudit_labels])
+                        global_idle = _op.ComposedOp([_op.EmbeddedOp(self.state_space, (qlbl,), gate)
+                                                      for qlbl in self.processor_spec.qudit_labels])
                         self.operation_blks['layers'][_Lbl('{auto_global_idle}')] = global_idle
-                        global_idle_layer_label = layer_rules.global_idle_layer_label = _Lbl('{auto_global_idle}')
+                        global_idle_layer_label = self._layer_rules.global_idle_layer_label = _Lbl('{auto_global_idle}')
             else:
                 gate = None  # this is set to something useful in the "elif independent_gates" block below
 
@@ -233,7 +253,7 @@ def __init__(self, processor_spec, gatedict, prep_layers=None, povm_layers=None,
                 # when just ordering doesn't align (e.g. Gcnot:1:0 on 2-qudits needs to embed)
                 allowed_sslbls_fn = resolved_avail if callable(resolved_avail) else None
                 gate_nQudits = self.processor_spec.gate_num_qudits(gateName)
-                embedded_op = _opfactory.EmbeddingOpFactory(state_space, base_gate,
+                embedded_op = _opfactory.EmbeddingOpFactory(self.state_space, base_gate,
                                                             num_target_labels=gate_nQudits,
                                                             allowed_sslbls_fn=allowed_sslbls_fn)
                 self.factories['layers'][_Lbl(gateName)] = embedded_op
@@ -282,16 +302,16 @@ def __init__(self, processor_spec, gatedict, prep_layers=None, povm_layers=None,
                     # into inds (except in the special case inds[0] == '*' where we make an EmbeddingOpFactory)
                     try:
                         if gate_is_factory:
-                            if inds is None or inds == tuple(qudit_labels):  # then no need to embed
+                            if inds is None or inds == tuple(self.processor_spec.qudit_labels):  # then no need to embed
                                 embedded_op = base_gate
                             else:
-                                embedded_op = _opfactory.EmbeddedOpFactory(state_space, inds, base_gate)
+                                embedded_op = _opfactory.EmbeddedOpFactory(self.state_space, inds, base_gate)
                             self.factories['layers'][_Lbl(gateName, inds)] = embedded_op
                         else:
-                            if inds is None or inds == tuple(qudit_labels):  # then no need to embed
+                            if inds is None or inds == tuple(self.processor_spec.qudit_labels):  # then no need to embed
                                 embedded_op = base_gate
                             else:
-                                embedded_op = _op.EmbeddedOp(state_space, inds, base_gate)
+                                embedded_op = _op.EmbeddedOp(self.state_space, inds, base_gate)
                             self.operation_blks['layers'][_Lbl(gateName, inds)] = embedded_op
 
                             # If a 1Q idle gate (factories not supported yet) then turn this into a global idle
@@ -301,21 +321,63 @@ def __init__(self, processor_spec, gatedict, prep_layers=None, povm_layers=None,
 
                     except Exception as e:
                         if on_construction_error == 'warn':
-                            _warnings.warn("Failed to embed %s gate. Dropping it." % str(_Lbl(gateName, inds)))
-                        if on_construction_error in ('warn', 'ignore'): continue
-                        else: raise e
+                            _warnings.warn(f"Failed to embed {str(_Lbl(gateName, inds))} gate. Dropping it.")
+                        if on_construction_error in ('warn', 'ignore'): 
+                            continue
+                        else: 
+                            raise e
 
                 if len(singleQ_idle_layer_labels) > 0:
                     if implicit_idle_mode in ('add_global', 'only_global') and global_idle_layer_label is None:
                         # then create a global idle based on 1Q idle gates
                         global_idle = _op.ComposedOp([self.operation_blks['layers'][lbl]
                                                       for lbl in singleQ_idle_layer_labels.values()])
-                        global_idle_layer_label = layer_rules.global_idle_layer_label = _Lbl('{auto_global_idle}')
+                        global_idle_layer_label = self._layer_rules.global_idle_layer_label = _Lbl('{auto_global_idle}')
                         self.operation_blks['layers'][_Lbl('{auto_global_idle}')] = global_idle
                     elif implicit_idle_mode == 'pad_1Q':
-                        layer_rules.single_qubit_idle_layer_labels = singleQ_idle_layer_labels
+                        self._layer_rules.single_qubit_idle_layer_labels = singleQ_idle_layer_labels
 
-        self._clean_paramvec()
+    def _init_instrument_layers(self, independent_ops, mm_instdict, on_construction_error):
+        for instname in self.processor_spec.instrument_names:
+            inst_spec = self.processor_spec.instrument_specifiers[instname]
+            resolved_avail = self.processor_spec.resolved_availability(instname)
+
+            if not independent_ops:
+                inst = mm_instdict.get(instname, None)
+                if inst is not None:
+                    self.instrument_blks['gates'][_Lbl(instname)] = inst
+            else:
+                inst = None
+        
+        if callable(resolved_avail) or resolved_avail == '*':
+            raise NotImplementedError() #TODO: implement this case.
+        else:
+            for inds in resolved_avail:
+                if _Lbl(instname, inds) in mm_instdict and inds is not None:
+                    base_inst = mm_instdict[_Lbl(instname, inds)]
+                    self.instrument_blks['gates'][_Lbl(instname, inds)] = base_inst
+                elif independent_ops:
+                    base_inst = mm_instdict[instname].copy()
+                    self.instrument_blks['gates'][_Lbl(instname, inds)] = base_inst
+                else:
+                     base_inst = inst
+                
+                try:
+                    if inds is None or inds == tuple(self.processor_spec.qudit_labels):
+                        embedded_inst = base_inst
+                    else:
+                        embedded_instrument_members = dict()
+                        for key, member in base_inst.items():
+                            embedded_instrument_members[key] = _op.EmbeddedOp(self.state_space, inds, member)
+                        embedded_inst = _instruments.Instrument(embedded_instrument_members)
+                    self.instrument_blks['layers'][_Lbl(instname, inds)] = embedded_inst
+                except Exception as e:
+                    if on_construction_error == 'warn':
+                            _warnings.warn(f"Failed to embed {str(_Lbl(instname, inds))} gate. Dropping it.")
+                    if on_construction_error in ('warn', 'ignore'): 
+                        continue
+                    else: 
+                        raise e
 
     def create_processor_spec(self):
         import copy as _copy
@@ -354,6 +416,8 @@ def _from_nice_serialization(cls, state):
                                                          modelmembers.get('operation_blks|gates', []))
         mdl.operation_blks['layers'] = _OrderedMemberDict(mdl, None, None, flags,
                                                           modelmembers.get('operation_blks|layers', []))
+        mdl.instrument_blks['gates'] = _OrderedMemberDict(mdl, None, None, flags,
+                                                           modelmembers.get('instrument_blks|gates', []))
         mdl.instrument_blks['layers'] = _OrderedMemberDict(mdl, None, None, flags,
                                                            modelmembers.get('instrument_blks|layers', []))
         mdl.factories['gates'] = _OrderedMemberDict(mdl, None, None, flags, modelmembers.get('factories|gates', []))
diff --git a/pygsti/models/modelconstruction.py b/pygsti/models/modelconstruction.py
index 42c4763d2..25fdef610 100644
--- a/pygsti/models/modelconstruction.py
+++ b/pygsti/models/modelconstruction.py
@@ -736,19 +736,57 @@ def create_explicit_model(processor_spec, custom_gates=None,
                           lindblad_parameterization='auto',
                           evotype="default", simulator="auto",
                           ideal_gate_type='auto', ideal_spam_type='computational',
-                          embed_gates=False, basis='pp'):
+                          embed_gates=False, basis='pp',
+                          ideal_instrument_type='auto'):
+    """
+    Create a new `ExplicitOpModel` based on a processor specification.
+
+    Parameters
+    ----------
+    processor_spec
+
+    custom_gates
+
+    depolarization_strengths
+
+    stochastic_error_probs
+
+    lindblad_error_coeffs
+
+    lindblad_paramterization
+
+    evotype
+
+    simulator
+
+    ideal_gate_type
+
+    ideal_spam_type
+
+    embed_gates
+
+    basis
+
+    ideal_instrument_type
+
+
+    Returns
+    -------
+    ExplicitOpModel
+    """    
+
 
     modelnoise = _build_modelnoise_from_args(depolarization_strengths, stochastic_error_probs, lindblad_error_coeffs,
                                              depolarization_parameterization, stochastic_parameterization,
                                              lindblad_parameterization, allow_nonlocal=True)
 
     return _create_explicit_model(processor_spec, modelnoise, custom_gates, evotype,
-                                  simulator, ideal_gate_type, ideal_spam_type, ideal_spam_type, embed_gates, basis)
+                                  simulator, ideal_gate_type, ideal_spam_type, ideal_spam_type, embed_gates, basis, ideal_instrument_type)
 
 
 def _create_explicit_model(processor_spec, modelnoise, custom_gates=None, evotype="default", simulator="auto",
                            ideal_gate_type='auto', ideal_prep_type='auto', ideal_povm_type='auto',
-                           embed_gates=False, basis='pp'):
+                           embed_gates=False, basis='pp', ideal_instrument_type='auto'):
     qudit_labels = processor_spec.qudit_labels
     state_space = _statespace.QubitSpace(qudit_labels) if all([udim == 2 for udim in processor_spec.qudit_udims]) \
         else _statespace.QuditSpace(qudit_labels, processor_spec.qudit_udims)
@@ -766,7 +804,11 @@ def _create_explicit_model(processor_spec, modelnoise, custom_gates=None, evotyp
         ideal_prep_type = _state.state_type_from_op_type(ideal_gate_type)
     if ideal_povm_type == "auto":
         ideal_povm_type = _povm.povm_type_from_op_type(ideal_gate_type)
+    if ideal_instrument_type == 'auto':
+        ideal_instrument_type = _instrument.instrument_type_from_op_type(ideal_gate_type)
 
+    #-------------Gates-------------------#
+    #-------------------------------------#
     def _embed_unitary(statespace, target_labels, unitary):
         dummyop = _op.EmbeddedOp(statespace, target_labels,
                                  _op.StaticUnitaryOp(unitary, basis=basis, evotype="statevec_slow"))
@@ -871,128 +913,36 @@ def _embed_unitary(statespace, target_labels, unitary):
                     layer = _op.ComposedOp([ideal_gate, noiseop]) if (noiseop is not None) else ideal_gate
                     ret.operations[key] = layer
 
-    # Instruments:
+    #-----------Instruments-----------#:
+    #---------------------------------#
+    def _embed_instrument_op(statespace, target_labels, instrument_op):
+        dummyop = _op.EmbeddedOp(statespace, target_labels,
+                                 _op.StaticArbitraryOp(instrument_op, basis=basis))
+        return dummyop.to_dense("HilbertSchmidt")
+    
     for instrument_name in processor_spec.instrument_names:
         instrument_spec = processor_spec.instrument_specifier(instrument_name)
-
-        #FUTURE: allow instruments to be embedded
-        #resolved_avail = processor_spec.resolved_availability(instrument_name)
-        resolved_avail = [None]  # all instrument (so far) act on all the qudits
+        resolved_avail = processor_spec.resolved_availability(instrument_name)
+        dense_instrument_spec ={k: _bt.change_basis(mx, from_basis='std', to_basis=basis) 
+                                for k, mx in instrument_spec.to_dense_spec().items()}
 
         # resolved_avail is a list/tuple of available sslbls for the current gate/factory
         for inds in resolved_avail:  # inds are target qudit labels
-            key = _label.Label(instrument_name, inds)
-
-            if isinstance(instrument_spec, str):
-                if instrument_spec == "Iz": #TODO: Create a set of standard instruments the same way we handle gates.
-                    #NOTE: this is very inefficient currently - there should be a better way of
-                    # creating an Iz instrument in the FUTURE
-                    inst_members = {}
-                    if not all([udim == 2 for udim in processor_spec.qudit_udims]):
-                        raise NotImplementedError("'Iz' instrument can only be constructed on a space of *qubits*")
-                    for ekey, effect_vec in _povm.ComputationalBasisPOVM(nqubits=len(qudit_labels), evotype=evotype,
-                                                                         state_space=state_space).items():
-                        E = effect_vec.to_dense("HilbertSchmidt").reshape((state_space.dim, 1))
-                        inst_members[ekey] = _np.dot(E, E.T)  # (effect vector is a column vector)
-                    ideal_instrument = _instrument.Instrument(inst_members)
-                else:
-                    raise ValueError("Unrecognized instrument spec '%s'" % instrument_spec)
-
-            elif isinstance(instrument_spec, dict):
-
-                def _spec_to_densevec(spec, is_prep):
-                    num_qudits = len(qudit_labels)
-                    if isinstance(spec, str):
-                        if spec.isdigit():  # all([l in ('0', '1') for l in spec]): for qubits
-                            bydigit_index = spec
-                            assert (len(bydigit_index) == num_qudits), \
-                                "Wrong number of qudits in '%s': expected %d" % (spec, num_qudits)
-                            #Map a (possibly possibly mixed-base) string qudit state specifiers to a
-                            #integer mapping into the corresponding standard basis state.
-                            v = _np.zeros(state_space.udim)
-                            inc = _np.flip(_np.cumprod(list(reversed(processor_spec.qudit_udims[1:] + (1,)))))
-                            index = _np.dot(inc, list(map(int, bydigit_index)))
-                            v[index] = 1.0
-                        elif (not is_prep) and spec.startswith("E_") and spec[len('E_'):].isdigit():
-                            bydigit_index = spec[len('E_'):]
-                            assert (len(bydigit_index) == num_qudits), \
-                                "Wrong number of qudits in '%s': expected %d" % (spec, num_qudits)
-                            v = _np.zeros(state_space.udim)
-                            inc = _np.flip(_np.cumprod(list(reversed(processor_spec.qudit_udims[1:] + (1,)))))
-                            index = _np.dot(inc, list(map(int, bydigit_index)))
-                            v[index] = 1.0
-                        elif (not is_prep) and spec.startswith("E") and spec[len('E'):].isdigit():
-                            index = int(spec[len('E'):])
-                            assert (0 <= index < state_space.udim), \
-                                "Index in '%s' out of bounds for state space with udim %d" % (spec, state_space.udim)
-                            v = _np.zeros(state_space.udim); v[index] = 1.0
-                        elif is_prep and spec.startswith("rho_") and spec[len('rho_'):].isdigit():
-                            bydigit_index = spec[len('rho_'):]
-                            assert (len(bydigit_index) == num_qudits), \
-                                "Wrong number of qudits in '%s': expected %d" % (spec, num_qudits)
-                            v = _np.zeros(state_space.udim)
-                            inc = _np.flip(_np.cumprod(list(reversed(processor_spec.qudit_udims[1:] + (1,)))))
-                            index = _np.dot(inc, list(map(int, bydigit_index)))
-                            v[index] = 1.0
-                        elif is_prep and spec.startswith("rho") and spec[len('rho'):].isdigit():
-                            index = int(effect_spec[len('rho'):])
-                            assert (0 <= index < state_space.udim), \
-                                "Index in '%s' out of bounds for state space with udim %d" % (spec, state_space.udim)
-                            v = _np.zeros(state_space.udim); v[index] = 1.0
-                        else:
-                            raise ValueError("Unrecognized instrument member spec '%s'" % spec)
-                    elif isinstance(spec, _np.ndarray):
-                        assert (len(spec) == state_space.udim), \
-                            "Expected length-%d (not %d!) array to specify a state of %s" % (
-                                state_space.udim, len(spec), str(state_space))
-                        v = spec
-                    else:
-                        raise ValueError("Invalid effect or state prep spec: %s" % str(spec))
-
-                    return _bt.change_basis(_ot.state_to_dmvec(v), 'std', basis)
-
-                # elements are key, list-of-2-tuple pairs or numpy array 
-                inst_members = {}
-                for k, inst_effect_spec in instrument_spec.items():
-                    #one option is to specify the full dense instrument effect as a numpy array.
-                    if isinstance(inst_effect_spec, _np.ndarray):
-                        inst_members[k] = inst_effect_spec.copy()
-                        continue
-                    member = None
-                    if isinstance(inst_effect_spec, tuple):
-                        inst_effect_spec = [inst_effect_spec]
-                    elif isinstance(inst_effect_spec, list):
-                        #elements should be 2-tuples corresponding to effect specs and prep effects respectively.
-                        for (effect_spec, prep_spec) in inst_effect_spec:
-                            effect_vec = _spec_to_densevec(effect_spec, is_prep=False)
-                            prep_vec = _spec_to_densevec(prep_spec, is_prep=True)
-                            if member is None:
-                                member = _np.outer(effect_vec, prep_vec)
-                            else:
-                                member += _np.outer(effect_vec, prep_vec)
-                    else:
-                        raise ValueError('Unsupported instrument effect specification. See documentation of `QuditProcessorSpec` or `QubitProcessorSpec` for supported formats.')
-
-                    assert (member is not None), \
-                        "You must provide at least one rank-1 specifier for each instrument member!"
-                    inst_members[k] = member
-                ideal_instrument = _instrument.Instrument(inst_members)
-            else:
-                raise ValueError("Invalid instrument spec: %s" % str(instrument_spec))
+            inst_label = _label.Label(instrument_name, inds)
 
             if inds is None or inds == tuple(qudit_labels):  # then no need to embed
-                #ideal_gate = _op.create_from_unitary_mx(gate_unitary, ideal_gate_type, 'pp',
-                #                                        None, evotype, state_space)
-                pass  # ideal_instrument already created
+                ideal_instrument = _instrument.Instrument(dense_instrument_spec)
             else:
-                raise NotImplementedError("Embedded Instruments aren't supported yet")
-                # FUTURE: embed ideal_instrument onto qudits given by layer key (?)
+                embedded_instrument_members = dict()
+                for key, member in dense_instrument_spec.items():
+                    embedded_instrument_members[key] = _embed_instrument_op(state_space, inds, member)
+                ideal_instrument = _instrument.Instrument(embedded_instrument_members)
 
             #TODO: once we can compose instruments, compose with noise op here
             #noiseop = modelnoise.create_errormap(key, evotype, state_space, target_labels=inds)
             #layer = _op.ComposedOp([ideal_gate, noiseop]) if (noiseop is not None) else ideal_gate
-            layer = ideal_instrument
-            ret.instruments[key] = layer
+            ret.instruments[inst_label] = _instrument.convert(ideal_instrument, to_type=ideal_instrument_type,
+                                                              basis=basis) 
 
     # SPAM:
     local_noise = False; independent_gates = True; independent_spam = True
@@ -1006,12 +956,13 @@ def _spec_to_densevec(spec, is_prep):
 
     modelnoise.warn_about_zero_counters()
 
-    if ideal_gate_type == "full" and ideal_prep_type == "full" and ideal_povm_type == "full":
+    if ideal_gate_type == "full" and ideal_prep_type == "full" and ideal_povm_type == "full" and ideal_instrument_type=='full':
         ret.default_gauge_group = _gg.FullGaugeGroup(ret.state_space, basis, evotype)
-    elif (ideal_gate_type in ("full TP", "TP") and ideal_prep_type in ("full TP", "TP")
-          and ideal_povm_type in ("full TP", "TP")):
+    elif (ideal_gate_type in ("full TP", "TP", 'GLND') and ideal_prep_type in ("full TP", "TP", 'GLND')
+          and ideal_povm_type in ("full TP", "TP", 'GLND') and ideal_instrument_type in ('full TP', 'TP')):
         ret.default_gauge_group = _gg.TPGaugeGroup(ret.state_space, basis, evotype)
-    elif ideal_gate_type == "CPTP" and ideal_prep_type == "CPTP" and ideal_povm_type == "CPTP":
+    elif (ideal_gate_type in ("CPTP", "CPTPLND") and ideal_prep_type in ("CPTP", 'CPTPLND')  
+         and ideal_povm_type in ("CPTP", 'CPTPLND') and ideal_instrument_type in ('full TP', 'TP')):
         ret.default_gauge_group = _gg.UnitaryGaugeGroup(ret.state_space, basis, evotype)
     else:
         ret.default_gauge_group = _gg.TrivialGaugeGroup(ret.state_space)
diff --git a/pygsti/processors/processorspec.py b/pygsti/processors/processorspec.py
index 514dd2cd8..4f4a5ee60 100644
--- a/pygsti/processors/processorspec.py
+++ b/pygsti/processors/processorspec.py
@@ -15,6 +15,7 @@
 import collections as _collections
 import warnings as _warnings
 from functools import lru_cache
+from math import log
 
 
 from pygsti.tools import internalgates as _itgs
@@ -23,6 +24,7 @@
 from pygsti.tools import basistools as _bt
 from pygsti.baseobjs import qubitgraph as _qgraph
 from pygsti.baseobjs.label import Label as _Lbl
+from pygsti.baseobjs.statespace import QubitSpace as _QubitSpace
 from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable
 from pygsti.modelmembers.operations import LinearOperator as _LinearOp
 from pygsti.modelmembers.operations import FullArbitraryOp as _FullOp
@@ -187,10 +189,16 @@ def __init__(self, qudit_labels, qudit_udims, gate_names, nonstd_gate_unitaries=
             Any additional information that should be attached to this processor spec.
         """
         num_qudits = len(qudit_labels)
-        assert(len(qudit_udims) == num_qudits), "length of `qudit_labels` must equal that of `qubit_udims`!"
+        assert(len(qudit_udims) == num_qudits), "length of `qudit_labels` must equal that of `qudit_udims`!"
         assert(not (len(qudit_labels) > 1 and availability is None and geometry is None)), \
             "For multi-qudit processors you must specify either the geometry or the availability!"
 
+        first_udim = qudit_udims[0]
+        if not all([udim == first_udim for udim in qudit_udims]):
+            msg = 'Warning: automatic availability resolution for operations does not currently work for processor specifications'+\
+                  ' with mixed qudit dimensions. Please manually specify the availability for each operation.'
+            _warnings.warn(msg)
+
         if nonstd_gate_unitaries is None: nonstd_gate_unitaries = {}
 
         #Store inputs for adding models later
@@ -238,7 +246,18 @@ def __init__(self, qudit_labels, qudit_udims, gate_names, nonstd_gate_unitaries=
                 raise ValueError(
                     str(gname) + " is not a valid 'standard' gate name, it must be given in `nonstd_gate_unitaries`")
 
-        # Note: do *not* store the complex vectors defining states, POVM effects, and instrument members
+        #store the ideal instrument effects for each of the instruments.
+        std_instruments = _itgs.standard_instruments()
+        self.instrument_specifiers = dict()
+        for instname in instrument_names:
+            if nonstd_instruments is not None and instname in nonstd_instruments:
+                self.instrument_specifiers[instname] = nonstd_instruments[instname]
+            elif instname in std_instruments:
+                self.instrument_specifiers[instname] = std_instruments[instname]
+            else:
+                raise ValueError(f"{str(instname)} is not a valid 'standard' gate name, it must be given in `nonstd_instruments`")        
+
+        # Note: do *not* store the complex vectors defining states, POVM effects
         # as these can be large n-qudit vectors.
 
         # Set self.qudit_graph  (note QubitGraph is not actually qubit specific -- works fine for labelled qudits too)
@@ -257,7 +276,8 @@ def __init__(self, qudit_labels, qudit_udims, gate_names, nonstd_gate_unitaries=
         # Set availability
         if availability is None: 
             availability = {}
-        self.availability = {gatenm: availability.get(gatenm, 'all-edges') for gatenm in self.gate_names}
+        self.availability = {op_name: availability.get(op_name, 'all-edges') 
+                            for op_name in _itertools.chain(self.gate_names,self.instrument_names)}
         # if _Lbl(gatenm).sslbls is not None NEEDED?
 
         self.compiled_from = None  # could hold (QuditProcessorSpec, compilations) tuple if not None
@@ -309,7 +329,7 @@ def _serialize_instrument(obj):
 
         nonstd_preps = {k: _serialize_state(obj) for k, obj in self.nonstd_preps.items()}
         nonstd_povms = {k: _serialize_povm(obj) for k, obj in self.nonstd_povms.items()}
-        nonstd_instruments = {':'.join(map(str, k)): _serialize_instrument(obj) for k, obj in self.nonstd_instruments.items()}
+        nonstd_instruments = {k: obj.to_nice_serialization() for k, obj in self.nonstd_instruments.items()}
 
         state.update({'qudit_labels': list(self.qudit_labels),
                       'qudit_udims': list(self.qudit_udims),
@@ -384,7 +404,7 @@ def _unserialize_instrument(obj):
 
         nonstd_preps = {k: _unserialize_state(obj) for k, obj in state.get('nonstd_preps', {}).items()}
         nonstd_povms = {k: _unserialize_povm(obj) for k, obj in state.get('nonstd_povms', {}).items()}
-        nonstd_instruments = {tuple(k.split(':')): _unserialize_instrument(obj) for k, obj in state.get('nonstd_instruments', {}).items()}
+        nonstd_instruments = {k: InstrumentSpec.from_nice_serialization(obj) for k, obj in state.get('nonstd_instruments', {}).items()}
 
         return nonstd_gate_unitaries, nonstd_preps, nonstd_povms, nonstd_instruments
 
@@ -478,11 +498,10 @@ def instrument_specifier(self, name):
         -------
         str or dict
         """
-        if name in self.nonstd_instruments:
-            return self.nonstd_instruments[name]
+        if name in self.instrument_specifiers:
+            return self.instrument_specifiers[name]
         else:
-            # assert(is_standard_instrument_name(name)) TODO
-            return name
+            raise ValueError(f'{name} not found in `nonstd_instruments` and is not a built-in instrument name.')
 
     @property
     def primitive_op_labels(self):
@@ -494,23 +513,36 @@ def primitive_op_labels(self):
             ret.extend([_Lbl(gn, sslbls) for sslbls in avail])
         return tuple(ret)
 
-    def gate_num_qudits(self, gate_name):
+    def op_num_qudits(self, op_name):
         """
-        The number of qudits that a given gate acts upon.
+        The number of qudits that a given gate or instrument acts upon.
 
         Parameters
         ----------
-        gate_name : str
-            The name of the gate.
+        op_name : str
+            The name of the gate or instrument.
 
         Returns
         -------
         int
         """
-        unitary = self.gate_unitaries[gate_name]
-        if unitary is None: return len(self.qudit_labels)  # unitary=None => identity on all qudits
-        if isinstance(unitary, (int, _np.int64)): return unitary  # unitary=int => identity in n qudits
-        return int(round(_np.log2(unitary.shape[0])))  # possibly factory *function* SHAPE (unitary may be callable)
+        first_udim = self.qudit_udims[0]
+        if not all([udim == first_udim for udim in self.qudit_udims]):
+            msg = 'Automatic resolution of the number of qudits an for operation acts upon does not currently work' \
+                  'for processor specifications'+\
+                  ' with mixed qudit dimensions. Please manually specify the availability for each operation.'
+            raise ValueError(msg)
+
+        if op_name[0] == 'G':
+            unitary = self.gate_unitaries[op_name]
+            if isinstance(unitary, (int, _np.int64)): 
+                return unitary  # unitary=int => identity in n qudits
+            return int(round(log(unitary.shape[0], first_udim)))  # possibly factory *function* SHAPE (unitary may be callable)
+        elif op_name[0] == 'I':
+            instrument_spec = self.instrument_specifiers[op_name]
+            return instrument_spec.num_qudits  
+        else:
+            raise ValueError("Operation names should begin with either 'G' for gates or 'I' for instruments.")
 
     def rename_gate_inplace(self, existing_gate_name, new_gate_name):
         """
@@ -539,24 +571,24 @@ def rename(nm):
         self.gate_unitaries = {rename(k): v for k, v in self.gate_unitaries.items()}
         self.availability = {rename(k): v for k, v in self.availability.items()}
 
-    def resolved_availability(self, gate_name, tuple_or_function="auto"):
+    def resolved_availability(self, op_name, tuple_or_function="auto"):
         """
-        The availability of a given gate, resolved as either a tuple of sslbl-tuples or a function.
+        The availability of a given gate or instrument, resolved as either a tuple of sslbl-tuples or a function.
 
         This function does more than just access the `availability` attribute, as this may
-        hold special values like `"all-edges"`.  It takes the value of `self.availability[gate_name]`
+        hold special values like `"all-edges"`.  It takes the value of `self.availability[op_name]`
         and resolves and converts it into the desired format: either a tuple of state-space labels
         or a function with a single state-space-labels-tuple argument.
 
         Parameters
         ----------
-        gate_name : str
-            The gate name to get the availability of.
+        op_name : str
+            The gate or instrument name to get the availability of.
 
         tuple_or_function : {'tuple', 'function', 'auto'}
             The type of object to return.  `'tuple'` means a tuple of state space label tuples,
             e.g. `((0,1), (1,2))`.  `'function'` means a function that takes a single state
-            space label tuple argument and returns `True` or `False` to indicate whether the gate
+            space label tuple argument and returns `True` or `False` to indicate whether the gate or instrument
             is available on the given target labels.  If `'auto'` is given, then either a tuple or
             function is returned - whichever is more computationally convenient.
 
@@ -565,15 +597,25 @@ def resolved_availability(self, gate_name, tuple_or_function="auto"):
         tuple or function
         """
         assert(tuple_or_function in ('tuple', 'function', 'auto'))
-        avail_entry = self.availability.get(gate_name, 'all-edges')
-        gate_nqudits = self.gate_num_qudits(gate_name)
-        return self._resolve_availability(avail_entry, gate_nqudits, tuple_or_function)
-
-    def _resolve_availability(self, avail_entry, gate_nqudits, tuple_or_function="auto"):
+        avail_entry = self.availability.get(op_name, 'all-edges')
+        try:  
+            op_nqudits = self.op_num_qudits(op_name)
+        except ValueError:
+            #String values for availability all require logic we don't
+            #have for cases with mixed qudit dimension, which are the most likely
+            #reasons the above `op_num_qudits` case would fail.
+            msg = 'Automatic availability resolution for operations does not currently work for processor specifications'+\
+                  ' with mixed qudit dimensions. Please manually specify the availability for each operation'+\
+                  ' using the list of tuples format (a list of tuples of qudit targets).'
+            assert not isinstance(avail_entry, str) and not callable(avail_entry), msg        
+
+        return self._resolve_availability(avail_entry, op_nqudits, tuple_or_function)
+
+    def _resolve_availability(self, avail_entry, op_nqudits, tuple_or_function="auto"):
 
         if callable(avail_entry):  # a boolean function(sslbls)
             if tuple_or_function == "tuple":
-                return tuple([sslbls for sslbls in _itertools.permutations(self.qudit_labels, gate_nqudits)
+                return tuple([sslbls for sslbls in _itertools.permutations(self.qudit_labels, op_nqudits)
                               if avail_entry(sslbls)])
             return avail_entry  # "auto" also comes here
 
@@ -582,18 +624,18 @@ def _resolve_availability(self, avail_entry, gate_nqudits, tuple_or_function="au
                 def _f(sslbls):
                     return set(sslbls).issubset(self.qudit_labels) and tuple(sslbls) == tuple(sorted(sslbls))
                 return _f
-            return tuple(_itertools.combinations(self.qudit_labels, gate_nqudits))  # "auto" also comes here
+            return tuple(_itertools.combinations(self.qudit_labels, op_nqudits))  # "auto" also comes here
 
         elif avail_entry == 'all-permutations':
             if tuple_or_function == "function":
                 def _f(sslbls):
                     return set(sslbls).issubset(self.qudit_labels)
                 return _f
-            return tuple(_itertools.permutations(self.qudit_labels, gate_nqudits))  # "auto" also comes here
+            return tuple(_itertools.permutations(self.qudit_labels, op_nqudits))  # "auto" also comes here
 
         elif avail_entry == 'all-edges':
-            assert(gate_nqudits in (1, 2)), \
-                "I don't know how to place a %d-qudit gate on graph edges yet" % gate_nqudits
+            assert(op_nqudits in (1, 2)), \
+                "I don't know how to place a %d-qudit gate on graph edges yet" % op_nqudits
             if tuple_or_function == "function":
                 def _f(sslbls):
                     if len(sslbls) == 1: return True
@@ -602,8 +644,8 @@ def _f(sslbls):
                 return _f
 
             # "auto" also comes here:
-            if gate_nqudits == 1: return tuple([(i,) for i in self.qudit_labels])
-            elif gate_nqudits == 2: return tuple(self.qudit_graph.edges(double_for_undirected=True))
+            if op_nqudits == 1: return tuple([(i,) for i in self.qudit_labels])
+            elif op_nqudits == 2: return tuple(self.qudit_graph.edges(double_for_undirected=True))
             else: raise NotImplementedError()
 
         elif avail_entry in ('arbitrary', '*'):  # indicates user supplied factory determines allowed sslbls
@@ -824,7 +866,76 @@ def global_idle_layer_label(self):
                 return _Lbl(gn, self.qudit_labels)
         return None
 
+    def _inst_spec_to_dense_elems(self, instrument_spec):
+        """Helper method for converting an instrument spec into a dictionary of dense elements."""
+        # elements are key, list-of-2-tuple pairs or numpy array 
+        inst_members = {}
+        for k, inst_effect_spec in instrument_spec.items():
+            #one option is to specify the full dense instrument effect as a numpy array.
+            if isinstance(inst_effect_spec, _np.ndarray):
+                inst_members[k] = inst_effect_spec.copy()
+                continue
+            member = None
+            if isinstance(inst_effect_spec, tuple):
+                inst_effect_spec = [inst_effect_spec]
+            elif isinstance(inst_effect_spec, list):
+                #elements should be 2-tuples corresponding to effect specs and prep effects respectively.
+                for (effect_spec, prep_spec) in inst_effect_spec:
+                    effect_vec = _spec_to_densevec(effect_spec, is_prep=False)
+                    prep_vec = _spec_to_densevec(prep_spec, is_prep=True)
+                    if member is None:
+                        member = _np.outer(effect_vec, prep_vec)
+                    else:
+                        member += _np.outer(effect_vec, prep_vec)
+            else:
+                raise ValueError('Unsupported instrument effect specification. See documentation of `QuditProcessorSpec` or `QubitProcessorSpec` for supported formats.')
+
+            assert (member is not None), \
+                "You must provide at least one rank-1 specifier for each instrument member!"
+            inst_members[k] = member
+        return inst_members
+    
+    def _num_qudits_for_inst_spec(self, instrument_spec):
+        """A helper function which identifies the number of qudits an instrument specification acts on"""
+        first_udim = self.qudit_udims[0]
+        if not all([udim == first_udim for udim in self.qudit_udims]):
+            msg = 'Automatic resolution of the number of qudits an for operation acts upon does not currently work' \
+                  'for processor specifications'+\
+                  ' with mixed qudit dimensions. Please manually specify the availability for each operation.'
+            raise ValueError(msg)
+        # elements are key, list-of-2-tuple pairs or numpy array
+        # get the first instrument effect spec for the instrument spec.
+        first_inst_effect_spec = instrument_spec[next(iter(instrument_spec))]  
+
+        if isinstance(first_inst_effect_spec, _np.ndarray):
+            return int(round(log(_np.sqrt(first_inst_effect_spec.shape[0]), first_udim)))
+            
+        if isinstance(inst_effect_spec, tuple):
+            inst_effect_spec = [inst_effect_spec]
+
+        if isinstance(inst_effect_spec, list):
+            #elements should be 2-tuples corresponding to effect specs and prep effects respectively.
+            for (effect_spec, _) in inst_effect_spec:
+                effect_vec = _spec_to_densevec(effect_spec, is_prep=False)
+                if member is None:
+                    member = _np.outer(effect_vec, prep_vec)
+                else:
+                    member += _np.outer(effect_vec, prep_vec)
+        else:
+            raise ValueError('Unsupported instrument effect specification. See documentation of `QuditProcessorSpec` or `QubitProcessorSpec` for supported formats.')
+
+        assert (member is not None), \
+            "You must provide at least one rank-1 specifier for each instrument member!"
+        inst_members[k] = member
+        return inst_members
+
+
 
+
+def _num_qubits_for_spec(spec, qudit_labels, qudit_udims, state_space, is_prep):
+    """Helper function used in getting the number of qudits an instrument acts upon
+       from an instrument specification."""
+       
 class QubitProcessorSpec(QuditProcessorSpec):
     """
     The device specification for a one or more qudit quantum computer.
@@ -1367,3 +1478,240 @@ def compute_2Q_connectivity(self):
 
         return _qgraph.QubitGraph(qubit_labels, twoQ_connectivity)
     
+
+class InstrumentSpec(_NicelySerializable):
+    """
+    Class for storing a specifier for an ideal quantum instrument for use as part of
+    a `QuditProcessorSpec` or `QubitProcessorSpec`. This class enables compact storage
+    of data associated with ideal quantum instruments as well as methods for producing
+    dense or sparse representations of these instruments.
+    """
+
+    @classmethod
+    def cast(cls, obj):
+        """
+        Casts obj to an instance of `InstrumentSpec`, if possible.
+        """
+        pass
+
+
+    def __init__(self, name, inst_effect_labels, inst_effect_specs, state_space):
+        """
+        Create an instance of a specifier for ideal quantum instruments.
+        
+        Parameters
+        ----------
+        name : str
+            Name for the specified quantum instrument.
+
+        inst_effect_labels : list of str
+            A list of labels for each of the instrument effect elements.
+        
+        inst_effect_specs : list of lists of 2-tuples, or lists of ndarrays
+            The values this list give the specifications for each of the corresponding instrument elements labeled
+            in `inst_effect_labels`. The values of this list can be the following specifiers:
+
+            - numpy array: A numpy array corresponding to the dense representation of the instrument effect.
+            - lists of 2-tuples: Each tuple in this list of 2-tuples is such that the first element corresponds to a POVM effect
+              specifier (see `nonstd_povms` for supported options), and the second element is a state preparation specifier
+              (see `nonstd_preps` for supported options). These specifiers are used to construct appropriate effect and preparation
+              representations which are then have their outer product taken. This is done for each 2-tuple, and the outer products are
+              then summed to get the overall instrument effect.
+
+        state_space : `StateSpace`
+            The state space upon which this instrument acts.
+        """
+
+        self.name = name
+        self.instrument_effect_labels = list(inst_effect_labels)
+        self.instrument_spec = {key:val for key,val in zip(inst_effect_labels, inst_effect_specs)}
+        self.state_space = state_space
+
+        super().__init__()
+
+    @property
+    def num_qudits(self):
+        """
+        Number of qudits upon which this instrument acts.        
+        """
+        return self.state_space.num_qudits
+    
+    def to_dense_spec(self):
+        """
+        Method for converting this instrument spec into a dictionary of dense elements.
+        
+        Returns
+        -------
+        dict of np.ndarray
+            A dictionary whose keys are labels for instrument effects, and whose
+            values are numpy arrays with the dense representations of the instrument
+            effects specified in `self.instrument_spec`, in the standard/computational basis.        
+            
+        """
+        # elements are key, list-of-2-tuple pairs or numpy array 
+        inst_members = {}
+        for k, inst_effect_spec in self.instrument_spec.items():
+            #one option is to specify the full dense instrument effect as a numpy array.
+            if isinstance(inst_effect_spec, _np.ndarray):
+                inst_members[k] = inst_effect_spec.copy()
+                continue
+            member = None
+            if isinstance(inst_effect_spec, tuple):
+                inst_effect_spec = [inst_effect_spec]
+            if isinstance(inst_effect_spec, list):
+                #elements should be 2-tuples corresponding to effect specs and prep effects respectively.
+                for (effect_spec, prep_spec) in inst_effect_spec:
+                    effect_vec = self._inst_effect_spec_to_densevec(effect_spec, is_prep=False)
+                    prep_vec = self._inst_effect_spec_to_densevec(prep_spec, is_prep=True)
+                    if member is None:
+                        member = _np.outer(effect_vec, prep_vec)
+                    else:
+                        member += _np.outer(effect_vec, prep_vec)
+            else:
+                raise ValueError('Unsupported instrument effect specification. See documentation of `QuditProcessorSpec` or `QubitProcessorSpec` for supported formats.')
+
+            assert (member is not None), \
+                "You must provide at least one rank-1 specifier for each instrument member!"
+            inst_members[k] = member
+        return inst_members
+
+    def _inst_effect_spec_to_densevec(self, inst_effect_spec, is_prep):
+        """
+        Helper method used in converting parts of an instrument effect specification into the corresponding
+        dense vectors.
+
+        Parameters
+        ----------
+        inst_effect_spec : str or np.ndarray
+            Either a string specifier or a numpy array specifying part of instrument effect specifier.
+        
+        is_prep : bool
+            Flag indicating whether the input specification component is for an effect vector or a state
+            preparation.
+        """
+        if isinstance(inst_effect_spec, str):
+            if inst_effect_spec.isdigit():  # all([l in ('0', '1') for l in spec]): for qubits
+                bydigit_index = inst_effect_spec
+                assert (len(bydigit_index) == self.num_qudits), \
+                    "Wrong number of qudits in '%s': expected %d" % (inst_effect_spec, self.num_qudits)
+                #Map a (possibly possibly mixed-base) string qudit state specifiers to a
+                #integer mapping into the corresponding standard basis state.
+                v = _np.zeros(self.state_space.udim)
+                inc = _np.flip(_np.cumprod(list(reversed(self.state_space.qudit_udims[1:] + (1,)))))
+                index = _np.dot(inc, list(map(int, bydigit_index)))
+                v[index] = 1.0
+            elif (not is_prep) and inst_effect_spec.startswith("E_") and inst_effect_spec[len('E_'):].isdigit():
+                bydigit_index = inst_effect_spec[len('E_'):]
+                assert (len(bydigit_index) == self.num_qudits), \
+                    "Wrong number of qudits in '%s': expected %d" % (inst_effect_spec, self.num_qudits)
+                v = _np.zeros(self.state_space.udim)
+                inc = _np.flip(_np.cumprod(list(reversed(self.state_space.qudit_udims[1:] + (1,)))))
+                index = _np.dot(inc, list(map(int, bydigit_index)))
+                v[index] = 1.0
+            elif (not is_prep) and inst_effect_spec.startswith("E") and inst_effect_spec[len('E'):].isdigit():
+                index = int(inst_effect_spec[len('E'):])
+                assert (0 <= index < self.state_space.udim), \
+                    "Index in '%s' out of bounds for state space with udim %d" % (inst_effect_spec, self.state_space.udim)
+                v = _np.zeros(self.state_space.udim); v[index] = 1.0
+            elif is_prep and inst_effect_spec.startswith("rho_") and inst_effect_spec[len('rho_'):].isdigit():
+                bydigit_index = inst_effect_spec[len('rho_'):]
+                assert (len(bydigit_index) == self.num_qudits), \
+                    "Wrong number of qudits in '%s': expected %d" % (inst_effect_spec, self.num_qudits)
+                v = _np.zeros(self.state_space.udim)
+                inc = _np.flip(_np.cumprod(list(reversed(self.qudit_udims[1:] + (1,)))))
+                index = _np.dot(inc, list(map(int, bydigit_index)))
+                v[index] = 1.0
+            elif is_prep and inst_effect_spec.startswith("rho") and inst_effect_spec[len('rho'):].isdigit():
+                index = int(inst_effect_spec[len('rho'):])
+                assert (0 <= index < self.state_space.udim), \
+                    "Index in '%s' out of bounds for state space with udim %d" % (inst_effect_spec, self.state_space.udim)
+                v = _np.zeros(self.state_space.udim); v[index] = 1.0
+            else:
+                raise ValueError("Unrecognized instrument member spec '%s'" % inst_effect_spec)
+        elif isinstance(inst_effect_spec, _np.ndarray):
+            assert (len(inst_effect_spec) == self.state_space.udim), \
+                "Expected length-%d (not %d!) array to specify a state of %s" % (
+                    self.state_space.udim, len(inst_effect_spec), str(self.state_space))
+            v = inst_effect_spec
+        else:
+            raise ValueError("Invalid effect or state prep spec: %s" % str(inst_effect_spec))
+
+        return _ot.state_to_dmvec(v)
+    
+    def _to_nice_serialization(self):
+        state = super()._to_nice_serialization()
+
+        def _serialize_state(obj):
+            return (obj.to_nice_serialization() if isinstance(obj, _NicelySerializable)
+                    else (obj if isinstance(obj, str) else self._encodemx(obj)))
+
+        # NicelySerializable is commented out while ModelMembers inherit from it but do not implement
+        # a non-base to_nice_serialization() method
+        def _serialize_instrument_member(obj):
+            if isinstance(obj, str): 
+                return obj
+            #if isinstance(obj, _NicelySerializable): return obj.to_nice_serialization()
+            if isinstance(obj, tuple):
+                obj = [obj]
+            if isinstance(obj, list):
+                assert(all([isinstance(rank1op_spec, (list, tuple)) and len(rank1op_spec) == 2
+                           for rank1op_spec in obj]))
+                return [(_serialize_state(espec), _serialize_state(rspec)) for espec, rspec in obj]
+            if isinstance(obj, _np.ndarray):
+                return [self._encodemx(obj)]
+            raise ValueError("Cannot serialize Instrument member specifier of type %s!" % str(type(obj)))
+
+        state.update({
+            'name': self.name,
+            'instrument_effect_labels': self.instrument_effect_labels,
+            'instrument_effect_specs': [_serialize_instrument_member(val) for val in self.instrument_spec.values()],
+            'state_space' : self.state_space._to_nice_serialization()
+        })
+
+        return state
+    
+    @classmethod
+    def _from_nice_serialization(cls, state):
+
+        from pygsti.baseobjs.statespace import QubitSpace as _QubitSpace, ExplicitStateSpace as _ExplicitStateSpace, QuditSpace as _QuditSpace
+        #HACK: Replace this with call to _state_class
+        if 'unitary_space_dimensions' in state['state_space']:
+            state_space = _ExplicitStateSpace._from_nice_serialization(state['state_space'])
+        elif 'qudit_labels' in state['state_space']:
+            state_space = _QuditSpace._from_nice_serialization(state['state_space'])
+        elif 'qubit_labels' in state['state_space']:
+            state_space = _QubitSpace._from_nice_serialization(state['state_space'])
+        else:
+            raise ValueError('Cannot deserialize state space object.')
+
+        def _unserialize_state(obj):
+            if isinstance(obj, str): return obj
+            elif isinstance(obj, dict) and "module" in obj:  # then a NicelySerializable object
+                return _NicelySerializable.from_nice_serialization(obj)
+            else:  # assume a matrix encoding of some sort (could be list or dict)
+                return cls._decodemx(obj)
+
+        def _unserialize_instrument_member(obj):
+            if isinstance(obj, str): return obj
+            elif isinstance(obj, dict) and "module" in obj:  # then a NicelySerializable object
+                return _NicelySerializable.from_nice_serialization(obj)
+            elif isinstance(obj, list):
+                if obj and isinstance(obj[0], tuple):
+                    return [(_unserialize_state(espec), _unserialize_state(rspec)) for espec, rspec in obj]
+                else: #assume this is a serialized ndarray
+                    return _unserialize_state(obj[0])
+            raise ValueError("Cannot unserialize Instrument member specifier of type %s!" % str(type(obj)))
+
+        instrument_effect_specs = [_unserialize_instrument_member(inst_effect_spec) for inst_effect_spec in state['instrument_effect_specs']]
+        instrument_effect_labels = state['instrument_effect_labels']
+        name = state['name']
+        return cls(name, instrument_effect_labels, instrument_effect_specs, state_space) 
+    
+    
+
+
+
+
+
+
+
diff --git a/pygsti/tools/internalgates.py b/pygsti/tools/internalgates.py
index 6439cd3ce..e197b39c5 100644
--- a/pygsti/tools/internalgates.py
+++ b/pygsti/tools/internalgates.py
@@ -340,6 +340,29 @@ def u_op(exp):
 
     return std_unitaries
 
+def standard_instruments():
+    """
+    Construct a dictionary of standard instrument elements for commonly used quantum instruments.
+
+    Returns:
+    --------
+    std_instruments: dict of dicts
+        Returns a dictionary whose keys are strings corresponding to recognized built-in instruments.
+        The values of the dictionary are themselves dictionaries whose keys are instrument effect labels
+        and whose values are numpy arrays corresponding to the ideal instrument elements in the pauli-product
+        basis.
+    """
+    from pygsti.processors.processorspec import InstrumentSpec as _InstrumentSpec
+    from pygsti.baseobjs.statespace import QubitSpace as _QubitSpace
+    std_instruments = dict()
+    #Iz
+    Iz_inst_members = {}
+    Iz_effect_specs = [[(_np.array([1,0]), _np.array([1,0]))], [(_np.array([0,1]), _np.array([0,1]))]]
+    Iz_effect_lbls = ['p0', 'p1'] 
+    std_instruments['Iz'] = _InstrumentSpec(name='Iz', inst_effect_labels=Iz_effect_lbls, inst_effect_specs=Iz_effect_specs,
+                                            state_space= _QubitSpace(1))
+    #TODO: n-qubit Iz instruments to match old conventions.
+    return std_instruments
 
 def unitary_to_standard_gatename(unitary, up_to_phase = False, return_phase = False):
     """